探索人类肠道免疫细胞图谱

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树。

生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。

https://www.gutcellatlas.org/

下载h5ad格式的文件，我们一看应该用scanpy来分析啊，所以我们：

import scanpy as sc
import pandas as pd
import numpy as np  
import matplotlib.mlab as mlab  
import matplotlib.pyplot as plt 
import seaborn as sns

adata = sc.read_h5ad("C:\\Users\\Administrator\\Desktop\\Colon_cell_atlas.h5ad")
adata

看看数据格式，大部分的分析已经为我们做好了：

AnnData object with n_obs × n_vars = 41650 × 18927 
    obs: 'donor', 'region', 'n_genes', 'percent_mito', 'n_counts', 'cell_type', 'BCR_ChainCombination', 'BCR_SEQUENCE_ID', 'IGH_V_CALL_GENOTYPED', 'IGH_D_CALL', 'IGH_J_CALL', 'BCR_ISOTYPE', 'BCR_CLONE', 'BCR_PANDONOR_CLONE', 'IGL_C_Gene', 'IGL_VDJ_Gene', 'IGH_READS', 'IGH_UMIS', 'IGH_MU_FREQ', 'IGK_READS', 'IGK_UMIS', 'IGL_READS', 'IGL_UMIS', 'TCR_v_gene', 'TCR_d_gene', 'TCR_j_gene', 'TCR_c_gene', 'TCR_cdr3', 'TCR_cdr3_nt', 'TCR_reads', 'TCR_umis', 'TCR_Clone', 'percent_mito1'
    var: 'gene_ids-Pan_T7935487', 'feature_types-Pan_T7935487', 'gene_ids-Pan_T7935488', 'feature_types-Pan_T7935488', 'gene_ids-Pan_T7935489', 'feature_types-Pan_T7935489', 'gene_ids-Pan_T7935494', 'feature_types-Pan_T7935494', 'gene_ids-4861STDY7208411', 'feature_types-4861STDY7208411', 'gene_ids-4861STDY7208413', 'feature_types-4861STDY7208413', 'gene_ids-4861STDY7208410', 'feature_types-4861STDY7208410', 'gene_ids-4861STDY7208412', 'feature_types-4861STDY7208412', 'gene_ids-4861STDY7208414', 'feature_types-4861STDY7208414', 'gene_ids-4861STDY7208415', 'feature_types-4861STDY7208415', 'gene_ids-4861STDY7208416', 'feature_types-4861STDY7208416', 'gene_ids-4861STDY7208417', 'feature_types-4861STDY7208417', 'gene_ids-4861STDY7135911', 'feature_types-4861STDY7135911', 'gene_ids-4861STDY7135912', 'feature_types-4861STDY7135912', 'gene_ids-4861STDY7135913', 'feature_types-4861STDY7135913', 'gene_ids-4861STDY7135914', 'feature_types-4861STDY7135914', 'gene_ids-4861STDY7135915', 'feature_types-4861STDY7135915', 'gene_ids-4861STDY7135916', 'feature_types-4861STDY7135916', 'gene_ids-4861STDY7135917', 'feature_types-4861STDY7135917', 'gene_ids-4861STDY7135918', 'feature_types-4861STDY7135918', 'gene_ids-Human_colon_16S7255675', 'feature_types-Human_colon_16S7255675', 'gene_ids-Human_colon_16S7255676', 'feature_types-Human_colon_16S7255676', 'gene_ids-Human_colon_16S7255677', 'feature_types-Human_colon_16S7255677', 'gene_ids-Human_colon_16S7255678', 'feature_types-Human_colon_16S7255678', 'gene_ids-Human_colon_16S7255679', 'feature_types-Human_colon_16S7255679', 'gene_ids-Human_colon_16S7255680', 'feature_types-Human_colon_16S7255680', 'gene_ids-Human_colon_16S7255681', 'feature_types-Human_colon_16S7255681', 'gene_ids-Human_colon_16S7255682', 'feature_types-Human_colon_16S7255682', 'gene_ids-Human_colon_16S8001863', 'feature_types-Human_colon_16S8001863', 'gene_ids-Human_colon_16S8001867', 'feature_types-Human_colon_16S8001867', 'gene_ids-Human_colon_16S8001871', 'feature_types-Human_colon_16S8001871', 'gene_ids-Human_colon_16S8000484', 'feature_types-Human_colon_16S8000484', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'ISOTYPE_colors', 'TCR_Clone_colors', 'TCR_cdr3_colors', 'TCR_chain_colors', 'annotation3_colors', 'annotation4_colors', 'donor_colors', 'leiden', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'region_colors', 'cell_type_colors'
    obsm: 'X_pca', 'X_umap'

