bioconductor系列教程之二简单的Affymetrix 外显子(exon)分析流程中（质量控制）

通过上一节的介绍，我们了解了外显子分析的基本流程。这一节，我们关注一下质量分析。

其实直到目前为止，还没有什么很好的关于Affymetrix exon array的质量控制的分析手段。人们还是沿用了之前的一些统计分析方法，具体的来讲，就是看原图(image)，看RLE和NUSE图像（参见：bioconductor系列教程之一分析基因芯片中（质量控制））。在实现上，仍然可以使用之前同样的R语句。

> setwd("~/Documents/exonArray/ExonarraysMCF7andMCF10Adata_cel")
> library(affy)
> library(simpleaffy)
> pd<-read.AnnotatedDataFrame("filelist.txt",header=T,sep="\t",row.names=1)
> pData(pd)
                              group
0306_YH148_EX_exMCF7_r1.CEL       a
0406_YH156_EX_exMCF7_r2.CEL       a
0406_YH157_EX_exMCF7_r3.CEL       a
0306_YH149_EX_exMCF10A_r1.CEL     b
0406_YH158_EX_exMCF10A_r2.CEL     b
0406_YH159_EX_exMCF10A_r3.CEL     b
> data<-ReadAffy(filenames=rownames(pData(pd)),phenoData=pd,verbose=T)
1 reading 0306_YH148_EX_exMCF7_r1.CEL ...instantiating an AffyBatch (intensity a 6553600x6 matrix)...done.
Reading in : 0306_YH148_EX_exMCF7_r1.CEL
Reading in : 0406_YH156_EX_exMCF7_r2.CEL
Reading in : 0406_YH157_EX_exMCF7_r3.CEL
Reading in : 0306_YH149_EX_exMCF10A_r1.CEL
Reading in : 0406_YH158_EX_exMCF10A_r2.CEL
Reading in : 0406_YH159_EX_exMCF10A_r3.CEL
> library("exon.pmcdf")
> data@cdfName<-"exon.pmcdf"
> Pset<-fitPLM(data)
> library(RColorBrewer)
> colors<-brewer.pal(12,"Set3")
> pdf("RLE.pdf")
> Mbox(Pset,ylim=c(-1,1),col=colors,main="RLE")
> dev.off()
> pdf("NUSE.pdf")
> boxplot(Pset,ylim=c(0.9, 1.25),col=colors,main="NUSE")
> dev.off()

但是，使用这样的语句，即使你的机器非常的强大也会消耗很长的时间，就我的机器而言：2.66GHz Intel Core i7处理器，8GB 1067MHz DDR3内存，就处理上次教程一半的数据，就花掉大约两个半小时的时间来运行上面几行代码。

很明显，即使强大的处理能力，对于外显子分析来讲，依然是一个耗时的过程。反正是耗时，我们不如使用支持依机器配置的aroma.affymetrix库。上次我们有提到过这个库，它使用大量的IO操作来降低对内存的要求，你只需要有1.5GB的内存就足以应对外显子的分析了，还是很强大的。但是，时间上，你需要更强的耐心。

首先安装aroma.affymetrix库文件：

> source("http://aroma-project.org/hbLite.R")
> hbInstall("aroma.affymetrix")

而后就是文件的摆放了。aroma-project要求你必须把文件按照指定的格式摆放在指定的文件夹内。我们现在就来研究一下如何放置文件。首先进入工作目录，而后在目录内新建两个文件夹，一个是：annotationData/chipTypes/，一个是：rawData。在chipTypes目录下新建以芯片名为名的目录，比如:HuEx-1_0-st-v2，在rawData下新建以实验名为名的目录，比如qiuworld，而后再在qiuworld下新建与刚刚提到的芯片名一致的目录名HuEx-1_0-st-v2。具体总结如下：

