DALS025-批次效應(yīng)01-什么是批次效應(yīng)

title: DALS025-批次效應(yīng)01-什么是批次效應(yīng)
date: 2019-08-25 12:0:00
type: "tags"
tags:

  • 批次效應(yīng)
    categories:
  • Data Analysis for the life sciences

前言

這一部分是《Data Analysis for the life sciences》的第10章批次效應(yīng)的第1小節(jié)，這一部分的主要內(nèi)容涉及批次效應(yīng)(Batch Effects)的介紹，有關(guān)批次效應(yīng)的前言部分Rmarkdown文檔參見作者的Github。

什么是批次效應(yīng)

高通量研究中一個(gè)經(jīng)常被忽視的問題就是批次效應(yīng)(batch effects)，批次效應(yīng)受到當(dāng)時(shí)檢測的實(shí)驗(yàn)室條件、試劑批次和人員差異的影響。當(dāng)批次效應(yīng)與我們目標(biāo)結(jié)果混淆并愜以不正確的結(jié)果時(shí)，這就成了一個(gè)主要問題。在這一章中，我們將詳細(xì)地描述批次效應(yīng)：對于批次效應(yīng)如何檢測、解釋、建模和調(diào)整。

批次效應(yīng)是基因組學(xué)研究中面臨的最大挑戰(zhàn)，尤其是在精確醫(yī)學(xué)這個(gè)背景下。大多數(shù)情況（但并非全部）下，高通量技術(shù)已經(jīng)被報(bào)道存在著一種形式或另外一種形式的批次效應(yīng)[Leek et al. (2010) Nature Reviews Genetics 11, 733-739]。但是，批次效應(yīng)并非基因組學(xué)所特有的。實(shí)際上， 在1972年一篇文獻(xiàn)中，Mj Youden就在對物理常數(shù)的經(jīng)驗(yàn)估計(jì)的背景下提到了批次效應(yīng)。他指出，物理常數(shù)“存在著主觀估計(jì)的特征”，以及如何在不同實(shí)驗(yàn)室之間變化的。例如，在下面的Table1中，Youden就展示了來自于不同實(shí)驗(yàn)室又的天文單的估計(jì)值。這些報(bào)告包括對離差(spread)（現(xiàn)在我們稱之為置信區(qū)間）的估計(jì)。

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

從不同實(shí)驗(yàn)室之間的可變性以及報(bào)道的界限（可以理解為置信區(qū)間）來看，都不能解釋這種可變性，這就清楚地表明不同實(shí)驗(yàn)室之間，而非某個(gè)實(shí)驗(yàn)室內(nèi)部存在著某些效應(yīng)。這種類型的變異就是我們所謂的批次效應(yīng)。請注意，有些實(shí)驗(yàn)室報(bào)告了兩個(gè)估計(jì)值（紫色和橙色），這樣我們就在相同的實(shí)驗(yàn)室也看到了批次效應(yīng)。
我們可以使用統(tǒng)計(jì)學(xué)符號來精確地描述這個(gè)問題。我們使用下面的公式來表示這些檢測值：
Y_{i,j} = \mu + \varepsilon_{i,j}, j=1,\dots,N

其中，Y_{i,j} 表示第 i 個(gè)實(shí)驗(yàn)室的第 j 個(gè)檢測值，而 \mu 表示真實(shí)的物理常數(shù)，\varepsilon_{i,j} 表示獨(dú)立的檢測誤差。為了解釋可變性，我們引入了 \varepsilon_{i,j} ，隨后我們根據(jù)數(shù)據(jù)來計(jì)算標(biāo)準(zhǔn)誤。就像本書前面提到的那樣，我們使用 N 值均值來估計(jì)物理物理常數(shù)，如下所示：
\bar{Y}_i = \frac{1}{N} \sum_{i=1}^{N} Y_{i,j}

再來構(gòu)建一個(gè)置信區(qū)間：
\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

但是，這個(gè)置信區(qū)間太小，它無法覆蓋批次效應(yīng)的變異，下面是一個(gè)更加合適的模型：
Y_{i,j} = \mu + \gamma_i + \varepsilon_{i,j}, j=1, \dots, N

