本節(jié)概覽：
1.DESeq2打月、 edgeR揭朝、limma的使用
2.三類差異分析軟件的結(jié)果比較——相關(guān)性酸员、韋恩圖
3.選取差異基因繪制火山圖和熱圖

承接前期文章：RNA-seq入門實(shí)戰(zhàn)（四）：差異分析前的準(zhǔn)備——數(shù)據(jù)檢查

一、DESeq2续室、 edgeR、limma的使用

強(qiáng)烈建議查看官方說(shuō)明書進(jìn)行這三種差異分析的學(xué)習(xí)谒养，鏈接在文章末尾給出挺狰。
注意，這三個(gè)包都需要輸入counts進(jìn)行分析，不能用tpm她渴、fpkm等歸一化后的數(shù)據(jù)达址。
正式分析前先進(jìn)行目錄設(shè)置、實(shí)驗(yàn)組和對(duì)照組的指定：

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

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

##設(shè)定 實(shí)驗(yàn)組exp / 對(duì)照組ctr
exp="primed"
ctr="naive"

1. DESeq2

DESeq2是目前最常用的差異分析R包趁耗。除了可以導(dǎo)入counts外沉唠，如果上游使用salmon，DESeq2官方還給出了直接導(dǎo)入tximport生成的txi對(duì)象的方法苛败。counts與txi的獲取見(jiàn)RNA-seq入門的簡(jiǎn)單實(shí)戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣

library(DESeq2)
library("BiocParallel") #啟用多核計(jì)算

##構(gòu)建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 满葛，過(guò)濾低表達(dá)基因
dds <- dds[keep,] 
dds <- DESeq(dds,quiet = F) 
res <- results(dds,contrast=c("group_list", exp, ctr))  #指定提取為exp/ctr結(jié)果
resOrdered <- res[order(res$padj),]  #order根據(jù)padj從小到大排序結(jié)果
tempDEG <- as.data.frame(resOrdered)
DEG_DEseq2 <- na.omit(tempDEG)

2. edgeR

使用EdgeR時(shí)注意選擇合適的分析算法，官方建議bulk RNA-seq選擇quasi-likelihood(QL) F-test tests罢屈，scRNA-seq 或是沒(méi)有重復(fù)樣品的數(shù)據(jù)選用 likelihood ratio test嘀韧。

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

#分組矩陣design構(gòu)建
group <- factor(group_list)
group <- relevel(group,ctr)     #將對(duì)照組的因子設(shè)置為1
design <- model.matrix(~0+group)
rownames(design) <- colnames(counts)
colnames(design) <- levels(group)

## 表達(dá)矩陣DGEList構(gòu)建與過(guò)濾低表達(dá)基因
dge <- DGEList(counts=counts, group=group)
keep.exprs <- filterByExpr(dge) #自動(dòng)篩選過(guò)濾低表達(dá)基因
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)) #統(tǒng)計(jì)檢驗(yàn)#注意比對(duì)順序:實(shí)驗(yàn)-1 /對(duì)照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))  #比對(duì):順序?qū)嶒?yàn)/對(duì)照，已設(shè)對(duì)照為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進(jìn)行差異分析有兩種方法 ：limma-trend和voom缠捌，在樣本測(cè)序深度相差不大時(shí)兩種方法差距不大锄贷，而測(cè)序深度相差大時(shí)voom更有優(yōu)勢(shì)，因此我們一般都選擇voom的方法進(jìn)行差異分析曼月。Limma-voom vs limma-trend (bioconductor.org)

library(limma)
library(edgeR)

#分組矩陣design構(gòu)建
design <- model.matrix(~0+factor(group_list)) #構(gòu)建分組矩陣
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(counts)

## 表達(dá)矩陣DGEList構(gòu)建與過(guò)濾低表達(dá)基因
dge <- DGEList(counts=counts) 
keep.exprs <- filterByExpr(dge,design=design) #過(guò)濾低表達(dá)基因
dge <- dge[keep.exprs,,keep.lib.sizes=FALSE] 
dge <- calcNormFactors(dge)  #歸一化基因表達(dá)分布,得到的歸一化系數(shù)被用作文庫(kù)大小的縮放系數(shù)
cont.matrix <- makeContrasts(contrasts=paste0(exp,'-',ctr), #比對(duì)順序?qū)嶒?yàn)/對(duì)照
                             levels = design)

## DE分析  limma-trend（logCPM,有相似文庫(kù)大幸耆础）  or  voom（文庫(kù)大小差異大）
# de <- cpm(dge, log=TRUE, prior.count=3)  #如選擇logCPM，則eBayes設(shè)trend=TRUE
de <- voom(dge,design,plot=TRUE, normalize="quantile")
fit1 <- lmFit(de, design)               #線性擬合
fit2 <- contrasts.fit(fit1,cont.matrix) #統(tǒng)計(jì)檢驗(yàn)
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數(shù)據(jù)

#### 保存DEG數(shù)據(jù) ####
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種差異分析結(jié)果比較