<current working directory>
  +- annotationData/
  |  +- chipTypes/
  |     +- <chipTypeA>/ <-- must match exactly the name of the CDF file (fullname minus tags)
  |     |  +- CDF file(s) and other annotation (possibly subdirectories)
  |     |
  |     +- <chipTypeB>/ <-- must match exactly the name of the CDF file (fullname minus tags)
  |        +- CDF file(s) and other annotation (possibly subdirectories)
  |        ...
  |
  +- rawData/
  |  +- <dataSet1>/
  |  |  +- <chipTypeA>/ <-- must match exactly a chip type folder under annotationData/
  |  |     +- CEL files
  |  |
  |  +- <dataSet2>/
  |  |  +- <chipTypeB>/ <-- must match exactly a chip type folder under annotationData/
  |  |     +- CEL files
  |  |
  |  +- <dataSet3>/
  |  |  +- <chipTypeA>/ <-- must match exactly a chip type folder under annotationData/
  |  |     +- CEL files
  |  ...

接下来，去下载CDF文件至annotationData/chipTypes/HuEx-1_0-st-v2目录下。

我们读入数据：

> setwd("~/Documents/exonArray/qiuworld")
> library(aroma.affymetrix)
> verbose<-Arguments$getVerbose(-8,timestamp=TRUE)
> chipType<-"HuEx-1_0-st-v2"
> cdf<-AffymetrixCdfFile$byChipType(chipType,tags="core,A20071112,EP")
> print(cdf)
AffymetrixCdfFile:
Path: annotationData/chipTypes/HuEx-1_0-st-v2
Filename: HuEx-1_0-st-v2,core,A20071112,EP.cdf
Filesize: 38.25MB
Chip type: HuEx-1_0-st-v2,core,A20071112,EP
RAM: 0.00MB
File format: v4 (binary; XDA)
Dimension: 2560x2560
Number of cells: 6553600
Number of units: 18708
Cells per unit: 350.31
Number of QC units: 1
> cs<-AffymetrixCelSet$byName("qiuworld",cdf=cdf)
> print(cs)
AffymetrixCelSet:
Name: qiuworld
Tags: 
Path: rawData/qiuworld/HuEx-1_0-st-v2
Platform: Affymetrix
Chip type: HuEx-1_0-st-v2,core,A20071112,EP
Number of arrays: 6
Names: 0306_YH148_EX_exMCF7_r1, 0306_YH149_EX_exMCF10A_r1, ..., 0406_YH159_EX_exMCF10A_r3
Time period: 2006-03-24 04:32:49 -- 2006-04-21 13:04:27
Total file size: 378.20MB
RAM: 0.01MB

去噪和标准化

我们知道对于普通的基因芯片的表达定量过程其实有四个步骤：背景较正，如：bgcorrect.method=”rma”,标准化数据，如：normalize.method=”quantiles”, PM较正，如：pmcorrect.method=”pmonly”,以及表达定量，如：summary.method=”medianpolish”。因为外显子芯片没有MM探针，所以不存在PM较正的过程。

这里我们使用RMA算法来进行背景较正。

> bc <- RmaBackgroundCorrection(cs,tag="core")
> csBC <- process(bc,verbose=verbose)

上面的命令中有tag=”core”字样，这是个什么意思呢？当你回到工作目录，你就会看到新建了一个名为probeData的目录，其下有目录，名为qiuworld,core。这里的core就是我们命名的标签。你可以使用多种不能的背景较正和CDF文件来生成不同的probeData文件，其中为了区别彼此，这个tag就非常有用了。

接下来，我们用四分位标准化方法来对数据进行标准化。

> qn <- QuantileNormalization(csBC, typesToUpdate="pm")
> print(qn)
QuantileNormalization:
Data set: qiuworld
Input tags: core
User tags: *
Asterisk ('*') tags: QN
Output tags: core,QN
Number of files: 6 (378.20MB)
Platform: Affymetrix
Chip type: HuEx-1_0-st-v2,core,A20071112,EP
Algorithm parameters: (subsetToUpdate: NULL, typesToUpdate: chr "pm", subsetToAvg: NULL, typesToAvg: chr "pm", .targetDistribution: NULL)
Output path: probeData/qiuworld,core,QN/HuEx-1_0-st-v2
Is done: FALSE
RAM: 0.00MB
> csN <- process(qn,verbose=verbose)
> print(csN)
AffymetrixCelSet:
Name: qiuworld
Tags: core,QN
Path: probeData/qiuworld,core,QN/HuEx-1_0-st-v2
Platform: Affymetrix
Chip type: HuEx-1_0-st-v2,core,A20071112,EP
Number of arrays: 6
Names: 0306_YH148_EX_exMCF7_r1, 0306_YH149_EX_exMCF10A_r1, ..., 0406_YH159_EX_exMCF10A_r3
Time period: 2006-03-24 04:32:49 -- 2006-04-21 13:04:27
Total file size: 378.20MB
RAM: 0.01MB

