简单使用DESeq做差异分析

DESeq这个R包主要针对count data,其数据来源可以是RNA-Seq或者其他高通量测序数据。类似地,对于CHIP-Seq数据或者质谱肽段数据也是使用的。

由于DESeq是一个R包,因此使用它需要一点点R基础语法。

  1. 首先需要读入一个数据框,列代表每个sample,行代表每个gene

    database_all <- read.table(file = "readcount", sep = "\t", header = T, row.names = 1)
    database <- database_all[,1:6]
    

    这里主要对于两两比较的数据,因此我取了数据的前6列,分别是两组样品,每组3个生物学重复

  2. 设定分组信息,也就是样本分组的名称

    type <- factor(c(rep("LC_1",3), rep("LC_2",3)))
    

    我这里是样品1是LC_1,样品2是LC_2

  3. 由于DESeq包要求接下来的count data必须要整数型,因此我们需要对数据进行取整,然后将数据database和分组信息type读入到cds对象中

    database <- round(as.matrix(database))
    cds <- newCountDataSet(database,type)
    
  4. 接下里对于不同类型的数据要进行不同的处理,可以粗略分为有生物学重复数据、有部分生物学重复数据以及无生物学重复数据

    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”即可

  5. 最后就是查看符合阈值的差异基因有多少个即可,然后将结果输出到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文档,还有其他对应对因素分析的使用方法。