采用DESeq2对表达量进行PCA和聚类分析-CFANZ编程社区

欢迎关注”生信修炼手册”!

得到基因/转录本的表达量之后，通常会通过以下三种类型的图表来检验和分析生物学样本和实验设计间关系。

1. 样本的聚类树

利用所有样本的表达量数据，对样本进行聚类。理论上如果样本和实验操作都没有问题，那么属于同一组的生物学重复样本会聚到一起。示意图如下

采用DESeq2对表达量进行PCA和聚类分析_数据

上图中，样本的名称用组别代替，可以看到，同一条件的样本聚在了一起。

2. PCA图

通过主成分分析进行降维，在二维或者三维平面上展示样本点的分布，根据点的位置，也可以看出属于同一组的样本是否在一起，不同组之间的样本有没有明显分开，示意如下

采用DESeq2对表达量进行PCA和聚类分析_数据_02
从图中可以看到，不同条件的样本区分的很明显，而生物学重复之间距离较近，表明生物学重复的一致性和不同分组的差异性较好。

3. 热图

相比样本的聚类树，热图包含了更多的信息，比如可以直观的展示不同分组间表达量的差异，也是常见的可视化手段之一，示意如下

采用DESeq2对表达量进行PCA和聚类分析_聚类_03

只要有样本的表达量矩阵，DESeq2可以轻松的画出以上3种图表。但是我们应该选择原始的表达量矩阵，还是归一化之后的表达量矩阵来画呢？或者有没有其他的选择呢？

输入的矩阵不同，得出的结论也会不同。由于基因的表达水平在不同样本间本身就存在一定的差异，所以无论是采用原始的还是归一化之后的表达量矩阵，效果都不理想。针对这一问题，DESeq2提出了两种count值的转换算法，rlog和VST转换。

1. rlog 转换

rlog 转换的用法如下

rld <- rlog(dds)

2. VST 转换

用法如下

vsd <- vst(dds)

两种转换本质上是在降低生物学重复之间的差异，使得样本聚类和PCA分析的效果更好。转换之后的表达量数据可以采用assay函数进行提取，代码如下

> head(assay(rld)[, 1:2])
       sample1   sample2
gene1 2.049029 1.6828707
gene2 8.151262 6.8552583
gene3 0.818971 0.2964686
gene4 5.340361 4.4766682
gene5 6.316175 6.8345783
gene6 2.157821 1.9264385

对于raw count定量表格，建议采用rlog或者VST转换之后的数据去进行PCA和聚类分析，效果会更好。

利用DESeq2提供的示例数据pasilla，分别用原始的count, 归一化之后的count, rlog, vst 转换的count 进行PCA分析，代码如下

dds <- estimateSizeFactors(dds)
raw <- SummarizedExperiment(counts(dds, normalized=FALSE),
                                colData=colData(dds))
nor <- SummarizedExperiment(counts(dds, normalized=TRUE),
                                colData=colData(dds))
vsd <- vst(dds)
rld <- rlog(dds)
pdf("PCA.pdf")
plotPCA( DESeqTransform(raw), intgroup=c("condition", "type") )
plotPCA( DESeqTransform(nor), intgroup=c("condition", "type") )
plotPCA(vsd, intgroup=c("condition", "type"))
plotPCA(rld, intgroup=c("condition", "type"))
dev.off()