接下来我们就需要知道基因水平的表达定量了。拿到这个，可以用于差异表达基因的确定。

> plmTr <- ExonRmaPlm(csN, mergeGroups=TRUE)
> print(plmTr)
ExonRmaPlm:
Data set: qiuworld
Chip type: HuEx-1_0-st-v2,core,A20071112,EP
Input tags: core,QN
Output tags: core,QN,RMA,merged
Parameters: (probeModel: chr "pm"; shift: num 0; flavor: chr "affyPLM"; treatNAsAs: chr "weights"; mergeGroups: logi TRUE).
Path: plmData/qiuworld,core,QN,RMA,merged/HuEx-1_0-st-v2
RAM: 0.00MB

当然，我们也可以这个时候也获得外显子水平的表达定量，可以用于alternative splicing分析。

> plmEx <- ExonRmaPlm(csN, mergeGroups=FALSE)
> print(plmEx)
ExonRmaPlm:
Data set: qiuworld
Chip type: HuEx-1_0-st-v2,core,A20071112,EP
Input tags: core,QN
Output tags: core,QN,RMA
Parameters: (probeModel: chr "pm"; shift: num 0; flavor: chr "affyPLM"; treatNAsAs: chr "weights"; mergeGroups: logi FALSE).
Path: plmData/qiuworld,core,QN,RMA/HuEx-1_0-st-v2
RAM: 0.00MB

质量分析

而后我们用线性模型拟合(PLM)所有的数据，并用于质量控制分析。

> fit(plmTr,verbose=verbose)
> rs <-calculateResiduals(plmTr,verbose=verbose)
> ae <- ArrayExplorer(rs)
> setColorMaps(ae,c("log2,log2neg,rainbow","log2,log2pos,rainbow"))
> process(ae,interleaved="auto",verbose=verbose)
> display(ae)
> qamTr <- QualityAssessmentModel(plmTr)
> plotNuse(qamTr)
> plotRle(qamTr)

质量分析也只能是这么多了……还有很多工作可以做的。

接下来我们就很快的再过一点，比如说我们想看一下表达定量的结果，

#先是基因水平的
> cesTr <- getChipEffectSet(plmTr)
> trFit <- extractDataFrame(cesTr, units=1:3, addNames=TRUE)
#再是外显子水平的
> cesEx <- getChipEffectSet(plmEx)
> exFit <- extractDataFrame(cesEx, units=1:3, addNames=TRUE)

OK,而后的工作就是计算FIRMA score了。基本步骤是：

> firma<-FirmaModel(plmTr)
> fit(firma,verbose=verbose)
> fsRes<-extractDataFrame(fs, units=NULL, addNames=TRUE)
> fsLog2<-log2(fsRes[,6:ncol(fsRes)])
> fsLog2<-cbind(fsRes[,1:5], fsLog2)
> write.csv(fsLog2,"FIRMA.score.core.csv", row.names=F)

在介绍FIRMA的文章中指出，rijk = yijk – ci -pk = (实测的探针信号强度) – (基于PLM模型拟合后信号强度的期待值)。所以如果FIRMA的值如果越小的一个负值，就越有可能说明那个exon被splice out了，反之，如果这个FIRMA的值越大，说明其被splice in的可能性就越大。

但是，因为整个过程并没有对探针进行筛选，所以，单独看FIRMA的值并不足以确定它就是一个spilce event.还需要依照affymetrix白皮书中介绍的步骤去除一些假阳性的值，以提高准确率。

16 thoughts on “bioconductor系列教程之二简单的Affymetrix 外显子(exon)分析流程中（质量控制）”

grassule说道：

