DALS025-批次效应01-什么是批次效应

title: DALS025-批次效应01-什么是批次效应
date: 2019-08-25 12:0:00
type: "tags"
tags:

  • 批次效应
    categories:
  • Data Analysis for the life sciences

前言

这一部分是《Data Analysis for the life sciences》的第10章批次效应的第1小节，这一部分的主要内容涉及批次效应(Batch Effects)的介绍，有关批次效应的前言部分Rmarkdown文档参见作者的Github。

什么是批次效应

高通量研究中一个经常被忽视的问题就是批次效应(batch effects)，批次效应受到当时检测的实验室条件、试剂批次和人员差异的影响。当批次效应与我们目标结果混淆并惬以不正确的结果时，这就成了一个主要问题。在这一章中，我们将详细地描述批次效应：对于批次效应如何检测、解释、建模和调整。

批次效应是基因组学研究中面临的最大挑战，尤其是在精确医学这个背景下。大多数情况（但并非全部）下，高通量技术已经被报道存在着一种形式或另外一种形式的批次效应[Leek et al. (2010) Nature Reviews Genetics 11, 733-739]。但是，批次效应并非基因组学所特有的。实际上， 在1972年一篇文献中，Mj Youden就在对物理常数的经验估计的背景下提到了批次效应。他指出，物理常数“存在着主观估计的特征”，以及如何在不同实验室之间变化的。例如，在下面的Table1中，Youden就展示了来自于不同实验室又的天文单的估计值。这些报告包括对离差(spread)（现在我们称之为置信区间）的估计。

library(rafalib)
library(downloader)
##Download the data from
url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/astronomicalunit.csv"
filename <- tempfile() 
if (!file.exists(filename)) download(url, destfile=filename)

dat <- read.csv(filename)
year <-  jitter(dat[,2]) ##add jitter so points are not on top of each other

##Use color to denote the labs that reported more than one measurement
labs <- as.character(dat[,1]) ##what lab did it
labs[ !labs%in%c("Jodrell Bank","Spencer Jones")] <- "Others"
labs <- factor(labs, levels=c("Others","Spencer Jones","Jodrell Bank"))
cols=as.numeric(labs)

current <- 92.956039 ##this is the current estimate in millions of mph

mypar()
plot(year, dat[,3], ylim=c(min(dat[,4]),max(dat[,5])), pch=16, col=cols, 
     xlab="Year",ylab="Astronomical unit (millions of miles)")
for(i in 1:nrow(dat))
  lines(c(year[i],year[i]),c(dat[i,4],dat[i,5]),col=cols[i],lwd=3)
legend("topright", legend=levels(labs), col=seq_along( labs ) ,cex=0.75, lty=1,pch=16)
abline(h=current,lty=2)
text(1905,current,"Current estimate",pos=3)
image

从不同实验室之间的可变性以及报道的界限（可以理解为置信区间）来看，都不能解释这种可变性，这就清楚地表明不同实验室之间，而非某个实验室内部存在着某些效应。这种类型的变异就是我们所谓的批次效应。请注意，有些实验室报告了两个估计值（紫色和橙色），这样我们就在相同的实验室也看到了批次效应。
我们可以使用统计学符号来精确地描述这个问题。我们使用下面的公式来表示这些检测值：
Y_{i,j} = \mu + \varepsilon_{i,j}, j=1,\dots,N

其中，Y_{i,j} 表示第 i 个实验室的第 j 个检测值，而 \mu 表示真实的物理常数，\varepsilon_{i,j} 表示独立的检测误差。为了解释可变性，我们引入了 \varepsilon_{i,j} ，随后我们根据数据来计算标准误。就像本书前面提到的那样，我们使用 N 值均值来估计物理物理常数，如下所示：
\bar{Y}_i = \frac{1}{N} \sum_{i=1}^{N} Y_{i,j}

再来构建一个置信区间：
\bar{Y}_i \pm 2 s_i / \sqrt{N} \mbox{ with } s_i^2= \frac{1}{N-1} \sum_{i=1}^N (Y_{i,j} - \bar{Y}_i)^2

但是，这个置信区间太小，它无法覆盖批次效应的变异，下面是一个更加合适的模型：
Y_{i,j} = \mu + \gamma_i + \varepsilon_{i,j}, j=1, \dots, N

其中 \gamma_i 表示一个实验室特定的偏差或批次效应(batch effect)。从图片上我们可以明显看出来，实验室之间 \gamma 的变化要远大于一个实验室内部 \varepsilon 的变化。用统计学的术语来描述这个问题就是，\mu 和 \gamma 无法被识别。我们可以估计 \mu_i+\gamma_i ，但是无法区分开。我们可以将 \gamma 视为一个随机变量。在这个案例中，每个实验室都有一个错误项 \gamma_i ，这一项在贯穿于该实验室的所有检测值，在每个检测值中，这一项是相同的，但是实验室与实验室间的这一项则不同。因此，这个问题就可以用以下方程表示：
s_i / \sqrt{N} \mbox{ with } s_i^2= \frac{1}{N-1} \sum_{i=1}^N (Y_{ij} - \bar{Y}_i)^2

