寫這篇文章一部分原因是填2年前的一個(gè)坑 轉(zhuǎn)錄組入門（7）：差異表達(dá)分析. 另一部分原因是GQ最近又在搞一波無(wú)重復(fù)的差異表達(dá)分析, 所以專門去學(xué)了edgeR
我個(gè)人是不太推薦沒(méi)有重復(fù)的差異表達(dá)分析，畢竟統(tǒng)計(jì)學(xué)上的p值是為了證明兩個(gè)樣本的差異是真實(shí)存在而不是抽樣誤差導(dǎo)致, 但是你單個(gè)樣本如何計(jì)算變異呢?
因此每當(dāng)別人提問(wèn)的時(shí)候, 我個(gè)人的建議就是定性看看倍數(shù)變化吧. 但是如果真的強(qiáng)行要算p值, 其實(shí)也不是不行, edgeR就是一種選擇.
環(huán)境準(zhǔn)備
我們需要安裝兩個(gè)R包，一個(gè)是edgeR, 一個(gè)是airway. 其中airway是一個(gè)數(shù)據(jù)集包, 功能就是提供一個(gè)用于分析的數(shù)據(jù)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!requireNamespace("edgeR", quietly = TRUE))
BiocManager::install("edgeR")
if (!requireNamespace("airway", quietly = TRUE))
BiocManager::install("airway")
加載R包
library(edgeR)
library(airway)
構(gòu)建DGEList
DGEList是edgeR分析流程中必須的對(duì)象. 構(gòu)建該對(duì)象需要提供兩類信息: 表達(dá)量矩陣和分組信息.
為了方便大家重復(fù),我們這里的數(shù)據(jù)來(lái)自于airway. 對(duì)于你自己的數(shù)據(jù), 可以用read.table等函數(shù)進(jìn)行導(dǎo)入.
data("airway")
expr_matrix <- assay(airway)
meta_info <- colData(airway)
expr_matrix 是一個(gè) 64102 個(gè)基因和8個(gè)樣本的矩陣.meta_info 里存放的是樣本的元信息, 記錄樣本的處理, 來(lái)源等信息. 我們這里就用一部分?jǐn)?shù)據(jù), 也就是前兩列構(gòu)建DGEList對(duì)象
counts <- expr_matrix[,1:2]
group <- 1:2
y <- DGEList(counts=counts, group = group)
數(shù)據(jù)過(guò)濾
由于原來(lái)的表達(dá)量矩陣基因數(shù)太大, 可能存在某些基因根本沒(méi)有表達(dá), 因此需要預(yù)先過(guò)濾
keep <- rowSums(cpm(y)>1) >= 1
y <- y[keep, , keep.lib.sizes=FALSE]
這部分代碼的意思指的是保留在至少在一個(gè)樣本里有表達(dá)的基因(CPM > 1)。 基因數(shù)就從原來(lái)的64102降到14756
標(biāo)準(zhǔn)化
考慮到測(cè)序深度不同, 我們需要對(duì)其進(jìn)行標(biāo)準(zhǔn)化, 避免文庫(kù)大小不同導(dǎo)致的分析誤差.
edgeR里默認(rèn)采用TMM(trimmed mean of M-values) 對(duì)配對(duì)樣本進(jìn)行標(biāo)準(zhǔn)化，用到的函數(shù)是calcNormFactors
y <- calcNormFactors(y)
差異表達(dá)分析
不同差異表達(dá)分析工具的目標(biāo)就是預(yù)測(cè)出dispersion(離散值), 有了離散值就能夠計(jì)算p值. 那么dispersion怎么計(jì)算呢? edgeR給了幾個(gè)方法
方法一: 根據(jù)經(jīng)驗(yàn)給定一個(gè)值(BCV, square-root-dispersion). edgeR給的建議是, 如果你是人類數(shù)據(jù), 且實(shí)驗(yàn)做的很好(無(wú)過(guò)多的其他因素影響), 設(shè)置為0.4, 如果是遺傳上相似的模式物種, 設(shè)置為0.1, 如果是技術(shù)重復(fù), 那么設(shè)置為0.01
這里用的數(shù)據(jù)是人類, 此處設(shè)置為 0.4
y_bcv <- y
bcv <- 0.4
et <- exactTest(y_bcv, dispersion = bcv ^ 2)
我們用decideTestsDGE看下有多少基因上調(diào), 多少基因下調(diào). 設(shè)置p.value=0.05
gene1 <- decideTestsDGE(et, p.value = 0.05, lfc = 0)
summary(gene1)
2-1
Down 4
NotSig 14733
Up 19
差異基因少的可憐, 只有4+19個(gè)。我們可以嘗試調(diào)整下BCV
y_bcv <- y
bcv <- 0.2
et2 <- exactTest(y_bcv, dispersion = bcv ^ 2)
gene2 <- decideTestsDGE(et2, p.value = 0.05, lfc = 0)
summary(gene2)
2-1
Down 159
NotSig 14380
Up 217
這個(gè)時(shí)候的差異基因上升到了159 + 217個(gè).
方法2: 根據(jù)已知一些不會(huì)發(fā)生改變的基因推測(cè)dispersion. 假設(shè)你已經(jīng)知道了一些基因是不會(huì)發(fā)生變化，例如管家基因，那么我們就可以通過(guò)它們來(lái)預(yù)測(cè)dispersion.
先復(fù)制原來(lái)對(duì)象的一份拷貝.
y1 <- y
y1$samples$group <- 1
然后你需要提供一個(gè)變量，housekeeping存放著已知不改變的基因名，例如管家基因中，然后進(jìn)行dispersion估計(jì)
這里需要一個(gè)housekeeping的向量, 來(lái)自教程的最后部分
y0 <- estimateDisp(y1[housekeeping,], trend="none", tagwise=FALSE)
最后加入到原來(lái)的數(shù)據(jù)集中進(jìn)行分析。
y$common.dispersion <- y0$common.dispersion
design <- model.matrix(~group)
fit <- glmFit(y, design)
lrt <- glmLRT(fit)
這個(gè)時(shí)候我們?cè)倏聪虏町惢? 是341 + 388
gene3 <- decideTestsDGE(lrt, p.value = 0.05, lfc = 0)
summary(gene3)
group
Down 341
NotSig 14027
Up 388
說(shuō)實(shí)話, 好像和預(yù)先設(shè)置的沒(méi)啥區(qū)別. 但是我們用韋恩圖比較下
library(gplots)
venn(list(gene1=names(gene1@.Data[gene1@.Data == 1,]),
gene2=names(gene1@.Data[gene2@.Data == 1,]),
gene3=names(gene1@.Data[gene3@.Data == 1,])
))
venn(list(gene1=names(gene1@.Data[gene1@.Data == -1,]),
gene2=names(gene2@.Data[gene2@.Data == -1,]),
gene3=names(gene3@.Data[gene3@.Data == -1,])
))
從中可以發(fā)現(xiàn)根據(jù)已知不變基因預(yù)測(cè)dispersion得到的差異基因最多。 還可以畫個(gè)火山圖看下
library(ggplot2)
df <- lrt$table
ggplot(df, aes(x=logFC, y=-log10(df$PValue)) ) +
geom_point() +
ylab("-log10(p value)") +
theme_bw()
和有重復(fù)表達(dá)進(jìn)行比較
由于我們的數(shù)據(jù)集原來(lái)是有重復(fù)的，那么我們就可以比較下無(wú)重復(fù)和有重復(fù)之間會(huì)相差多少基因
group <- meta_info$dex
y_rep <- DGEList(counts=expr_matrix, group = group)
keep <- rowSums(cpm(y_rep)>1) >= 5
y_rep <- y_rep[keep, , keep.lib.sizes=FALSE]
y_rep <- calcNormFactors(y_rep)
design <- model.matrix(~group)
y_rep <- estimateDisp(y_rep, design = design)
fit <- glmQLFit(y_rep, design)
qlf.2vs1 <- glmQLFTest(fit, coef=2)
我們看下差異基因的數(shù)目, 是1050 + 869, 明顯多于之前.
res <- decideTestsDGE(qlf.2vs1, p.value = 0.05, lfc = 0)
summary(res)
將我們有重復(fù)的前100差異基因和無(wú)重復(fù)的前100差異基因進(jìn)行比較
c1 <- row.names(topTags(et2,n = 100 ))
c2 <- row.names(topTags(qlf.2vs1, n = 100 ))
intersect_gene <- intersect(c1,c2)
我們將這些交集基因標(biāo)記在有重復(fù)的火山圖上
library(ggplot2)
library(ggrepel)
data <- qlf.2vs1$table
data$significant <- as.factor(data$PValue<0.05 & abs(data$logFC) > 1)
data2 <- data[intersect_gene,]
data2$geneID <- rownames(data2)
p <- ggplot(data=data, aes(x=logFC, y =-log10(PValue),color=significant)) +
geom_point(alpha=0.8, size=1.2)+
scale_color_manual(values =c("black","red"))+
labs(title="Volcanoplot", x="log2 (fold change)",y="-log10 (q-value)")+
theme(plot.title = element_text(hjust = 0.4))+
geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+
#theme(legend.position='none')
theme_bw()+
theme(panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
geom_text(data=data2, aes(label=geneID),col="red",alpha = 1) +
geom_text_repel(data=data2, aes(label=geneID),col="black",alpha = 0.8)
print(p)
好消息是無(wú)重復(fù)情況下的確能找到一些明顯差異表達(dá)的基因，但是壞消息是改變不怎么明顯的重要基因可能就會(huì)因?yàn)槟阍O(shè)置一個(gè)閾值被篩選掉。因此無(wú)重復(fù)的分析還是能做的，就是閾值需要放寬些。
下面的代碼是用來(lái)提取一些沒(méi)有顯著變化的基因作為之前說(shuō)的housekeeping
nosig <- names(res@.Data[res@.Data == 0,])
gene_name <- sample(nosig,500)
# y1來(lái)自上面無(wú)重復(fù)差異表達(dá)代碼
y1_gene <- rownames(y1)
housekeeping <- which(y1_gene %in% gene_name)
