基于轉(zhuǎn)錄組測(cè)序獲得的定量表達(dá)值但指，識(shí)別差異表達(dá)變化的基因或其它非編碼RNA分子谎柄，實(shí)際上方法還是非常多的。但就目前來看朋鞍，DESeq2和edgeR是出現(xiàn)頻率最高的兩種方法了已添。

DESeq2已經(jīng)在上一篇文章中作了簡(jiǎn)介，本篇繼續(xù)展示R包edgeR的差異基因分析流程滥酥。類似DESeq2更舞，edgeR作為被廣泛使用的R包，文獻(xiàn)中經(jīng)常能看到它的身影坎吻，如下舉例缆蝉。

相關(guān)文獻(xiàn)描述

1 安裝edgeR

edgeR可直接使用Bioconductor安裝，還是非常簡(jiǎn)單的禾怠。

#Bioconductor 安裝 edgeR
#install.packages('BiocManager')  #需要首先安裝 BiocManager返奉，如果尚未安裝請(qǐng)先執(zhí)行該步
BiocManager::install('edgeR')

2 使用edgeR鑒定差異表達(dá)基因

edgeR使用經(jīng)驗(yàn)貝葉斯估計(jì)和基于負(fù)二項(xiàng)模型的精確檢驗(yàn)來確定差異基因，通過在基因之間來調(diào)節(jié)跨基因的過度離散程度吗氏，使用類似于Fisher精確檢驗(yàn)但適應(yīng)過度分散數(shù)據(jù)的精確檢驗(yàn)用于評(píng)估每個(gè)基因的差異表達(dá)。

以下是edgeR分析差異表達(dá)基因的一般過程雷逆。

2.1 準(zhǔn)備數(shù)據(jù)和文件讀取

首先準(zhǔn)備基因表達(dá)值矩陣弦讽。
本文的所有測(cè)試數(shù)據(jù)和R代碼，可在文末獲取

• “control_treat.count.txt”，是6個(gè)測(cè)序樣本的基因表達(dá)值矩陣往产，包括3個(gè)處理組（treat）和3個(gè)對(duì)照組（control）被碗，需注意將對(duì)照組在前，處理組在后仿村。第1列是基因名稱锐朴，注意不能有重復(fù)值。

  
  
  示例文件樣式蔼囊，行是基因焚志，列是樣本，內(nèi)容為基因表達(dá)值
  
  

  然后將準(zhǔn)備好的基因表達(dá)值矩陣讀到R中畏鼓，并預(yù)定義樣本分組信息酱酬。

#讀取基因表達(dá)矩陣
targets <- read.delim('control_treat.count.txt', row.names = 1, sep = '\t', check.names = FALSE)

#指定分組，注意要保證表達(dá)矩陣中的樣本順序和這里的分組順序是一一對(duì)應(yīng)的
#對(duì)照組在前云矫，處理組在后
group <- rep(c('control', 'treat'), each = 3)

2.2 表達(dá)值預(yù)處理和標(biāo)準(zhǔn)化

在執(zhí)行差異表達(dá)基因分析前膳沽，將輸入的基因表達(dá)矩陣和分組信息構(gòu)建edgeR的DGEList對(duì)象，便于儲(chǔ)存數(shù)據(jù)和中間運(yùn)算让禀；并對(duì)表達(dá)值標(biāo)準(zhǔn)化挑社，以消除由于樣品制備或建庫測(cè)序過程中帶來的影響。

同時(shí)巡揍，需要手動(dòng)過濾一些低表達(dá)量的基因滔灶。可能會(huì)有同學(xué)疑問吼肥，為什么上一篇DESeq2中無需手動(dòng)過濾呢录平？這是因?yàn)镈ESeq2設(shè)計(jì)之初就將自動(dòng)過濾步驟封裝在函數(shù)中，可以實(shí)現(xiàn)對(duì)低表達(dá)基因的自動(dòng)處理缀皱；而edgeR則沒有斗这，因此需額外操作一下。推薦根據(jù)CPM（count-per-million啤斗，每百萬堿基中目標(biāo)基因的read count值）值進(jìn)行過濾表箭，例如使用CPM值為1作為標(biāo)準(zhǔn)，即當(dāng)某個(gè)基因在read count最低的樣本（文庫）中的count值大于(read count最低的樣品count總數(shù)/1000000)钮莲，則保留免钻。

library(edgeR)

#數(shù)據(jù)預(yù)處理
#（1）構(gòu)建 DGEList 對(duì)象
dgelist <- DGEList(counts = targets, group = group)

#（2）過濾 low count 數(shù)據(jù)，例如 CPM 標(biāo)準(zhǔn)化（推薦）
keep <- rowSums(cpm(dgelist) > 1 ) >= 2
dgelist <- dgelist[keep, , keep.lib.sizes = FALSE]

#（3）標(biāo)準(zhǔn)化崔拥，以 TMM 標(biāo)準(zhǔn)化為例
dgelist_norm <- calcNormFactors(dgelist, method = 'TMM')

2.3 計(jì)算差異倍數(shù)列表

整體過程分為兩步极舔。

• 一是估計(jì)離散度。根據(jù)試驗(yàn)設(shè)計(jì)的樣本分組链瓦，通過加權(quán)似然經(jīng)驗(yàn)貝葉斯模型拆魏，估算每個(gè)基因表達(dá)的離散度盯桦，其代表了基因表達(dá)值在個(gè)體間的差異。
• 二是模型擬合渤刃。edgeR包中提供了多種計(jì)算差異基因的算法模型拥峦，之間略有不同，如下展示了兩種卖子。

