肿瘤样品的单细胞需要提取上皮细胞继续细分

肿瘤样品的单细胞需要提取上皮细胞继续细分 by 单细胞天地

单细胞图谱时代早就过去了，不再是随便选取一个物种，一种组织或者一种疾病，挑选几个样品做单细胞，简单的说明清楚其单细胞亚群组成比例和生物学意义就足够的时机了。但是图谱既然是被做烂了的思路，就应该是有系统性梳理的价值，比如肿瘤样品的单细胞就应该是首先是按照如下所示的标记基因进行第一次分群：

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讲：

本次我们介绍的重点是对上皮细胞进行降维聚类分群，以及单细胞高级分析，文章里面的主要的图表我们分开介绍：

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转录因子分析

转录因子分析和细胞通讯分析属于单细胞数据分析里面的高级分析了，主要是对计算机资源消耗会比较大，我也多次分享过细节教程：

如下所示：

转录因子分析

看起来已经是非常完美了，单独提取上皮细胞并且进行了细分亚群，搞清楚每个亚群的生物学功能富集以及其特异性的转录因子激活。因为是肿瘤研究，也可以加上一个画龙点睛的生存分析从上面那么多的转录因子里面挑选出有统计学显著区分生存的。

生存分析是画龙点睛：

我在生信技能树多次分享过生存分析的细节；

如下所示的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单细胞转录组成本摆在那里，参考我们的《明码标价》专栏里面的单细胞内容

仅仅是测数据每个样品就100G以上，对计算资源的消耗也很大。还等什么呢，赶快扫描下面二维码添加微信咨询吧！

ixxmu / mp_duty