所以我们可以做一些可视化的工作：

sc.pl.highest_expr_genes(adata, n_top=20,palette="dark" )

adata.obs['percent_mito1'] = adata.obs['percent_mito']*10000
adata.obs.melt(id_vars=['n_counts','percent_mito'], var_name='myVariable', value_name='myValue')
sns.set(style="ticks")
g = sns.FacetGrid(adata.obs.melt(id_vars=['n_counts','percent_mito'], var_name='myVariable', value_name='myValue'),col = 'myVariable')
g.map(plt.scatter,'n_counts',"myValue", alpha=.7)
g.add_legend();

sc.pl.umap(adata, color=['cell_type'])

图例画到图上看看：

sc.pl.umap(adata, color=['cell_type'],legend_loc='on data',legend_fontsize=7)

看看各个细胞类型的比例：

se = pd.Series(adata.obs['cell_type'])
countDict = dict(se.value_counts())
proportitionDict = dict(se.value_counts(normalize=True))
plt.bar(proportitionDict.keys(),proportitionDict.values())
plt.xticks(rotation=90)

plt.pie(proportitionDict.values(),labels =proportitionDict.keys(), autopct='%1.1f%%', shadow=False,startangle=150)

看看每个群的差异基因,其实作者已经计算好了，h5ad里面：

sc.get.rank_genes_groups_df(adata, group="Th1")
Out[*]: 
       scores  names  logfoldchanges          pvals      pvals_adj
0   74.605721  ANXA1        2.031630   0.000000e+00   0.000000e+00
1   71.477692  KLRB1        2.275039   0.000000e+00   0.000000e+00
2   65.698090   IL7R        1.729141   0.000000e+00   0.000000e+00
3   60.053253   IL32        1.783656   0.000000e+00   0.000000e+00
4   57.090515    ID2        1.872997   0.000000e+00   0.000000e+00
..        ...    ...             ...            ...            ...
95  31.323795    FOS        0.580307  1.350933e-197  2.635990e-195
96  31.206936  RPL30        0.353758  1.141130e-188  1.945780e-186
97  31.120506  SYTL3        1.492344  1.389941e-196  2.657314e-194
98  30.957434  RPL34        0.288940  5.245242e-186  8.785548e-184
99  30.864616  RPS14        0.303797  4.713286e-185  7.690377e-183

小提琴图：

sc.pl.violin(adata, ['ANXA1', 'KLRB1', 'PPBP'], groupby='cell_type',rotation=90)

点图：

ax = sc.pl.dotplot(adata, ['ANXA1', 'KLRB1', 'PPBP'], groupby='cell_type')

堆积小提琴图：

ax = sc.pl.stacked_violin(adata, ['ANXA1', 'KLRB1', 'PPBP'], groupby='cell_type', rotation=90)

当然，scanpy的分析都是可以做的，关键是找到有意思的生物学解释。其实，这个数据中大量的免疫组库的信息，我们并没有去利用。

adata.obs.columns
Out[59]: 
Index(['donor', 'region', 'n_genes', 'percent_mito', 'n_counts', 'cell_type',
       'BCR_ChainCombination', 'BCR_SEQUENCE_ID', 'IGH_V_CALL_GENOTYPED',
       'IGH_D_CALL', 'IGH_J_CALL', 'BCR_ISOTYPE', 'BCR_CLONE',
       'BCR_PANDONOR_CLONE', 'IGL_C_Gene', 'IGL_VDJ_Gene', 'IGH_READS',
       'IGH_UMIS', 'IGH_MU_FREQ', 'IGK_READS', 'IGK_UMIS', 'IGL_READS',
       'IGL_UMIS', 'TCR_v_gene', 'TCR_d_gene', 'TCR_j_gene', 'TCR_c_gene',
       'TCR_cdr3', 'TCR_cdr3_nt', 'TCR_reads', 'TCR_umis', 'TCR_Clone',
       'percent_mito1'],
      dtype='object')