由于本次樣品兩組間差異十分顯著炎辨，差異基因很多，因此篩選閾值調(diào)整為：FoldChang=10聪姿，padj=0.001碴萧。一般情況下選擇FoldChang=1.5~4，padj<=0.05即可末购，根據(jù)樣本情況而定破喻。
下面查看三種差異分析結(jié)果的相關(guān)性和差異基因的重疊情況。

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

#查看相關(guān)性盟榴、一致性
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){#根據(jù)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設(shè)置為按數(shù)量比例繪圖
plot(Vennplot1, doWeights = T)
plot(Vennplot2, doWeights = T)
plot(Vennplot3, doWeights = T)
dev.off()

從以下圖中可以看出，三種方法所得結(jié)果相關(guān)性很高曹货，都在0.99以上咆繁，顯著差異基因的重疊也很高。

三顶籽、選取差異基因繪制火山圖和熱圖

以下示范選取DESeq2差異分析結(jié)果進(jìn)行繪制, 篩選閾值設(shè)置為：FoldChang=10玩般，padj=0.001

1. 火山圖的繪制

library(ggplot2)
library(pheatmap)
##篩選條件設(shè)置
log2FC_cutoff = log2(10)
padj_cutoff = 0.001
##選取差異分析結(jié)果
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)) +
  #點(diǎn)和背景
  geom_point(alpha=0.4, size=1) +
  theme_classic()+ #無(wú)網(wǎng)格線
  #坐標(biāo)軸
  xlab("log2 ( FoldChange )") + 
  ylab("-log10 ( P.adjust )") +
  #標(biāo)題文本
  ggtitle( title ) +
  #分區(qū)顏色                  
  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) +
  #圖例標(biāo)題間距等設(shè)置
  theme(plot.title = element_text(hjust = 0.5), 
        plot.margin=unit(c(2,2,2,2),'lines'), #上右下左
        legend.title = element_blank(), #不顯示圖例標(biāo)題
        legend.position="right")  #圖例位置

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

2. 熱圖的繪制

##選擇要展示基因表達(dá)量的數(shù)據(jù)
# dat <- log2(edgeR::cpm(counts)+1) 
dat <- log2(tpm+1)
# dat <- read.table("../2.check/Deseq2_rld.txt"); colnames(dat) <- rownames(gl)  #R讀取數(shù)據(jù)列名可能會(huì)出錯(cuò)，需要重新對(duì)應(yīng)一下

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上下調(diào)基因名
        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')

到此就完成了基因差異分析的基本過(guò)程礼饱，得到了不同分組間的差異基因相關(guān)信息坏为，接下來(lái)就要對(duì)差異基因進(jìn)行富集分析啦

參考資料
三種R包的官方說(shuō)明書：
RNA-seq workflow: gene-level exploratory analysis and differential expression (bioconductor.org)
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數(shù)據(jù)分析易如反掌 (bioconductor.org)
部分代碼修改參考自：
GitHub - jmzeng1314/GEO
【生信技能樹】轉(zhuǎn)錄組測(cè)序數(shù)據(jù)分析_嗶哩嗶哩_bilibili
【生信技能樹】GEO數(shù)據(jù)庫(kù)挖掘_嗶哩嗶哩_bilibili

RNA-seq實(shí)戰(zhàn)系列文章：
RNA-seq入門實(shí)戰(zhàn)（零）：RNA-seq流程前的準(zhǔn)備——Linux與R的環(huán)境創(chuàng)建
RNA-seq入門實(shí)戰(zhàn)（一）：上游數(shù)據(jù)下載匀伏、格式轉(zhuǎn)化和質(zhì)控清洗
RNA-seq入門實(shí)戰(zhàn)（二）：上游數(shù)據(jù)的比對(duì)計(jì)數(shù)——Hisat2+ featureCounts 與 Salmon
RNA-seq入門實(shí)戰(zhàn)（三）：從featureCounts與Salmon輸出文件獲取counts矩陣
RNA-seq入門實(shí)戰(zhàn)（四）：差異分析前的準(zhǔn)備——數(shù)據(jù)檢查
RNA-seq入門實(shí)戰(zhàn)（五）：差異分析——DESeq2 edgeR limma的使用與比較
RNA-seq入門實(shí)戰(zhàn)（六）：GO洒忧、KEGG富集分析與enrichplot超全可視化攻略
RNA-seq入門實(shí)戰(zhàn)（七）：GSEA——基因集富集分析
RNA-seq入門實(shí)戰(zhàn)（八）：GSVA——基因集變異分析
RNA-seq入門實(shí)戰(zhàn)（九）：PPI蛋白互作網(wǎng)絡(luò)構(gòu)建（上）——STRING數(shù)據(jù)庫(kù)的使用
RNA-seq入門實(shí)戰(zhàn)（十）：PPI蛋白互作網(wǎng)絡(luò)構(gòu)建（下）——Cytoscape軟件的使用
RNA-seq入門實(shí)戰(zhàn)（十一）：WGCNA加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析——關(guān)聯(lián)基因模塊與表型