2014-11-14 8:31 上午

博主你好，我最近学习使用bioconductor分析芯片数据，你的博客我是通过你的红宝书 R语言与bioconductor生物信息学应用上面的介绍来到这里的，本书讲的不错，不过在运行中有时会有各种各样的古怪的问题，通过阅读你的博客部分问题得到了解决，你的博客确实给了我很大的帮助，确实非常感谢。下面有几个问题：
1、我想学习一下你的 bioconductor全自动处理基因芯片CEL文件R程序这篇文章，不知能不能告诉我密码。
2、想向你问一下，我想你可以再自己的网站上发一些视频讲座，使用屏幕录像专家之类的软件，这样既可以看的清楚，也可以将某些问题讲的更清楚，你同时将R 的命令发布到同期的视频上，估计会引起很大的关注。

回复
1. admin说道：
  
  2014-11-17 9:44 上午
  
  很好建议。但是我时间和精力有限，可能短期内都无法实现这一目标了。
  
  回复
波波说道：

2014-12-29 5:41 上午

老师您好！
我也是想学一下您的“ bioconductor全自动处理基因芯片CEL文件R程序”这篇文章，不知道老师能不能给发一下密码。真心感谢老师！！！

回复
1. admin说道：
  
  2014-12-31 11:11 上午
  
  已发送，请查收。
  
  回复
lsf说道：

2015-01-21 1:12 上午

老师：
你好。我最近也在学习用R语言来处理cel文件，想看下您的“ bioconductor全自动处理基因芯片CEL文件R程序”这篇文章，不知道老师能不能给发一下密码。谢谢。

回复
1. admin说道：
  
  2015-01-22 9:16 上午
  
  已发送，请查收。
  
  回复
lsf说道：

2015-01-30 6:30 上午

博主：
麻烦问下你，就是在数据中有两个表头时，我使用read.table和readlines时候，都是读取第一个表头的数据，怎么样才能这两个表头的数据一起读取呢？注：表头有一列是相同的，貌似跟主键差不多。。。

回复
1. admin说道：
  
  2015-02-06 10:04 上午
  
  你的意思是表头占了两行对吧？如果是我的话，我会分开读表头的两行，以及数据部分，然后再将表头处理之后附给数据。
  
  回复
hufanglq说道：

2015-02-12 12:57 上午

博主，您好！当我按照您的方式尝试aroma.affymetrix，用fit(plmTr, verbose = verbose)拟合的时候，出现了错误：
20150211 17:02:51| Number of chunks: 7 (2705 units/chunk)
20150211 17:02:51| Reading, extracting, and writing units…
20150211 17:02:51| Chunk #1 of 7 (2705 units)
20150211 17:02:51| Reading CDF list structure…
20150211 17:02:51| Reading CDF list structure…done
=> RAM: 75MB
错误于(…) : 3个参数给'(‘，但它只需要1个
此外: 警告信息：
In (ns == 0L) & (nchar(unitNames) == 0L) :
长的对象长度不是短的对象长度的整倍数

不知道是怎么回事，麻烦您帮忙看看。

回复
1. admin说道：
  
  2015-02-12 11:52 上午
  
  需要更多的信息来查错。
  
  回复
Gang D说道：

2015-06-05 1:47 下午

亲，你好，想看下您的“ bioconductor全自动处理基因芯片CEL文件R程序”这篇文章，求密码，跪谢！

回复
1. admin说道：
  
  2015-06-08 12:15 下午
  
  已发送，请查收。
  
  回复
zhenzh说道：

2016-10-23 6:27 上午

老师您好，我想看下您的“bioconductor全自动处理基因芯片CEL文件R程序”这篇文章，请问下老师您能发下密码吗？谢谢！

回复
1. admin说道：
  
  2016-10-24 11:47 上午
  
  已发送，请查收。
  
  回复
yangtao说道：

2016-11-02 1:09 上午

老师您好，我想学习您的“bioconductor全自动处理基因芯片CEL文件R程序”这篇文章，希望老师能发我一下密码，麻烦了！

回复
1. admin说道：
  
  2016-11-02 8:16 上午
  
  已发送，请查收。
  
  回复