使用edgeR进行无重复差异表达分析

写这篇文章一部分原因是填2年前的一个坑 转录组入门（7）：差异表达分析. 另一部分原因是GQ最近又在搞一波无重复的差异表达分析, 所以专门去学了edgeR

我个人是不太推荐没有重复的差异表达分析，毕竟统计学上的p值是为了证明两个样本的差异是真实存在而不是抽样误差导致, 但是你单个样本如何计算变异呢?

因此每当别人提问的时候, 我个人的建议就是定性看看倍数变化吧. 但是如果真的强行要算p值, 其实也不是不行, edgeR就是一种选择.

环境准备

我们需要安装两个R包，一个是edgeR, 一个是airway. 其中airway是一个数据集包, 功能就是提供一个用于分析的数据

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)

构建DGEList

DGEList是edgeR分析流程中必须的对象. 构建该对象需要提供两类信息: 表达量矩阵和分组信息.

为了方便大家重复,我们这里的数据来自于airway. 对于你自己的数据, 可以用read.table等函数进行导入.

data("airway")
expr_matrix <- assay(airway)
meta_info <- colData(airway)

expr_matrix 是一个 64102 个基因和8个样本的矩阵.meta_info 里存放的是样本的元信息, 记录样本的处理, 来源等信息. 我们这里就用一部分数据, 也就是前两列构建DGEList对象

counts <- expr_matrix[,1:2]
group <- 1:2
y <- DGEList(counts=counts, group = group)

数据过滤

由于原来的表达量矩阵基因数太大, 可能存在某些基因根本没有表达, 因此需要预先过滤

keep <- rowSums(cpm(y)>1) >= 1
y <- y[keep, , keep.lib.sizes=FALSE]

这部分代码的意思指的是保留在至少在一个样本里有表达的基因(CPM > 1)。 基因数就从原来的64102降到14756

标准化

考虑到测序深度不同, 我们需要对其进行标准化, 避免文库大小不同导致的分析误差.

edgeR里默认采用TMM(trimmed mean of M-values) 对配对样本进行标准化，用到的函数是calcNormFactors

y <- calcNormFactors(y)

差异表达分析

不同差异表达分析工具的目标就是预测出dispersion(离散值), 有了离散值就能够计算p值. 那么dispersion怎么计算呢? edgeR给了几个方法

方法一: 根据经验给定一个值(BCV, square-root-dispersion). edgeR给的建议是, 如果你是人类数据, 且实验做的很好(无过多的其他因素影响), 设置为0.4, 如果是遗传上相似的模式物种, 设置为0.1, 如果是技术重复, 那么设置为0.01

这里用的数据是人类, 此处设置为 0.4

y_bcv <- y
bcv <- 0.4
et <- exactTest(y_bcv, dispersion = bcv ^ 2)

我们用decideTestsDGE看下有多少基因上调, 多少基因下调. 设置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个。我们可以尝试调整下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

这个时候的差异基因上升到了159 + 217个.

方法2: 根据已知一些不会发生改变的基因推测dispersion. 假设你已经知道了一些基因是不会发生变化，例如管家基因，那么我们就可以通过它们来预测dispersion.

先复制原来对象的一份拷贝.

y1 <- y
y1$samples$group <- 1

然后你需要提供一个变量，housekeeping存放着已知不改变的基因名，例如管家基因中，然后进行dispersion估计

这里需要一个housekeeping的向量, 来自教程的最后部分

y0 <- estimateDisp(y1[housekeeping,], trend="none", tagwise=FALSE)

最后加入到原来的数据集中进行分析。

y$common.dispersion <- y0$common.dispersion
design <-  model.matrix(~group)
fit <- glmFit(y, design)
lrt <- glmLRT(fit)

这个时候我们再看下差异基因, 是341 + 388

gene3 <- decideTestsDGE(lrt, p.value = 0.05, lfc = 0)
summary(gene3)

       group
Down     341
NotSig 14027
Up       388

说实话, 好像和预先设置的没啥区别. 但是我们用韦恩图比较下

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,])
          ))

下调基因

从中可以发现根据已知不变基因预测dispersion得到的差异基因最多。 还可以画个火山图看下

library(ggplot2)
df <- lrt$table
ggplot(df, aes(x=logFC, y=-log10(df$PValue)) ) +
  geom_point() +
  ylab("-log10(p value)") +
  theme_bw()

火山图

和有重复表达进行比较

由于我们的数据集原来是有重复的，那么我们就可以比较下无重复和有重复之间会相差多少基因

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)

我们看下差异基因的数目, 是1050 + 869, 明显多于之前.

res <- decideTestsDGE(qlf.2vs1, p.value = 0.05, lfc = 0)
summary(res)

将我们有重复的前100差异基因和无重复的前100差异基因进行比较

c1 <- row.names(topTags(et2,n = 100 ))
c2 <- row.names(topTags(qlf.2vs1, n = 100 ))
intersect_gene <- intersect(c1,c2)

我们将这些交集基因标记在有重复的火山图上

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)

火山图2

好消息是无重复情况下的确能找到一些明显差异表达的基因，但是坏消息是改变不怎么明显的重要基因可能就会因为你设置一个阈值被筛选掉。因此无重复的分析还是能做的，就是阈值需要放宽些。

下面的代码是用来提取一些没有显著变化的基因作为之前说的housekeeping

nosig <- names(res@.Data[res@.Data == 0,])
gene_name <- sample(nosig,500)
# y1来自上面无重复差异表达代码
y1_gene <- rownames(y1)
housekeeping <- which(y1_gene %in% gene_name)