这是对标准差的低估，因为它没有解释由 \gamma 导致的实验室内的相关性（这一句不懂，原文如下：

Under this interpretation the problem is that:is an underestimate of the standard error since it does not account for the within lab correlation induced by \gamma.

如果我们假设 \gamma=0 ，那么利用来自几个实验室的数据，我们实际上可以估计出 \gamma 。或者说我们可以将它们视为随机效应，简单地将它们当成一个新的估计值，并使用所有的检测值来计算其标准差。以下是我们将报告中的均值值视为随机观察值的置信区间：

avg <- mean(dat[,3])
se <- sd(dat[,3]) / sqrt(nrow(dat))
cat("95% confidence interval is: [",avg-1.96*se,",", avg+1.96*se,"]")
cat("which does include the current estimate is:",current)

结果如下所示：

> cat("95% confidence interval is: [",avg-1.96*se,",", avg+1.96*se,"]")
95% confidence interval is: [ 92.8727 , 92.98542 ]> cat("which does include the current estimate is:",current)
which does include the current estimate is: 92.95604

Youden的论文还包括最近关于光速以及策略常数的批次效应估计的案例。在这一章里，我们只展示高通量生物数据中批次效应的广泛性和复杂性。

混杂

书中有一个术语，即confounding，这里译为混杂。这一部分内容的Rmarkdown文档可以参考作者的Github。

当批次效应与我们的目标结果混杂时，就会导致严重的𨤹。这里我们描述一下混杂(confounding)，以及它与我们数据解释的关系。

我们从这本书或者说从任何其它数据分析课程中尝到最重要的思想之一就是“相关不等于因果”。这句话多数情况就是真的，一个常见的案例就是混杂(confounding)。简单地讲，当我们观察到 X 和 Y 之间存在着相关(correlation)或关联(association)时，往往就会存在着混杂，但严格来说，这是因为 X 和 Y 都依赖于一个无关的变量 Z 。这里我们会描述一个Simposon悖论，这是一个基于一个著名的法律案件的案例，接着，我们还会提到一个在高通量生物学研究中的一个混杂案例。

案例之Simpson悖论

加州大学伯克分校(UCB)1973年的入学数据显示，在录取的学生中，44%是男性，30%是女性，显著男性更多，以下是这些数据：

library(dagdata)
data(admissions)
admissions$total=admissions$Percent*admissions$Number/100
##percent men get in
sum(admissions$total[admissions$Gender==1]/sum(admissions$Number[admissions$Gender==1]))
##percent women get in
sum(admissions$total[admissions$Gender==0]/sum(admissions$Number[admissions$Gender==0]))

结果如下所示：

> sum(admissions$total[admissions$Gender==1]/sum(admissions$Number[admissions$Gender==1]))
[1] 0.4451951
> ##percent women get in
> sum(admissions$total[admissions$Gender==0]/sum(admissions$Number[admissions$Gender==0]))
[1] 0.3033351

卡方检验的结果明确拒绝零假设，即录取的性别是独立的，两者不受影响（也就是说，男女录取公平），如下所示：

##make a 2 x 2 table
index = admissions$Gender==1
men = admissions[index,]
women = admissions[!index,]
menYes = sum(men$Number*men$Percent/100)
menNo = sum(men$Number*(1-men$Percent/100))
womenYes = sum(women$Number*women$Percent/100)
womenNo = sum(women$Number*(1-women$Percent/100))
tab = matrix(c(menYes,womenYes,menNo,womenNo),2,2)
print(chisq.test(tab)$p.val)

p值如下所示：

> print(chisq.test(tab)$p.val)
[1] 9.139492e-22

但经过更仔细的观察会发现一个自相矛盾的结果，以下是按专业划分后，不同性别的录取百分比：

y=cbind(admissions[1:6,c(1,3)],admissions[7:12,3])
colnames(y)[2:3]=c("Male","Female")
y

结果如下所示：

> y
  Major Male Female
1     A   62     82
2     B   63     68
3     C   37     34
4     D   33     35
5     E   28     24
6     F    6      7

从结果中我们可以发现，不同专业之间没有性别偏见。

但是，我们在前面使用了卡方检验发现，入学和性别之间存在着某种关系。然而，当我们的数据按照不同专业进行分组时，这种依赖性似乎消失了，这是怎么一回事呢？

这就是Simpson悖论的一个案例。

