摘要
測(cè)序產(chǎn)生的數(shù)據(jù)內(nèi)容為，每個(gè)樣本中，每個(gè)基因分配到多少個(gè)測(cè)序片段桩蓉，各種類型的測(cè)序（RNA-seq，CHIP-Seq抓艳，HiC）產(chǎn)生的數(shù)據(jù)類似触机。RNA-seq數(shù)據(jù)分析的一個(gè)基本任務(wù)就是檢測(cè)差異性表達(dá)的基因。其中一個(gè)重要的問(wèn)題就是對(duì)不同條件下產(chǎn)生的系統(tǒng)學(xué)差異進(jìn)行量化和統(tǒng)計(jì)學(xué)推斷玷或。DESeq2使用負(fù)二項(xiàng)式廣義線性模型來(lái)檢測(cè)基因表達(dá)的差異性儡首，對(duì)分散和LFC的估計(jì)包含數(shù)據(jù)的先驗(yàn)分布。這篇短文介紹了包的使用和工作流程偏友。
1 標(biāo)準(zhǔn)流程
簡(jiǎ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á)矩陣蔬胯，通過(guò)給定樣本信息，因子公式設(shè)計(jì)來(lái)進(jìn)行差異分析位他，最后可生成普通差異分析結(jié)果和收縮后的差異分析結(jié)果氛濒，詳細(xì)的后續(xù)。
如何獲取幫助
Bioconductor
導(dǎo)入的數(shù)據(jù)
- 使用未標(biāo)準(zhǔn)化的數(shù)據(jù)
DESeq2內(nèi)部會(huì)根據(jù)樣本大小對(duì)counts進(jìn)行調(diào)整鹅髓，自帶標(biāo)準(zhǔn)化過(guò)程舞竿。 - DESeqDataSet
其在DESeq2中是一種類型，在代碼中常用dds來(lái)表示窿冯，其實(shí)例對(duì)象用來(lái)存儲(chǔ)counts和中間估計(jì)量骗奖。中間估計(jì)量中就包括跨基因收集的信息。這個(gè)對(duì)象必須包含design formula醒串，用來(lái)構(gòu)建模型的變量（分組信息）执桌。這個(gè)對(duì)象可以通過(guò)四種上游途徑來(lái)構(gòu)建：轉(zhuǎn)錄豐度文件；counts矩陣芜赌；htseq-count文件仰挣；SummarizedExperiment對(duì)象。這里只記錄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對(duì)象
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樣本信息順序要一致
- 附加信息操作一般用不上
- 前過(guò)濾
代碼：
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
注：
并不是必須颓芭，不影響計(jì)算結(jié)果，這樣做的優(yōu)點(diǎn)是dds對(duì)象占用的內(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ù)操作前指定馍惹，如果沒(méi)有指定躺率，采用的是默認(rèn)
- levels對(duì)應(yīng)的分別為玛界，分母，分子
- ref為指定的對(duì)照組悼吱，即分母
- 當(dāng)公式中的某個(gè)變量對(duì)應(yīng)的樣本沒(méi)有的時(shí)候慎框，可以通過(guò)dropleveles移除
- 折疊重復(fù)
DESeq2提供collapseReplicates函數(shù)進(jìn)行去重復(fù)（非生物平行），詳見(jiàn)手冊(cè)后添。
差異性分析
代碼：
dds <- DESeq(dds)
res <- results(dds)
res <- results(dds, contrast=c("condition","treated","untreated"))
差異性分析的計(jì)算和估計(jì)過(guò)程整合到了一個(gè)函數(shù)DESeq中笨枯，其中的具體細(xì)節(jié)步驟后續(xù)會(huì)有介紹。對(duì)DESeq分析產(chǎn)生的結(jié)果遇西，通過(guò)results函數(shù)生成結(jié)果表馅精。
注：
- 可以通過(guò)contrast參數(shù)設(shè)置生成結(jié)果表的特定因子
- LFC收縮
代碼：
resultsNames(dds)
resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
就是將DESeq函數(shù)處理后生成的的dds對(duì)象傳遞給lfcShrink函數(shù)即可，參數(shù)后續(xù)粱檀。
- 任務(wù)并行
代碼：
library("BiocParallel")
register(MulticoreParam(4))
就是利用一個(gè)包洲敢，開(kāi)啟多線程。
- 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ǎn)單總結(jié)函數(shù)茄蚯，來(lái)得到想要的初步結(jié)果压彭。
- 獨(dú)立假設(shè)權(quán)重
DESeq2中p值的矯正使用的是IHW包，具體原理見(jiàn)對(duì)應(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))
效果好的話壮不，收縮的圖形不會(huì)出現(xiàn)左寬右窄。
- 多種收縮方法
#查看你coef的順序序號(hào)
resultsNames(dds)
#通過(guò)序號(hào)代替全稱皱碘；三種收縮方法
resNorm <- lfcShrink(dds, coef=2, type="normal")
resAsh <- lfcShrink(dds, coef=2, type="ashr")
resApe <- lfcShrink(dds, coef=2, type="apeglm")
目前DESeq2提供了三種可選的方法询一，不同LFC擬合的分布來(lái)進(jìn)行后續(xù)的收縮操作，具體細(xì)節(jié)參考對(duì)應(yīng)方法的具體文檔尸执。
- counts繪圖
代碼：
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")
繪制特定基因的counts圖家凯，可以通過(guò)設(shè)置return=T，來(lái)返回一個(gè)用于ggplot可以進(jìn)一步設(shè)置具體參數(shù)的data.frame對(duì)象如失。
- 結(jié)果的更多信息
代碼：
mcols(res)$description
可以通過(guò)mcol函數(shù)來(lái)查看結(jié)果表中各個(gè)變量的含義绊诲。其中變量LFC就是變化的倍數(shù)，存在NA的可能原因有：
- 一個(gè)基因在所有樣本中的表達(dá)值都為0
- 一個(gè)基因的某個(gè)樣本的表達(dá)值為離群值（通過(guò)Cook距離檢測(cè)）
- 一個(gè)基因counts標(biāo)準(zhǔn)化后數(shù)值過(guò)低
- 結(jié)果的導(dǎo)出和豐富的可視化
ReportingTools褪贵，regionReport掂之，Glimma，pcaExplorer等工具可以生成一個(gè)報(bào)表形式的結(jié)果脆丁。 - 導(dǎo)出結(jié)果到csv文件
利用R基礎(chǔ)的write.csv就可以
多因素設(shè)計(jì)
用的時(shí)候再說(shuō)
2 數(shù)據(jù)轉(zhuǎn)換及可視化
- count數(shù)據(jù)的轉(zhuǎn)換
- 數(shù)據(jù)質(zhì)量的檢測(cè)
3 流程中的可變步驟
4 DESeq2理論基礎(chǔ)
5 常見(jiàn)問(wèn)題
參考資料
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
