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

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

其实直到目前为止,还没有什么很好的关于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)

NUSE

NUSE


RLE

RLE

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

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

#先是基因水平的
> 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)分析流程中(质量控制)

  1. Reply grassule 11月 14,2014 8:31 上午

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

  2. Reply 波波 12月 29,2014 5:41 上午

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

  3. Reply lsf 1月 21,2015 1:12 上午

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

  4. Reply lsf 1月 30,2015 6:30 上午

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

    • Reply admin 2月 6,2015 10:04 上午

      你的意思是表头占了两行对吧?如果是我的话,我会分开读表头的两行,以及数据部分,然后再将表头处理之后附给数据。

  5. Reply hufanglq 2月 12,2015 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) :
    长的对象长度不是短的对象长度的整倍数

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

  6. Reply Gang D 6月 5,2015 1:47 下午

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

  7. Reply zhenzh 10月 23,2016 6:27 上午

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

  8. Reply yangtao 11月 2,2016 1:09 上午

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

Leave a Reply

  

  

  

%d 博主赞过: