简单使用DESeq做差异分析

DESeq这个R包主要针对count data，其数据来源可以是RNA-Seq或者其他高通量测序数据。类似地，对于CHIP-Seq数据或者质谱肽段数据也是使用的。由于DESeq是一个R包，因此使用它需要一点点R基础语法。

首先需要读入一个数据框，列代表每个sample，行代表每个gene
```
 database_all <- read.table(file = "readcount", sep = "\t", header = T, row.names = 1)
 database <- database_all[,1:6]
```
这里主要对于两两比较的数据，因此我取了数据的前6列，分别是两组样品，每组3个生物学重复
设定分组信息，也就是样本分组的名称
```
 type <- factor(c(rep("LC_1",3), rep("LC_2",3)))
```
我这里是样品1是LC_1，样品2是LC_2
由于DESeq包要求接下来的count data必须要整数型，因此我们需要对数据进行取整,然后将数据database和分组信息type读入到cds对象中
```
 database <- round(as.matrix(database))
 cds <- newCountDataSet(database,type)
```
接下里对于不同类型的数据要进行不同的处理，可以粗略分为有生物学重复数据、有部分生物学重复数据以及无生物学重复数据

4.1 有生物学重复
```
 cds <- estimateSizeFactors(cds)
 cds <- estimateDispersions(cds)
 res <- nbinomTest(cds,"LC_1","LC_2")
```
其通过estimate the dispersion并对count data进行标准化，然后得到每个gene做T test检验

4.2 对于部分样品有生物学重复
```
 cds <- estimateSizeFactors(cds)
 cds <- estimateDispersions(cds)
 res <- nbinomTest(cds,"LC_1","LC_2")
```
其步骤跟有上述的一样的，DESeq会根据有生物学重复的样品来estimate the dispersion，当然要保证unreplicated condition does not have larger variation than the replicated one

4.3 对于没有生物学重复
```
 cds <- estimateDispersions(cds, method="blind", sharingMode="fit-only" )
 res <- nbinomTest(cds,"LC_1","LC_2")
```
注意参数method="blind" 和 sharingMode="fit-only"即可

最后就是查看符合阈值的差异基因有多少个即可，然后将结果输出到csv文件中方便查看

 table(res$padj <0.05)
 res <- res[order(res$padj),]
 sum(res$padj<=0.01,na.rm = T)
 write.csv(resdata,file = "LC_1_vs_LC_2_DESeq.csv")

Summary

DESeq在前几年的文章中经常被使用，但是现在有了其升级版DESeq2，后者相比前者对于犯第一类错误卡的并不是那么严格了，所以在同样的padj的阈值下，筛选到的差异基因的数目也会多一点。

并且上述中，我都只使用nbinomTest来获得由显著差异的gene，其实DESeq包还提供了其他检验方法，比如nbinomGLMTest函数，具体可以查看DESeq文档，还有其他对应对因素分析的使用方法。

本文出自于http://www.bioinfo-scrounger.com转载请注明出处