简单的转录组差异基因表达分析 -- DESeq2

Posted 2023-03-09

参考技术A

经典的转录组差异分析通常会使用到三个工具 limma/voom , edgeR 和 DESeq2 。今天我们就通过一个小规模的转录组测序数据来演示 DESeq2 的简单流程。

对于 DESeq2 的分析流程而言，我们需要输入的数据包括：

​ 下面就以 mobData 中的数据为例简单介绍 DESeq2 的分析流程

​ 由于 mobData 中的行名没有提供基因的ID，我们也不是为了探究生物学问题，就以 mobData 的行数作为其ID

DESeqDataSet 是 DESeq2 流程中储存read counts和中间统计分析数据的对象，之后的分析都建立在该对象之上进行。

​ 在进行差异分析之前，需要对样本数据的表达矩阵进行预处理，包括：

​ 通过PCA结果来看各组样本分组情况还是不错的，但hclust的聚类结果反映的分组就略微有点混杂了，可能要聚类计算的距离函数选用不当有关。

使用 DESeq() 函数进行差异分析时，该函数干了以下三件事：

​ counts() 可以提取 DESeq object 中的表达矩阵，而 results() 可以提取差异分析的结果，其中包括了：

样本间的均值, log2 fold changes, standard errors, test statistics, p-values and adjusted p-values.

​ 使用 results() 函数时需要指明进行比较的样本，这里用 contrast=c("group_list","MM","WW") 提取 MM 组和 WW 组进行差异分析的结果。如果想要比较 WM 组和 WW 组，只要改变 contrast=c("group_list","WM","WW") 即可。

​ 检查结果中是否包含 NA 值

​ 这里 padj 中有1142个 NA 值是因为使用 results() 提取差异分析结果时，大于 alpha 值（这里是0.1）的矫正后p-value都会被当做是 NA 。因此，我们将这些 padj 值都设为 1

排序后以 log2FoldChange 绝对值大于1， padj 小于0.05为条件筛选显著的差异表达基因

​ 至此，便筛选出了 217个 在 MM 组和 WW 组之间的显著差异表达基因。至于后续的可视化分析则是因课题而异了，等以后有空了再补坑吧！

​ 有同学可能注意到，虽然我们的样本有多个组，但在差异分析时进行的还是pairwise的分析，为什么我们不可以三个组一起分析呢？学过ANOVA的同学应该都知道，ANOVA就是可以应对这种多组差异分析的情况。但要注意的是，ANOVA只可以告诉我们对于某个基因在这三组中是否存在差异，想要找出是哪一组有其他组别有差异还是需要进行pairwise t-test之类的分析。所以，在这里我们两组两组地进行分析正是出于这个考虑，并且有更方便我们解释差异分析的结果，说明在A组基因的表达量相对于B组的是上调还是下调。另外，本文的差异分析还是处于单因子水平（只有一个变量），至于多因子的差异分析以后研究透了再和大家进行分享。

完。