我们上面进行的卡方检验表明，入学和性别之间存在依赖关系。然而，当数据按专业分组时，这种依赖性似乎消失了。到底怎么回事？现在我们绘制一个能显示申请专业的人，以及最终进入这个专业读书的人的百分比的图形，用于说明男女在录取比例方面的问题，如下所示：

y=cbind(admissions[1:6,5],admissions[7:12,5])
y=sweep(y,2,colSums(y),"/")*100
x=rowMeans(cbind(admissions[1:6,3],admissions[7:12,3]))
library(rafalib)
mypar()
matplot(x,y,xlab="percent that gets in the major",
        ylab="percent that applies to major",
        col=c("blue","red"),cex=1.5)
legend("topleft",c("Male","Female"),col=c("blue","red"),pch=c("1","2"),box.lty=0)
image

从图片上我们可以发现，男生更想申请那些“容易”的专业。也就是说，男生这个变量与“容易”专业这个变量发生了混杂。

混杂的图形解释

在这一部分里，我们将混杂图形化。在下面的图形中，每个字母表示一个学生。被录取的人用绿色表示，字母表示专业。在第1张图中，所有相同性别的学生都被放在一起，我们可以发现，男生中绿色的比较更大，如下所示：

###make data for plot
library(rafalib)
mypar()
CEX=0.5
NC <- 70
tmp=rowSums(tab)
FNC <- round(NC*tmp[2]/tmp[1])
SCALE <- 1
makematrix<-function(x,n,addx=0,addy=0){
  m<-ceiling(length(x)/n)
  expand.grid(1:n+addx,addy+1:m)[seq(along=x),] 
}
males<- sapply(1:6,function(i){
  tot=admissions[i,2]*SCALE
  p=admissions[i,3]/100
  x=rep(c(0,1),round(tot*c(1-p,p)))
})
allmales<-Reduce(c,males)
females<- sapply(7:12,function(i){
  tot=admissions[i,2]*SCALE
  p=admissions[i,3]/100
  rep(c(0,1),round(tot*c(1-p,p)))
})
allfemales<-Reduce(c,females)
mypar(1,1)
malepoints <- makematrix(allmales,NC)
femalepoints <- makematrix(allfemales,FNC,NC+NC/10)
NR <- max(c(malepoints[,2],femalepoints[,2]))
plot(0,type="n",xlim=c(min(malepoints[,1]),max(femalepoints[,1])),ylim=c(0,NR),xaxt="n",yaxt="n",xlab="",ylab="")
PCH=LETTERS[rep(1:6,sapply(males,length))]
o<-order(-allmales)
points(malepoints,col=2-allmales[o],pch=PCH[o],cex=CEX)
PCH=LETTERS[rep(1:6,sapply(females,length))]
o<-order(-allfemales)
points(femalepoints,col=2-allfemales[o],pch=PCH[o],cex=CEX)
abline(v=NC+NC/20)
axis(side=3,c(NC/2,NC+NC/2),c("Male","Female"),tick=FALSE)
image

现在我们按照专业对不同性别的学生进行分级。这里我们要注意，大多数被录取的学生（绿色）来源于两个容易的专业，即A和B，如下所示：

mypar()
malepoints <- vector("list",length(males))
femalepoints <- vector("list",length(males))
N<- length(males)
 
ADDY <- vector("numeric",N+1)
for(i in 1:N){
  malepoints[[i]] <- makematrix(males[[i]],NC,0,ADDY[i])
  femalepoints[[i]] <- makematrix(females[[i]],FNC,NC+NC/10,ADDY[i])
   ADDY[i+1] <- max(malepoints[[i]][,2],femalepoints[[i]][,2])+1
}
plot(0,type="n",
     xlim=c( min(sapply(malepoints,function(x)min(x[,1]))),max(sapply(femalepoints,function(x)max(x[,1])))),
  ylim=c(0,max(sapply(femalepoints,function(x)max(x[,2])))),xaxt="n",yaxt="n",xlab="",ylab="")
          
for(i in 1:N){
  points(malepoints[[i]],col=2+sort(-males[[i]]),pch=LETTERS[i],cex=CEX)
  points(femalepoints[[i]],col=2+sort(-females[[i]]),pch=LETTERS[i],cex=CEX)
  if(i>1) abline(h=ADDY[i])
  }
abline(v=NC+NC/20)
axis(side=3,c(NC/2,NC+FNC/2),c("Male","Female"),tick=FALSE)
axis(side=2,ADDY[-1]/2+ADDY[-length(ADDY)]/2,LETTERS[1:N],tick=FALSE,las=1)
image

分层后的均值

在下图中，我们可以看到，我们根据专业设定条件或分层后，然后再看差异，也就是说，当我们控制了混杂项后，这就混杂效应就消失了，如下所示：

