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

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了简单的转录组差异基因表达分析 -- DESeq2相关的知识,希望对你有一定的参考价值。

参考技术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组的是上调还是下调。另外,本文的差异分析还是处于单因子水平(只有一个变量),至于多因子的差异分析以后研究透了再和大家进行分享。

完。

以上是关于简单的转录组差异基因表达分析 -- DESeq2的主要内容,如果未能解决你的问题,请参考以下文章

转录组差异分析流程三大R包比较

转录组入门(7):差异表达分析

转录组数据标准化--Normalization

Corset轻松搞定无参转录组差异基因(转载)

R语言DESeq2基因差异表达分析

怎么分析关注的功能基因集在转录组结果中表现如何?