RNA-seq入门实战（五）：差异分析——DESeq2 edgeR limma的使用与比较

连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中几个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！

他前面的分享是：

下面是他对我们b站转录组视频课程的详细笔记

本节概览：

1.DESeq2、 edgeR、limma的使用

2.三类差异分析软件的结果比较——相关性、韦恩图

3.选取差异基因绘制火山图和热图

一、DESeq2、 edgeR、limma的使用

强烈建议查看官方说明书进行这三种差异分析的学习，链接在文章末尾给出。
注意，这三个包都需要输入counts进行分析，不能用tpm、fpkm等归一化后的数据。

承接上节 RNA-seq入门实战（三）：在R里面整理表达量counts矩阵 和  RNA-seq入门实战（二）：上游数据的比对计数——Hisat2+ featureCounts 与 Salmon  正式分析前先进行目录设置、实验组和对照组的指定：

rm(list = ls())  
options(stringsAsFactors = F)

setwd("C:/Users/Lenovo/Desktop/test")
load(file = '1.counts.Rdata')
dir.create("3.DEG")
setwd("3.DEG")

##设定 实验组exp / 对照组ctr
exp="primed"
ctr="naive"

1. DESeq2

DESeq2是目前最常用的差异分析R包。除了可以导入counts外，如果上游使用salmon，DESeq2官方还给出了直接导入tximport生成的txi对象的方法。counts与txi的获取见  RNA-seq入门实战（三）：在R里面整理表达量counts矩阵 和  RNA-seq入门实战（二）：上游数据的比对计数——Hisat2+ featureCounts 与 Salmon   

library(DESeq2)
library("BiocParallel") #启用多核计算

##构建dds DESeqDataSet
if(T){
    dds <- DESeqDataSetFromMatrix(countData = counts,
                                  colData = gl,
                                  design = ~ group_list)
    }
if(F){  #若上游为salmon
    dds <- DESeqDataSetFromTximport(txi, 
                                    colData = gl,
                                    design = ~ group_list)
    }


dds$group_list <- relevel(dds$group_list, ref = ctr)   #指定 control group

keep <- rowSums(counts(dds)) >= 1.5*ncol(counts)  #Pre-filtering ，过滤低表达基因
dds <- dds[keep,] 
dds <- DESeq(dds,quiet = F) 
res <- results(dds,contrast=c("group_list", exp, ctr))  #指定提取为exp/ctr结果
resOrdered <- res[order(res$padj),]  #order根据padj从小到大排序结果
tempDEG <- as.data.frame(resOrdered)
DEG_DEseq2 <- na.omit(tempDEG)

2. edgeR

使用EdgeR时注意选择合适的分析算法，官方建议bulk RNA-seq选择quasi-likelihood(QL) F-test tests，scRNA-seq 或是没有重复样品的数据选用 likelihood ratio test。

library(edgeR)  #install.packages("statmod")
library(statmod)

#分组矩阵design构建
group <- factor(group_list)
group <- relevel(group,ctr)     #将对照组的因子设置为1
design <- model.matrix(~0+group)
rownames(design) <- colnames(counts)
colnames(design) <- levels(group)

## 表达矩阵DGEList构建与过滤低表达基因
dge <- DGEList(counts=counts, group=group)
keep.exprs <- filterByExpr(dge) #自动筛选过滤低表达基因
dge <- dge[keep.exprs,,keep.lib.sizes=FALSE] 
dge <- calcNormFactors(dge, method = 'TMM') #归一化因子用于 normalizes the library sizes
dge <- estimateDisp(dge, design, robust=T) 

## To perform quasi-likelihood(QL) F-test tests:  bulk RNA-seq 
fit <- glmQLFit(dge, design, robust=T)  #拟合模型 
lt <- glmQLFTest(fit, contrast=c(-1,1)) #统计检验#注意比对顺序:实验-1 /对照1

## To perform likelihood ratio test：scRNA-seq and no replicates data
# fit <- glmFit(dge, design, robust=T)
# lt <- glmLRT(fit, contrast=c(-1,1))  #比对:顺序实验/对照，已设对照为1

tempDEG <- topTags(lt, n = Inf) #sort by PValue, n is the number of genes/tags to return
tempDEG <- as.data.frame(tempDEG)
DEG_edgeR <- na.omit(tempDEG)

3. limma

limma进行差异分析有两种方法 ：limma-trend和voom，在样本测序深度相差不大时两种方法差距不大，而测序深度相差大时voom更有优势，因此我们一般都选择voom的方法进行差异分析。Limma-voom vs limma-trend (bioconductor.org)

