收到我們生信入門班的一個(gè)學(xué)員提問：使用edgeR算法做bulk轉(zhuǎn)錄組的差異表達(dá)分析，發(fā)現(xiàn)有一些基因看logFC值明明是下調(diào)，但是使用cpm值看象箱線圖的時(shí)候卻顯示基因是上調(diào)？為什么呢？
這是一個(gè)非常經(jīng)典的轉(zhuǎn)錄組數(shù)據(jù)分析里面有意思的案例，下面就來看看具體細(xì)節(jié)。
此外，我們技能樹最新一期生信入門就在11月3號(hào)開課，還沒上車的可以看一看瞄一瞄，詳細(xì)介紹頁面可以點(diǎn)擊：生信入門&數(shù)據(jù)挖掘線上直播課11月班，完全適合0基礎(chǔ)小白。
快來添加小助手報(bào)名咨詢吧：
簡(jiǎn)單做一個(gè)差異分析
首先使用edgeR算法做一個(gè)差異表達(dá)分析，這里你可以使用任何一個(gè)count矩陣作為數(shù)據(jù)輸入，再加上一個(gè)樣本分組信息，比如我們前面用的：
count轉(zhuǎn)TPM/FPKM實(shí)戰(zhàn)（GSE229904）
去下載：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE229904
我這里輸入的數(shù)據(jù)是一個(gè)隨便的具有四個(gè)樣本的數(shù)據(jù)，兩個(gè)case，兩個(gè)control：
rm(list = ls())
options(stringsAsFactors = F)
# 加載包
library(edgeR)
library(ggplot2)
## 讀取 count矩陣
head(count)
# 分組信息
group_list <- c("control","control","treat","treat")
group_list <- factor(group_list,levels = c("treat","control"))
group_list
table(group_list)
輸入數(shù)據(jù)都準(zhǔn)備好了，走edgeR算法：
# 構(gòu)建線性模型。0代表x線性模型的截距為0
design <- model.matrix(~0+group_list)
rownames(design) <- colnames(count)
colnames(design) <- levels(factor(group_list))
design
# 構(gòu)建edgeR的DGEList對(duì)象
DEG <- DGEList(counts=count, group=factor(group_list))
DEG <- calcNormFactors(DEG)
# 計(jì)算線性模型的參數(shù)
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)
# 擬合線性模型
fit <- glmFit(DEG, design)
# 進(jìn)行差異分析
lrt <- glmLRT(fit, contrast=c(1,-1)) #實(shí)驗(yàn)組在前
# 提取過濾差異分析結(jié)果
DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG),adjust.method = "BH"))
head(DEG_edgeR)
# 篩選上下調(diào)，設(shè)定閾值
fc_cutoff <- 1.2 # 1.2 1.5 2
pvalue <- 0.05 # 0.05, 0.01, 0.001
DEG_edgeR$regulated <- "normal"#給它一個(gè)初始值
DEG_edgeR$regulated[ DEG_edgeR$logFC > log2(fc_cutoff) & DEG_edgeR$PValue < pvalue ] <- "up"
DEG_edgeR$regulated[ DEG_edgeR$logFC < -log2(fc_cutoff) & DEG_edgeR$PValue < pvalue ] <- "down"
table(DEG_edgeR$regulated)
## 添加一列g(shù)ene symbol
library(org.Hs.eg.db)
library(clusterProfiler)
id2symbol <- bitr(rownames(DEG_edgeR), fromType = "ENSEMBL", toType = "SYMBOL", OrgDb = org.Hs.eg.db)
head(id2symbol)
DEG_edgeR <- cbind(GeneID=rownames(DEG_edgeR),DEG_edgeR)
DEG_edgeR_symbol <- merge(id2symbol,DEG_edgeR, by.x="ENSEMBL",by.y="GeneID",all.y=T)
head(DEG_edgeR_symbol)
差異結(jié)果如下：
挑一個(gè)基因看差異結(jié)果和表達(dá)值比較
這里用的上課的方法，使用簡(jiǎn)單的cpm進(jìn)行的標(biāo)準(zhǔn)化值來進(jìn)行同一個(gè)基因在不同樣本間的比較：
##====== 檢查是否上下調(diào)設(shè)置錯(cuò)了
# 挑選一個(gè)基因看其上下調(diào)方向
# ENSG00000037474 NSUN2
library(ggplot2)
cpm <- cpm(count) # count標(biāo)準(zhǔn)化
head(cpm)
exp <- c(t(cpm[match("ENSG00000037474",rownames(cpm)),]))
test <- data.frame(value=exp, group=group_list)
test
ggplot(data=test,aes(x=group,y=value,fill=group)) +
geom_boxplot() +
ggtitle("ENSG00000037474 NSUN2")
# 在差異結(jié)果中的方向
DEG_edgeR_symbol[DEG_edgeR_symbol$ENSEMBL=="ENSG00000037474", ]
可以看到 cpm 值顯示的是：相對(duì)對(duì)照組，在treat里面上調(diào)，但是差異結(jié)果里面的log2FC<0，卻是下調(diào)？
為什么呢？
為什么會(huì)出現(xiàn)上面的情況：算法下調(diào)，但是表達(dá)值上調(diào)？
其實(shí)主要在于edger使用的標(biāo)準(zhǔn)化方法為TMM算法，而我們?cè)诶L制箱線圖的時(shí)候卻用的 cpm 標(biāo)準(zhǔn)化，我后面試了 tmp 算法標(biāo)準(zhǔn)化，也會(huì)出現(xiàn)與 edger 不一致的情況，參考：count轉(zhuǎn)TPM/FPKM實(shí)戰(zhàn)（GSE229904）
edger 算法標(biāo)準(zhǔn)化 TMM 后的表達(dá)矩陣提?。?/span>
library(edgeR)
y <- DGEList(counts = count)
y <- calcNormFactors(y)
counts.tmm.edger <- round(cpm(y), 4)
head(counts.tmm.edger)
head(counts.tmm.edger["ENSG00000037474",c("control-1","control-1","treat-1","treat-2")])
# 對(duì)照組
mean(counts.tmm.edger["ENSG00000037474",c("control-1","control-1")])
# [1] 50.2231
# 實(shí)驗(yàn)組
mean(counts.tmm.edger["ENSG00000037474",c("treat-1","treat-2")])
# [1] 48.92365
這個(gè)時(shí)候要是想展示關(guān)鍵基因的表達(dá)值箱線圖，可以選擇算法內(nèi)部的標(biāo)準(zhǔn)化方法~
今天分享到這。
