单细胞学习

R包Seurat instruction2 主成分分析

2019-04-17  本文已影响0人  阿糖胞苷_SYSU

2019-2-1 来自Seurat官网https://satijalab.org/seurat/

主成分分析

1.进行降维分析

用上一步缩放的数据进行PCA分析。

默认输入object@var.genes中的基因,可以定义为使用 pc.genes。用上一步发现的高度可变基因进行降维可以改善效果。

但是UMI数据,特别是去除了技术变异,当用更大的基因集合去跑时,PCA的结果基本一致,包括整个转录组。

具体command:

pbmc <- RunPCA(object = pbmc, pc.genes = pbmc@var.genes, do.print = TRUE, pcs.print = 1:5,

    genes.print = 5)

##可视化PCA结果,方法有:PrintPCA,VizPCA,PCAPlot和PCHeatmap

PrintPCA(object = pbmc, pcs.print = 1:5, genes.print = 5, use.full = FALSE)

VizPCA(object = pbmc, pcs.use = 1:2)

PCAPlot(object = pbmc, dim.1 = 1, dim.2 = 2)

#ProjectPCA 可以基于基因与计算得出的成分的关联性进行评分(包括未参与PCA的基因)。可以找到与细胞异质性强相关但未能被可变基因筛选出来的细胞标志物。通过设置use.full=T可以执行上述功能。

pbmc <- ProjectPCA(object = pbmc, do.print = FALSE)

##PCHeatmap比较特别,可以针对异质性来源进行简单探索,这对决定用哪个PC进行下游分析很有用。根据PCA评分(PCA scores)对细胞和基因进行排序。

##将cell.use设置为能够将“最极端”的细胞划分为两端的数值,可以显著提高大数据库的速度。这属于监督分析,是探索关联基因集的有力工具。

PCHeatmap(object = pbmc, pc.use = 1, cells.use = 500, do.balanced = TRUE, label.columns = FALSE)

PCHeatmap(object = pbmc, pc.use = 1:12, cells.use = 500, do.balanced = TRUE,  label.columns = FALSE, use.full = FALSE)

2 决定最显著的主成分

为了克服技术噪音,Seurat基于PCA scores进行细胞聚类。每个PC代表一个包含关联基因的信息的基因集(metagene)。纳入多少PC进入下游分析是重要的一步。

http://www.cell.com/abstract/S0092-8674(15)00549-8

## 用jackStraw程序进行重采样的test,随机转化了一个子集(默认1%),重新进行PCA,构建了一个gene scores的“空白分布”(null distribution),重复进行。

##将富集了低P值基因的PC作为“显著”PC。这一程序耗时较长,可以PCElbowPlot()减少运行时间。

pbmc <- JackStraw(object = pbmc, num.replicate = 100, display.progress = FALSE)

##JackStrawPlot提供可视化的工具,和统一标准分布(虚线表示)比较每个PC的P值分布情况,显著的PC会呈现低P值基因的富集(虚线上方的实线表示)。

JackStrawPlot(object = pbmc, PCs = 1:12)

JackStrawPlot

#另一种特别的方式是利用PCElbowPlot绘制主成分的标准差的散点图,将曲线的拐点作为cutoff

值。

PCElPlot

PC的选择,确定数据集真实维度,对Seurat是重要的一步,也是最有挑战性的一步。建议以下3种方法:

1.更加监督性,通过探索PC寻找异质性的相关来源,可用于联合GSEA分析。

2.基于随机空白模型进行test,对于大数据集来说比较费时,且可能不能呈现比较好的PCcutoff值。

3.常用的启发式算法,可以快速计算。

例子中得到的结果基本一致,但在PC7-10之间选择哪个作为cutoff需要判断。

最终选择了JackStraw的结果,因为PC热图的结果很好的解释了PC的意义。尽管cutoff小范围的变动并不影响结果,但强烈建议对所选PC进行探索。

3 补充:

1 载入数据library(Seurat)

pbmc <- readRDS(file = "~/Projects/datasets/pbmc3k_final.rds")

2.探索新的降维结构

上一篇下一篇

猜你喜欢

热点阅读