library(limma)
library(edgeR)

#分组矩阵design构建
design <- model.matrix(~0+factor(group_list)) #构建分组矩阵
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(counts)

## 表达矩阵DGEList构建与过滤低表达基因
dge <- DGEList(counts=counts) 
keep.exprs <- filterByExpr(dge,design=design) #过滤低表达基因
dge <- dge[keep.exprs,,keep.lib.sizes=FALSE] 
dge <- calcNormFactors(dge)  #归一化基因表达分布,得到的归一化系数被用作文库大小的缩放系数
cont.matrix <- makeContrasts(contrasts=paste0(exp,'-',ctr), #比对顺序实验/对照
                             levels = design)

## DE分析  limma-trend（logCPM,有相似文库大小）  or  voom（文库大小差异大）
# de <- cpm(dge, log=TRUE, prior.count=3)  #如选择logCPM，则eBayes设trend=TRUE
de <- voom(dge,design,plot=TRUE, normalize="quantile")
fit1 <- lmFit(de, design)               #线性拟合
fit2 <- contrasts.fit(fit1,cont.matrix) #统计检验
efit <- eBayes(fit2, trend=F)  #Apply empirical Bayes smoothing to the standard errors

tempDEG <- topTable(efit, coef=paste0(exp,'-',ctr), n=Inf)  #padj值从小到大排列
DEG_limma_voom  <- na.omit(tempDEG)

4. 保存DEG数据

#### 保存DEG数据 ####
save(DEG_limma_voom,DEG_DEseq2,DEG_edgeR,
       file ='test_DEG_results.Rdata') 
  
  #### 合并三种DEG文件交集至csv
  allg <- intersect(rownames(DEG_limma_voom),rownames(DEG_edgeR))#取交集
  allg <- intersect(allg,rownames(DEG_DEseq2))
  ALL_DEG <- cbind(DEG_limma_voom[allg,c(1,4,5)],
                                 DEG_edgeR[allg,c(1,4,5)],
                                 DEG_DEseq2[allg,c(2,5,6)]) 
  colnames(ALL_DEG)
  colnames(ALL_DEG) <- c('limma_log2FC','limma_pvalue','limma_padj',
                         'edgeR_log2FC','edgeR_pvalue','edgeR_padj',
                         'DEseq2_log2FC','DEseq2_pvalue','DEseq2_padj')
  write.csv(ALL_DEG,file = 'test_DEG_results.csv')

---

二、3种差异分析结果比较

由于本次样品两组间差异十分显著，差异基因很多，因此筛选阈值调整为：FoldChang=10，padj=0.001。一般情况下选择FoldChang=1.5~4，padj<=0.05即可，根据样本情况而定。
下面查看三种差异分析结果的相关性和差异基因的重叠情况。

#### 比较一下三种DEG方法结果 ####
load(list.files(pattern = 'DEG_results.Rdata'))
a <- read.csv(file = list.files(pattern = 'DEG_results.csv'),header = T,row.names = 1)

#查看相关性、一致性
colnames(a)
print(cor(a[,c(1,4,7)])) 

#查看显著差异基因重叠性，绘制韦恩图
#BiocManager::install("RBGL") #安装依赖包
#install.packages("Vennerable", repos="http://R-Forge.R-project.org") #安装Vennerable包
library(Vennerable) 

log2FC_cutoff=log2(10);  p_cutoff=0.001   #筛选显著差异基因比较
if(T){#根据Padj筛选
DEseq2_deg <- rownames(DEG_DEseq2[with(DEG_DEseq2,abs(log2FoldChange)>log2FC_cutoff & padj<p_cutoff),])
edgeR_deg <- rownames(DEG_edgeR[with(DEG_edgeR,abs(logFC)>log2FC_cutoff & FDR<p_cutoff),])
limma_deg <- rownames(DEG_limma_voom[with(DEG_limma_voom,abs(logFC)>log2FC_cutoff & adj.P.Val<p_cutoff),])
}

mylist <- list(DEseq2=DEseq2_deg, edgeR=edgeR_deg, limma=limma_deg)
str(mylist)
Vennplot <- Venn(mylist)
Vennplot1 <- Vennplot[,c('DEseq2','edgeR')]
Vennplot2 <- Vennplot[,c('DEseq2','limma')]
Vennplot3 <- Vennplot[,c('limma','edgeR')]

pdf(file = paste0('3DEG_Vennplot_lg2FC',log2FC_cutoff,'.pdf'))
plot(Vennplot, doWeights = F)
plot(Vennplot, doWeights = T) #doWeights=T设置为按数量比例绘图
plot(Vennplot1, doWeights = T)
plot(Vennplot2, doWeights = T)
plot(Vennplot3, doWeights = T)
dev.off()