其中 \gamma_i 表示一個(gè)實(shí)驗(yàn)室特定的偏差或批次效應(yīng)(batch effect)。從圖片上我們可以明顯看出來，實(shí)驗(yàn)室之間 \gamma 的變化要遠(yuǎn)大于一個(gè)實(shí)驗(yàn)室內(nèi)部 \varepsilon 的變化。用統(tǒng)計(jì)學(xué)的術(shù)語來描述這個(gè)問題就是，\mu 和 \gamma 無法被識別。我們可以估計(jì) \mu_i+\gamma_i ，但是無法區(qū)分開。我們可以將 \gamma 視為一個(gè)隨機(jī)變量。在這個(gè)案例中，每個(gè)實(shí)驗(yàn)室都有一個(gè)錯(cuò)誤項(xiàng) \gamma_i ，這一項(xiàng)在貫穿于該實(shí)驗(yàn)室的所有檢測值，在每個(gè)檢測值中，這一項(xiàng)是相同的，但是實(shí)驗(yàn)室與實(shí)驗(yàn)室間的這一項(xiàng)則不同。因此，這個(gè)問題就可以用以下方程表示：
s_i / \sqrt{N} \mbox{ with } s_i^2= \frac{1}{N-1} \sum_{i=1}^N (Y_{ij} - \bar{Y}_i)^2

這是對標(biāo)準(zhǔn)差的低估，因?yàn)樗鼪]有解釋由 \gamma 導(dǎo)致的實(shí)驗(yàn)室內(nèi)的相關(guān)性（這一句不懂，原文如下：

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.

如果我們假設(shè) \gamma=0 ，那么利用來自幾個(gè)實(shí)驗(yàn)室的數(shù)據(jù)，我們實(shí)際上可以估計(jì)出 \gamma ?；蛘哒f我們可以將它們視為隨機(jī)效應(yīng)，簡單地將它們當(dāng)成一個(gè)新的估計(jì)值，并使用所有的檢測值來計(jì)算其標(biāo)準(zhǔn)差。以下是我們將報(bào)告中的均值值視為隨機(jī)觀察值的置信區(qū)間：

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)

結(jié)果如下所示：

> 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的論文還包括最近關(guān)于光速以及策略常數(shù)的批次效應(yīng)估計(jì)的案例。在這一章里，我們只展示高通量生物數(shù)據(jù)中批次效應(yīng)的廣泛性和復(fù)雜性。

混雜

書中有一個(gè)術(shù)語，即confounding，這里譯為混雜。這一部分內(nèi)容的Rmarkdown文檔可以參考作者的Github。

當(dāng)批次效應(yīng)與我們的目標(biāo)結(jié)果混雜時(shí)，就會導(dǎo)致嚴(yán)重的??。這里我們描述一下混雜(confounding)，以及它與我們數(shù)據(jù)解釋的關(guān)系。

我們從這本書或者說從任何其它數(shù)據(jù)分析課程中嘗到最重要的思想之一就是“相關(guān)不等于因果”。這句話多數(shù)情況就是真的，一個(gè)常見的案例就是混雜(confounding)。簡單地講，當(dāng)我們觀察到 X 和 Y 之間存在著相關(guān)(correlation)或關(guān)聯(lián)(association)時(shí)，往往就會存在著混雜，但嚴(yán)格來說，這是因?yàn)?X 和 Y 都依賴于一個(gè)無關(guān)的變量 Z 。這里我們會描述一個(gè)Simposon悖論，這是一個(gè)基于一個(gè)著名的法律案件的案例，接著，我們還會提到一個(gè)在高通量生物學(xué)研究中的一個(gè)混雜案例。

案例之Simpson悖論

加州大學(xué)伯克分校(UCB)1973年的入學(xué)數(shù)據(jù)顯示，在錄取的學(xué)生中，44%是男性，30%是女性，顯著男性更多，以下是這些數(shù)據(jù)：

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

結(jié)果如下所示：

> 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

卡方檢驗(yàn)的結(jié)果明確拒絕零假設(shè)，即錄取的性別是獨(dú)立的，兩者不受影響（也就是說，男女錄取公平），如下所示：

##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

但經(jīng)過更仔細(xì)的觀察會發(fā)現(xiàn)一個(gè)自相矛盾的結(jié)果，以下是按專業(yè)劃分后，不同性別的錄取百分比：

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

結(jié)果如下所示：

