python/salus/cmplatform/fig1.py

   1 #!/usr/bin/env python3
   2
   3 import sys
   4 import os
   5 from typing import NamedTuple
   6
   7 PlatformTuple = ('Illumina', 'Salus')
   8 SamplesDict = {
   9     'mbrain': {
  10         'sid' : 'mbrain',
  11         'sub' : 'Mouse Brain Sptial',
  12         'type': 'visium',
  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         'prefix'   : '/share/result/spatial/data/BoAo_sp',
  24         'suffixOut': dict.fromkeys(PlatformTuple,"outs"),
  25         'suffixMtx': 'filtered_feature_bc_matrix',
  26         'platforms': {PlatformTuple[0]:'illumina', PlatformTuple[1]: 'salus'},
  27         'pattern': ('prefix', 'platformV', 'sid', 'suffixOutV', 'suffixMtx')
  28     },
  29     'human': {
  30         'sid' : 'human',
  31         'sub' : 'Human Single Cell',
  32         'type': 'mobivision',
  33         'prefix'   : '/share/result/spatial/data/MoZhuo_sc/FX20230913',
  34         'suffixOut': {PlatformTuple[0]: 'out/R22045213-220914-LYY-S11-R03-220914-LYY-S11-R03_combined_outs',
  35                       PlatformTuple[1]: 'out_subset/20221124-LYY-S09-R03_AGGCAGAA_fastq_outs'},
  36         'suffixMtx': 'filtered_cell_gene_matrix',
  37         'platforms': {PlatformTuple[0]:'illumina', PlatformTuple[1]: 'sailu'},
  38         'pattern': ('prefix', 'platformV', 'suffixOutV', 'suffixMtx')
  39     }
  40 }
  41
  42 def checkModules() -> None:
  43     import importlib.metadata
  44     from packaging import version
  45     pkgname = "squidpy"
  46     min_ver = "1.2.3"
  47     got_ver = importlib.metadata.version(pkgname)
  48     if version.parse(got_ver) < version.parse(min_ver):
  49         raise importlib.VersionConflict(f"{pkgname}>={min_ver} is needed, but found {pkgname}=={got_ver}")
  50
  51 if __name__ == "__main__":
  52     if len(sys.argv) > 1:
  53         thisID = sys.argv[1]
  54         if thisID not in SamplesDict:
  55             print(f"[x]sid can only be {SamplesDict.keys()}", file=sys.stderr)
  56             exit(1)
  57     else:
  58         thisID = 'mbrain'
  59     print(sys.argv, file=sys.stderr)
  60     print(f"[i]{thisID}")
  61     sys.stdout.flush()
  62     checkModules()
  63
  64 import matplotlib; matplotlib.use("module://mplcairo.base")
  65 from matplotlib import pyplot as plt
  66 import mplcairo
  67
  68 plt.rcParams['figure.figsize'] = (6.0, 6.0) # set default size of plots
  69 font = {'family' : 'STIX Two Text',
  70         #'size'   : 22,
  71         'weight' : 'normal'}
  72 matplotlib.rc('font', **font)
  73
  74 import numpy as np
  75 import pandas as pd
  76 import fast_matrix_market
  77 import anndata as ad
  78 import scanpy as sc
  79 import squidpy as sq
  80 import seaborn as sns
  81 import scipy
  82
  83 def main() -> None:
  84
  85     class scDatItem(NamedTuple):
  86         name: str
  87         bgRaw: tuple[int,int]
  88         bgFlt: tuple[int,int]
  89         annDat: ad.AnnData
  90
  91         def __repr__(self) -> str:
  92             return f'[sc:{self.name}, Raw_BC*Gene={self.bgRaw[0]}x{self.bgRaw[1]}, NonZero_BC*Gene={self.bgFlt[0]}x{self.bgFlt[1]} ({self.annDat.n_obs}x{self.annDat.n_vars})]'
  93
  94     scDat = []
  95     nfoDict = SamplesDict[thisID]
  96     print("[i]Start.", file=sys.stderr)
  97     for platform in PlatformTuple:
  98         nfoDict['platformK']  = platform
  99         nfoDict['platformV']  = nfoDict['platforms'][platform]
 100         nfoDict['suffixOutV'] = nfoDict['suffixOut'][platform]
 101         mtxPath = os.path.join( *[nfoDict[v] for v in nfoDict['pattern']] )
 102         print(f"[i]Reading {mtxPath}", file=sys.stderr)
 103         adata=sc.read_10x_mtx(mtxPath, var_names='gene_symbols', make_unique=True, gex_only=True)
 104         adata.var_names_make_unique()  # this is necessary if using `var_names='gene_symbols'` in `sc.read_10x_mtx`
 105         nnRaw = (adata.n_obs,adata.n_vars)
 106         adata.var['mt'] = adata.var_names.str.startswith('MT-') | adata.var_names.str.startswith('mt-')
 107         sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=True, inplace=True)
 108         sc.pp.filter_cells(adata, min_genes=1)
 109         sc.pp.filter_genes(adata, min_cells=1)
 110         nnFlt = (adata.n_obs,adata.n_vars)
 111         scDat.append(scDatItem(platform,nnRaw,nnFlt,adata))
 112
 113     print("\n".join(map(str,scDat)))
 114
 115     with pd.option_context("mode.copy_on_write", True):
 116         obsmbi = scDat[0].annDat.obs[['n_genes_by_counts', 'total_counts']].copy(deep=False)
 117         obsmbs = scDat[1].annDat.obs[['n_genes_by_counts', 'total_counts']].copy(deep=False)
 118         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()
 119         p2df = obsmbi.join(obsmbs,lsuffix='_'+scDat[0].name,rsuffix='_'+scDat[1].name,how='inner').replace([np.inf, -np.inf, 0], np.nan).dropna()
 120         p3tuple = (frozenset(scDat[0].annDat.var_names), frozenset(scDat[1].annDat.var_names))
 121
 122     print("[i]Begin fig A. 1D", file=sys.stderr)
 123     custom_params = {"axes.spines.right": False, "axes.spines.top": False}
 124     sns.set_theme(style="ticks", rc=custom_params, font="STIX Two Text")
 125     figA=sns.JointGrid(data=p1df, x="total_counts", y="n_genes_by_counts", hue='Platform', dropna=True)
 126     #figA.plot(sns.scatterplot, sns.histplot, alpha=.7, edgecolor=".2", linewidth=.5)
 127     figA.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 128     figA.plot_marginals(sns.histplot, kde=True, alpha=.618)
 129     figA.figure.suptitle(f"Gene to UMI plot - {nfoDict['sub']}")
 130     figA.set_axis_labels(xlabel='UMIs per Barcode', ylabel='Genes per Barcode')
 131     figA.savefig(f"1D_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'Gene to UMI plot', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 132
 133     print("[i]Begin fig B. 1E", file=sys.stderr)
 134     figB=sns.JointGrid(data=p2df, x="total_counts_Illumina", y="total_counts_Salus", dropna=True)
 135     figB.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 136     figB.plot_marginals(sns.histplot, kde=True, alpha=.618)
 137     figB.figure.suptitle(f"UMI per Barcode Counts Comparing - {nfoDict['sub']}")
 138     figB.set_axis_labels(xlabel='UMI Counts from Illumina', ylabel='UMI Counts from Salus')
 139     figB.savefig(f"1E_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'UMI per Barcode Counts Comparing', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 140
 141     print("[i]Begin fig . 1F", file=sys.stderr)
 142     from matplotlib_venn import venn2
 143     plt.figure(figsize=(4,4))
 144     plt.title("Sample Venn diagram")
 145     p3intersection = p3tuple[0] & p3tuple[1]
 146     p3veen = (p3tuple[0]-p3intersection, p3tuple[1]-p3intersection, p3intersection)
 147     GenesA = scDat[0].annDat.var.loc[p3veen[0]-p3veen[2]]
 148     GenesB = scDat[1].annDat.var.loc[p3veen[1]-p3veen[2]]
 149     GenesC = scDat[0].annDat.var.loc[p3veen[2]]
 150     p3vd=venn2(subsets=tuple(map(len,p3veen)), set_labels=(scDat[0].name, scDat[1].name))
 151     plt.savefig(f"1F_Genes_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'Veen of Genes', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 152     GenesA.to_csv(f"1F_Genes_{nfoDict['sid']}_{scDat[0].name}_only.csv",encoding='utf-8')
 153     GenesB.to_csv(f"1F_Genes_{nfoDict['sid']}_{scDat[1].name}_only.csv",encoding='utf-8')
 154     GenesC.to_csv(f"1F_Genes_{nfoDict['sid']}_intersection.csv.zst",encoding='utf-8',compression={'method': 'zstd', 'level': 9, 'write_checksum': True})
 155
 156     print("[i]Begin fig C. 2A", file=sys.stderr)
 157     # https://www.kaggle.com/code/lizabogdan/top-correlated-genes?scriptVersionId=109838203&cellId=21
 158     p4xdf = scDat[0].annDat.to_df()
 159     p4ydf = scDat[1].annDat.to_df()
 160     p4corraw = p4xdf.corrwith(p4ydf,axis=1)
 161     p4corr = p4corraw.dropna()
 162     plt.figure(figsize=(6,4))
 163     plt.title("Pearson correlation")
 164     figC=sns.histplot(p4corr,stat='percent',binwidth=0.01)
 165     plt.savefig(f"2A_Correlation_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'Pearson correlation', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 166
 167     print("[i]Begin fig D. 2B", file=sys.stderr)
 168     var_names = scDat[0].annDat.var_names.intersection(scDat[1].annDat.var_names)
 169     xadata = scDat[0].annDat[:, var_names]
 170     yadata = scDat[1].annDat[:, var_names]
 171     xdf=getOBSMdf(xadata)
 172     ydf=getOBSMdf(yadata)
 173     #p4df = xdf.assign(Platform=scDat[0].name).join(ydf.assign(Platform=scDat[1].name),lsuffix='_'+scDat[0].name,rsuffix='_'+scDat[1].name,how='inner')
 174     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()
 175     figD=sns.JointGrid(data=p4df, x="P1", y="P2", hue='Platform', dropna=True)
 176     figD.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 177     figD.plot_marginals(sns.histplot, kde=True, alpha=.618)
 178     figD.figure.suptitle(f"PCA - {nfoDict['sub']}")
 179     figD.set_axis_labels(xlabel='PC1', ylabel='PC2')
 180     figD.savefig(f"2B_PCA_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'PCA', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 181     print("[i]Begin fig E. 2C", file=sys.stderr)
 182     xdf=getOBSMdf(xadata,'X_umap')
 183     ydf=getOBSMdf(yadata,'X_umap')
 184     p5df = 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()
 185     figE=sns.JointGrid(data=p5df, x="P1", y="P2", hue='Platform', dropna=True)
 186     figE.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 187     figE.plot_marginals(sns.histplot, kde=True, alpha=.618)
 188     figE.figure.suptitle(f"UMAP - {nfoDict['sub']}")
 189     figE.set_axis_labels(xlabel='UMAP1', ylabel='UMAP2')
 190     figE.savefig(f"2C_UMAP_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'UMAP', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 191     print("[i]Begin fig E. 2Cn", file=sys.stderr)
 192     xdf=getOBSMdf(xadata,'X_draw_graph_fa')
 193     ydf=getOBSMdf(yadata,'X_draw_graph_fa')
 194     p5df = 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()
 195     figE=sns.JointGrid(data=p5df, x="P1", y="P2", hue='Platform', dropna=True)
 196     figE.plot_joint(sns.scatterplot, s=12.7, alpha=.6)
 197     figE.plot_marginals(sns.histplot, kde=True, alpha=.618)
 198     figE.figure.suptitle(f"ForceAtlas2 - {nfoDict['sub']}")
 199     figE.set_axis_labels(xlabel='FA1', ylabel='FA2')
 200     figE.savefig(f"2C_ForceAtlas2_{nfoDict['sid']}.pdf", transparent=True, dpi=300, metadata={'Title': 'ForceAtlas2', 'Subject': f"{nfoDict['sub']} Data", 'Author': 'HU Xuesong'})
 201
 202
 203 def getOBSMdf(anndata, obsmkey='X_pca') -> pd.DataFrame:
 204     if not obsmkey in anndata.obsm:
 205         if obsmkey=='X_pca':
 206             sc.tl.pca(anndata,zero_center=True)
 207         elif obsmkey=='X_umap':
 208             if not 'neighbors' in anndata.uns:
 209                 if not 'X_pca' in anndata.obsm:
 210                     sc.pp.pca(anndata,zero_center=True)
 211                 sc.pp.neighbors(anndata)
 212             sc.tl.umap(anndata)
 213         elif obsmkey=='X_draw_graph_fa':
 214             if not 'neighbors' in anndata.uns:
 215                 if not 'X_pca' in anndata.obsm:
 216                     sc.pp.pca(anndata,zero_center=True)
 217                 sc.pp.neighbors(anndata)
 218             sc.tl.draw_graph(anndata)
 219     data=anndata.obsm[obsmkey][0:,0:2]
 220     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])])
 221     return df
 222
 223 if __name__ == "__main__":
 224     main()  # time (./fig1.py human; ./fig1.py mbrain ; ./fig1.py mkidney ) | tee plot.log
 225
 226 '''
 227 x1 = np.random.randn(1000)
 228 y1 = np.random.randn(1000)
 229 x2 = np.random.randn(1000) * 5
 230 y2 = np.random.randn(1000)
 231 fig, ax = plt.subplots()
 232 # The figure and axes background must be made transparent.
 233 fig.patch.set(alpha=0)
 234 ax.patch.set(alpha=0)
 235 pc1 = ax.scatter(x1, y1, c='b', edgecolors='none')
 236 pc2 = ax.scatter(x2, y2, c='r', edgecolors='none')
 237 mplcairo.operator_t.ADD.patch_artist(pc2)  # Use additive blending.
 238 plt.show()
 239
 240 1、N和Q<5比率大于4%
 241 2、Q平均值小于20
 242 3、Q<20和purity<0.6的比率大于18%
 243
 244 import patchworklib as pw
 245 #from blend_modes import addition
 246
 247 ToDo:
 248  * Try layers of annData.
 249 '''