从以下图中可以看出，三种方法所得结果相关性很高，都在0.99以上，显著差异基因的重叠也很高。（所以一般来说大家无需纠结使用哪种方法，都是认可的）

---

三、选取差异基因绘制火山图和热图

以下示范选取DESeq2差异分析结果进行绘制, 筛选阈值设置为：FoldChang=10，padj=0.001

1. 火山图的绘制

library(ggplot2)
library(pheatmap)
##筛选条件设置
log2FC_cutoff = log2(10)
padj_cutoff = 0.001
##选取差异分析结果
need_DEG <- DEG_DEseq2[,c(2,6)] #选取log2FoldChange, padj信息
colnames(need_DEG) <- c('log2FoldChange','padj') 

need_DEG$significance  <- as.factor(ifelse(need_DEG$padj < padj_cutoff & abs(need_DEG$log2FoldChange) > log2FC_cutoff,
                                           ifelse(need_DEG$log2FoldChange > log2FC_cutoff ,'UP','DOWN'),'NOT'))

title <- paste0(' Up :  ',nrow(need_DEG[need_DEG$significance =='UP',]) ,
                     '\n Down : ',nrow(need_DEG[need_DEG$significance =='DOWN',]),
                     '\n FoldChange >= ',round(2^log2FC_cutoff,3))

g <- ggplot(data=need_DEG, 
            aes(x=log2FoldChange, y=-log10(padj), 
                color=significance)) +
  #点和背景
  geom_point(alpha=0.4, size=1) +
  theme_classic()+ #无网格线
  #坐标轴
  xlab("log2 ( FoldChange )") + 
  ylab("-log10 ( P.adjust )") +
  #标题文本
  ggtitle( title ) +
  #分区颜色                  
  scale_colour_manual(values = c('blue','grey','red'))+ 
  #辅助线
  geom_vline(xintercept = c(-log2FC_cutoff,log2FC_cutoff),lty=4,col="grey",lwd=0.8) +
  geom_hline(yintercept = -log10(padj_cutoff),lty=4,col="grey",lwd=0.8) +
  #图例标题间距等设置
  theme(plot.title = element_text(hjust = 0.5), 
        plot.margin=unit(c(2,2,2,2),'lines'), #上右下左
        legend.title = element_blank(), #不显示图例标题
        legend.position="right")  #图例位置

ggsave(g,filename = 'volcano_padj.pdf',width =8,height =7.5)

2. 热图的绘制

##选择要展示基因表达量的数据
# dat <- log2(edgeR::cpm(counts)+1) 
dat <- log2(tpm+1)
# dat <- read.table("../2.check/Deseq2_rld.txt"); colnames(dat) <- rownames(gl)  #R读取数据列名可能会出错，需要重新对应一下

gene_up <- rownames(need_DEG[with(need_DEG,log2FoldChange>log2FC_cutoff & padj<padj_cutoff),])
gene_down <- rownames(need_DEG[with(need_DEG,log2FoldChange< -log2FC_cutoff & padj<padj_cutoff),])
cg <- c(head(gene_up, 50),   #取前50 padj上下调基因名
        head(gene_down, 50))
cg <- na.omit(match(cg,rownames(dat))) 

pheatmap::pheatmap(dat[cg,],
                   show_colnames =T,
                   show_rownames = F,
                   #scale = "row",
                   fontsize = 7 ,
                   cluster_cols = T,
                   annotation_col=gl,
                   filename = 'heatmap_top50up&down_DEG.pdf')

---

到此就完成了基因差异分析的基本过程，得到了不同分组间的差异基因相关信息，接下来就要对差异基因进行富集分析啦。

目前简单的差异分析流程，基本上转录组测序技术和芯片技术拿到的表达量矩阵后续分析大同小异，公众号推文在：

参考资料

三种R包的官方说明书：
Analyzing RNA-seq data with DESeq2 (bioconductor.org)

edgeR: differential analysis of sequence read count data User's Guide (bioconductor.org)

limma usersguide.pdf (bioconductor.org)

bioconductor的worlk flow：
使用limma、Glimma和edgeR，RNA-seq数据分析易如反掌 (bioconductor.org)

部分代码修改参考自：
GitHub - jmzeng1314/GEO

【生信技能树】转录组测序数据分析_哔哩哔哩_bilibili

【生信技能树】GEO数据库挖掘_哔哩哔哩_bilibili

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年5~6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年5~6月场，你的生物信息学第一课