Merge pull request #8 from heberlr/development

Updates for the data analysis.

Author: heberlr

.gitignore

@@ -1,9 +1,13 @@
 *.o
 *.exe
 *.pyc
+*.log
 /output_*
-config/tmp.xml
-*.bmp
+/output/*.svg
+/output/*.xml
+/output/*.mat
+/output/*.bmp
+/output/*.dat
 /Results/
 *.pickle
-*.ipynb_checkpoints*
+*.ipynb_checkpoints*

Data_Analysis/DataSet.py

@@ -3,7 +3,14 @@
 import pickle
 import wget, os
 import torchvision.transforms as transforms
-from scipy.stats import ttest_ind
+from scipy.stats import ttest_ind, mannwhitneyu
+from fitter import Fitter
+
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib.ticker import ScalarFormatter
+from matplotlib_venn import venn3_unweighted
+import seaborn as sns
 
 NameParameters = ["macrophage_max_recruitment_rate","macrophage_recruitment_min_signal","macrophage_recruitment_saturation_signal","DC_max_recruitment_rate","DC_recruitment_min_signal","DC_recruitment_saturation_signal","DC_leave_prob","TC_death_rate","T_Cell_Recruitment","DM_decay"]
 NC_Patients_samples = np.array([515, 723, 1810, 2539, 2723, 2958, 3160, 5379, 6040, 7642]) # (No control)
@@ -40,16 +47,20 @@ def Loading_AditionalDataset():
     df = loading_pkl_file(local_path)
     return df
 
-def Loading_UMAP():
+def Loading_UMAP_traj():
     # Download the UMAP of trajectories
     local_path = 'Clust_UMAP_15_01_2.pickle' # trajectories of all simulations
     link = 'https://figshare.com/ndownloader/files/38481611'
+    download_file(local_path, link)
+    # return loading_pkl_file(local_path)
+    return pd.read_pickle(local_path)
+
+def Loading_UMAP_avg():
     local_path2 = 'Patient_UMAP_15_01_2.pickle' # trajectory averages of all simulations
     link2 = 'https://figshare.com/ndownloader/files/38898579'
-    download_file(local_path, link)
     download_file(local_path2, link2)
-    return loading_pkl_file(local_path), loading_pkl_file(local_path2)
-
+    return loading_pkl_file(local_path2)
+    
 def ZscoreImagesToTensor(dataframe): # Normalize images (zscore) and convert to tensor (pytorch pattern)
     images = np.asarray(dataframe['imagens'].tolist(), dtype="float64") # shape (1e4, 5, 5, 3)
     images_norm = np.zeros(shape = (images.shape[0],images.shape[3],images.shape[1],images.shape[2]), dtype="float64") # shape (1e4, 3, 5, 5)
@@ -116,22 +127,45 @@ def Normalize_Parameters(df_input, Pat_df_3classes):
 
 def T_test(df_input_min_max_scaled):
     df_SC_NC_min_max_scaled = df_input_min_max_scaled.loc[(df_input_min_max_scaled["label"] == 'NC') | (df_input_min_max_scaled["label"] == 'SC')]
-    Parameters_data = df_SC_NC_min_max_scaled.to_numpy()
-    # remove sample ID and label
-    NC_Parameters = Parameters_data[Parameters_data[:,-1] == 'NC'][:,1:-1]
-    SC_Parameters = Parameters_data[Parameters_data[:,-1] == 'SC'][:,1:-1]
+    Parameters_data = df_SC_NC_min_max_scaled
     T_test = []
-    for i in range(len(NameParameters)):
-        T_test.append( ttest_ind(NC_Parameters[:,i],SC_Parameters[:,i]) )
-    dict_tTest = {NameParameters[i]: T_test[i][0] for i in range(len(NameParameters))}
-    dict_pValue = {NameParameters[i]: T_test[i][1] for i in range(len(NameParameters))}
+    for parName in NameParameters:
+        T_test.append( ttest_ind(Parameters_data.loc[Parameters_data['label'] == 'NC'][parName], Parameters_data.loc[Parameters_data['label'] == 'SC'][parName] ) )
+    dict_tTest = {NameParameters[i]: T_test[i].statistic for i in range(len(NameParameters))}
+    dict_pValue = {NameParameters[i]: T_test[i].pvalue for i in range(len(NameParameters))}
     return pd.DataFrame([dict_tTest, dict_pValue],index=['t-stat', 'p-value'])
+    
+def MannWhitneyU_test(df_input_min_max_scaled):
+    df_SC_NC_min_max_scaled = df_input_min_max_scaled.loc[(df_input_min_max_scaled["label"] == 'NC') | (df_input_min_max_scaled["label"] == 'SC')]
+    Parameters_data = df_SC_NC_min_max_scaled
+    MWU_test = []
+    for parName in NameParameters:
+        MWU_test.append( mannwhitneyu(Parameters_data.loc[Parameters_data['label'] == 'NC'][parName], Parameters_data.loc[Parameters_data['label'] == 'SC'][parName], method="asymptotic" ) )
+    dict_MWUTest = {NameParameters[i]: MWU_test[i].statistic for i in range(len(NameParameters))}
+    dict_pValue = {NameParameters[i]: MWU_test[i].pvalue for i in range(len(NameParameters))}
+    return pd.DataFrame([dict_MWUTest, dict_pValue],index=['MWU-stat', 'p-value'])
+
+def Fitter_dist(df_input_min_max_scaled):
+    dic_Fit = {}
+    fileName = 'FitterPars.pickle'
+    if not os.path.exists(fileName):
+        df_SC_NC_min_max_scaled = df_input_min_max_scaled.loc[(df_input_min_max_scaled["label"] == 'NC') | (df_input_min_max_scaled["label"] == 'SC')]
+        data_plot = pd.melt(df_SC_NC_min_max_scaled, id_vars=['sample','label'], value_vars=NameParameters)
+        for parName in ['DC_max_recruitment_rate','macrophage_max_recruitment_rate']:
+            f_SC = Fitter(data_plot.loc[(data_plot['label'] == 'SC') & (data_plot['variable'] == parName)]['value'].to_numpy()); f_SC.fit()
+            f_NC = Fitter(data_plot.loc[(data_plot['label'] == 'NC') & (data_plot['variable'] == parName)]['value'].to_numpy()); f_NC.fit()
+            dic_Fit[parName+'SC'] = f_SC
+            dic_Fit[parName+'NC'] = f_NC
+        with open(fileName, 'wb') as handle:
+            pickle.dump(dic_Fit, handle, protocol=pickle.HIGHEST_PROTOCOL)
+    else:
+        dic_Fit = loading_pkl_file(fileName)
+    return dic_Fit
 
 def PatientsAnalysis():
     df = Loading_AditionalDataset()
     Mean_liveCancer_NC = []; Std_liveCancer_NC = []; ParMac_NC = []; ParDC_NC = []; Patient_NC = []; PatientLabel_NC = []
     Mean_liveCancer_MC = []; Std_liveCancer_MC = []; ParMac_MC = []; ParDC_MC = []; Patient_MC = []; PatientLabel_MC = []
-
     for sample in df['sample'].unique():
         if ( df.loc[(df['sample'] == sample) & (df['replicate'] == 0)]['sample_ref'].to_list() in  NC_Patients_samples): # (No control patients)
             Mean_liveCancer_NC.append(np.array(df.loc[df['sample'] == sample]['cancer_cell_day10']).mean())

Data_Analysis/Movie.py

@@ -8,7 +8,7 @@ def annotate_axes(ax, text, fontsize=18):
 def PlotFrame(popcells, imagecell, frameID, FigName=None):
   # ax0: Tumor microenvironment
   # ax1: Legend
-  # ax2: Melanoma and healthy tissue population
+  # ax2: Cancer and healthy tissue population
   # ax3: Lymph Node population
   # ax4: Immune cells population
   fig = plt.figure(layout=None, facecolor='0.9',figsize=(12,12))
@@ -25,7 +25,7 @@ def PlotFrame(popcells, imagecell, frameID, FigName=None):
   ax1.axes.yaxis.set_visible(False)
   # Population plots
   colours = {'live_lung': 'b', 'live_cancer': 'y', 'inac_DC': '#810F7C', 'act_DC': '#ff1493', 'inac_Mac': '#238B45', 'act_Mac': '#C0FF00', 'exas_Mac': '#A8DD76', 'hyper_MAC': '#A8DDB5', 'live_CD4': 'orange','live_CD8': 'red', 'LN_DC':'#016ca4', 'LN_TC':'#cd5000', 'LN_TH1':'#ffbb7b', 'LN_TH2':'#595959', 'LN_TCt':'#ababab', 'LN_Tht':'#a1c8f0'}
-  # Melanoma and Lung
+  # Cancer and Lung
   ax2 = fig.add_subplot(gs[-1, 0])
   ax2.plot(popcells[:,0]/1440.0, popcells[:,7], color=colours['live_lung'], label='live_lung')
   ax2.plot(popcells[:,0]/1440.0, popcells[:,14], color=colours['live_cancer'], label='live_cancer')

Data_Analysis/PlotResults.ipynb

@@ -1,14 +1,8 @@
-import matplotlib.pyplot as plt
-import matplotlib.gridspec as gridspec
-from matplotlib_venn import venn3_unweighted
-import seaborn as sns
-import numpy as np
-import pandas as pd
 from DataSet import *
 
 colours = {'NC':'#CA8151','MC':'#FFE1B3','SC':'#307CA1','NC+MC':'#999999','NC+SC':'#797979','MC+SC':'#898989','NC+MC+SC':'#696969','Mixed':'#595959', 'Mixed/MC': '#494949'}
 
-def Plot_Patients(Sim_df_lastFrame,Pat_df_7classes,Pat_df_3classes, Title=False, FigName=None):
+def plot_Patients(Sim_df_lastFrame,Pat_df_7classes,Pat_df_3classes, Title=False, FigName=None):
   fig, (ax1,ax2,ax3) = plt.subplots(1,3,figsize=(16, 6))
   Sim_sizes = Sim_df_lastFrame['label'].value_counts()
   Pat_7classes_sizes =  Pat_df_7classes['label'].value_counts()
@@ -30,59 +24,71 @@ def Plot_Patients(Sim_df_lastFrame,Pat_df_7classes,Pat_df_3classes, Title=False,
   if (FigName): plt.savefig(FigName)
   else: plt.show()
 
-def Plot_Parameters(df_input_min_max_scaled, FigName=None):
+def plot_Parameters(df_input_min_max_scaled, FigName=None):
   df_SC_NC_min_max_scaled = df_input_min_max_scaled.loc[(df_input_min_max_scaled["label"] == 'NC') | (df_input_min_max_scaled["label"] == 'SC')]
   data_plot = pd.melt(df_SC_NC_min_max_scaled, id_vars=['sample','label'], value_vars=NameParameters)
   plt.figure(figsize=(16, 6))
   ax = sns.boxplot(x="variable", y="value", hue="label", data=data_plot, palette=colours, saturation=1)
-  ax.set_xticklabels([r"$r_{recruit}[M]$",r"$\rho_{min}[M]$",r"$\rho_{sat}[M]$",r"$r_{recruit}[D]$",r"$\rho_{min}[D]$",r"$\rho_{sat}[D]$",r"$P_{leave}[D]$",r"$\delta_C$",r"$\kappa_{T}$",r"$\delta_{DM}$"],fontsize=14)
+  ax.set_xticklabels([r"$r_{recruit}[MP]$",r"$\rho_{min}[MP]$",r"$\rho_{sat}[MP]$",r"$r_{recruit}[DC]$",r"$\rho_{min}[DC]$",r"$\rho_{sat}[DC]$",r"$r_{leave}$",r"$\delta_C$",r"$\kappa_{T}$",r"$\delta_{DM}$"],fontsize=14)
   ax.set(xlabel=None)
   ax.set_ylabel('Normalized parameters',fontsize=14)
   ax.legend(loc='upper right', bbox_to_anchor=(1.1, 1),fontsize=14)
   if (FigName): plt.savefig(FigName)
   else: plt.show()
 
-def Plot_PatientAnalysis(PlotPatientPoint = False, FigName=None):
+def plot_PatientAnalysis(PlotPatientPoint = False, FigName=None):
     df_PatientVariation = PatientsAnalysis()
-    fig = plt.figure(figsize=(12, 8))
-    G = gridspec.GridSpec(4, 6, width_ratios=[1,1,1,1,1,0.3])
+    fig = plt.figure(figsize=(12, 10))
+    G = gridspec.GridSpec(4, 2, width_ratios=[1,0.05],height_ratios=[0.01,1,0.01,1])
+    # Patient classified as marginal control
+    ax_title_1 = plt.subplot(G[0, 0])
+    ax_title_1.axis('off')
+    ax_title_1.set_title('Ten virtual patients initially classified as marginal control (MC)')
+    G0 = gridspec.GridSpecFromSubplotSpec(2, 5, subplot_spec=G[1,0])
+    # Patient classified as no control
+    ax_title_2 = plt.subplot(G[2, 0])
+    ax_title_2.axis('off')
+    ax_title_2.set_title('Ten virtual patients initially classified as no control (NC)')
+    G1 = gridspec.GridSpecFromSubplotSpec(2, 5, subplot_spec=G[3,0])
     # Axes color bar
-    axcb = plt.subplot(G[:, 5])
+    axcb = plt.subplot(G[:, 1])
     NewColours = {'NC': colours['NC'], 'Mixed/MC': colours['Mixed/MC'], 'SC': colours['SC']}
     customPalette = sns.color_palette(NewColours.values(),len(NewColours))
 
     for ind, patientID in enumerate(MC_Patients_samples):
-        ax = plt.subplot(G[int(ind/5), ind%5])
+        ax = plt.subplot(G0[int(ind/5), ind%5])
         df = df_PatientVariation.loc[df_PatientVariation['patientID'] == patientID]
+        df["label"] = pd.to_numeric(df["label"])
         df_wide = df.pivot_table( index='DC_recruit', columns='Mac_recruit', values='label')
         g1 = sns.heatmap(df_wide,cmap=customPalette,vmin=0,vmax=2,cbar=False,ax=ax,linewidths=.5)#,annot=True)
         g1.invert_yaxis()
         g1.set(xlabel = (''), ylabel = (''), xticks = ([]), yticks = ([]))
-        g1.set_title('Patient '+str(patientID),color="#4783af")
+        g1.set_title('Patient '+str(patientID))#,color="#4783af")
         if ( ind%5 == 0): g1.set_yticks(ticks=[0.01,1.5,2.99],labels=['0','0.5','1'])
         if (PlotPatientPoint):
             MAC_rec_value = df_input.loc[df_input['sample'] == patientID]['macrophage_max_recruitment_rate'].values[0]
             DC_rec_value = df_input.loc[df_input['sample'] == patientID]['DC_max_recruitment_rate'].values[0]
             ax.scatter(3*MAC_rec_value/8e-9,3*DC_rec_value/4e-9, color=colours['MC'], edgecolors='k')
 
     for ind, patientID  in enumerate(NC_Patients_samples):
-        ax = plt.subplot(G[2+int((ind)/5), ind%5])
+        ax = plt.subplot(G1[int((ind)/5), ind%5])
         df = df_PatientVariation.loc[df_PatientVariation['patientID'] == patientID]
+        df["label"] = pd.to_numeric(df["label"])
         df_wide = df.pivot_table( index='DC_recruit', columns='Mac_recruit', values='label')
         if (ind == 9): g1 = sns.heatmap(df_wide,cmap=customPalette,vmin=0,vmax=2,ax=ax,linewidths=.5,cbar_ax=axcb)#,annot=True)
         else: g1 = sns.heatmap(df_wide,cmap=customPalette,vmin=0,vmax=2,cbar=False,ax=ax,linewidths=.5)#,annot=True)
         g1.invert_yaxis()
         g1.set(xlabel = (''), ylabel = (''), xticks = ([]), yticks = ([]))
-        g1.set_title('Patient '+str(patientID),color="#ff7f7f")
+        g1.set_title('Patient '+str(patientID))#,color="#ff7f7f")
         if ( int(ind/5) == 1): g1.set_xticks(ticks=[0.02,1.5,2.98],labels=['0','0.5','1'])
         if ( ind%5 == 0): g1.set_yticks(ticks=[0.01,1.5,2.99],labels=['0','0.5','1'])
         if (PlotPatientPoint):
             MAC_rec_value = df_input.loc[df_input['sample'] == patientID]['macrophage_max_recruitment_rate'].values[0]
             DC_rec_value = df_input.loc[df_input['sample'] == patientID]['DC_max_recruitment_rate'].values[0]
             ax.scatter(3*MAC_rec_value/8e-9,3*DC_rec_value/4e-9, color=colours['NC'], edgecolors='k')
 
-    fig.supxlabel(r'$r_{recruit}[M]$')
-    fig.supylabel(r'$r_{recruit}[D]$')
+    fig.supxlabel(r'$r_{recruit}[MP]$')
+    fig.supylabel(r'$r_{recruit}[DC]$')
     # modify colorbar:
     colorbar = ax.collections[0].colorbar
     r = colorbar.vmax - colorbar.vmin
@@ -93,11 +99,22 @@ def Plot_PatientAnalysis(PlotPatientPoint = False, FigName=None):
     if (FigName): plt.savefig(FigName)
     else: plt.show()
 
-def Plot_UMAP(Sim_df_lastFrame, FigName=None):
-    fig = plt.figure(figsize=(8, 7))
-    df_UMAP, df = Loading_UMAP()
+def plot_UMAP(Sim_df_lastFrame, FigName=None, NumSamples=None):
+    fig, ax = plt.subplots(nrows=1, ncols=1,figsize=(8, 7))
+    df_UMAP = Loading_UMAP_traj()
+    print(df_UMAP)
     df_UMAP['label'] = Sim_df_lastFrame['label']
-    ax = sns.scatterplot(data=df_UMAP,x='UMAP1',y='UMAP2',hue='label',palette=[colours['NC'],colours['MC'],colours['SC']], hue_order=['NC','MC','SC'])
+    if NumSamples:
+        cluster1 = df_UMAP.loc[ (df_UMAP['UMAP1'] < 10) & (df_UMAP['label'] == 'NC') ].sample(n=NumSamples, random_state=42).index.to_numpy()
+        cluster2 = df_UMAP.loc[ (df_UMAP['UMAP1'] >= 10) & (df_UMAP['label'] == 'NC')].sample(n=NumSamples, random_state=42).index.to_numpy()
+        df_UMAP['label'] = ['NC1' if temp_s in cluster1 else  temp_label for temp_s, temp_label  in zip(df_UMAP.index.to_numpy(),df_UMAP['label'])]
+        df_UMAP['label'] = ['NC2' if temp_s in cluster2 else  temp_label for temp_s, temp_label  in zip(df_UMAP.index.to_numpy(),df_UMAP['label'])]
+        print(df_UMAP.loc[df_UMAP['label'] == 'NC1'] )
+        print(df_UMAP.loc[df_UMAP['label'] == 'NC2'] )
+        hue_order = ['NC','MC','SC','NC1','NC2']
+        sns.scatterplot(data=df_UMAP.sort_values('label', key=np.vectorize(hue_order.index)),x='UMAP1',y='UMAP2',hue='label',palette=[[0.79216, 0.50588, 0.31765, 0.1],[1.0, 0.88235, 0.70196, 0.1],[0.18823, 0.48627, 0.63137, 0.1],[0, 0, 0, 1.0],[1.0, 0.0, 0.0, 1.0]], hue_order=hue_order, ax=ax)
+    else:
+        sns.scatterplot(data=df_UMAP,x='UMAP1',y='UMAP2',hue='label',palette=[colours['NC'],colours['MC'],colours['SC']], hue_order=['NC','MC','SC'], ax=ax)
     ax.legend(markerscale=2.0, fontsize=14)
     ax.set_xlabel('UMAP 1', fontsize=14)
     ax.set_ylabel('UMAP 2', fontsize=14)
@@ -106,11 +123,11 @@ def Plot_UMAP(Sim_df_lastFrame, FigName=None):
     if (FigName): plt.savefig(FigName)
     else: plt.show()
 
-def Plot_UMAP_MeanPatients(Pat_df_7classes, FigName=None):
-    fig = plt.figure(figsize=(8, 7))
-    df,df_UMAP = Loading_UMAP()
+def plot_UMAP_MeanPatients(Pat_df_7classes, FigName=None):
+    fig, ax = plt.subplots(nrows=1, ncols=1,figsize=(8, 7))
+    df_UMAP = Loading_UMAP_avg()
     df_UMAP['label'] = Pat_df_7classes['label']
-    ax = sns.scatterplot(data=df_UMAP,x='UMAP1',y='UMAP2',hue='label',palette=[colours['NC'],colours['MC'],colours['SC'],colours['NC+MC'],colours['NC+SC'],colours['MC+SC'],colours['NC+MC+SC']], hue_order=['NC','MC','SC','NC+MC','NC+SC','MC+SC','NC+MC+SC'])
+    sns.scatterplot(data=df_UMAP,x='UMAP1',y='UMAP2',hue='label',palette=[colours['NC'],colours['MC'],colours['SC'],colours['NC+MC'],colours['NC+SC'],colours['MC+SC'],colours['NC+MC+SC']], hue_order=['NC','MC','SC','NC+MC','NC+SC','MC+SC','NC+MC+SC'], ax=ax)
     ax.legend(markerscale=2.0, fontsize=14)
     ax.set_xlabel('UMAP 1', fontsize=14)
     ax.set_ylabel('UMAP 2', fontsize=14)
@@ -122,12 +139,19 @@ def Plot_UMAP_MeanPatients(Pat_df_7classes, FigName=None):
 if __name__ == '__main__':
     df_input, df_output = Loading_dataset()
     Sim_df_lastFrame = Classifier_Simulations(df_output)
+    print("DataFrame info (UMAP trajectories):")
+    print(Sim_df_lastFrame.info())
+    # plot_UMAP(Sim_df_lastFrame,FigName='Figure5_B.svg') # Plot UMAP of the trajectories
+    plot_UMAP(Sim_df_lastFrame,NumSamples=4,FigName='Figure5_B_analysis.png')
+    exit()
+    print("DataFrame info (UMAP trajectories): ",Sim_df_lastFrame.info())
     Pat_df_7classes, Pat_df_3classes = Classifier_Patients(Sim_df_lastFrame)
-    Plot_UMAP(Sim_df_lastFrame,FigName='Figure4_B.svg') # Plot UMAP of the trajectories
-    Plot_UMAP_MeanPatients(Pat_df_7classes,FigName='Figure4_D.svg') # Plot UMAP of the trajectory averages
-    Plot_Patients(Sim_df_lastFrame,Pat_df_7classes,Pat_df_3classes, FigName='Figure4_AC.svg') # Patient statistics
+    plot_UMAP_MeanPatients(Pat_df_7classes,FigName='Figure5_D.svg') # Plot UMAP of the trajectory averages
+    print("DataFrame info (UMAP trajectory averages): ",Pat_df_7classes.info())
+    exit()
+    plot_Patients(Sim_df_lastFrame,Pat_df_7classes,Pat_df_3classes, FigName='Figure5_AC.svg') # Patient statistics
     # Plot the quartiles of patient features
     df_input_min_max_scaled = Normalize_Parameters(df_input, Pat_df_3classes)
-    Plot_Parameters(df_input_min_max_scaled, FigName='Figure6_A.svg')
+    plot_Parameters(df_input_min_max_scaled, FigName='Figure4_A.svg')
     # Plot the analysis of 20 patients
-    Plot_PatientAnalysis(FigName='Figure6_B.svg')
+    plot_PatientAnalysis(FigName='Figure4_B.svg')