RNA-seq差异基因分析

生信菜鸟屋

接着上文的RNA-seq数据分析后，我们得到了一个基因表达矩阵，其中包括了原始count，fpkm，TPM。

随后我们需要选择一款差异分析的软件，目前的话主流的软件主要有三个limma，edger，DESeq2。如果是有生物学重复的情况下一般选择DESeq2（该文献被引用了3000多次，是一款非常权威的差异分析工具）。如果没有生物学重复的话就是用edger；虽然看样子没有生物学重复好像也可以做差异分析，但是个人是非常不推荐，生物学重复尽量在三个以上。

接下来简单介绍一下DESeq2的原理。首先在早期的时候，做差异分析的时候使用的是泊松分布做的（其实就是观察到的数据与期望值相比有没有差异），但是泊松分布有一个特点就是数据的平均数与方差是相等的。然后测序的基因表达矩阵的均值与方差是不等的，因此使用泊松分布是存在一定的问题的。 于是科学家们发现负二项分布可以很好的拟合基因表达数据，因此DESeq2就选择了该模型来进行差异分析。那么具体是怎么做的呢？

首先，DESeq2输入的是原始count数据；然后使用中位数标准化的方式计算样本间的缩放因子s，随后使用负二项分布对原始count进行拟合计算出其均值u，之后假设样本中某个基因的真实浓度为q，它们之间满足u=sq，这样就可计算出该基因在这个样本中真实的浓度是多少了；最后使用广义线性模型来计算该基因在处理间的差异表达。在这里是主要介绍了一下DESeq2做差异分析的过程。其实其内部是十分复杂的，例如处理样本数太少的时候会共享基因的表达量来进行压缩等过程。

接下来直接上代码实操，首先DESeq2的输入是原始count：

Mycount <-read.csv (“readcount. csv”)

数据格式见下图：

随后构建meta信息，就是告诉软件什么样本是什么处理：

condition <- factor(c(rep("control",2),rep("treat",2)), levels = c("control","treat"))
colData <- data.frame(row.names=colnames(mycounts), condition)

该meta信息也可以使用excle进行编辑后读取，就是下面这样子。

接下来就是构建DESeq2所需要的dds对象：

dds <- DESeqDataSetFromMatrix(countData = mycounts, colData = colData, design= ~ batch + condition)

注意这里的batch是一些可能对结果产生影响的协变量。例如研究病人和正常人的时候，性别就是一个协变量。

之后进行标准化及差异分析：

dds <- DESeq(dds)

最后就是提取结果了：

res <- results(dds)

之后有了这个差异结果后，大家就可以用来做后续的分析了，例如经典的火山图等。

微信扫一扫分享文章