最近,总有小伙伴和我说:“我跟着Seurat的官网跟了一遍教程,但是处理自己的数据的时候,效果总是不让人满意,有没有什么办法?”今天就来试试用调整参数的方式来看看能不能让Seurat的UMAP可视化变得好一些。
在开始之前,先一起回顾一下用Seurat进行单细胞分析流程吧~
Seurat分析流程
质量控制
在拿到一套单细胞数据时,要先对细胞进行筛选,一般情况下,需要根据以下三种情况来进行筛选:
1、基因含量过少的细胞(测序不可靠)
2、基因含量过多的细胞(剔除双细胞或多细胞)
3、线粒体含量过高的细胞(这些很可能是破裂的细胞)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
经过过滤之后,剩余的细胞,才是单细胞分析真正要处理的细胞。
标准化
为使不同细胞之间的基因表达水平具有可比性,必须进行标准化处理。即将每个细胞的基因表达量归一化(Normalization)到其总表达量,然后乘以一个缩放因子scale.factor。最后,对所得的表达水平进行对数转换,使表达值更符合正态分布。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
挑选高变基因
由于并不是所有基因能够在识别不同细胞时具有相同的贡献度和信息量,因此需要将这些具有识别作用的特征基因挑选出来,以用于进一步的探索。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
选取的高变基因的数量依赖于nfeatures参数,一般情况下,这个数值在1000~5000之间,具体哪个数值合适,需要经过多次迭代尝试不同的值后检验结果而决定。
数据缩放
在挑选出高变基因后,需要对数据进行缩放处理,以此来让这些基因对细胞的贡献处于同一水平,并在这一步当中剔除一些不需要的基因。在这一步当中,根据不同的细胞,features和vars.to.regress参数是可变的,不过在这里不多做讲述。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
主成分分析——PCA线性降维
通过这一步,将进行缩放后的基因分成了指定数量的主成分,将具有同质性的基因汇聚到一个主成分中(PC)。经过主成分分析之后,大多数具有贡献的基因都被集中到前几个主成分中,剩下的大多是没有意义的噪音。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), npcs = 50)
在这行代码中,npcs就是指定分出多少个主成分的数量的参数。至于选取多少个主成分进行分析,官方给出的建议是通过绘制肘部图(Elbowplot),选取其中的肘部来决定。
plot1 <- ElbowPlot(pbmc)
通过这种方式,将有意义的基因挑选出来,剔除不必要的噪音之后,进行进一步的探索。
细胞聚类
pbmc <- FindNeighbors(pbmc, dims = 1:30)
pbmc <- FindClusters(pbmc, resolution = 0.5)
通过这二行代码进行细胞聚类,就可以将差异较小的细胞归类到一起。
UMAP非线性降维分析
通过UMAP方法,将挑选出来的主成分进行非线性降维,使得结果可视化。除了UMAP外,还有一种常用的方法是tSNE,在这里以讨论UMAP为主,在Seurat中,进行UMAP非线性降维的代码如下。
pbmc <- RunUMAP(pbmc, dims = 1:10, n.neighbors = 30L, min.dist = 0.3)
plot1 <- DimPlot(pbmc, reduction = "umap", label = TRUE)
通过以上代码,就可以在选定主成分后进行UMAP降维,并通过图表进行可视化。在进行可视化之后,可以通过绘制一些已知细胞类型的经典标志基因的特征图,来进行基因注释。
通过以上的流程和代码,使用seurat教程中使用的参数处理了一套数据,可以得到以下的UMAP降维图。
可以看到第0,1,2,3聚类之间并没有分的很开,同质性并不好。有没有什么办法让这张可视化的图变得更加好看一些,具有区分性一些呢?基于之前提到的代码和其中的参数,让我们来进行一些讨论。
高变基因的挑选对UMAP投影的影响
按照单细胞分析流程的顺序,首先要做的就是挑选高变基因。在不改变其他参数时,挑选不同的高变基因(nfeatures),查看对这套数据集能够产生什么样的影响。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
可以看到,随着高变基因的取值变化,可视化的UMAP图也发生了变化。当nfeatures为1000时,细胞完全交织在一起,随着nfeatures变多、簇与簇的交界也变得更清晰,但当n=5000时,数量过多的特征基因会淹没原本的特征信号。经过多次迭代调整后,最终决定选取nfeatures=4000进行进一步探索。
主成分分析——挑选的主成分越多越好吗?
在确定了高变基因选取为4000之后,在进行主成分分析时,需要给定一个npcs参数,来决定划分出多少个主成分,而这个参数的不同,会导致每一个主成分的内容都发生改变。而在进行UMAP降维分析前,需要通过如下的一个肘图来挑选具体需要多少个主成分来进行聚类。当信息接近饱和时,可以从肘图当中看到一个肘部(拐点),在下图当中大约是横坐标为15的点,因此取前15个主成分进行UMAP。(注:在上一个小节中用于UMAP的主成分是前10个,仅用于迭代找出合适的nfeatures)
在这里,有人可能会产生疑惑:肘部的标准是什么?
根据Seurat官方的描述,没有肘部的绝对定义,如下图所示,PC10~20中的任何一个点都可以视为一个肘部,排位靠前的PC含有需要的基因信息,而靠后的PC通常含有大量的噪音,具体需要选择哪个点实际上需要经过多次尝试,不过一般来说,当信息量趋近饱和时,实际效果并不会相差太远。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc), npcs = 50)
pbmc <- RunUMAP(pbmc, dims = 1:15, n.neighbors = 30L, min.dist = 0.3)
plot1 <- DimPlot(pbmc, reduction = "umap", label = TRUE)
和之前相比产生了非常明显的变化,在之前由于信息不饱和,导致聚类分的并不是那么令人满意,比如在1:10时的第三和第五聚类,直接被分出了2个部分。而在取的PC增多后,这个问题得到了一定程度的解决。那么,如果在PCA时分出更多的主成分,或者是划定更多的主成分进行UMAP分析,是否会对结果产生影响呢?
结论是几乎没有区别!说明无论npcs的参数选取了多少,有意义的特征基因都已经被分到了前几个PC中,当信息量趋近于饱和时,进行UMAP非线性降维可视化的结果并没有区别。
在关注一下选取的PC的范围。实际上,关于肘部并没有一个明确的定义,在后几个PC中也许会含有一些有趣的基因,但是选取的PC过多又会让信息量过大,淹没前几个PC中有趣的基因。如何选取合适的PC范围也是需要技巧。在已经确认npcs不会影响UMAP的情况下,这里取npcs参数的值为50,分别查看取前5,20,50个PC后进行UMAP降维后的结果。
可以看出,在取不同范围的主成分(PC)后,UMAP的投影也发生了较大的变化。在dims=5时,由于取的主成分过少,导致取出来的结果根本没有参考价值。而当使用的PC多到一定程度后,对细胞的实际分群影响便不大了,说明信息已经趋近于饱和(dims=15和dims=20区别不大)。那是不是意味着取的主成分越多便越好呢?并不是的。当取了所有主成分后,发现细胞反而分不开了。大量的噪音反而影响了对细胞的聚类,因此取的主成分并不是越多越好。而在这套数据集中,进行UMAP非线性降维的主成分取前15~25个主成分都是比较合适的,在后续的分析中,经过多次尝试,npcs取前15个。
min.dist和n.neighbour,让图变得赏心悦目
在选取好主成分(PC)的范围后,经过一定的处理,最后来到UMAP非线性降维。在这一步当中,min.dist决定了点之间的最小距离,而较大的n.neighbour值会导致更多的全局结构得以保留,但会损失详细的局部结构。这二项参数并不改变每个成分或聚类的内容,因此不会对注释产生影响,但是会对可视化产生巨大的影响。
pbmc <- RunUMAP(pbmc,dims = 1:15, n.neighbors = 30L, min.dist = 0.3)
plot1 <- DimPlot(pbmc, reduction = "umap", label = TRUE)
可以看出,n.neighbour越大,分群越发紧密,更加接近一个水滴的形状,保留的整体性越强;而n.neighbour越小,细胞聚类之间就越发的松散,展示的细节就更多。根据多次尝试,当n.neighbour=15L时,可以达到比较好的视觉效果。
而对于min.dist而言,可以调整点与点之间的最小距离,改变聚类之间的紧密程度。在固定n.neighbour=15L的情况下分别尝试不同的min.dist,结果如下所示:
通过对比,可以看到细胞投影的整体形状没有改变,但是点之间的距离发生了变化。经过多次尝试,最终还是选取min.dist=0.3作为最终的选择。
通过上述流程进行参数调整,确定了四个比较重要的参数,再和原先的图进行对比,可以看到视觉效果好了许多:
那么,通过这种方式,就可以让UMAP可视化的效果在迭代中不断地提升,不知道对老师您是否有帮助呢?如果觉得不错的话,给我们点一个赞吧!感谢您的支持和鼓励!
注意!!!!一般来讲,在合理范围内调整参数并不会影响基因注释结果,但仍然存在特殊情况,请务必在保证注释结果不变的情况下进行参数的调整!
嘿嘿这里偷偷藏了一点常用代码,用+号相连就可以使用啦:
ggtitle()#为图表添加标题
theme(plot.title = element_text(hjust = 0.5))#调整标题的位置,现在的参数是居中
theme(legend.position = "none") #设置图例为无
theme(aspect.ratio = 1)#将图表输出为正方形
theme(plot.title = element_text(family = " "))#将图标标题字体设置为对应字体
theme(plot.title = element_text(size = "9"))#字号为9磅
theme(plot.title = element_text(color = "red"))#字体设置为红色
theme(plot.title = element_text(face = "bold"))#字体加粗
DimPlot(pbmc, reduction = "umap") #一些可用的调色用参数:group.by = "cell_type", features = "gene_name", 点大小与透明度 pt.size = 0.5, alpha = 0.8 聚类标签label = TRUE, 自动调整标签位置repel = TRUE
希望能对各位老师作图有所帮助哦!
微生信助力高分文章,用户320000+,谷歌学术7100+
转载本文请联系原作者获取授权,同时请注明本文来自陈明杰科学网博客。
链接地址:https://wap.sciencenet.cn/blog-707141-1502067.html?mobile=1
收藏