1_Process_data.py

#!/usr/bin/env python
# coding: utf-8

# In[1]: 0_process_data


import copy as cp
import itertools as itt
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import os
import pandas as pd
import pickle
import scipy.optimize as spo
import scipy.special as spsp
import scipy.stats as sps
import seaborn as sns

from matplotlib.ticker import FixedLocator

mpl.rcParams['axes.titlesize'] = 'xx-large'
mpl.rcParams['axes.labelsize'] = 'xx-large'
mpl.rcParams['xtick.labelsize'] = 'x-large'
mpl.rcParams['ytick.labelsize'] = 'x-large'
mpl.rcParams['xtick.direction'] = 'out'
mpl.rcParams['ytick.direction'] = 'out'
mpl.rcParams['legend.frameon'] = True
mpl.rcParams['legend.framealpha'] = 0.5
mpl.rcParams['legend.fontsize'] = 'large'

mydir = os.path.dirname(os.getcwd())
path_proj = os.path.join(mydir, 'MultiGen_analysis')
path_sc = os.path.join(path_proj, 'csv', 'Single_cell')
path_pickle = os.path.join(path_proj, 'pickled_data')
path_plot = os.path.join(path_proj, 'figures')


# In[2]:


#(USER-DEFINED) cell_cols is a dictionary that assign to each cell type a color for the plots. 
#cell_class_exp_time is a function that determines a cell's class upon its markers. The user define thresholds from Gating_matrix, which are imported in the next cell.

#The following snippets are for the analysis of primary cells, from human HSPCs

cell_cols = {'HSC':'seagreen', 'MPP':'Limegreen', 'LMPP':'Red', 'CMP':'Deepskyblue','GMP':'Violet', 'MEP':'Royalblue', 'CD34-':'aqua'}
sns.palplot(cell_cols.values())

def cell_class_exp_time(df, thr_cd34, thr_cd38, thr_cd90, thr_cd45l, thr_cd45h, thr_cd123): #the order of the thresholds in the Gating Matrix.csv needs to be retained
    if df['CD34'] > thr_cd34:
        if df['CD38'] <= thr_cd38:
            if df['CD90'] > thr_cd90:
                return {'class':'HSC', 'color':cell_cols['HSC'], 'rank':1}
            else:
                if df['CD45'] > thr_cd45h:
                    return {'class':'LMPP', 'color':cell_cols['LMPP'], 'rank':2}
                else:
                    return {'class':'MPP', 'color':cell_cols['MPP'], 'rank':3}
        else:
            if df['CD123'] > thr_cd123:
                if df['CD45'] > thr_cd45l:
                    return {'class':'GMP', 'color':cell_cols['GMP'], 'rank':4}
                else:
                    return {'class':'CMP', 'color':cell_cols['CMP'], 'rank':5}
            else:
                return {'class':'MEP', 'color':cell_cols['MEP'], 'rank':6}
    
    else:
        return {'class':'CD34-', 'color':cell_cols['CD34-'], 'rank':7}


#The following snippets are for a simpler analysis, using only CD34 and CD38 to distinguish more immature and mature progenitors

#cell_cols = {'CD34+CD38-':'seagreen', 'CD34+CD38+':'Limegreen', 'CD34-CD38+':'hotpink', 'CD34-CD38-':'midnightblue'}
#sns.palplot(cell_cols.values())

#def cell_class_exp_time(df, thr_cd34, thr_MHCII, thr_cd38, thr_PDL1, thr_cKit):#, thr_cd10, thr_cd45l, thr_cd45h, thr_cd90, thr_cd123=-1):
#    if df['CD34'] > thr_cd34:
#        if df['CD38'] <= thr_cd38:
#            return {'class':'CD34+CD38-', 'color':cell_cols['CD34+CD38-'], 'rank':1}
#        else:
#            return {'class':'CD34+CD38+', 'color':cell_cols['CD34+CD38+'], 'rank':2}
               
#    else:
#         if df['CD38'] > thr_cd38:
#             return {'class':'CD34-CD38+', 'color':cell_cols['CD34-CD38+'], 'rank':3}
#         else:
#             return {'class':'CD34-CD38-', 'color':cell_cols['CD34-CD38-'], 'rank':4}
#
#    return {'class':'NC', 'color':'Black', 'rank':100}#'NC'


#The following snippets are using CD34 only

#cell_cols = {'CD34+':'seagreen', 'CD34-':'hotpink'}
#sns.palplot(cell_cols.values())

