本文转载自:presentlife的简书、括囊无誉的简书
建议直接跳转：康君爱上了蕊酱的简书

摘要

测序产生的数据内容为，每个样本中，每个基因分配到多少个测序片段，各种类型的测序（RNA-seq，CHIP-Seq，HiC）产生的数据类似。RNA-seq数据分析的一个基本任务就是检测差异性表达的基因。其中一个重要的问题就是对不同条件下产生的系统学差异进行量化和统计学推断。DESeq2使用负二项式广义线性模型来检测基因表达的差异性，对分散和LFC的估计包含数据的先验分布。这篇短文介绍了包的使用和工作流程。

Deseq2的基本流程

#准备count tables
cnts <- matrix(rnbinom(n=1000, mu=100, size=1/0.5),ncol = 10) #构建满足随机负二项分布的一组数
cond <- factor(rep(1:2,each=5))

#创建DEseqDataSet
dds2 <- DESeqDataSetFromMatrix(cnts,DataFrame(cond), ~cond)
#前过滤
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

#标准分析
dds2 <- DESeq(dds2)
res <- results(dds2)
##另一种检验，LRT，似然率检验
#ddsLRT <- DESeq(dds2,test="LRT",reduced=~1)
#resLRT <- results(ddsLRT)

#moderate lfc
resultsNames(dds2)
resLFC <- lfcShrink(dds2,coef = 2,type = "apeglm")

#提取结果
summary()、subset()

已有基因表达矩阵，通过给定样本信息，因子公式设计来进行差异分析，最后可生成普通差异分析结果和收缩后的差异分析结果，详细的后续。

示例

1、加载DESeq2并设置样品信息

library(DESeq2) # 加载DESeq2包
countData <- raw_count_filt2 # 表达矩阵
condition <- factor(c("control","control","KD","KD")) # 定义condition
colData <- data.frame(row.names=colnames(countData), condition) # 样品信息矩阵

2、构建dds矩阵并进行差异分析

dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition ) # 构建dds矩阵
head(dds) # 查看dds矩阵的前6行
dds2 <- DESeq(dds) # 对dds进行Normalize
resultsNames（dds) # 查看结果的名称
res <- results(dds) # 使用results()函数获取结果，并赋值给res
head (res, n=5) # 查看res矩阵的前5行
mcols(res,use.names= TRUE) # 查看res矩阵每一列的含义
summary(res) # 对res矩阵进行总结

3、提取差异分析结果

table(res$padj<0.05) # 取padj小于0.05的数据，得到743行
res <- res[order(res$padj),] # 按照padj的大小将res重新排列
diff_gene_deseq2 <- subset(res,padj < 0.05 & (log2FoldChange >1 | log2FoldChange < -1)) # 获取padj小于0.05，表达倍数取以2为对数后绝对值大于1的差异表达基因，赋值给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") # 将结果写入control_vs_akap95.csv文件

示例2

1、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)

2、前过滤

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

3、标注因子水平

dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
dds$condition <- relevel(dds$condition, ref = "untreated")
dds$condition <- droplevels(dds$condition)

注：在进行后续操作前指定，如果没有指定，采用的是默认

levels对应的分别为，分母，分子

ref为指定的对照组，即分母

当公式中的某个变量对应的样本没有的时候，可以通过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）任务并行

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）独立假设权重

library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult

6、探索并导出结果

1）MA-plot

plotMA(res, ylim=c(-2,2))
plotMA(resLFC, ylim=c(-2,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")

3）counts绘图

plotCounts(dds, gene=which.min(res$padj), intgroup="condition")

4）结果的更多信息

mcols(res)$description

可以通过mcol函数来查看结果表中各个变量的含义。其中变量LFC就是变化的倍数，存在NA的可能原因有：

1、一个基因在所有样本中的表达值都为0
2、一个基因的某个样本的表达值为离群值（通过Cook距离检测）
3、一个基因counts标准化后数值过低
后续：结果的导出和丰富的可视化
ReportingTools，regionReport，Glimma，pcaExplorer等工具可以生成一个报表形式的结果。