Bioconductor:DESeq2
摘要
测序产生的数据内容为，每个样本中，每个基因分配到多少个测序片段，各种类型的测序（RNA-seq，CHIP-Seq，HiC）产生的数据类似。RNA-seq数据分析的一个基本任务就是检测差异性表达的基因。其中一个重要的问题就是对不同条件下产生的系统学差异进行量化和统计学推断。DESeq2使用负二项式广义线性模型来检测基因表达的差异性，对分散和LFC的估计包含数据的先验分布。这篇短文介绍了包的使用和工作流程。
1 标准流程
简单例子
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, name="condition_trt_vs_untrt")
res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")
已有基因表达矩阵，通过给定样本信息，因子公式设计来进行差异分析，最后可生成普通差异分析结果和收缩后的差异分析结果，详细的后续。
如何获取帮助
Bioconductor
导入的数据
- 使用未标准化的数据
DESeq2内部会根据样本大小对counts进行调整，自带标准化过程。 - DESeqDataSet
其在DESeq2中是一种类型，在代码中常用dds来表示，其实例对象用来存储counts和中间估计量。中间估计量中就包括跨基因收集的信息。这个对象必须包含design formula，用来构建模型的变量（分组信息）。这个对象可以通过四种上游途径来构建：转录丰度文件；counts矩阵；htseq-count文件；SummarizedExperiment对象。这里只记录counts矩阵方法。 - count 矩阵导入
代码：
#使用pasilla包中附带的数据
library("pasilla")
#读入表达矩阵，样本注释
pasCts <- system.file("extdata",
"pasilla_gene_counts.tsv",
package="pasilla", mustWork=TRUE)
pasAnno <- system.file("extdata",
"pasilla_sample_annotation.csv",
package="pasilla", mustWork=TRUE)
cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))
coldata <- read.csv(pasAnno, row.names=1)
coldata <- coldata[,c("condition","type")]
#样本注释修改，排序
rownames(coldata) <- sub("fb", "", rownames(coldata))
cts <- cts[, rownames(coldata)]
#加载DESeq2包，构建dds对象
library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ condition)
#为结果附加描述信息
featureData <- data.frame(gene=rownames(cts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
注：
- counts数值为整数
- 列名顺序和colData样本信息顺序要一致
- 附加信息操作一般用不上
- 前过滤
代码：
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
注：
并不是必须，不影响计算结果，这样做的优点是dds对象占用的内存小点，后续的计算耗时小点。
- 标注因子水平
代码：
dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
dds$condition <- relevel(dds$condition, ref = "untreated")
dds$condition <- droplevels(dds$condition)
注：
- 在进行后续操作前指定，如果没有指定，采用的是默认
- levels对应的分别为，分母，分子
- ref为指定的对照组，即分母
- 当公式中的某个变量对应的样本没有的时候，可以通过dropleveles移除
- 折叠重复
DESeq2提供collapseReplicates函数进行去重复（非生物平行），详见手册。
差异性分析
代码：
dds <- DESeq(dds)
res <- results(dds)
res <- results(dds, contrast=c("condition","treated","untreated"))
差异性分析的计算和估计过程整合到了一个函数DESeq中，其中的具体细节步骤后续会有介绍。对DESeq分析产生的结果，通过results函数生成结果表。
注：
- 可以通过contrast参数设置生成结果表的特定因子
- LFC收缩
代码：
resultsNames(dds)
resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
就是将DESeq函数处理后生成的的dds对象传递给lfcShrink函数即可，参数后续。
- 任务并行
代码：
library("BiocParallel")
register(MulticoreParam(4))
就是利用一个包，开启多线程。
- p值和矫正p值
代码：
resOrdered <- res[order(res$pvalue),]
summary(res)
sum(res$padj < 0.1, na.rm=TRUE)
res05 <- results(dds, alpha=0.05)
sum(res05$padj < 0.05, na.rm=TRUE)
可以利用p值和矫正p值集合一些简单总结函数，来得到想要的初步结果。
- 独立假设权重
DESeq2中p值的矫正使用的是IHW包，具体原理见对应文献及文档。
代码：
library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult
探索并导出结果
- MA-plot
代码：
plotMA(res, ylim=c(-2,2))
plotMA(resLFC, ylim=c(-2,2))
效果好的话，收缩的图形不会出现左宽右窄。
- 多种收缩方法
#查看你coef的顺序序号
resultsNames(dds)
#通过序号代替全称；三种收缩方法
resNorm <- lfcShrink(dds, coef=2, type="normal")
resAsh <- lfcShrink(dds, coef=2, type="ashr")
resApe <- lfcShrink(dds, coef=2, type="apeglm")
目前DESeq2提供了三种可选的方法，不同LFC拟合的分布来进行后续的收缩操作，具体细节参考对应方法的具体文档。
- counts绘图
代码：
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")
绘制特定基因的counts图，可以通过设置return=T，来返回一个用于ggplot可以进一步设置具体参数的data.frame对象。
- 结果的更多信息
代码：
mcols(res)$description
可以通过mcol函数来查看结果表中各个变量的含义。其中变量LFC就是变化的倍数，存在NA的可能原因有：
- 一个基因在所有样本中的表达值都为0
- 一个基因的某个样本的表达值为离群值（通过Cook距离检测）
- 一个基因counts标准化后数值过低
- 结果的导出和丰富的可视化
ReportingTools，regionReport，Glimma，pcaExplorer等工具可以生成一个报表形式的结果。 - 导出结果到csv文件
利用R基础的write.csv就可以
多因素设计
用的时候再说
2 数据转换及可视化
- count数据的转换
- 数据质量的检测
3 流程中的可变步骤
4 DESeq2理论基础
5 常见问题
参考资料
Michael I. Love, Simon Anders, and Wolfgang Huber.2018."Analyzing RNA-seq data with DESeq2".http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
