如何采用DESeq2对表达量进行PCA和聚类分析

如何采用DESeq2对表达量进行PCA和聚类分析，相信很多没有经验的人对此束手无策，为此本文总结了问题出现的原因和解决方法，通过这篇文章希望你能解决这个问题。

站在用户的角度思考问题，与客户深入沟通，找到石嘴山网站设计与石嘴山网站推广的解决方案，凭借多年的经验，让设计与互联网技术结合，创造个性化、用户体验好的作品，建站类型包括：网站设计、成都网站制作、企业官网、英文网站、手机端网站、网站推广、申请域名、雅安服务器托管、企业邮箱。业务覆盖石嘴山地区。

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

1. 样本的聚类树

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

如何采用DESeq2对表达量进行PCA和聚类分析

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

2. PCA图

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

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

3. 热图

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

如何采用DESeq2对表达量进行PCA和聚类分析

只要有样本的表达量矩阵，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()

raw count 的结果如下

如何采用DESeq2对表达量进行PCA和聚类分析
归一化之后count结果如下