PH525 series - Discovering Batch

本章节的主要内容是如何使用探索性数据分析（EDA）检测批次效应，使用的示例数据依旧是基因表达数据。

• 数据准备

library(rafalib)
library(Biobase)
library(GSE5859) ##Available from GitHub 
data(GSE5859)

"##去除实际为同一样本的不同数据(如果一对样本的相关性为1，
那么这必定意味着这两个数据是同一样本的数据以不同名字被上传至了公共数据库)"

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

"##创建去趋势数据，并提取日期与人种信息"
y <- exprs(e)-rowMeans(exprs(e))
dates <- pData(e)$date
eth <- pData(e)$ethnicity

"##在原始的样本信息中未提供性别信息，但是我们可以通过Y染色体上的基因表达量的中位数来判断是男是女"

library(hgfocus.db) ##install from Bioconductor
annot <- select(hgfocus.db, keys=featureNames(e), keytype="PROBEID",
columns=c("CHR")) 
annot <- annot[match(featureNames(e),annot$PROBEID),]
annot$CHR <- ifelse(is.na(annot$CHR),NA,paste0("chr",annot$CHR)) 
##计算Y染色体上基因表达量的平均数
chryexp <- colMeans(y[which(annot$CHR=="chrY"),])
mypar()
hist(chryexp)

202001161413.png

从上图中我们可以看到，男性女性的区别是很明显的，所以我们可以通过chryexp的正负来判断性别：

sex <- factor(ifelse(chryexp<0,"F","M"))

• 方差解释图

作者首先讲了一个概念：结构，大意就是说如果样本之间实际上是互相独立的，那么数据会呈现某种状态，那么术语“结构”（structure）指的就是我们所看到的数据的状态对其的一种偏离。翻译的可能有点不准，所以还是把原文贴上吧：

Here we are using the term structure to refer to the deviation from what one would see if the samples were in fact independent from each other.

决定我们需要多少主成分去描述这种结构的一种简单的探索方法就是绘制方差解释图，如果数据是互相独立的，那么对于主成分的方差解释长这样：

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

202001161521.png

但实际上我们的数据长这样子：

plot(s$d^2/sum(s$d^2))

202001161527.png

从上图中可以看出，至少有20多个点发生了偏离，那么造成这种偏离的原因是什么？人种？性别？日期还是其他什么？

• MDS 图

结合之前所学知识，我们可以通过绘制MDS图来解答这些问题，我们可以通过给不同变量的样本上色的有段探索变量和主成分之间的关系。比如我们用种族变量涂色：

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)

202001161540.png

从图中可以看出，第一主成分和种族之间存在着非常清晰的相关性，但同样的，橘黄色数据点内部存在聚集趋势较为明显的子集。由于种族和数据预处理的时间其实是相关的：

year = factor(format(dates,"%y"))
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

所以这次我们用时间（年份）去涂色，看一看日期是不是造成这种变异的主要原因：

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)

202001161552.png

可以看到，年份这一变量同样与第一主成分非常相关。所以究竟是什么变量造成的这一现象的？

• 为主成分绘制箱型图

因为月份个数比较多，通过涂色去探索数据会比较复杂（不好上色），所以在这里作者通过“月份”这一变量对数据进行分层然后为每一主成分绘制箱型图：

month <- format(dates,"%y%m")
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) }

202001161552.png

可以看到月份这一变量与第一主成分有着很强的相关性，然后再来看一下月份的方差分析结果：

month <- format(dates,"%y%m")
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")

202001161659.png

可以看到，前20个主成分与月份之间的相关性都是比较强的，然后我们还可以用F检验比较一下月份内与月份间的方差:

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

202001161712.png

阅读原文请戳