2.2转录组 | 差异表达分析（DESeq,edgeR,limm

3种不同软件用于差异表达分析：DESeq2,edgeR,limma

套路总结：

①导入read count， 保存为专门的对象用于后续分析

②原始数据过滤，根据标准化read count 或者 raw count 作为筛选标准

③raw read count 标准化

④通过各种算法（如经验贝叶斯，EM）预测dispersion离散值

⑤广义线性模型拟合数据

⑥差异分析，也就是统计检验部分

---

导入数据

该部分和2.1转录组 | 差异表达分析（DESeq）内容相同

rm(list = ls())
options(stringsAsFactors = T)
library(DESeq2)
control <- read.table("SRR3589956_matrix.count",
                      sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table("SRR3589957_matrix.count",
                   sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table("SRR3589958_matrix.count",
                   sep="\t",col.names = c("gene_id","rep2"))
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL

delta_mean <- abs(mean(raw_count_filt$rep1) - mean(raw_count_filt$rep2))

sampleNum <- length(raw_count_filt$control)
sampleMean <- mean(raw_count_filt$control)
control2 <- integer(sampleNum)

for (i in 1:sampleNum){
  if(raw_count_filt$control[i] < sampleMean){
    control2[i] <- raw_count_filt$control[i] + abs(raw_count_filt$rep1[i] - raw_count_filt$rep2[i])
  }
  else{
    control2[i] <- raw_count_filt$control[i] + rpois(1,delta_mean)
  }
}


raw_count_filt$control2 <- control2

1.使用DESeq2进行差异基因分析（这部分与2.1转录组 | 差异表达分析（DESeq）

内容相同）
DESeq2要求的数据是raw count， 没必要进行FPKM/TPM/RPFKM/TMM标准化

#构建DESeq2所需的DESeqDataSet对象
library(DESeq2)
countData <- raw_count_filt[,2:5]    #countData为基因表达矩阵
condition <- factor(c("control","KD","KD","control"))     #condition为处理条件
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )

#过滤low_count的数据
nrow(dds)
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds)  #看下图，过滤了将近一半的数据

#使用DESeq进行差异表达分析
dds <- DESeq(dds)        #出现如下红色字符，表示运行成功

#查看结果
res <- results(dds)
mcols(res, use.names = TRUE)
summary(res)

##提取差异基因分析的结果
table(res$padj<0.05)     #计算p值小于0.05的基因的个数
res_deseq <- res[order(res$padj),]   #根据res的padj值进行排序，并赋值给res_deseq
diff_gene_deseq2 <- subset(res_deseq, padj<0.05 & (log2FoldChange > 1 | log2FoldChange < -1))    #将res_deseq过滤：要求是p值小于0.05且log2FoldChange在1~-1范围内的

diff_gene_deseq2 <- row.names(diff_gene_deseq2)
res_diff_data <- merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
write.csv(res_diff_data,file = "11_30_humen_data.csv",row.names = F)

2.使用edgeR进行差异基因分析

（1）构建DGElist对象
edgeR使用DGEList函数读取count matrix数据， 即提供一个现成的matrix数据（ 也可以利用 tximport 或 DESeqDataSetFromHTSeq 读取单独的文件，然后传递给DGEList）

library(edgeR)
group <- factor(c("control","KD","KD","control"))
# genelist <- DGEList(counts=raw_count_filt[,2:5], group = group)
counts=raw_count_filt[,2:5]
genelist <- DGEList(counts, group = group)

（2）过滤 low_counts数据
根据经验，基因至少在某一些文库的count超过10~15才被认为是表达。（但也可以根据实际过滤的情况，做相应的调整）。①可以简单粗暴的对每个基因的raw count进行比较，②但建议先用CPM标准化后再比较，避免了文库大小的影响。

# 简单粗暴的方法
keep <- rowSums(genelist$count) > 50
# 利用CPM标准化
keep <- rowSums(cpm(genelist) > 0.5 ) >=2     #这里0.5（即阈值）等于10。
table(keep)
genelist.filted <- genelist[keep, ,keep.lib.sizes=FALSE]    #keep.lib.size=FALSE表示重新计算文库大小。

（3）根据组成偏好(composition bias)标准化

#edgeR的calcNormFactors函数使用TMM算法对DGEList标准化
genelist.norm <- calcNormFactors(genelist.filted)

注：大部分的mRNA-Seq数据分析用TMM标准化就行了，但是也有例外，比如说single-cell RNA-Seq(Lun, Bach, and Marioni 2016), 还有就是global differential expression， 基因组一半以上的基因都是差异表达的，请尽力避免，(D. Wu et al. 2013)， 不然就需要用到内参进行标准化了(Risso et al. 2014).

（4）实验设计矩阵(Design matrix)
edgeR的线性模型和差异表达分析需要定义一个实验设计矩阵，类似于DESeq2中的design参数。 很直白的就能发现是1vs0

design <- model.matrix(~0+group)
colnames(design) <- levels(group)
design

image.png

（5）估计离散值（Dispersion）
edgeR对每个基因都估测了一个经验贝叶斯稳健离散值、一个公共离散值、一个趋势离散值

genelist.Disp <- estimateDisp(genelist.norm, design, robust = TRUE)
plotBCV(genelist.Disp)

plotBCV

通过quasi-likelihood (QL)拟合NB模型，用于解释生物学和技术性导致的基因特异性变异

fit <- glmQLFit(genelist.Disp, design, robust=TRUE)
head(fit$coefficients)

