摘要
測序產(chǎn)生的數(shù)據(jù)內(nèi)容為，每個(gè)樣本中，每個(gè)基因分配到多少個(gè)測序片段，各種類型的測序（RNA-seq，CHIP-Seq，HiC）產(chǎn)生的數(shù)據(jù)類似。RNA-seq數(shù)據(jù)分析的一個(gè)基本任務(wù)就是檢測差異性表達(dá)的基因。其中一個(gè)重要的問題就是對不同條件下產(chǎn)生的系統(tǒng)學(xué)差異進(jìn)行量化和統(tǒng)計(jì)學(xué)推斷。DESeq2使用負(fù)二項(xiàng)式廣義線性模型來檢測基因表達(dá)的差異性，對分散和LFC的估計(jì)包含數(shù)據(jù)的先驗(yàn)分布。這篇短文介紹了包的使用和工作流程。
1 標(biāo)準(zhǔn)流程
簡單例子
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")
已有基因表達(dá)矩陣，通過給定樣本信息，因子公式設(shè)計(jì)來進(jìn)行差異分析，最后可生成普通差異分析結(jié)果和收縮后的差異分析結(jié)果，詳細(xì)的后續(xù)。
如何獲取幫助
Bioconductor
導(dǎo)入的數(shù)據(jù)
- 使用未標(biāo)準(zhǔn)化的數(shù)據(jù)
DESeq2內(nèi)部會根據(jù)樣本大小對counts進(jìn)行調(diào)整，自帶標(biāo)準(zhǔn)化過程。 - DESeqDataSet
其在DESeq2中是一種類型，在代碼中常用dds來表示，其實(shí)例對象用來存儲counts和中間估計(jì)量。中間估計(jì)量中就包括跨基因收集的信息。這個(gè)對象必須包含design formula，用來構(gòu)建模型的變量（分組信息）。這個(gè)對象可以通過四種上游途徑來構(gòu)建：轉(zhuǎn)錄豐度文件；counts矩陣；htseq-count文件；SummarizedExperiment對象。這里只記錄counts矩陣方法。 - count 矩陣導(dǎo)入
代碼：
#使用pasilla包中附帶的數(shù)據(jù)
library("pasilla")
#讀入表達(dá)矩陣，樣本注釋
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包，構(gòu)建dds對象
library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ condition)
#為結(jié)果附加描述信息
featureData <- data.frame(gene=rownames(cts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
注：
- counts數(shù)值為整數(shù)
- 列名順序和colData樣本信息順序要一致
- 附加信息操作一般用不上
- 前過濾
代碼：
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
注：
并不是必須，不影響計(jì)算結(jié)果，這樣做的優(yōu)點(diǎn)是dds對象占用的內(nèi)存小點(diǎn)，后續(xù)的計(jì)算耗時(shí)小點(diǎn)。
- 標(biāo)注因子水平
代碼：
dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
dds$condition <- relevel(dds$condition, ref = "untreated")
dds$condition <- droplevels(dds$condition)
注：
- 在進(jìn)行后續(xù)操作前指定，如果沒有指定，采用的是默認(rèn)
- levels對應(yīng)的分別為，分母，分子
- ref為指定的對照組，即分母
- 當(dāng)公式中的某個(gè)變量對應(yīng)的樣本沒有的時(shí)候，可以通過dropleveles移除
- 折疊重復(fù)
DESeq2提供collapseReplicates函數(shù)進(jìn)行去重復(fù)（非生物平行），詳見手冊。
差異性分析
代碼：
dds <- DESeq(dds)
res <- results(dds)
res <- results(dds, contrast=c("condition","treated","untreated"))
差異性分析的計(jì)算和估計(jì)過程整合到了一個(gè)函數(shù)DESeq中，其中的具體細(xì)節(jié)步驟后續(xù)會有介紹。對DESeq分析產(chǎn)生的結(jié)果，通過results函數(shù)生成結(jié)果表。
注：
- 可以通過contrast參數(shù)設(shè)置生成結(jié)果表的特定因子
- LFC收縮
代碼：
resultsNames(dds)
resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
就是將DESeq函數(shù)處理后生成的的dds對象傳遞給lfcShrink函數(shù)即可，參數(shù)后續(xù)。
- 任務(wù)并行
代碼：
library("BiocParallel")
register(MulticoreParam(4))
就是利用一個(gè)包，開啟多線程。
- 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值集合一些簡單總結(jié)函數(shù)，來得到想要的初步結(jié)果。
- 獨(dú)立假設(shè)權(quán)重
DESeq2中p值的矯正使用的是IHW包，具體原理見對應(yīng)文獻(xiàn)及文檔。
代碼：
library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult
探索并導(dǎo)出結(jié)果
- MA-plot
代碼：
plotMA(res, ylim=c(-2,2))
plotMA(resLFC, ylim=c(-2,2))
效果好的話，收縮的圖形不會出現(xiàn)左寬右窄。
- 多種收縮方法
#查看你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擬合的分布來進(jìn)行后續(xù)的收縮操作，具體細(xì)節(jié)參考對應(yīng)方法的具體文檔。
- counts繪圖
代碼：
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")
繪制特定基因的counts圖，可以通過設(shè)置return=T，來返回一個(gè)用于ggplot可以進(jìn)一步設(shè)置具體參數(shù)的data.frame對象。
- 結(jié)果的更多信息
代碼：
mcols(res)$description
可以通過mcol函數(shù)來查看結(jié)果表中各個(gè)變量的含義。其中變量LFC就是變化的倍數(shù)，存在NA的可能原因有：
- 一個(gè)基因在所有樣本中的表達(dá)值都為0
- 一個(gè)基因的某個(gè)樣本的表達(dá)值為離群值（通過Cook距離檢測）
- 一個(gè)基因counts標(biāo)準(zhǔn)化后數(shù)值過低
- 結(jié)果的導(dǎo)出和豐富的可視化
ReportingTools，regionReport，Glimma，pcaExplorer等工具可以生成一個(gè)報(bào)表形式的結(jié)果。 - 導(dǎo)出結(jié)果到csv文件
利用R基礎(chǔ)的write.csv就可以
多因素設(shè)計(jì)
用的時(shí)候再說
2 數(shù)據(jù)轉(zhuǎn)換及可視化
- count數(shù)據(jù)的轉(zhuǎn)換
- 數(shù)據(jù)質(zhì)量的檢測
3 流程中的可變步驟
4 DESeq2理論基礎(chǔ)
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
