title: DALS026-批次效應(yīng)02-發(fā)現(xiàn)批次效應(yīng)
date: 2019-08-26 12:0:00
type: "tags"
tags:
- 批次效應(yīng)
categories: - Data Analysis for the life sciences
前言
這一部分是《Data Analysis for the life sciences》的第10章批次效應(yīng)的第2小節(jié)，這一部分的主要內(nèi)容涉及批次效應(yīng)(Batch Effects)的發(fā)現(xiàn)與處理，有關(guān)這一部分Rmarkdown文檔參見(jiàn)作者的Github。
前面我們已經(jīng)介紹過(guò)了PCA，我們?cè)趯?shí)際運(yùn)用過(guò)程中就可以使用PCA來(lái)進(jìn)行探索性數(shù)據(jù)分析了。為了說(shuō)明這一點(diǎn)，我們將會(huì)使用一個(gè)實(shí)際中的數(shù)據(jù)集，出于我們教學(xué)的目的，這個(gè)數(shù)據(jù)還是原始數(shù)據(jù)，并沒(méi)有經(jīng)過(guò)數(shù)據(jù)清洗。我們首先從公共數(shù)據(jù)庫(kù)中下載這個(gè)原始數(shù)據(jù)，接著你要做的就是預(yù)處理這些數(shù)據(jù)，并使用Bioconductor中的一個(gè)包進(jìn)行后續(xù)的數(shù)據(jù)分析。
基因表達(dá)數(shù)據(jù)
第一步，下載數(shù)據(jù)，如下所示：
library(rafalib)
library(Biobase)
# devtools::install_github("genomicsclass/GSE5859")
library(GSE5859) ##Available from GitHub
data(GSE5859)
注：這里書(shū)本中有所出入，
第二步，數(shù)據(jù)探索。這里先從探索樣本相關(guān)性矩陣開(kāi)始，我們注意到有一對(duì)樣本的相關(guān)性是1.這就問(wèn)題時(shí)著有一個(gè)樣本被上傳到了公共數(shù)據(jù)庫(kù)2次，但是上傳后的名稱不同，下面的代碼將會(huì)剔除這個(gè)樣本，如下所示：
cors <- cor(exprs(e))
Pairs=which(abs(cors)>0.9999,arr.ind=TRUE)
out = Pairs[which(Pairs[,1]<Pairs[,2]),,drop=FALSE]
if(length(out[,2])>0) e=e[,-out[2]]
還可以移除控制組的探針名，如下所示：
out <- grep("AFFX",featureNames(e))
e <- e[-out,]
第三步，處理數(shù)據(jù)?，F(xiàn)在我們創(chuàng)建一個(gè)去趨勢(shì)(detrended)基因表達(dá)數(shù)據(jù)矩陣，并從樣本注釋表中提取日期與結(jié)果，如下所示：
y <- exprs(e)-rowMeans(exprs(e))
dates <- pData(e)$date
eth <- pData(e)$ethnicity
原始數(shù)據(jù)集中不包括樣本信息中的性別。但是，我們?yōu)榱顺鲇诮虒W(xué)的目的，我們使用下面的數(shù)據(jù)來(lái)研究一下這個(gè)性別問(wèn)題，也就是說(shuō)，在下面的代碼中，我們會(huì)展示一下如何預(yù)測(cè)每個(gè)樣本的性別。這種計(jì)算的基本思想就是觀察Y染色體中基因的中位數(shù)基因表達(dá)水平。男性的這個(gè)數(shù)字應(yīng)該更高。為此，我們需要上傳一個(gè)注釋包，用于提供這個(gè)實(shí)驗(yàn)中所用平臺(tái)的一些特征信息，如下所示：
annotation(e)
結(jié)果如下所示：
> annotation(e)
[1] "hgfocus"
此時(shí)，我們需要下載并安裝hgfocus.db包，然后提取染色體位置的信息，如下所示：
map2gene <- mapIds(hgfocus.db, keys=featureNames(e),
column="ENTREZID", keytype="PROBEID",
multiVals="first")
library(Homo.sapiens)
map2chr <- mapIds(Homo.sapiens, keys=map2gene,
column="TXCHROM", keytype="ENTREZID",
multiVals="first")
chryexp <- colMeans(y[which(unlist(map2chr)=="chrY"),])
如果我們創(chuàng)建Y染色體中位數(shù)基因表達(dá)水平的直方圖，我們就能清楚地看到差異，如下所示：
mypar()
hist(chryexp)
現(xiàn)在我們預(yù)測(cè)一下性別，如下所示：
sex <- factor(ifelse(chryexp<0,"F","M"))
計(jì)算主成分
現(xiàn)在我們計(jì)算一下主成分，如下所示：
s <- svd(y)
dim(s$v)
結(jié)果如下所示：
> s <- svd(y)
> dim(s$v)
[1] 207 207
我們也可以使用prcomp()函數(shù)來(lái)創(chuàng)建一個(gè)主成分對(duì)象。
變異解釋
第一步就是解釋由結(jié)構(gòu)(structure)誘導(dǎo)的樣本的相關(guān)性程度是什么樣的，如下所示：
library(RColorBrewer)
cols=colorRampPalette(rev(brewer.pal(11,"RdBu")))(100)
image ( cor(y) ,col=cols,zlim=c(-1,1))
上圖顯示了相關(guān)性。每個(gè)小單元絡(luò)的i與j代表了樣本i和樣本j的相關(guān)性。紅色表示高，白色表示0，藍(lán)色表示低。
這時(shí)我們使用了術(shù)語(yǔ)結(jié)構(gòu)(structure)，結(jié)構(gòu)是指，如果樣本實(shí)際上是相互獨(dú)立時(shí)我們所能發(fā)現(xiàn)的偏差。從上面圖片上我們可以明顯看到，樣本實(shí)際上有分組，也就是說(shuō)，不同組內(nèi)的樣本相關(guān)性比組外的更強(qiáng)。
我們需要生成一個(gè)簡(jiǎn)單的探索性圖來(lái)確定我們?多少主成分來(lái)描述這種結(jié)構(gòu)，這就是方差解釋圖。 如果數(shù)據(jù)是獨(dú)立的，那么方差的解釋就如下所示：
y0 <- matrix( rnorm( nrow(y)*ncol(y) ) , nrow(y), ncol(y) )
d0 <- svd(y0)$d
plot(d0^2/sum(d0^2),ylim=c(0,.25))
事實(shí)上，我們看到的方差解釋圖是下面的這個(gè)樣子：
plot(s$d^2/sum(s$d^2))
至少有20個(gè)PC高于我們對(duì)獨(dú)立數(shù)據(jù)的預(yù)期。我們下一步就是嘗試用檢測(cè)的變量來(lái)解釋這些PC。這些PC是由種族(ethnicity)，性別(sex)還是時(shí)間(date)或其它的因素驅(qū)動(dòng)的嗎？
MDS圖
如前面一樣，我們可以先用MDS圖來(lái)探索一下數(shù)據(jù)，回答上述提到的問(wèn)題。探索感興趣的變量和PC之間的關(guān)系的一種方法就是使用顏色來(lái)表示這些變量，例如，以下是使用表示了顏色來(lái)標(biāo)記種族(ethnicity)這個(gè)變量，通過(guò)前兩個(gè)PC來(lái)展示的數(shù)據(jù)結(jié)果：
cols = as.numeric(eth)
mypar()
plot(s$v[,1],s$v[,2],col=cols,pch=16,
xlab="PC1",ylab="PC2")
legend("bottomleft",levels(eth),col=seq(along=levels(eth)),pch=16)
從上圖來(lái)看，第1 PC與種族(ethnicity)有著很強(qiáng)的聯(lián)系。但是，我們也會(huì)看到一些橙色的點(diǎn)有的形成了亞簇(subslusters)。我們從以前的預(yù)分析中也知道，種族(ethnicity)和預(yù)處理時(shí)間(preprocessing date)有相關(guān)性，如下所示：
year = factor(format(dates,"%y"))
table(year,eth)
結(jié)果如下所示：
> table(year,eth)
eth
year ASN CEU HAN
02 0 32 0
03 0 54 0
04 0 13 0
05 80 3 0
06 2 0 23
因此，我們通過(guò)與前面相同的圖來(lái)研究一下，日期作為主要變異來(lái)源的可能性，在下面的圖形中，我們使用顏色來(lái)表示年，如下所示：
cols = as.numeric(year)
mypar()
plot(s$v[,1],s$v[,2],col=cols,pch=16,
xlab="PC1",ylab="PC2")
legend("bottomleft",levels(year),col=seq(along=levels(year)),pch=16)
從上面圖形我們看到，日期（以年為單位）也與第1PC非常相關(guān)。那么到底是哪個(gè)變量促進(jìn)了這種情況呢？由于存在著很高的混雜(confounding)效應(yīng)，現(xiàn)在還不清楚是哪個(gè)因素起了重要作用。但是，在解決這個(gè)問(wèn)題方面，我們還會(huì)進(jìn)一步探索數(shù)據(jù)。
PC箱線圖
樣本之間的相關(guān)性的圖形可以展示出一個(gè)復(fù)雜的結(jié)構(gòu)，它似乎有5個(gè)以的因素（其中一個(gè)是年份）參與其中。很明顯，這種復(fù)雜程度遠(yuǎn)非種族(ethnicity)這一個(gè)因素能解釋的。我們此時(shí)還探索一下月份之間的相關(guān)性，如下所示：
month <- format(dates,"%y%m")
length( unique(month))
variable <- as.numeric(month)
mypar(2,2)
for(i in 1:4){
boxplot(split(s$v[,i],variable),las=2,range=0)
stripchart(split(s$v[,i],variable),add=TRUE,vertical=TRUE,pch=1,cex=.5,col=1)
}
一共有21個(gè)月份，下圖是按照不同月份劃分場(chǎng)次層后，根據(jù)第1PC到第4PC來(lái)繪制的箱線圖，如下所示：
從上圖我們可以看到，月份與第1PC有著非常強(qiáng)烈的相關(guān)性，即使按照種族(ethnicity)和其它因素來(lái)劃分層次，也是如此。我們要知道，2002-2004之間處理的樣本都來(lái)自于同一種族群體。在這樣的情況下，我們可以使用方差分析來(lái)查看一下PC與哪個(gè)月份相關(guān)，如下所示：
corr <- sapply(1:ncol(s$v),function(i){
fit <- lm(s$v[,i]~as.factor(month))
return( summary(fit)$adj.r.squared )
})
mypar()
plot(seq(along=corr), corr, xlab="PC")
從上圖我們可以看到，第1 PC的相關(guān)性非常強(qiáng)，而對(duì)于前20左右的PC，相關(guān)性較強(qiáng)。我們還可以計(jì)算一下月份內(nèi)與月份之間的F統(tǒng)計(jì)量，如下所示：
Fstats<- sapply(1:ncol(s$v),function(i){
fit <- lm(s$v[,i]~as.factor(month))
Fstat <- summary(aov(fit))[[1]][1,4]
return(Fstat)
})
mypar()
plot(seq(along=Fstats),sqrt(Fstats))
p <- length(unique(month))
abline(h=sqrt(qf(0.995,p-1,ncol(s$v)-1)))
上圖展示的是，方差分析的F平方根，用于解釋PC與月份的關(guān)系。
至此為止，我們就了解了如何使用PCA聯(lián)合EDA來(lái)做為一個(gè)強(qiáng)大的功能檢測(cè)并理解批次效應(yīng)的。在后面的部分里，我們將會(huì)了解如何使用PC用于因子分析(factor analysis)，用于改進(jìn)模型估計(jì)。
練習(xí)
P431
