通过上一节的介绍,我们了解了外显子分析的基本流程。这一节,我们关注一下质量分析。
其实直到目前为止,还没有什么很好的关于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白皮书中介绍的步骤去除一些假阳性的值,以提高准确率。
博主你好,我最近学习使用bioconductor分析芯片数据,你的博客我是通过你的红宝书 R语言与bioconductor生物信息学应用上面的介绍来到这里的,本书讲的不错,不过在运行中有时会有各种各样的古怪的问题,通过阅读你的博客部分问题得到了解决,你的博客确实给了我很大的帮助,确实非常感谢。下面有几个问题:
1、 我想学习一下你的 bioconductor全自动处理基因芯片CEL文件R程序 这篇文章,不知能不能告诉我密码。
2、想向你问一下,我想你可以再自己的网站上发一些视频讲座,使用屏幕录像专家之类的软件,这样既可以看的清楚,也可以将某些问题讲的更清楚,你同时将R 的命令发布到同期的视频上,估计会引起很大的关注。
很好建议。但是我时间和精力有限,可能短期内都无法实现这一目标了。
老师您好!
我也是想学一下您的“ bioconductor全自动处理基因芯片CEL文件R程序”这篇文章,不知道老师能不能给发一下密码。真心感谢老师!!!
已发送,请查收。
老师:
你好。我最近也在学习用R语言来处理cel文件,想看下您的“ bioconductor全自动处理基因芯片CEL文件R程序”这篇文章,不知道老师能不能给发一下密码。谢谢。
已发送,请查收。
博主:
麻烦问下你,就是在数据中有两个表头时,我使用read.table和readlines时候,都是读取第一个表头的数据,怎么样才能这两个表头的数据一起读取呢?注:表头有一列是相同的,貌似跟主键差不多。。。
你的意思是表头占了两行对吧?如果是我的话,我会分开读表头的两行,以及数据部分,然后再将表头处理之后附给数据。
博主,您好!当我按照您的方式尝试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) :
长的对象长度不是短的对象长度的整倍数
不知道是怎么回事,麻烦您帮忙看看。
需要更多的信息来查错。
亲,你好,想看下您的“ bioconductor全自动处理基因芯片CEL文件R程序”这篇文章,求密码,跪谢!
已发送,请查收。
老师您好,我想看下您的“bioconductor全自动处理基因芯片CEL文件R程序”这篇文章,请问下老师您能发下密码吗?谢谢!
已发送,请查收。
老师您好,我想学习您的“bioconductor全自动处理基因芯片CEL文件R程序”这篇文章,希望老师能发我一下密码,麻烦了!
已发送,请查收。