#差異表達(dá)基因分析
#首先根據(jù)分組信息構(gòu)建試驗(yàn)設(shè)計(jì)矩陣略号，分組信息中一定要是對(duì)照組在前，處理組在后
design <- model.matrix(~group)

#（1）估算基因表達(dá)值的離散度
dge <- estimateDisp(dgelist_norm, design, robust = TRUE)

#（2）模型擬合洋闽，edgeR 提供了多種擬合算法
#負(fù)二項(xiàng)廣義對(duì)數(shù)線性模型
fit <- glmFit(dge, design, robust = TRUE)
lrt <- topTags(glmLRT(fit), n = nrow(dgelist$counts))

write.table(lrt, 'control_treat.glmLRT.txt', sep = '\t', col.names = NA, quote = FALSE)

#擬似然負(fù)二項(xiàng)廣義對(duì)數(shù)線性模型
fit <- glmQLFit(dge, design, robust = TRUE)
lrt <- topTags(glmQLFTest(fit), n = nrow(dgelist$counts))

write.table(lrt, 'control_treat.glmQLFit.txt', sep = '\t', col.names = NA, quote = FALSE)

最后輸出了各基因的差異表達(dá)倍數(shù)表玄柠。

以負(fù)二項(xiàng)廣義對(duì)數(shù)線性模型的結(jié)果為例，包含了基因id喊递、log2轉(zhuǎn)化后的差異倍數(shù)（Fold Change）值随闪、顯著性p值以及校正后p值（默認(rèn)FDR校正）等主要信息。

差異表達(dá)基因分析結(jié)果

2.4 篩選差異表達(dá)基因

通過輸出的差異表達(dá)倍數(shù)表骚勘，即可自定義閾值篩選差異表達(dá)基因了铐伴。

例如這里以上述輸出的負(fù)二項(xiàng)廣義對(duì)數(shù)線性模型的結(jié)果為例，將它再次讀入到R中俏讹，并根據(jù)|log2FC|≥1以及FDR<0.01定義顯著變化的基因当宴，以及通過“up”和“down”分別區(qū)分上、下調(diào)的基因泽疆。

##篩選差異表達(dá)基因
#讀取上述輸出的差異倍數(shù)計(jì)算結(jié)果
gene_diff <- read.delim('control_treat.glmLRT.txt', row.names = 1, sep = '\t', check.names = FALSE)

#首先對(duì)表格排個(gè)序户矢，按 FDR 值升序排序，相同 FDR 值下繼續(xù)按 log2FC 降序排序
gene_diff <- gene_diff[order(gene_diff$FDR, gene_diff$logFC, decreasing = c(FALSE, TRUE)), ]

#log2FC≥1 & FDR<0.01 標(biāo)識(shí) up殉疼，代表顯著上調(diào)的基因
#log2FC≤-1 & FDR<0.01 標(biāo)識(shí) down梯浪，代表顯著下調(diào)的基因
#其余標(biāo)識(shí) none，代表非差異的基因
gene_diff[which(gene_diff$logFC >= 1 & gene_diff$FDR < 0.01),'sig'] <- 'up'
gene_diff[which(gene_diff$logFC <= -1 & gene_diff$FDR < 0.01),'sig'] <- 'down'
gene_diff[which(abs(gene_diff$logFC) <= 1 | gene_diff$FDR >= 0.01),'sig'] <- 'none'

#輸出選擇的差異基因總表
gene_diff_select <- subset(gene_diff, sig %in% c('up', 'down'))
write.table(gene_diff_select, file = 'control_treat.glmQLFit.select.txt', sep = '\t', col.names = NA, quote = FALSE)

#根據(jù) up 和 down 分開輸出
gene_diff_up <- subset(gene_diff, sig == 'up')
gene_diff_down <- subset(gene_diff, sig == 'down')

write.table(gene_diff_up, file = 'control_treat.glmQLFit.up.txt', sep = '\t', col.names = NA, quote = FALSE)
write.table(gene_diff_down, file = 'control_treat.glmQLFit.down.txt', sep = '\t', col.names = NA, quote = FALSE)

標(biāo)識(shí)顯著上下調(diào)的基因瓢娜，添加在最后一列

2.5 差異表達(dá)火山圖示例

已經(jīng)識(shí)別了顯著差異表達(dá)的基因，最后期望通過火山圖將它們表示出來褒颈，火山圖是文獻(xiàn)中常見統(tǒng)計(jì)圖表之一柒巫。

ggplot2為此提供了優(yōu)秀的作圖方案，參考以下示例谷丸。

##ggplot2 差異火山圖
library(ggplot2)

#默認(rèn)情況下堡掏，橫軸展示 log2FoldChange，縱軸展示 -log10 轉(zhuǎn)化后的 FDR
p <- ggplot(data = gene_diff, aes(x = logFC, y = -log10(FDR), color = sig)) +
geom_point(size = 1) +  #繪制散點(diǎn)圖
scale_color_manual(values = c('red', 'gray', 'green'), limits = c('up', 'none', 'down')) +  #自定義點(diǎn)的顏色
labs(x = 'log2 Fold Change', y = '-log10 FDR', title = 'control vs treat', color = '') +  #坐標(biāo)軸標(biāo)題
theme(plot.title = element_text(hjust = 0.5, size = 14), panel.grid = element_blank(), #背景色淤井、網(wǎng)格線布疼、圖例等主題修改
    panel.background = element_rect(color = 'black', fill = 'transparent'), 
    legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-1, 1), lty = 3, color = 'black') +  #添加閾值線
geom_hline(yintercept = 2, lty = 3, color = 'black') +
xlim(-12, 12) + ylim(0, 35)  #定義刻度邊界

p

差異表達(dá)基因火山圖