python/salus/cmplatform/nfig1.py

   1 #!/usr/bin/env python3
   2 import sys
   3 import os
   4 from typing import NamedTuple
   5
   6 PlatformTuple = ('Illumina', 'Salus')
   7 SamplesDict = {
   8     'mbrain': {
   9         'sid' : 'mbrain',
  10         'sub' : 'Mouse Brain Sptial',
  11         'type': 'visium',
  12         'fltPct'   : 99.5,
  13         'prefix'   : '/share/result/spatial/data/BoAo_sp',
  14         'suffixOut': dict.fromkeys(PlatformTuple,"outs"),
  15         'suffixMtx': 'filtered_feature_bc_matrix',
  16         'platforms': {PlatformTuple[0]:'illumina', PlatformTuple[1]: 'salus'},
  17         'pattern': ('prefix', 'platformV', 'sid', 'suffixOutV', 'suffixMtx')
  18     },
  19     'mkidney': {
  20         'sid' : 'mkidney',
  21         'sub' : 'Mouse Kindey Sptial',
  22         'type': 'visium',
  23         'fltPct'   : 99.5,
  24         'prefix'   : '/share/result/spatial/data/BoAo_sp',
  25         'suffixOut': dict.fromkeys(PlatformTuple,"outs"),
  26         'suffixMtx': 'filtered_feature_bc_matrix',
  27         'platforms': {PlatformTuple[0]:'illumina', PlatformTuple[1]: 'salus'},
  28         'pattern': ('prefix', 'platformV', 'sid', 'suffixOutV', 'suffixMtx')
  29     },
  30     'human': {
  31         'sid' : 'human',
  32         'sub' : 'Human Single Cell',
  33         'type': 'mobivision',
  34         'fltPct'   : 85,
  35         'prefix'   : '/share/result/spatial/data/MoZhuo_sc/FX20230913',
  36         'suffixOut': {PlatformTuple[0]: 'out/R22045213-220914-LYY-S11-R03-220914-LYY-S11-R03_combined_outs',
  37                       PlatformTuple[1]: 'out_subset/20221124-LYY-S09-R03_AGGCAGAA_fastq_outs'},
  38         'suffixMtx': 'filtered_cell_gene_matrix',
  39         'platforms': {PlatformTuple[0]:'illumina', PlatformTuple[1]: 'sailu'},
  40         'pattern': ('prefix', 'platformV', 'suffixOutV', 'suffixMtx')
  41     }
  42 }
  43
  44 thisID = 'mbrain'
  45 if __name__ == "__main__":
  46     if len(sys.argv) > 1:
  47         thisID = sys.argv[1]
  48         if thisID not in SamplesDict:
  49             print(f"[x]sid can only be {SamplesDict.keys()}", file=sys.stderr)
  50             exit(1)
  51     print(sys.argv, file=sys.stderr)
  52     print(f"[i]{thisID}")
  53     sys.stdout.flush()
  54 nfoDict = SamplesDict[thisID]
  55
  56 import matplotlib; matplotlib.use("module://mplcairo.base")
  57 from matplotlib import pyplot as plt
  58 import mplcairo
  59
  60 plt.rcParams['figure.figsize'] = (6.0, 6.0) # set default size of plots
  61 plt.rcParams['figure.dpi'] = 300
  62 plt.rcParams['savefig.bbox'] = 'tight'
  63 plt.rcParams["savefig.transparent"] = True
  64 font = {'family' : 'STIX Two Text',
  65         #'size'   : 22,
  66         'weight' : 'normal'}
  67 matplotlib.rc('font', **font)
  68
  69 import numpy as np
  70 import pandas as pd
  71 import fast_matrix_market
  72 import anndata as ad
  73 import leidenalg
  74 import scanpy as sc
  75 sc._settings.ScanpyConfig.n_jobs = -1
  76 #import squidpy as sq
  77 import seaborn as sns
  78 import scipy
  79 import pynndescent
  80
  81 import warnings
  82 warnings.filterwarnings('ignore')
  83 from copy import deepcopy
  84
  85 class scDatItem(NamedTuple):
  86     name: str
  87     rawDat: ad.AnnData
  88     annDat: ad.AnnData
  89     def __repr__(self) -> str:
  90         return f'[sc:{self.name}, BC*Gene: Raw={self.rawDat.shape}, Filtered={self.annDat.shape}]'
  91
  92 def main() -> None:
  93     scDat = []
  94     #nfoDict = SamplesDict[thisID]
  95     print("[i]Start.", file=sys.stderr)
  96     for platform in PlatformTuple:
  97         nfoDict['platformK']  = platform
  98         nfoDict['platformV']  = nfoDict['platforms'][platform]
  99         nfoDict['suffixOutV'] = nfoDict['suffixOut'][platform]
 100         mtxPath = os.path.join( *[nfoDict[v] for v in nfoDict['pattern']] )
 101         print(f"[i]Reading {mtxPath}", file=sys.stderr)
 102         adata=sc.read_10x_mtx(mtxPath, var_names='gene_symbols', make_unique=True, gex_only=True)
 103         adata.var_names_make_unique()
 104         adata.var['mt'] = adata.var_names.str.startswith('MT-') | adata.var_names.str.startswith('mt-')
 105         rdata = deepcopy(adata)
 106         sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=True, inplace=True)
 107         sc.pp.filter_cells(adata, min_genes=1)
 108         sc.pp.filter_genes(adata, min_cells=1)
 109         scDat.append(scDatItem(platform,rdata,adata))
 110         rdata.write_h5ad(f"{nfoDict['sid']}_{platform}.raw.h5ad",compression='lzf')
 111     print("\n".join(map(str,scDat)))
 112     with pd.option_context("mode.copy_on_write", True):
 113         obsmbi = scDat[0].annDat.obs[['n_genes_by_counts', 'total_counts']].copy(deep=False)
 114         obsmbs = scDat[1].annDat.obs[['n_genes_by_counts', 'total_counts']].copy(deep=False)
 115         p1df = pd.concat([obsmbi.assign(Platform=scDat[0].name), obsmbs.assign(Platform=scDat[1].name)], ignore_index=True).replace([np.inf, -np.inf, 0], np.nan).dropna()
 116         p2df = obsmbi.join(obsmbs,lsuffix='_'+scDat[0].name,rsuffix='_'+scDat[1].name,how='inner').replace([np.inf, -np.inf, 0], np.nan).dropna()
 117         p3tuple = (frozenset(scDat[0].annDat.var_names), frozenset(scDat[1].annDat.var_names))
 118
 119     metapdf={'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'}
 120     print("[i]Begin fig A. 1D", file=sys.stderr)
 121     custom_params = {"axes.spines.right": False, "axes.spines.top": False}
 122     sns.set_theme(style="ticks", rc=custom_params, font="STIX Two Text")
 123     figA=sns.JointGrid(data=p1df, x="total_counts", y="n_genes_by_counts", hue='Platform', dropna=True)
 124     #figA.plot(sns.scatterplot, sns.histplot, alpha=.7, edgecolor=".2", linewidth=.5)
 125     figA.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 126     figA.plot_marginals(sns.histplot, kde=False, alpha=.618)
 127     figA.figure.suptitle(f"Gene to UMI plot - {nfoDict['sub']}")
 128     figA.set_axis_labels(xlabel='UMIs per Barcode', ylabel='Genes per Barcode')
 129     figA.savefig(f"1D_{nfoDict['sid']}.pdf", metadata={**metapdf, 'Title': 'Gene to UMI plot'})
 130
 131     print("[i]Begin fig B. 1E", file=sys.stderr)
 132     figB=sns.JointGrid(data=p2df, x="total_counts_Illumina", y="total_counts_Salus", dropna=True)
 133     figB.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 134     figB.plot_marginals(sns.histplot, kde=True, alpha=.618)
 135     figB.figure.suptitle(f"UMI per Barcode Counts Comparing - {nfoDict['sub']}")
 136     figB.set_axis_labels(xlabel='UMI Counts from Illumina', ylabel='UMI Counts from Salus')
 137     figB.savefig(f"1E_{nfoDict['sid']}.pdf", metadata={**metapdf, 'Title': 'UMI per Barcode Counts Comparing'})
 138
 139     print("[i]Begin fig . 1G", file=sys.stderr)
 140     from matplotlib_venn import venn2
 141     plt.figure(figsize=(4,4))
 142     plt.title(f"Genes Venn diagram - {nfoDict['sub']}")
 143     p3intersection = p3tuple[0] & p3tuple[1]
 144     p3veen = (p3tuple[0]-p3intersection, p3tuple[1]-p3intersection, p3intersection)
 145     GenesA = scDat[0].annDat.var.loc[p3veen[0]-p3veen[2]]
 146     GenesB = scDat[1].annDat.var.loc[p3veen[1]-p3veen[2]]
 147     GenesC = scDat[0].annDat.var.loc[p3veen[2]]
 148     p3vd=venn2(subsets=tuple(map(len,p3veen)), set_labels=(scDat[0].name, scDat[1].name))
 149     plt.savefig(f"1G_Genes_{nfoDict['sid']}.pdf", metadata={**metapdf, 'Title': 'Veen of Genes'})
 150     GenesA.to_csv(f"1G_Genes_{nfoDict['sid']}_{scDat[0].name}_only.csv",encoding='utf-8')
 151     GenesB.to_csv(f"1G_Genes_{nfoDict['sid']}_{scDat[1].name}_only.csv",encoding='utf-8')
 152     GenesC.to_csv(f"1G_Genes_{nfoDict['sid']}_intersection.csv.zst",encoding='utf-8',compression={'method': 'zstd', 'level': 9, 'write_checksum': True})
 153
 154     print("[i]Begin fig C. 2A", file=sys.stderr)
 155     # https://www.kaggle.com/code/lizabogdan/top-correlated-genes?scriptVersionId=109838203&cellId=21
 156     p4xdf = scDat[0].annDat.to_df()
 157     p4ydf = scDat[1].annDat.to_df()
 158     p4corraw = p4xdf.corrwith(p4ydf,axis=1)
 159     p4corr = p4corraw.dropna()
 160     plt.figure(figsize=(6,4))
 161     plt.title(f"Pearson correlation - {nfoDict['sub']}")
 162     figC=sns.histplot(p4corr,stat='count',binwidth=0.01)
 163     plt.savefig(f"2A_Correlation_{nfoDict['sid']}.pdf", metadata={**metapdf, 'Title': 'Pearson correlation'})
 164     '''
 165     print("[i]Begin fig D. 2B", file=sys.stderr)
 166     var_names = scDat[0].annDat.var_names.intersection(scDat[1].annDat.var_names)
 167     xadata = scDat[0].annDat[:, var_names]
 168     yadata = scDat[1].annDat[:, var_names]
 169     xdf=getOBSMdf(xadata)
 170     ydf=getOBSMdf(yadata)
 171     #p4df = xdf.assign(Platform=scDat[0].name).join(ydf.assign(Platform=scDat[1].name),lsuffix='_'+scDat[0].name,rsuffix='_'+scDat[1].name,how='inner')
 172     p4df = pd.concat([xdf.assign(Platform=scDat[0].name), ydf.assign(Platform=scDat[1].name)], ignore_index=True).replace([np.inf, -np.inf, 0], np.nan).dropna()
 173     figD=sns.JointGrid(data=p4df, x="P1", y="P2", hue='Platform', dropna=True)
 174     figD.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 175     figD.plot_marginals(sns.histplot, kde=True, alpha=.618)
 176     figD.figure.suptitle(f"PCA - {nfoDict['sub']}")
 177     figD.set_axis_labels(xlabel='PC1', ylabel='PC2')
 178     figD.savefig(f"2B_rawPCA_{nfoDict['sid']}.pdf", metadata={**metapdf, 'Title': 'PCA'})
 179     '''
 180     import scvi
 181     for IDlist in ([0],[1],[0,1]):
 182         rawList = [scDat[i].rawDat for i in IDlist]
 183         dataIDs = [scDat[i].name for i in IDlist]
 184         if len(rawList) == 1:
 185             adata = rawList[0]
 186             dataID = dataIDs[0]
 187         elif len(rawList) == 2:
 188             adata=ad.concat(rawList, label='Platform', keys=PlatformTuple, index_unique='-')
 189             dataID = 'Both'
 190         print(f"[i]Begin Tab 1. 1F Dropout rates - {dataID}. With scvi {scvi.__version__}", file=sys.stderr)
 191         adata.var['mt'] = adata.var_names.str.startswith('MT-') | adata.var_names.str.startswith('mt-')
 192         sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=True, inplace=True)
 193         if dataID == 'Both':
 194             scvi.data.poisson_gene_selection(adata,n_top_genes=8000,n_samples=10000,batch_key='Platform')
 195         else:
 196             scvi.data.poisson_gene_selection(adata,n_top_genes=8000,n_samples=10000)
 197         doDropOutPlot(dataID,adata)
 198         adata = None
 199
 200 def doDropOutPlot(dataID,adata) -> None:
 201     adata.var['mean_'] = np.array(adata.X.mean(0))[0]
 202     GenesM = adata.var.sort_values(by='prob_zero_enrichment_rank', ascending=False)
 203     GenesM.to_csv(f"1F_GenesDropout_{nfoDict['sid']}_{dataID}_PlatformAsBatch.csv.zst",encoding='utf-8',compression={'method': 'zstd', 'level': 9, 'write_checksum': True})
 204     print(f"[i]Begin Fig 1. 1F GenesM3DropSelected (added) - {dataID}", file=sys.stderr)
 205     highly_variable_df = adata.var.query('highly_variable')
 206     # Set up the figure and axes
 207     fig, ax = plt.subplots(figsize=(10, 6))
 208     # Create the scatter plot for the main points with color bar
 209     scatter = sns.scatterplot(x='mean_', y='observed_fraction_zeros', hue='prob_zero_enrichment', data=adata.var, palette='viridis', legend='brief')
 210     # Create the line plot for expected_fraction_zeros
 211     sns.lineplot(x='mean_', y='expected_fraction_zeros', data=adata.var, color='r', label='Expected Fraction Zeros')
 212     # Highlight highly variable points
 213     sns.scatterplot(x='mean_', y='observed_fraction_zeros', data=highly_variable_df, color='pink', marker='.', s=5, alpha=0.5)
 214     box_coords = adata.var.query('highly_variable').agg({'mean_': ['min', 'max'], 'observed_fraction_zeros': ['min', 'max']})
 215     # Draw a rectangle to cover highly variable points
 216     rect = plt.Rectangle(box_coords.loc['min'],
 217                          box_coords['mean_'].diff()['max'], box_coords['observed_fraction_zeros'].diff()['max'],
 218                          fill=None, edgecolor='blue', linewidth=2, alpha=0.5)
 219     ax.add_patch(rect)
 220     # Annotate right-top and left-bottom points
 221     fmt = '.4f'
 222     for mean_val, obs_frac_val in zip(box_coords['mean_'], box_coords['observed_fraction_zeros']):
 223         label = f'({mean_val:{fmt}},{obs_frac_val:{fmt}})'
 224         # Add padding to avoid overlapping with the rectangle
 225         bbox_props = dict(boxstyle="round,pad=0.3", fc="white", ec="white", lw=1, alpha=0.62)
 226         ax.text(mean_val, obs_frac_val, label, bbox=bbox_props)
 227     # Set x-axis to log scale
 228     ax.set_xscale('log')
 229     # Set plot title
 230     ax.set_title(f'Mean vs Observed Fraction Zeros - {nfoDict["sub"]} {dataID}')
 231     # Create a color bar for Prob Zero Enrichment
 232     cbar = fig.colorbar(scatter.get_children()[0], ax=ax, orientation='vertical', pad=0.1)
 233     cbar.set_label('Prob Zero Enrichment')
 234     plt.savefig(f"1F_GenesM3DropSelected_{nfoDict['sid']}_{dataID}_PlatformAsBatch.pdf", metadata={'Title': 'scvi.data.poisson_gene_selection', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 235     plt.close('all')
 236     print(f"[i]Begin Fig 1. 1F GenesDropoutHist (added) - {dataID}", file=sys.stderr)
 237     plt.figure(figsize=(6,4))
 238     plt.title(f"Gene DropRatio Histogram - {nfoDict['sub']} {dataID}")
 239     histplot = sns.histplot(adata.var, x='observed_fraction_zeros', bins=30, kde=False, hue='highly_variable', multiple="dodge", shrink=.8)
 240     bars_heights = [p.get_height() for p in histplot.patches if p.get_facecolor()[:3] == sns.color_palette()[1]]
 241     plt.ylim(0, max(bars_heights)*1.1)  # Adjust the margin as needed
 242     plt.savefig(f"1F_GenesDropoutHist_{nfoDict['sid']}_{dataID}_PlatformAsBatch.pdf", metadata={'Title': 'Gene DropRatio Histogram', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 243     plt.close('all')
 244
 245 def getOBSMdf(anndata, obsmkey='X_pca') -> pd.DataFrame:
 246     if not obsmkey in anndata.obsm:
 247         if obsmkey=='X_pca':
 248             sc.tl.pca(anndata,zero_center=True)
 249     data=anndata.obsm[obsmkey][0:,0:2]
 250     df=pd.DataFrame(data=data[0:,0:], index=[anndata.obs_names[i] for i in range(data.shape[0])], columns=['P'+str(1+i) for i in range(data.shape[1])])
 251     return df
 252
 253 if __name__ == "__main__":
 254     main()  # time (./nfig1.py human; ./nfig1.py mbrain ; ./nfig1.py mkidney ) | tee nplot.log
 255
 256 # pip install -U --force-reinstall lightning
 257 """
 258 pip3 install git+https://github.com/matplotlib/mplcairo
 259 pip3 install matplotlib_venn
 260 """