> 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

從結(jié)果中我們可以發(fā)現(xiàn)，不同專業(yè)之間沒有性別偏見。

但是，我們在前面使用了卡方檢驗(yàn)發(fā)現(xiàn)，入學(xué)和性別之間存在著某種關(guān)系。然而，當(dāng)我們的數(shù)據(jù)按照不同專業(yè)進(jìn)行分組時(shí)，這種依賴性似乎消失了，這是怎么一回事呢？

這就是Simpson悖論的一個(gè)案例。

我們上面進(jìn)行的卡方檢驗(yàn)表明，入學(xué)和性別之間存在依賴關(guān)系。然而，當(dāng)數(shù)據(jù)按專業(yè)分組時(shí)，這種依賴性似乎消失了。到底怎么回事？現(xiàn)在我們繪制一個(gè)能顯示申請專業(yè)的人，以及最終進(jìn)入這個(gè)專業(yè)讀書的人的百分比的圖形，用于說明男女在錄取比例方面的問題，如下所示：

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

從圖片上我們可以發(fā)現(xiàn)，男生更想申請那些“容易”的專業(yè)。也就是說，男生這個(gè)變量與“容易”專業(yè)這個(gè)變量發(fā)生了混雜。

混雜的圖形解釋

在這一部分里，我們將混雜圖形化。在下面的圖形中，每個(gè)字母表示一個(gè)學(xué)生。被錄取的人用綠色表示，字母表示專業(yè)。在第1張圖中，所有相同性別的學(xué)生都被放在一起，我們可以發(fā)現(xiàn)，男生中綠色的比較更大，如下所示：

###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

現(xiàn)在我們按照專業(yè)對不同性別的學(xué)生進(jìn)行分級。這里我們要注意，大多數(shù)被錄取的學(xué)生（綠色）來源于兩個(gè)容易的專業(yè)，即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

分層后的均值

在下圖中，我們可以看到，我們根據(jù)專業(yè)設(shè)定條件或分層后，然后再看差異，也就是說，當(dāng)我們控制了混雜項(xiàng)后，這就混雜效應(yīng)就消失了，如下所示：

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

事實(shí)上，在不同專業(yè)錄取的學(xué)生中，性別差異僅表示為，男生比女生高3.5%，如下所示：

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

結(jié)果如下所示：

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

棒球中的Simpson悖論

在棒球比賽的統(tǒng)計(jì)中我們也經(jīng)常看到Simpson悖論，現(xiàn)在我們來看一個(gè)有名的案例，在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)

這里的混雜項(xiàng)就是比賽場次，Jeter的比賽場次數(shù)目更多，但是結(jié)果卻是Justice擊球率更高。

混雜之高通量案例

為了描述我們在生物學(xué)中遇到的混雜問題，我們將會使用《Common genetic variants account for differences in gene expression among ethnic groups》這篇文獻(xiàn)中的數(shù)據(jù)集，這篇文獻(xiàn)指出，兩個(gè)不同種族的血液大約有50%的差異基因，現(xiàn)在我們來下載這個(gè)數(shù)據(jù)集：

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

我們使用Bioconductor中的函數(shù)exprs和pData就能提取基因的表達(dá)數(shù)據(jù)與樣本信息，如下所示：

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

需要注意，樣本按照不同的時(shí)間進(jìn)行了處理，如下所示：

Note that some samples were processed at different times.

head(sampleInfo)

結(jié)果如下所示：

> 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

日期是一個(gè)無關(guān)的變量，它理論上不影響這個(gè)數(shù)據(jù)集中基因的表達(dá)值，然而，正如我們在前面分析中看到的那樣，這個(gè)日期似乎也起了一些作用，因此我們在這里就來講一下這個(gè)日期的因素。
我們可以明顯看到，日期和種族這兩個(gè)因素幾乎完全混雜了：

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

結(jié)果如下所示：

> 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檢驗(yàn)以及生成的火山圖我們可以發(fā)現(xiàn)，不同種族之間似乎存在著數(shù)千個(gè)差異基因。但是，當(dāng)我們僅僅對2002年和2003年的CEU人群進(jìn)行比較發(fā)現(xiàn)，我們又發(fā)現(xiàn)了數(shù)千個(gè)差異基因：

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

練習(xí)

P419

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容