注意：①估计离散值这个步骤其实有许多estimate*Disp函数。当不存在实验设计矩阵(design matrix)的时候，estimateDisp 等价于 estimateCommonDisp 和 estimateTagwiseDisp 。而当给定实验设计矩阵(design matrix)时， estimateDisp 等价于 estimateGLMCommonDisp, estimateGLMTrendedDisp 和 estimateGLMTagwiseDisp。 其中tag与gene同义。
注意：② 这里的第三， 四， 五步对应的就是DESeq2的DESeq包含的2步，标准化和离散值估测。
（6） 差异表达检验之一
这一步主要构建比较矩阵，类似于DESeq2中的results函数的 contrast 参数。

cntr.vs.KD <- makeContrasts(control-KD, levels=design)
res <- glmQLFTest(fit, contrast=cntr.vs.KD)
ig.edger <- res$table[p.adjust(res$table$PValue, method = "BH") < 0.01, ]

注意： 这里用的是glmQLFTest而不是glmLRT是因为前面用了glmQLTFit进行拟合，所以需要用QL F-test进行检验。如果前面用的是glmFit，那么对应的就是glmLRT. 作者称QL F-test更加严格。多重试验矫正用的也是BH方法。

（7） 提取显著性差异的基因用作下游分析
plotMD图

topTags(res,n=10)
is.de <- decideTestsDGE(res)
summary(is.de)
plotMD(res, status=is.de, values=c(1,-1), col=c("red","blue"),
       legend="topright")

plotMD

（8）差异表达检验之二
前面“差异表达检验之一”找到的显著性差异的基因，没有考虑到效应值（即具体变化了多少倍）。我们也可以用找到的表达量变化比较大的基因，对应的函数是 glmTreat 。

tr <- glmTreat(fit, contrast=cntr.vs.KD, lfc=log2(1.5))
ig.edger <- res$table[p.adjust(res$table$PValue, method = "BH") < 0.01, ]
topTags(res,n=10)
is.de <- decideTestsDGE(res)
summary(is.de)
plotMD(res, status=is.de, values=c(1,-1), col=c("red","blue"),
       legend="topright")

plotMD

3.使用limma进行差异基因分析

Limma采用经验贝叶斯模型（ Empirical Bayesian model）让结果更稳健。
在处理RNA-Seq数据时，raw read count先被转成log2-counts-per-million (logCPM)，然后对mean-variance关系建模。建模有两种方法：①精确权重法，即“voom"；② 经验贝叶斯先验趋势，即 “limma-trend”。
（1）数据预处理
Limma使用edgeR的DGEList对象，并且过滤方法都是一致的，对应edgeR的（1）、（2）、（3）

library(edgeR)library(limma)group <- factor(c("control","KD","KD","control"))genelist <- DGEList(counts=raw_count_filt[,2:5], group = group)### filter base  use CPMkeep <- rowSums(cpm(genelist) > 0.5 ) >=2table(keep)genelist.filted <- genelist[keep, ,keep.lib.sizes=FALSE]### normalizaitionx <- calcNormFactors(genelist.filted, method = "TMM")

（2）差异表达分析之一：使用“limma-trend”

design <- model.matrix(~0+group)
colnames(design) <- levels(group)
logCPM <- cpm(genelist.norm, log=TRUE, prior.count=3)
fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE)
topTable(fit, coef=ncol(design))

（3）差异表达分析之二：使用“voom”

### DGE with voom
v <- voom(genelist.norm, design, plot=TRUE)
#v <- voom(counts, design, plot=TRUE)

image.png

fit <- lmFit(v, design)
fit <- eBayes(fit)
all <- topTable(fit, coef=ncol(design), number=10000)
sig.limma <- all[all$adj.P.Val < 0.01, ]
fit <- treat(fit, lfc=log2(1.2))
topTreat(fit, coef=ncol(design))

---

对以上不同软件（DESeq,edgeR,limma）做出的差异表达分析的结果进行比较

可以用维恩图，也可以用UpSetR

sig.deseq2 <- res_deseq
sig.edger <- ig.edger
sig.limma <- all[all$adj.P.Val < 0.01, ]

##韦恩图
library(VennDiagram)
deseq2_num <- row.names(sig.deseq2)
sedger_num <- row.names(sig.edger)
limma_num  <- row.names(sig.limma)
#三维韦恩图会输出得到3triple_Venn.tiff
venn.plot <- venn.diagram(
  x = list(
    DESeq2 = deseq2_num,
    edgeR = sedger_num,
    limma = limma_num
  ),
  filename = "3triple_Venn.tiff",
  col = "transparent",
  fill = c("red", "blue", "green"),
  alpha = 0.5,
  label.col = c("darkred", "white", "darkblue", "white",
                "white", "white", "darkgreen"),
  cex = 2.5,
  fontfamily = "serif",
  fontface = "bold",
  cat.default.pos = "text",
  cat.col = c("darkred", "darkblue", "darkgreen"),
  cat.cex = 2.5,
  cat.fontfamily = "serif",
  cat.dist = c(0.06, 0.06, 0.03),
  cat.pos = 0
)

##UpSetR
library(UpSetR)
input <- fromList(list(edgeR=rownames(sig.edger), DESeq2=rownames(sig.deseq2), limma=rownames(sig.limma)))
upset(input)

Venn UpSetR

edgeR和limma过程中，有些步骤没看明白。最后对3个软件找到的差异表达的基因数量构出的韦恩图，可以看到差异很大，尤其是DESeq，应该是我最后输入的3个数据集有问题。待修改~

---

参考：
01 转录组入门（7）：差异表达分析
02 转录组学习七（差异基因分析）