肿瘤样品的单细胞需要提取上皮细胞继续细分
单细胞图谱时代早就过去了,不再是随便选取一个物种,一种组织或者一种疾病,挑选几个样品做单细胞,简单的说明清楚其单细胞亚群组成比例和生物学意义就足够的时机了。但是图谱既然是被做烂了的思路,就应该是有系统性梳理的价值,比如肿瘤样品的单细胞就应该是首先是按照如下所示的标记基因进行第一次分群 :
immune (CD45+,PTPRC), epithelial/cancer (EpCAM+,EPCAM), stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)
然后每个亚群进行第二层次细分亚群,甚至第三层次,第四次分群,结构清晰明了。我们以上皮细胞亚群的 细分来举例说明每个分析点的工作量:
现在我们要介绍的是发表于2020的文章,标题 是:《Single-cell transcriptomic architecture and intercellular crosstalk of human intrahepatic cholangiocarcinoma》,数据集在;https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138709
是肝癌里面的 intrahepatic cholangiocarcinoma 这个细分疾病,本研究共5.6万个细胞,取样策略如下,
可以看到研究者的第一层次降维聚类分群并不是我们前面的提到的上皮,免疫和间质细胞,而是直接进行了第二层次分群,一般来说背诵下面的基因即可:
# T Cells (CD3D, CD3E, CD8A),
# B cells (CD19, CD79A, MS4A1 [CD20]),
# Plasma cells (IGHG1, MZB1, SDC1, CD79A),
# Monocytes and macrophages (CD68, CD163, CD14),
# NK Cells (FGFBP2, FCG3RA, CX3CR1),
# Photoreceptor cells (RCVRN),
# Fibroblasts (FGF7, MME),
# Endothelial cells (PECAM1, VWF).
# epi or tumor (EPCAM, KRT19, PROM1, ALDH1A1, CD24).
# immune (CD45+,PTPRC), epithelial/cancer (EpCAM+,EPCAM),
# stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)
library(ggplot2)
genes_to_check = c('PTPRC', 'CD3D', 'CD3E', 'CD4','CD8A',
'CD19', 'CD79A', 'MS4A1' ,
'IGHG1', 'MZB1', 'SDC1',
'CD68', 'CD163', 'CD14',
'TPSAB1' , 'TPSB2', # mast cells,
'RCVRN','FPR1' , 'ITGAM' ,
'C1QA', 'C1QB', # mac
'S100A9', 'S100A8', 'MMP19',# monocyte
'LAMP3', 'IDO1','IDO2',## DC3
'CD1E','CD1C', # DC2
'KLRB1','NCR1', # NK
'FGF7','MME', 'ACTA2', ## fibo
'DCN', 'LUM', 'GSN' , ## mouse PDAC fibo
'Amy1' , 'Amy2a2', # Acinar_cells
'PECAM1', 'VWF', ## endo
'EPCAM' , 'KRT19', 'PROM1', 'ALDH1A1' )
单细胞转录组数据分析的标准降维聚类分群,并且进行生物学注释后的结果。可以参考前面的例子:人人都能学会的单细胞聚类分群注释 ,我们演示了第一层次的分群。
如果你对单细胞数据分析还没有基础认知,可以看基础10讲:
01. 上游分析流程 02.课题多少个样品,测序数据量如何 03. 过滤不合格细胞和基因(数据质控很重要) 04. 过滤线粒体核糖体基因 05. 去除细胞效应和基因效应 06.单细胞转录组数据的降维聚类分群 07.单细胞转录组数据处理之细胞亚群注释 08.把拿到的亚群进行更细致的分群 09.单细胞转录组数据处理之细胞亚群比例比较
本次我们介绍的重点是对上皮细胞进行降维聚类分群,以及单细胞高级分析,文章里面的主要的图表我们分开介绍:
A) tSNE plots and CNV box plots for 6 distinct epithelial cell subclusters. B) Heatmap showing the top 10 DEGs (Wilcoxon test) in 6 epithelial cell subclusters. C) Differences in pathway activity (scored per cell by GSVA) in 4 malignant cell sub- clusters. D) Heatmap of the t-value for the area under the curve score of expression regulation by transcription factors, as estimated using SCENIC. E) Violin plots showing the expression of marker genes in distinct malignant subclusters. F) IHC staining showing the expression of CK19, TM4SF4, PROM1, CDH6, and SPINK1 in recurrent sample (ICC-20T) and adjacent liver tissue. G) Kaplan–Meier survival curve of TM4SF4, SPINK1, and KRT19 expression using the optimal group cut-off point in patients with ICC (from TCGA).
标准降维聚类分群是6个细胞亚群(上皮细胞的细分)
其中0,1,2,3是有拷贝数变异的癌细胞,而4和5是正常水平细胞,分别是cholangiocytes和hepatocytes,如下所示:
可以很明显的看到标号为4和5的亚群在癌旁部分,而且CNV程度超级低,所以是正常细胞 :
Subclusters 4 (cholangiocytes, top DEGs: FXYD2 and RHOB) Subclusters 5 (hepatocytes, top DEGs: ASGR1 and APOE)
而有拷贝数变异的癌细胞就异质性很大了,下面是有拷贝数变异的癌细胞的4个亚群的各自介绍 :
ubcluster 0 malignant cells were char- acterized by high expression levels of mesenchymal markers such as COL1A1, fibronectin, and IGFBP7, indicating epithelial- mesenchymal transition (EMT) characteristics. Subcluster 1 displayed high expression levels of the malignancy-promoting factors S100P and FABP5; subcluster 2 exhibited high expres- sion levels of the immune-associated genes CD74 and HLA-DRA; subcluster 3 malignant cells, from the recurrent patient, displayed high expression levels of SPINK1
各个上皮细胞亚群的特异性基因
如下所示,各自单细胞亚群的特异性基因分别是:
当然了,这样的每个亚群挑选自己生物学背景范围内的基因的策略,对纯粹的生物信息学工程师来说,非常的不友好, 因为大家能背诵的基因屈指可数。
gsva对各个细胞亚群进行生物学功能数据库注释
不过,我们生信工程师有大杀器,就是生物学功能数据库注释, 包括go和kegg的,通常是 使用 clusterProfiler 包进行 :
# 首先前面的降维聚类分群找到了 sce.markers 各个亚群的标记基因 :
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)
ids=bitr(sce.markers$gene,'SYMBOL','ENTREZID','org.Hs.eg.db')
sce.markers=merge(sce.markers,ids,by.x='gene',by.y='SYMBOL')
gcSample=split(sce.markers$ENTREZID, sce.markers$cluster)
gcSample # entrez id , compareCluster
xx <- compareCluster(gcSample, fun="enrichKEGG",
organism="hsa",
pvalueCutoff= 1,
qvalueCutoff=1)
p=dotplot(xx)
p+ theme(axis.text.x = element_text(angle = 45,
vjust = 0.5, hjust=0.5))
p
ggsave('compareCluster-sce.markers.pdf')
这里作者使用的gsva和gsea这样的策略,GSEA分析,我也多次讲解:
但实际上,绝大部分读者并没有去细看这个统计学原理,也不需要知道gsea分析的nes值如何计算,需要知道的是nes值的生物学意义。
而单细胞转录组数据的批量GSVA代码我在:借鉴escape包的一些可视化GSVA或者ssGSEA结果矩阵的方法 剧透过,其实跟 对单细胞表达矩阵做gsea分析的代码实现环节大同小异。无论是gsva还是gsea,本质上都是看某个生物学基因集的得分。
当然了,生物学功能数据库注释也有可能并不是很细致,还可以尝试SCENIC这样的转录因子分析,这样就定位到了具体的基因。
SCENIC转录因子分析
转录因子分析和细胞通讯分析属于单细胞数据分析里面的高级分析了,主要是对计算机资源消耗会比较大,我也多次分享过细节教程:
张泽民团队的单细胞研究把T细胞分的如此清楚 细胞通讯分析的背景知识 构建单细胞亚群网络(类似于细胞通讯分析) 细胞通讯分析结果的解读 SCENIC转录因子分析结果的解读 人人都能学会的单细胞聚类分群注释 对单细胞表达矩阵做gsea分析 单细胞转录因子分析之SCENIC流程
如下所示:
看起来已经是非常完美了,单独提取上皮细胞并且进行了细分亚群,搞清楚每个亚群的生物学功能富集以及其特异性的转录因子激活。因为是肿瘤研究,也可以加上一个画龙点睛的生存分析从上面那么多的转录因子里面挑选出有统计学显著区分生存的。
生存分析是画龙点睛:
我在生信技能树多次分享过生存分析的细节;
人人都可以学会生存分析(学徒数据挖掘) 学徒数据挖掘之谁说生存分析一定要按照表达量中位值或者平均值分组呢? 基因表达量高低分组的cox和连续变量cox回归计算的HR值差异太大? 学徒作业-两个基因突变联合看生存效应 TCGA数据库里面你的基因生存分析不显著那就TMA吧 对“不同数据来源的生存分析比较”的补充说明 批量cox生存分析结果也可以火山图可视化 既然可以看感兴趣基因的生存情况,当然就可以批量做完全部基因的生存分析 多测试几个数据集生存效应应该是可以找到统计学显著的! 我不相信kmplot这个网页工具的结果(生存分析免费做) 为什么不用TCGA数据库来看感兴趣基因的生存情况 200块的代码我的学徒免费送给你,GSVA和生存分析 集思广益-生存分析可以随心所欲根据表达量分组吗 生存分析时间点问题 寻找生存分析的最佳基因表达分组阈值 apply家族函数和for循环还是有区别的(批量生存分析出图bug) TCGA数据库生存分析的网页工具哪家强 KM生存曲线经logRNA检验后也可以计算HR值
如下所示的3个基因被生存分析挑出来了展现给大家。
这里面有一些缩略词:
CNV, copy number variation DEGs, differentially expressed genes GSVA, gene set variation analysis ICC, intrahepatic cholangiocarcinoma IHC, immunohistochemistry tSNE, t-distributed stochastic neighbor embedding. SCENIC,Single-Cell rEgulatory Network Inference and Clustering
因为10x单细胞转录组成本摆在那里,参考我们的《明码标价》专栏里面的单细胞内容