本文轉(zhuǎn)載自:presentlife的簡(jiǎn)書、括囊無(wú)譽(yù)的簡(jiǎn)書
建議直接跳轉(zhuǎn)：康君愛(ài)上了蕊醬的簡(jiǎn)書
摘要
測(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)分布堵泽。這篇短文介紹了包的使用和工作流程修己。
Deseq2的基本流程
#準(zhǔn)備count tables
cnts <- matrix(rnbinom(n=1000, mu=100, size=1/0.5),ncol = 10) #構(gòu)建滿足隨機(jī)負(fù)二項(xiàng)分布的一組數(shù)
cond <- factor(rep(1:2,each=5))
#創(chuàng)建DEseqDataSet
dds2 <- DESeqDataSetFromMatrix(cnts,DataFrame(cond), ~cond)
#前過(guò)濾
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
#標(biāo)準(zhǔn)分析
dds2 <- DESeq(dds2)
res <- results(dds2)
##另一種檢驗(yàn)，LRT迎罗，似然率檢驗(yàn)
#ddsLRT <- DESeq(dds2,test="LRT",reduced=~1)
#resLRT <- results(ddsLRT)
#moderate lfc
resultsNames(dds2)
resLFC <- lfcShrink(dds2,coef = 2,type = "apeglm")
#提取結(jié)果
summary()睬愤、subset()
已有基因表達(dá)矩陣，通過(guò)給定樣本信息纹安，因子公式設(shè)計(jì)來(lái)進(jìn)行差異分析尤辱，最后可生成普通差異分析結(jié)果和收縮后的差異分析結(jié)果，詳細(xì)的后續(xù)厢岂。
示例
1光督、加載DESeq2并設(shè)置樣品信息
library(DESeq2) # 加載DESeq2包
countData <- raw_count_filt2 # 表達(dá)矩陣
condition <- factor(c("control","control","KD","KD")) # 定義condition
colData <- data.frame(row.names=colnames(countData), condition) # 樣品信息矩陣
2、構(gòu)建dds矩陣并進(jìn)行差異分析
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition ) # 構(gòu)建dds矩陣
head(dds) # 查看dds矩陣的前6行
dds2 <- DESeq(dds) # 對(duì)dds進(jìn)行Normalize
resultsNames（dds) # 查看結(jié)果的名稱
res <- results(dds) # 使用results()函數(shù)獲取結(jié)果塔粒，并賦值給res
head (res, n=5) # 查看res矩陣的前5行
mcols(res,use.names= TRUE) # 查看res矩陣每一列的含義
summary(res) # 對(duì)res矩陣進(jìn)行總結(jié)
3结借、提取差異分析結(jié)果
table(res$padj<0.05) # 取padj小于0.05的數(shù)據(jù)，得到743行
res <- res[order(res$padj),] # 按照padj的大小將res重新排列
diff_gene_deseq2 <- subset(res,padj < 0.05 & (log2FoldChange >1 | log2FoldChange < -1)) # 獲取padj小于0.05卒茬，表達(dá)倍數(shù)取以2為對(duì)數(shù)后絕對(duì)值大于1的差異表達(dá)基因船老，賦值給diff_gene_deseq2
head (diff_gene_deseq2, n=5) # 查看diff_gene_deseq2矩陣的前5行
diff_gene_deseq2 <- row.names(diff_gene_deseq2) # 提取diff_gene_deseq2的行名
head (diff_gene_deseq2, n=5)
resdata <- merge (as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
head (resdata,n=5)
write.csv(resdata, file="control_vs_akap95.csv") # 將結(jié)果寫入control_vs_akap95.csv文件
示例2
1咖熟、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)
2、前過(guò)濾
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
3薪韩、標(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移除
4、差異性分析
dds <- DESeq(dds)
res <- results(dds)
res <- results(dds, contrast=c("condition","treated","untreated"))
5坎匿、可選操作
1）LFC收縮
resultsNames(dds)
resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
2）任務(wù)并行
library("BiocParallel")
register(MulticoreParam(4))
3）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)
4）獨(dú)立假設(shè)權(quán)重
library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult
6盾剩、探索并導(dǎo)出結(jié)果
1）MA-plot
plotMA(res, ylim=c(-2,2))
plotMA(resLFC, ylim=c(-2,2))
注：效果好的話，收縮的圖形不會(huì)出現(xiàn)左寬右窄替蔬。
2）多種收縮方法
#查看你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")
3）counts繪圖
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")
4）結(jié)果的更多信息
mcols(res)$description
可以通過(guò)mcol函數(shù)來(lái)查看結(jié)果表中各個(gè)變量的含義。其中變量LFC就是變化的倍數(shù)承桥，存在NA的可能原因有：
1驻粟、一個(gè)基因在所有樣本中的表達(dá)值都為0
2、一個(gè)基因的某個(gè)樣本的表達(dá)值為離群值（通過(guò)Cook距離檢測(cè)）
3凶异、一個(gè)基因counts標(biāo)準(zhǔn)化后數(shù)值過(guò)低
后續(xù)：結(jié)果的導(dǎo)出和豐富的可視化
ReportingTools蜀撑，regionReport，Glimma剩彬，pcaExplorer等工具可以生成一個(gè)報(bào)表形式的結(jié)果酷麦。