adata.obs['IGL_C_Gene']
Out[60]: 
index
AAACGGGAGGTGCAAC-1-Pan_T7935487              IGLC1
AAAGATGTCCTCAACC-1-Pan_T7935487               None
AAAGTAGTCTTACCGC-1-Pan_T7935487              IGLC2
AACCATGCAGATTGCT-1-Pan_T7935487              IGLC2
AACCATGTCCTGCCAT-1-Pan_T7935487               None

TTTGTCACAAGTTAAG-1-Human_colon_16S8000484    IGLC1
TTTGTCAGTACCGAGA-1-Human_colon_16S8000484      nan
TTTGTCATCAACACCA-1-Human_colon_16S8000484      nan
TTTGTCATCCCAACGG-1-Human_colon_16S8000484      nan
TTTGTCATCGGTGTCG-1-Human_colon_16S8000484      nan
Name: IGL_C_Gene, Length: 41650, dtype: category
Categories (7, object): [IGLC1, IGLC2, IGLC3, IGLC7, Multi, None, nan]

adata.obs['TCR_cdr3_nt']
Out[61]: 
index
AAACGGGAGGTGCAAC-1-Pan_T7935487                                                        nan
AAAGATGTCCTCAACC-1-Pan_T7935487                                                        nan
AAAGTAGTCTTACCGC-1-Pan_T7935487                                                        nan
AACCATGCAGATTGCT-1-Pan_T7935487                                                        nan
AACCATGTCCTGCCAT-1-Pan_T7935487                                                        nan

TTTGTCACAAGTTAAG-1-Human_colon_16S8000484                                              nan
TTTGTCAGTACCGAGA-1-Human_colon_16S8000484                                              nan
TTTGTCATCAACACCA-1-Human_colon_16S8000484    TGTGCCAGCAGTTACTTTGGGGGCGGCAATCAGCCCCAGCATTTT
TTTGTCATCCCAACGG-1-Human_colon_16S8000484       TGTGCTGTGAGGCCCCAGGATAGCAACTATCAGTTAATCTGG
TTTGTCATCGGTGTCG-1-Human_colon_16S8000484                                              nan
Name: TCR_cdr3_nt, Length: 41650, dtype: category
Categories (2423, object): [TGCAGCATGAAACGCCCCGGACTCACAGATACGCAGTATTTT, TGCAGCGCCAGGACTAGCGGGGGACACCAGTTCTTC, TGCAGCGCCCGGGGGGGGGACACCGGGGAGCTGTTTTTT, TGCAGCGCCGTCAGGGATGACGAAAAACTGTTTTTT, ..., TGTGTTGTGAGGGGTGGAGGAAGCCAAGGAAATCTCATCTTT, TGTGTTGTGAGTGACAGGGTCTCAGGAACCTACAAATACATCTTT, TGTTTTTATAACCAGGGAGGAAAGCTTATCTTC, nan]

数据是有的了，如何挖掘就要看讲故事的能力了。可以想象一下，免疫组库有哪些分析点，把他们用起来。

可以参考某公司的分析报告：

  • Sequence Alignment
  • CDR Abundancy Analysis
  • CDR3 Length Distribution
  • Top V(D)J Analysis
  • V Gene Usage Analysis
  • Clustering Analysis
  • Diversity Analysis

https://cdn2.hubspot.net/hubfs/3478602/Sell%20Sheet%20Collateral%20Library/NGS/NGS%20BI%20Reports/BCR%20Immuno-Profiling/Immunoglobulin%20Repertoire%20Sequencing%20Analysis%20Report.html

赞 (0)

相关推荐