如何采用DESeq2对表达量进行PCA和聚类分析
如何采用DESeq2对表达量进行PCA和聚类分析，相信很多没有经验的人对此束手无策，为此本文总结了问题出现的原因和解决方法，通过这篇文章希望你能解决这个问题。
得到基因/转录本的表达量之后，通常会通过以下三种类型的图表来检验和分析生物学样本和实验设计间关系。
1. 样本的聚类树
利用所有样本的表达量数据，对样本进行聚类。理论上如果样本和实验操作都没有问题，那么属于同一组的生物学重复样本会聚到一起。示意图如下
上图中，样本的名称用组别代替，可以看到，同一条件的样本聚在了一起。
2. PCA图
通过主成分分析进行降维，在二维或者三维平面上展示样本点的分布，根据点的位置，也可以看出属于同一组的样本是否在一起，不同组之间的样本有没有明显分开，示意如下
从图中可以看到，不同条件的样本区分的很明显，而生物学重复之间距离较近，表明生物学重复的一致性和不同分组的差异性较好。
3. 热图
相比样本的聚类树，热图包含了更多的信息，比如可以直观的展示不同分组间表达量的差异，也是常见的可视化手段之一，示意如下
只要有样本的表达量矩阵，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 sample2gene1 2.049029 1.6828707gene2 8.151262 6.8552583gene3 0.818971 0.2964686gene4 5.340361 4.4766682gene5 6.316175 6.8345783gene6 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 的结果如下
归一化之后count结果如下
VST转换之后的结果如下
rlog转换之后的结果如下
可以很明显看出，原始的count和归一化之后的count, 其PCA图是杂乱无序的，没什么明显规律，而VST和rlog转换之后，生物学重复之间更佳的接近，不同分组也区分的较为明显。
看完上述内容，你们掌握如何采用DESeq2对表达量进行PCA和聚类分析的方法了吗？如果还想学到更多技能或想了解更多相关内容，欢迎关注行业资讯频道，感谢各位的阅读！
