关于“数据的维度”(dims参数)的选择
完成PCA之后,我们获得了该数据集的所有主成分(PCs)信息,但是如何决定纳入多少个主成分进行下游分析呢?
主要参考以下方法:
热图
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
如上图所示,可以看出前15个主成分可以把细胞分成差异明显的两群,说明前15个主成分中含有的显著的差异基因更多,主成分也就更有意义,所以下游分析可以纳入前15个PCs。
碎石图 Elbow plot
ElbowPlot(pbmc)
通过碎石图可以看出每个PC对变异的贡献情况,从上图可以看出9~10PC以后逐渐趋于稳定(噪声主导),也就是说真实信号主要来自前10个左右的PCs,所以可以选择前10个PCs进行下游分析。
JackStraw法
随机置换一部分数据(默认为1%),然后重新 PCA,重复此过程。将包含较多低 P 值特征的主成分为「重要的」主成分。
JackStraw()
函数可以计算出每个主成分中各基因的P值,用于判断哪些主成分更具有统计学意义,ScoreJackStraw()
用于量化主成分的显著性强度,富含低P值基因较多的主成分更有统计学意义。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot()
函数可视化比较每个主成分的 p 值分布和均匀分布(虚线)。在这个例子中,在前 10 到 12 个主成分之后,主成分的重要性开始下降。
JackStrawPlot(pbmc, dims = 1:15)
JackStraw法相当于计算每个主成分的p值,根据p值选择显著性的PCs纳入下游分析,比较科学,但是当数据量比较大时,计算非常慢,并且可视化不够直观,个人推荐仅用于参考,不作为首选。
究极方法:循环
for(i in c(5,10,15,18,20,23,26,28,30,32,35,40)){
All.merge.singlets.new <- FindNeighbors(object = All.merge.singlets.new, dims = 1:i, verbose = T)
All.merge.singlets.new <- FindClusters(object = All.merge.singlets.new, resolution = 0.8, verbose = T)
All.merge.singlets.new <- RunUMAP(object = All.merge.singlets.new, dims = 1:i, verbose = T)
All.merge.singlets.new <- RunTSNE(object = All.merge.singlets.new, dims = 1:i, verbose = T,check_duplicates = FALSE)
gd1 <- DimPlot(object = All.merge.singlets.new, reduction = "umap", group.by = "orig.ident") #+ NoLegend()
gd2 <- DimPlot(object = All.merge.singlets.new, reduction = "umap",label = T) #+ NoLegend()
# gd3 <- DimPlot(object = Batch.merge.singlets, reduction = "umap", group.by = "celltype") #+ NoLegend()
gd4 <- DimPlot(object = All.merge.singlets.new, reduction = "tsne", group.by = "orig.ident") #+ NoLegend()
gd5 <- DimPlot(object = All.merge.singlets.new, reduction = "tsne",label = T) #+ NoLegend()
# gd6 <- DimPlot(object = Batch.merge.singlets, reduction = "tsne", group.by = "celltype") #+ NoLegend()
CombinePlots(plots = list(gd1,gd2,gd4,gd5), ncol = 2)
# CombinePlots(plots = list(gd1,gd2,gd3,gd4,gd5,gd6), ncol = 3)
ggsave(filename=paste('All_sample_UMap_tSNE_by_cluster_and_sample_DimPlot_pc',i,'.pdf',sep=''), width=20, height=13)
#saveRDS(dc, file = 'dc_tSNE_UMAP.rds')
}
通过循环将不同的dims值代入,进行细胞聚类分析,导出pdf查看每个dims取值时细胞分群情况及降维图,比如从dims=15
往后的分群变化都不大了,那么即可选择dims=15
。