python基本部分学完了之后，见专辑《python生信笔记 》，就需要开始进行大量的数据分析实战练习进行掌握了。今天来学习一篇顶刊杂志的单细胞数据处理，使用python！
note：这个数据主要来自最近开一个直播进行实战训练，详细情况以及进群方式见：python单细胞数据分析实战直播交流群
数据背景 今天来分析直播的那篇文献里面里面提到的 Adams, Kaminski et al. 数据，来自文献 10.1126/sciadv.aba1983 。文献标题为《Single-cell RNA-seq reveals ectopic and aberrant lung-resident cell populations in idiopathic pulmonary fibrosis》。发表：Sci Adv. 2020 Jul 8。
数据集为：adams_kaminski_2020_copd.h5ad，在文献里面的下载地址为 ：https://zenodo.org/records/7227571，在里面的 input_data.tar.xz
作者简单了解： 单位：耶鲁大学医学院和哈佛大学医学院布里格姆与妇女医院。
通讯为 耶鲁大学医学院的（肺部）医学教授Naftali Kaminski博士，是特发性肺纤维化（IPF）领域的先驱之一，一直致力于了解这种疾病的分子和细胞基础。通过单细胞测序，他的实验室建立了IPF细胞图谱，可帮助确定疾病特异性的转录图谱和疾病病理学中的功能。
图来自ref：https://mp.weixin.qq.com/s/7hKR9ZQE3RD0mmCAEUwl7g
本文通过分析32个特发性肺纤维化（IPF）样本、28个吸烟与非吸烟对照样本及18个慢性阻塞性肺疾病（COPD）样本中的312,928个细胞，构建了致死性间质性肺疾病——特发性肺纤维化（IPF）的单细胞图谱。
其他还有一些发现，我这里就想看看本文说的 一类未被报道过的异常基底样细胞群。
在IPF富集的上皮细胞中，我们鉴定出一类未被报道过的异常基底样细胞群。该类细胞同时表达基底上皮细胞标志物、间质标志物、衰老标志物及发育标志物，并定位于IPF肺组织中肌成纤维细胞灶边缘。
图谱构建 对 312,928 个细胞进行了分析，基于特异性标志物（图1C），鉴定出38种不同的细胞类型（图1B）。
数据：GSE136831。
可探索网址：www.ipfcellatlas.com
上皮细胞 肺上皮的亚群细分已经看过好几篇文献了，但是一直记不住细节，这里再来一遍！本文的亚群注释结果如下：
alveolar type 1 (AT1) 肺泡Ⅰ型细胞 alveolar type 2 (AT2) 肺泡Ⅱ型细胞 club cells 棒状细胞（clara cell，克拉拉细胞） pulmonary neuroendocrine cells 肺神经内分泌细胞 IPF 肺组织的上皮细胞组成特征表现为 ：气道源性上皮细胞比例显著增加（与对照组相比）basal cells 基底细胞、ciliated cells 纤毛细胞 和 goblet cells 杯状细胞，而肺泡上皮细胞数量则大幅减少（AT1和AT2细胞）。
上面这些细胞一次也没记住过，看看样子呢： AT1和AT2细胞 长相差这么多呢？还能相互转换？
(A) 气道上皮主要由纤毛细胞、杯状细胞和基底细胞组成，同时包含棒状细胞、离子细胞、神经内分泌细胞（NECs）、簇细胞和Deuterosomal细胞等更为罕见的细胞类型。
(B) 基底细胞作为气道上皮祖细胞的作用及其分化上皮细胞的谱系层级关系。图中标明了各细胞（亚）型最重要的细胞标志物。
图ref：https://www.frontiersin.org/journals/allergy/articles/10.3389/falgy.2021.787128/full
图ref：https://www.frontiersin.org/journals/cell-and-developmental-biology/articles/10.3389/fcell.2022.1022457/full
一类未被报道过的异常基底样细胞群 在上皮细胞中，作者鉴定出一类在转录水平上区别于以往已描述的肺内任何上皮细胞类型的细胞群，并将其命名为 异常基底样细胞 （图2C）。
除上皮标志物外，这些细胞的特点：
同时表达基底细胞标志物 TP63、KRT17、LAMB3 和LAMC2，但未表达其他已知基底标志物如KRT5与KRT15
它们还表达上皮-间质转化标志物（如VIM、CDH2、FN1、COL1A1、TNC、HMGA2）
衰老相关基因（包括CDKN1A、CDKN2A、CCND1、CCND2、MDM2、GDF15）（图2C右）
并在所有已检测细胞类型中显示出最高的标准化表达水平——包括已确定的IPF相关分子如基质金属肽酶7（MMP7）（13）、整合素αVβ6亚基（14）以及EPHB2（15）
转录因子基序富集分析： 这些细胞表现出高水平表达及下游激活证据的SOX9（图2C上）——该转录因子对远端气道发育（16）与修复（17）至关重要。
来源： 在鉴定出的483个异常基底样细胞中，448个来源于IPF肺组织，仅33个来源于COPD肺组织；该细胞群占IPF全部上皮细胞的3.5%，而在COPD中仅占1.1%（卡方检验P < 2.2 × 10⁻¹⁶）。
免疫组化： 使用p63、KRT17、HMGA2、COX2和p21作为标志物，证实了这些细胞在IPF肺组织中的存在并进行了定位（图2D；对照染色结果见图S5）。在IPF肺中，这些细胞始终定位于覆盖肌成纤维细胞灶的上皮层。
独立验证结果 ：重新分析Reyfman等人发表的IPF单细胞数据，鉴定出80个异常基底样细胞，其标志物谱与本文所述细胞具有相似特征（图S6）。我们观察到不同数据集间异常基底样细胞的关联性（Spearman相关系数rho=0.81）高于各数据集内部不同上皮细胞类型间的关联性（图2E），从而证实了本研究的结果。
数据：GSE136831 加载python模块： # %% # %load_ext autoreload # %autoreload 2 import scanpy as sc import numpy as np import pandas as pd import scipy.sparse as sp from scanpy_helpers.annotation import AnnotationHelper from nxfvars import nxfvars sc.set_figure_params(figsize=( 5 , 5 )) 读取h5ad数据： # %% adata = sc.read_h5ad(filename= "data/12_input_adatas/adams_kaminski_2020_copd.h5ad" ) adata 共312928个细胞。
一些其他属性查看： adata.obs adata.obs[ 'CellType_Category' ].value_counts() adata.obs[ 'Manuscript_Identity' ].value_counts() adata.obs[ 'Subclass_Cell_Identity' ].value_counts() adata.obs[ 'condition' ].value_counts() adata.obs[ 'sample' ].value_counts() adata.obs[ 'patient' ].value_counts() 给的信息还挺多：有详细的注释信息。
走一波降维聚类分群 走一波降维聚类分群，直接把已有注释信息放上来看看：
# 计算线粒体基因比例 adata.var[ 'mt' ]=adata.var_names.str.startswith( 'MT-' ) # 计算核糖体基因比例 adata.var[ 'ribo' ]=adata.var_names.str.match( '^RP[SL]' ) # 计算红血细胞基因比例 adata.var[ 'hb' ]=adata.var_names.str.match( '^HB[^(P)]' ) # 计算 sc.pp.calculate_qc_metrics(adata,qc_vars=[ 'mt' , 'ribo' , 'hb' ],percent_top=None,log1p=False,inplace=True) # 质控 adata.obs.head() # 每个基因在多少个细胞中表达 adata.var.head() # 可视化细胞的上述比例情况 # jitter：可以设置散点抖动的宽度 # rotation：可以调整x轴标签的旋转角度，样本多可以设置90度 # size：设置散点的大小 sc.pl.violin(adata,[ 'n_genes_by_counts' , 'total_counts' , 'pct_counts_mt' ],groupby= 'sample' ,jitter=0.4,rotation=90,multi_panel=True) # 对上面小提琴的优化 from matplotlib.pyplot import rc_context with rc_context({ "figure.figsize" : (9, 10)}): sc.pl.violin( adata, [ 'n_genes_by_counts' , 'total_counts' , 'pct_counts_mt' ], groupby= "sample" , jitter=0.4,rotation=90,multi_panel=True, #stripplot =False, # remove the internal dots #inner ="box", # adds a boxplot inside violins ) 作者已经做完了过滤等。
这里不再进行数据过滤。
降维聚类： # 首先将数据矩阵标准化（校正文库大小）： sc.pp.normalize_total(adata,target_sum= 1e4 ) # 对数据进行log sc.pp.log1p(adata) adata.raw = adata.copy() # 高变基因 sc.pp.highly_variable_genes(adata) # 归一化 sc.pp.scale(adata) # pca分析 sc.tl.pca(adata,use_highly_variable= True ) # harmony整合 sc.external.pp.harmony_integrate(adata, 'sample' ) # 聚类 sc.pp.neighbors(adata,use_rep= 'X_pca_harmony' ,n_pcs= 30 ) sc.tl.umap(adata) # Using the igraph implementation and a fixed number of iterations can be significantly faster, especially for larger datasets sc.tl.leiden(adata, flavor= "igraph" , n_iterations= 2 ,resolution= 0.5 ) 运行了28分钟。
好像scanpy现在可以支持GPU了？等我去看看。
画一下图看看：
# umap图 from matplotlib.pyplot import rc_context with rc_context({ "figure.figsize" : ( 8 , 7 )}): sc.pl.umap(adata, color=[ "leiden" ], legend_loc= "on data" , size= 7 ) # umap图 from matplotlib.pyplot import rc_context with rc_context({ "figure.figsize" : ( 8 , 7 )}): sc.pl.umap(adata, color=[ "condition" ], legend_loc= "on data" , size= 7 ) 原有的注释结果： # umap图 from matplotlib.pyplot import rc_context with rc_context({ "figure.figsize" : ( 8 , 7 )}): sc.pl.umap(adata, color=[ "CellType_Category" , 'Manuscript_Identity' ], legend_loc= "on data" , size= 7 ) 文献发现的那一坨细胞在哪里看到了没，我看到了！
marker基因图： sc.pl.umap(adata, color=[ "EPCAM" , "PTPRC" , "NKG7" , "CD3D" , "VWF" , "ACTA2" , "DCN" , "CD79A" , "CD68" ]) 基本上可以跟上面的注释对应上。
最后保存一下数据： ## 保存数据 adata.write( "01_qc/adams_kaminski_2020_copd_umap.h5ad" ) 今天分享到这里，可以把里面新发现的这群细胞提取出来分析看看！