#def cell_class_exp_time(df, thr_cd34, thr_MHCII, thr_cd38, thr_PDL1):#, thr_cd10, thr_cd45l, thr_cd45h, thr_cd90, thr_cd123=-1):
#    if df['CD34'] > thr_cd34:
#        return {'class':'CD34+', 'color':cell_cols['CD34+'], 'rank':1}
#    else:
#        return {'class':'CD34-', 'color':cell_cols['CD34-'], 'rank':2}
#
#    return {'class':'NC', 'color':'Black', 'rank':100}#'NC'

def cell_class(df):
    if df.Experiment+'_'+df.Culture_time in dct_thr.keys():
        return cell_class_exp_time(df, *(dct_thr[df.Experiment+'_'+df.Culture_time]))
    else:
        return {'class':'Experiment_or_time_not_found', 'color':'Black', 'rank':100}#'Experiment_or_time_not_found'    

df_gating = pd.read_excel(os.path.join(path_proj,'Gating_matrix.xlsx'), index_col=0)
dct_thr = {exp_time:df_gating[exp_time].values for exp_time in df_gating.columns}
dct_thr

# In[3]:


#(EXPERIMENT SPECIFIC) conds and cond_rule are used later to rename the culture conditions
def cond_rule(cond):
    conds = ['GT', 'Diff']
    cond_map = {'GT': conds[0], 'Diff': conds[1]}
    return cond_map.get(cond, 'NA')

#Impossible families are families that are including too many cells per generation (like three cells in Gen 1)     
    
def find_impossible_families(df):
    fams = np.unique(df.Family)
    lst_fam_vec = [np.unique(df[df.Family==fam].Generation, return_counts=True) for fam in fams]
    cohort_number = np.array([(el[1]/np.power(2.,el[0])).sum() for el in lst_fam_vec])
    if all(cohort_number<=1):
        print('No impossible families detected')
        return []
    else:
        impossible_families = fams[cohort_number>1]
        print('Impossible families:', impossible_families)
        return impossible_families


# In[4]:


#Import the single cell data from the flow cytometry experiment

lst_sc_files = [file for file in os.listdir(path_sc) if '.csv' in file and file!='Pooled_data.csv']
df_sc_lst = []
for file_name in lst_sc_files:
    print(f"Processing file: {file_name}")
    
#(EXPERIMENT SPECIFIC) as the filenames here ends with '.csv'
    experiment_name = file_name[:-4]

# Read file into DataFrame
    try:
        df = pd.read_csv(os.path.join(path_sc, file_name), sep=',', decimal='.')
    except Exception as e:
        print(f"Error reading {file_name}: {e}")
        continue
    
# Verify required columns exist
    required_columns = ['Condition', 'Well', 'Color', 'Culture_time', 'Original_cell']
    if not all(col in df.columns for col in required_columns):
        print(f"File {file_name} missing required columns. Skipping.")
        continue    

#(EXPERIMENT SPECIFIC) sep=',', decimal='.' are unusual
    df = pd.read_csv(os.path.join(path_sc, file_name), sep=',', decimal='.')
    print(df.columns)
    
#Here the user should decide how to format the imported data
    df = df.assign(
        Culture_condition=df.apply(lambda r: cond_rule(r.Condition), axis=1),
        Experiment=[experiment_name for k in range(len(df))]
    )
    df = df.assign(
        Well_experiment=lambda r:r.Well+'_'+r.Experiment,
        Family=lambda r: r.Well + r.Color + r.Culture_time + r.Original_cell + r.Experiment,
        Class=df.apply(func=lambda r:cell_class(r)['class'], axis=1),
        Cell_color=df.apply(func=lambda r:cell_class(r)['color'], axis=1),
        Cell_rank=df.apply(func=lambda r:cell_class(r)['rank'], axis=1)
    )
   
#Detect and remove impossible families
    impossible_fams = find_impossible_families(df)
    df = df[~(df.Family.isin(impossible_fams))]
    df_sc_lst.append(df)
 
# Concatenate and save the pooled data
df_pool = pd.concat(df_sc_lst, ignore_index=True, sort=True)
output_path = os.path.join(path_sc, 'Pooled_data.csv')
df_pool.to_csv(output_path, sep=',', decimal='.', index=False)
print(f"Pooled data saved to {output_path}")


# In[5]:
  
#7/len(np.unique(df_pool.Family)) * 100

np.unique(df_pool.Class, return_counts=True), np.unique(df_pool.Culture_time, return_counts=True)