y=cbind(admissions[1:6,3],admissions[7:12,3])
matplot(1:6,y,xaxt="n",xlab="major",ylab="percent",col=c("blue","red"),cex=1.5)
axis(1,1:6,LETTERS[1:6])
legend("topright",c("Male","Female"),col=c("blue","red"),pch=c("1","2"),
       box.lty=0)
image

事实上，在不同专业录取的学生中，性别差异仅表示为，男生比女生高3.5%，如下所示：

mean(y[,1]-y[,2])

结果如下所示：

> mean(y[,1]-y[,2])
[1] -3.5

棒球中的Simpson悖论

在棒球比赛的统计中我们也经常看到Simpson悖论，现在我们来看一个有名的案例，在1995年和1996年，David Justice的平均击球率高于Derek Jeter，但是Jeter的击球率却高于整体平均水平，如下所示：

1995 1996 Combined
Derek Jeter 12/48 (.250) 183/582 (.314) 195/630 (.310)
David Justice 104/411 (.253) 45/140 (.321) 149/551 (.270)

这里的混杂项就是比赛场次，Jeter的比赛场次数目更多，但是结果却是Justice击球率更高。

混杂之高通量案例

为了描述我们在生物学中遇到的混杂问题，我们将会使用《Common genetic variants account for differences in gene expression among ethnic groups》这篇文献中的数据集，这篇文献指出，两个不同种族的血液大约有50%的差异基因，现在我们来下载这个数据集：

library(Biobase) ##available from Bioconductor
library(genefilter) 
load("GSE5859.rda") ##available from github

我们使用Bioconductor中的函数exprs和pData就能提取基因的表达数据与样本信息，如下所示：

geneExpression = exprs(e)
sampleInfo = pData(e)

需要注意，样本按照不同的时间进行了处理，如下所示：

Note that some samples were processed at different times.

head(sampleInfo)

结果如下所示：

> head(sampleInfo)
  ethnicity       date        filename
1       CEU 2003-02-04 GSM25349.CEL.gz
2       CEU 2003-02-04 GSM25350.CEL.gz
3       CEU 2002-12-17 GSM25356.CEL.gz
4       CEU 2003-01-30 GSM25357.CEL.gz
5       CEU 2003-01-03 GSM25358.CEL.gz
6       CEU 2003-01-16 GSM25359.CEL.gz

日期是一个无关的变量，它理论上不影响这个数据集中基因的表达值，然而，正如我们在前面分析中看到的那样，这个日期似乎也起了一些作用，因此我们在这里就来讲一下这个日期的因素。
我们可以明显看到，日期和种族这两个因素几乎完全混杂了：

year = factor( format(sampleInfo$date,"%y") )
tab = table(year,sampleInfo$ethnicity)
print(tab)

结果如下所示：

> print(tab)
    
year ASN CEU HAN
  02   0  32   0
  03   0  54   0
  04   0  13   0
  05  80   3   0
  06   2   0  24

通过t检验以及生成的火山图我们可以发现，不同种族之间似乎存在着数千个差异基因。但是，当我们仅仅对2002年和2003年的CEU人群进行比较发现，我们又发现了数千个差异基因：

library(genefilter)
##remove control genes
out <- grep("AFFX",rownames(geneExpression))
eth <- sampleInfo$ethnicity
ind<- which(eth%in%c("CEU","ASN"))
res1 <- rowttests(geneExpression[-out,ind],droplevels(eth[ind]))
ind <- which(year%in%c("02","03") & eth=="CEU")
res2 <- rowttests(geneExpression[-out,ind],droplevels(year[ind]))
XLIM <- max(abs(c(res1$dm,res2$dm)))*c(-1,1)
YLIM <- range(-log10(c(res1$p,res2$p)))
mypar(1,2)
plot(res1$dm,-log10(res1$p),xlim=XLIM,ylim=YLIM,
     xlab="Effect size",ylab="-log10(p-value)",main="Populations")
plot(res2$dm,-log10(res2$p),xlim=XLIM,ylim=YLIM,
     xlab="Effect size",ylab="-log10(p-value)",main="2003 v 2002")
image

练习

P419

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 162,408评论 4赞 371
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 68,690评论 2赞 307
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 112,036评论 0赞 255
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 44,726评论 0赞 221
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 53,123评论 3赞 296
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 41,037评论 1赞 225
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 32,178评论 2赞 318
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 30,964评论 0赞 213
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 34,703评论 1赞 250
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,863评论 2赞 254
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,333评论 1赞 265
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 28,658评论 3赞 263
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 33,374评论 3赞 244
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,195评论 0赞 8
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,988评论 0赞 201
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 36,167评论 2赞 285
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 35,970评论 2赞 279

推荐阅读更多精彩内容