你以为纯粹的单细胞(系)至少也有3个亚群
考虑到咱们生信技能树粉丝对单细胞数据挖掘的需求，我开通了一个专栏《100个单细胞转录组数据降维聚类分群图表复现》，也亲自示范了几个，不过自己带娃，读博，时间精力有限，所以把剩余的90多个任务安排了学徒，实习生，学员。真的是太棒了，群策群力！
今天介绍的数据数据集**GSE150241**它是一个黑色素瘤细胞系A2058的10X单细胞转录组测序。
A2058 是一种已知对 BRAF 抑制具有抗性的黑色素瘤细胞系，而10X单细胞转录组测序这种方法使我们能够研究每个黑色素瘤细胞的转录情况，特别是那些表达 ABCB5 的细胞。
GSM4592552 A2058_scRNAseq
同样的，自己下载文件，然后复制粘贴本教程的全部代码，就可以复现出图表！
step1: 导入数据 构建seurat对象
下载好数据到工作目录下：
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)
library(clustree)
library(cowplot)
#读取单细胞数据，这里是h5文件
sce.all=CreateSeuratObject(Read10X_h5('GSM4592552_filtered_gene_bc_matrices_h5.h5'))
#计算线粒体基因比例
sce.all=PercentageFeatureSet(sce.all, pattern = "^MT-", col.name = "percent_mito")
#计算核糖体基因比例
sce.all=PercentageFeatureSet(sce.all, "^RP[SL]", col.name = "percent_ribo")
#计算红血细胞基因比例
sce.all=PercentageFeatureSet(sce.all, "^HB[^(P)]", col.name = "percent_hb")
#过滤
selected_mito <- WhichCells(sce.all, expression = percent_mito < 15)
selected_ribo <- WhichCells(sce.all, expression = percent_ribo > 3)
selected_hb <- WhichCells(sce.all, expression = percent_hb < 0.1)
sce.all.filt <- subset(sce.all, cells = selected_mito)
sce.all.filt <- subset(sce.all.filt, cells = selected_ribo)
sce.all.filt <- subset(sce.all.filt, cells = selected_hb)
step2：降维聚类分群
sce.all.filt = NormalizeData(sce.all.filt)
sce.all.filt = FindVariableFeatures(sce.all.filt)
sce.all.filt = ScaleData(sce.all.filt,
vars.to.regress = c("nFeature_RNA", "percent_mito"))
sce.all.filt = RunPCA(sce.all.filt, npcs = 20)
sce.all.filt = RunTSNE(sce.all.filt, npcs = 20)
sce.all.filt = RunUMAP(sce.all.filt, dims = 1:20)
# 聚类分析
sce.all.filt <- FindNeighbors(sce.all.filt, dims = 1:10)
for (res in c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8,1)) {
sce.all.filt <- FindClusters(sce.all.filt, graph.name = "RNA_snn", resolution = res, algorithm = 1)
}
clustree(sce.all.filt@meta.data, prefix = "RNA_snn_res.")
DimPlot(sce.all.filt, reduction = "tsne",group.by = "RNA_snn_res.0.2",label = T,label.box = T)+NoAxes()
可以看到这里有明显的4个亚群
根据一些常见的marker gene 我们没法很好的对单个细胞系的肿瘤细胞进行注释。
那我们来查看每个亚群的top10基因
Idents(sce.all.filt)="RNA_snn_res.0.2"
sce.all.markers <- FindAllMarkers(object = sce.all.filt, only.pos = TRUE, min.pct = 0.25,
thresh.use = 0.25)
top10 <- sce.all.markers %>% group_by(cluster) %>% top_n(10, avg_log2FC)
DoHeatmap(sce.all.filt,top10$gene,size=3)
从热图中可以看到，肿瘤细胞中也存在异质性。接着对每个亚群的top基因进行富集分析
sce.markers=sce.all.markers
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=0.05)
dotplot(xx)+ theme(axis.text.x = element_text(angle = 45,
vjust = 0.5, hjust=0.5))
前面我们提到过 ABCB5 是一种参与侵袭性和对 BRAF 抑制的抵抗力的细胞表面蛋白。应该是可以检查一下它和黑素细胞基因 (MITF)的表达量，看看有没有共表达或者互斥的现象。
cg=c('ABCB5','ALK','BRAF','MITF');cg
DotPlot(sce.all.filt,features = cg)
FeaturePlot(sce.all.filt,features = cg)
需要一定的生物学背景才能更好地解释这个黑色素瘤细胞系的异质性！另外，如果要更好地完成上面的图表，通常是需要掌握5个R包，分别是: scater,monocle,Seurat,scran,M3Drop，需要熟练掌握它们的对象，：一些单细胞转录组R包的对象 而且分析流程也大同小异：
step1: 创建对象 step2: 质量控制 step3: 表达量的标准化和归一化 step4: 去除干扰因素(多个样本整合) step5: 判断重要的基因 step6: 多种降维算法 step7: 可视化降维结果 step8: 多种聚类算法 step9: 聚类后找每个细胞亚群的标志基因 step10: 继续分类
单细胞转录组数据分析的标准降维聚类分群，并且进行生物学注释后的结果。可以参考前面的例子：人人都能学会的单细胞聚类分群注释 ，我们演示了第一层次的分群。
如果你对单细胞数据分析还没有基础认知，可以看基础10讲：
写在文末
咱们现在这个专栏《100个单细胞转录组数据降维聚类分群图表复现》分享的代码是到此为止，但是一般来说单细胞文章数据分析还有很多进阶图表制作，比如inferCNV看肿瘤拷贝数变异，monocle看拟时序等等。如果你也需要，可以加入我们这个专栏《100个单细胞转录组数据降维聚类分群图表复现》创作团队，获取进阶指引哦！见：急！计划招募100个单细胞爱好者，免费学全套单细胞降维聚类分群和生物学亚群注释