Closed ixxmu closed 4 years ago
微信公众号:生信小知识
关注可了解更多的生信教程及技巧。问题或建议,请公众号留言;
前言ReferenceRUVSeq文献介绍RUVSeq核心思想文章简介spike-in原理不足数据支持其他normalization方法与RUV之间的差异不同normalization对于DEG分析的影响最最重要的一点有spike-in内参时的RUVg简要教程step1 载入R包及数据step2 数据过滤step3 创建SeqExpressionSet对象step4 绘制PCA及箱线图检查数据step5 使用quantitle normalizationstep6 使用RUVg normalizationstep7 使用edgeR进行GLM建模分析DEG通用型数据去批次分析(重点)step1 载入R包及数据step2 数据过滤step3 创建SeqExpressionSet对象step4 使用quantitle normalizationstep5 常规edgeR找DEGstep6 找到Empirical control genesstep7 使用RUVg normalizationstep8 使用edgeR进行GLM建模分析DEG后记
RUVSeq包,全称是Remove Unwanted Variation from RNA-Seq Data,名字起的非常通俗易懂,让人一看就知道这个R包在干什么了。不过这个R包的介绍中,提到的是用来做normalization,通过normalization来消除batch effect。其实仔细想想,我们是可以把normalization理解为一种去除library size批次的方法,所以normalization和batch effect其实是在干一件事情,只不过batch的意义更宽泛罢了。
之前讲过关于SVA去批次的原理及方法:
作为业内大家同样认同的工具,有必要再学习下RUVSeq包的使用方法及原理。
https://bioconductor.org/packages/release/bioc/vignettes/RUVSeq/inst/doc/RUVSeq.pdf
https://www.nature.com/articles/nbt.2931
原文如下:
https://www.nature.com/articles/nbt.2931
假设我们有一个表达矩阵,有n个样本,j个基因。我们知道在做差异分析的时候,大家基本上都是利用limma/edgeR/DEseq把我们的表达矩阵拟合到一个GLM(广义线性模型)上去建模分析的。在这个模型中,我们会告诉程序我们的感兴趣的协变量和不感兴趣的协变量:
上面的公式中:
Y 是 n*j 的表达矩阵
W 是 n*k 的batch矩阵,代表着我们不想要的因素
α 是 k*j 的矩阵,代表着每个batch对于每个基因的影响大小
X 是 n*p 的矩阵,代表着我们感兴趣的因素,例如用药与不用药,肿瘤与正常(例如如果是正常和肿瘤,那么我们的设计矩阵就是一列是肿瘤,一列是正常)
β 是 p*j 的矩阵,代表这每个感兴趣的因素对于每个基因的影响大小
O 是 n*j 的补偿矩阵,如果没有可以设为0
有3种方法去估计batch因素:RUVg、RUVs和RUVr,考虑到为了方便以后查找,这里放上他们之间的差异和使用范围:
文章使用了2套数据:
SEQC的数据:有A,B两种不同的样本。每种样本有4个文库,每个文库在16个lane进行测序。但是缺乏生物学重复。
zebrafish数据:有3个treated和3个control zebrafish样本。
对于SEQC的数据:
圆和三角分别代表了样本在不同的lane(lane1和lane2)进行测序
蓝色和红色分别代表了不同的样本(A和B)
可以看到没有normalization和用UQ进行normalization后,在PC2和PC3中仍然会根据圆和三角(也就是不同的lane)进行分群聚类,说明存在这种批次效应,且normalization没有去除干净。
同样的对于zebrafish数据:
在PCA图中仍然没有很好的聚在一起。
当使用RUVg处理后:
可以看到,对于SEQC的数据(图a),除了PC1可以很好区分开A和B样本外,在其他PC中所有样本基本上没有什么差异。说明只保留下了生物学差异。
对于zebrafish数据(图d),在PCA图上,3个treated和3个control zebrafish聚类效果比之前的好很多。通过PC2我们可以区分开这2类细胞。
(为什么不是PC1??)
反正,经过了上述操作后,发现RUVg的效果还不错,虽然使用RUVg时是用计算的方法找到的所谓negative controls,也就是那些不受我们感兴趣分组因素影响的基因。那么作者就好奇了,如果我用大家公认的spike-in来做,岂不是更合理?因为spike-in是外源加入的,在不同分组中肯定是不受影响的!是肯定满足RUVg的前提假设的!于是,就做了后面的分析,但是结果和作者想的不一样了!哈哈哈哈哈哈哈哈
加入spike-in的好处就在于,这些control是肯定不受实验因素的影响,也就是所谓的绝对的非差异基因。
当基因的整体表达模式发生改变时,例如全部基因均上调或下调时,我们只能依靠spike-in来做normalization,具体原因见我之前写过的推文:
Lovén等人2012年在Cell上发表了一篇文章来说明他们使用spike-in对这种情况进行normalization,具体做法分为3步:
计数样本所有的细胞数
按照细胞数加入等比例的spike-in
根据spike-in进行normalization
这个方法依赖于一个前提假设:技术对于spike-in的影响和技术对于表达谱基因的影响是一致的。
但是如果建库过程中,建库技术对于spike-in和基因的影响不同,那么用spike-in来做normalization则是不可靠的了。
2013年有篇文章,mRNA enrichment protocols determine the quantification characteristics of external RNA spike-in controls in RNA-Seq studies,根据他们的研究,发现加入的ERCC spike-in测序后得到的reads比例在不同技术重复样本之间差异较大。
同样的,我们今天的主角,RUVSeq包,也是认为spike-in其实不准!
一个理想的spike-in应该满足下面2点才可以:
不受分组因素的影响
但会受到其他技术因素的影响
同样,为了研究不同浓度的spike-in测序后得到的count数之间的关系,他们进行了实验,拿到的数据支持:log-read count 和 log-concentration之间成线性相关。
于是用SEQC的数据去做拟合,找出log-read count 和 log-concentration之间线性斜率的系数,结果,有趣的事就出现了:
a图:理论上来说,这些系数应该是在1,但是现在他们却分了不同的群,A1-3A和B4的系数很像,而A4,B1-3则相差较大。说明不同的建库过程会影响spike-in counts。
b,c图:同样的,不同建库的数据中mapping到ERCC的reads比例也不同。而且似乎还会受到分组因素的影响,因为途中样本B总是比样本A更高,且treated组也总是比control组高。(很奇怪)
d图:Ctl.1和Ctl.5样本原始未normalization的count数据,横轴是在不同个体之间的平均log-count值,纵坐标是logFC值。红色的点是spike-in基因的情况,红色的线是spike-in的lowess拟合方程曲线。蓝色是密度分布情况。由于2个样本都没有受到分组因素的影响,所以spike-in的logFC应该=0才对。
其他使用spike-in进行normalization的方法假设:spike-in不受分组因素的影响,且其他技术因素对于所有基因和spike-in的影响是相同的。
而RUVg则仅仅假设:技术因素对于所有基因和spike-in都有影响,但是具体影响的大小W是要靠每个基因的系数α来决定的。
因为SEQC数据有qRT-PCR的数据,所以作者把qRT-PCR的结果作为金标准,对用不同normalization处理后得到的DEG进行分析,结果如a,b图:
a:纵坐标是用测序数据经过不同normalization处理后找到的DEG的logFC 与 qRT-PCR的logFC 之间的比值,如果是0说明两者相符。从图中可以看到RUV系列的效果较好。
b:ROC曲线
而对于zebrafish数据集,由于没有金标准,作者想到一个非常诡异的办法——看那些非DEG的分布。理论上来推断,非DEG的p值分布应该是一个均一分布。
简单来说就是看所有p>0.05的p值是不是一个均一分布,为了方便观察,我在c图上画了一条线,可以看到所有的RUV系列结果都是均一分布的,而其他方法得到的结果有明显的偏向性。
看到文章最后的discussion,我重要看到了一点我最关心的东西了:
就拿目前我的分析来看,所有的tumor来自2个数据库,但是所有的normal都来自另一个数据库,这时我们就不能使用RUVr和RUVs,但是!!!RUVg可以处理这种情况!!
鉴于我只需要用到RUVg,所以我就只看并整理RUVg部分的教程了。
# 载入R包及数据
if (T) {
rm(list = ls())
options(stringsAsFactors = F)
library(RUVSeq)
library(edgeR)
library(zebrafishRNASeq)
data(zfGenes)
}
# 数据过滤
if (T) {
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]# 20806
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]# 59个
}
# 创建SeqExpressionSet对象
if (T) {
x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
phenoData = data.frame(x, row.names=colnames(filtered)))
set@phenoData@data
set
}
# 绘制PCA及箱线图检查数据
if (T) {
library(RColorBrewer)
colors <- brewer.pal(3, "Set2")
colors
plotRLE(set, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(set, col=colors[x], cex=1.2)
}
可以看到无论是箱线图上,还是PCA图上,都说明了需要做normalization。
# 使用UQ normalization
if (T) {
set <- betweenLaneNormalization(set, which="upper")
plotRLE(set, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(set, col=colors[x], cex=1.2)
}
做完UQ normalization之后的数据:
虽然箱线图的结果好多了,但是在PCA上仍然聚类不对。
由于这里有spike-in基因,所以可以直接进行分析。如果没有spike-in基因,则需要我们先选出一些可能不受我们的分组因素影响的基因,例如:管家基因,传递给RUVg进行分析:
# 使用RUVg进行normalization
if (T) {
# 按照教程在uq_set上进行normalization
if (T) {
ruvg_uq_set <- RUVg(x=uq_set, cIdx=spikes, k=1)
pData(ruvg_uq_set)
# x W_1
# Ctl1 Ctl -0.04539413
# Ctl3 Ctl 0.50347642
# Ctl5 Ctl 0.40575319
# Trt9 Trt -0.30773479
# Trt11 Trt -0.68455406
# Trt13 Trt 0.12845337
plotRLE(ruvg_uq_set, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(ruvg_uq_set, col=colors[x], cex=1.2)
}
# 直接在原始set数据上进行normalization
if (T) {
ruvg_set <- RUVg(x=set, cIdx=spikes, k=1)
pData(ruvg_set)
# x W_1
# Ctl1 Ctl -0.0390975
# Ctl3 Ctl 0.4088907
# Ctl5 Ctl 0.4414666
# Trt9 Trt -0.2138095
# Trt11 Trt -0.7527070
# Trt13 Trt 0.1552567
plotRLE(ruvg_set, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(ruvg_set, col=colors[x], cex=1.2)
}
}
关于RUVg
函数几个参数的说明:
x
:SeqExpressionSet对象
cIdx
:用作negative controls的基因,可以是基因名,也可以是逻辑值
k
:可能存在的batch因素个数
isLog
:数据是否经过log转换,默认为FALSE
RUVg
函数的做了两件事:
在phenoData
上添加一列数据,是我们想要的批次效应W_1。
对数据进行了normalization,在assayData下的normalizdCounts中保存。
问题:这里我发现它做的RUVg函数是针对set对象的,而set是已经经过了UQnormalization之后的数据,于是我就探索了下如果不经过UQnormalization,直接以count进行RUVg,发现两者数据结果不一致!!
按照教程在uq_set上进行normalization:
直接在原始set数据上进行normalization:
目前感觉可能需要按照教程上的做法,先对数据进行UQ normalization,再用RUVg进行normalization吧。已经在Github上向作者提问,等待作者的回复吧。
关于edgeR的使用教程可以去看我很久之前写过的教程:
# 使用edgeR拟合广义线性模型GLM
# limma拟合的是线性模型lmfit()
if (T) {
# 设计矩阵
design <- model.matrix(~x + W_1, data=pData(ruvg_uq_set))
# 创建DEGList对象
y <- DGEList(counts=counts(ruvg_uq_set), group=x)
# UQ normalization
y <- calcNormFactors(y, method="upperquartile")
# 建模
if (T) {
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
}
# 提取结果
if (T) {
DEG_result <- topTags(lrt, n=nrow(y))$table
}
}
当然了,大家一般测序里很少会做spike-in内参,所以,这里就要我们自己找到一些非DEG的基因当作negative gene了。
# 载入R包及数据
if (T) {
rm(list = ls())
options(stringsAsFactors = F)
library(RUVSeq)
library(edgeR)
library(zebrafishRNASeq)
data(zfGenes)
}
# 数据过滤
if (T) {
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]# 20806
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]# 59个
}
# 创建SeqExpressionSet对象
if (T) {
x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
phenoData = data.frame(x, row.names=colnames(filtered)))
set@phenoData@data
set
}
# 使用UQ normalization
if (T) {
set <- betweenLaneNormalization(set, which="upper")
}
先不指定任何其他可能的影响因素,仅指定感兴趣的分组变量,找出存在的DEG:
# 先仅拟合感兴趣的分组变量
if (T) {
design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
}
if (T) {
# 去除TOP5000的差异基因,认为剩下的15865个基因不是DEG
# 这15865个非DEG当作negative genes进行后续RUVg的normalization
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
}
# RUVg进行normalization
if (T) {
set2 <- RUVg(uq_set, empirical, k=1)
pData(set2)
# x W_1
# Ctl1 Ctl -0.10879677
# Ctl3 Ctl 0.23066424
# Ctl5 Ctl 0.19926266
# Trt9 Trt 0.07672121
# Trt11 Trt -0.83540924
# Trt13 Trt 0.43755790
plotRLE(set2, outline=FALSE, ylim=c(-4, 4), col=colors[x])
plotPCA(set2, col=colors[x], cex=1.2)
}
从PCA的聚类图来看,似乎用非DEG计算得到的基因当作negative genes的效果比使用spike-in的效果还要好,因为在PCA结果上,似乎从PC2来看终于把Tr9从Ctl组区分开来了:
# 使用edgeR拟合广义线性模型GLM
if (T) {
# 设计矩阵
design <- model.matrix(~x + W_1, data=pData(set2))
# 创建DEGList对象
y <- DGEList(counts=counts(set2), group=x)
# UQ normalization
y <- calcNormFactors(y, method="upperquartile")
# 建模
if (T) {
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
}
# 提取结果
if (T) {
DEG_result <- topTags(lrt, n=nrow(y))$table
}
}
到此也算是把SVA和RUV 2个引用最多的去批次工具学完了,下面大家就把学到的东西运用到自己的数据上吧
https://mp.weixin.qq.com/s/riTzGYY-bERVgdNJ-EbTuA