bioconductor系列教程之三–Affymetrix 外显子(exon)分析综合 10

一,分析获取基因水平以及外显子水平的信号

使用R库oligo分析获取基因水平以及外显子水平的信号

library(oligo)
celFiles<-list.celfiles()#读取当前目录下所有的cel文件名
library(pd.huex.1.0.st.v2)#调取注释文件
affyRaw<-read.celfiles(celFiles)#读取CEL文件
dir.create("png")#创建一个目录,用于保存质量分析图片
for(i in 1:9){png(paste("png",i,"png",sep="."));image(affyRaw,i,col=gray((64:0)/64));dev.off();}#生成信号强度灰度图
pdf("png/hist.pdf")#生成柱状图
hist(affyRaw)
dev.off()
for(i in 1:9){png(paste("png/maplot",i,"png",sep=".")); MAplot(affyRaw,which=i,ylim=c(-1,1)); dev.off();}#生成原始信号MA图
png(paste("png/boxplot.png"))#生成原始信号箱线图
boxplot(affyRaw)
dev.off()
probesetSummaries<-rma(affyRaw, target="probeset")#计算外显子水平信号
geneSummaries<-rma(affyRaw,target="core")#计算基因水平信号
png("png/boxplot.rma.png")#生成RMA之后的箱线图
boxplot(geneSummaries,transfo=identity,main="After RMA")
dev.off()
pdf("png/hist.rma.pdf")#生成RMA之后的柱状图
hist(geneSummaries,transfo=identity,main="After RMA")
dev.off()
d.gene<-exprs(geneSummaries)#获得基因表达的数据
d.exon<-exprs(probesetSummaries)#获得外显子表达的数据

使用R库xps分析获取基因水平以及外显子水平的信号

library(xps)
dir<-"/path/to/cel/data"#保存工作目录,因为xps在读取数据时可能会改变工作目录,所以需要随时将其改回来。
setwd(dir)
celfiles<-dir(pattern=".cel",ignore.case=T)#查找工作目录下所有的CEL文件
libdir <- "/path/to/library/HuEx-1_0-st-v2"#设置芯片注释库文件的目录
scheme.huex10stv2 <- import.exon.scheme("Scheme_HuEx10stv2",
                                layoutfile=paste(libdir,"HuEx-1_0-st-v2.r2.clf",sep="/"),
                                schemefile=paste(libdir,"HuEx-1_0-st-v2.r2.pgf",sep="/"),
                                probeset=paste(libdir,"HuEx-1_0-st-v2.na31.hg19.probeset.csv",sep="/"),
                                transcript=paste(libdir,"HuEx-1_0-st-v2.na31.hg19.transcript.csv",sep="/"))#读入注释库。xps采用了root,会将数据写入一个以.root结尾的文件当中,以便将来再次调用。这里就会在工作目录下生成一个名为Scheme_HuEx10stv2.root的文件。
celnames<-gsub("\\(HuEx-1_0-st-v2\\).CEL","",celfiles)#将colname做一个修改
data <- import.data(scheme.huex10stv2,filename="tmpdt_Data",filedir=dir,celdir=dir,celfiles=celfiles,celnames=celnames)#使用之前准备好的注释文件,读入CEL数据。
data.rma.ps<-rma(data,"HuExonRMAcore_ps",background="antigenomic",normalize=T,option="probeset",exonlevel="core")#在exon水平上做RMA,需要注意的是,xps不象其它的包那样对rma的结果做log2的转变,所以如果需要使用这里的值,将来还需要做一步log2 transform.
xps.rma.ps<-validData(data.rma.ps)#在上一步中,数据同时被写入一个名为HuExonRMAcore_ps.root的文件,方便将来调用。
data.rma.tc <- rma(data,"HuExonRMAcore_tc",background="antigenomic",normalize=T,option="transcript",exonlevel="core")#在gene水平上做RMA
xps.rma.tc<-validData(data.rma.tc)#数据被validData调用后才能正常操作,否则会显示为空数据
d.gene<-log2(xps.rma.tc)#获得基因表达的数据
d.exon<-log2(xps.rma.ps)#获得外显子表达的数据

使用R库aroma.affymetrix分析获取基因水平以及外显子水平的信号
请一定先看一遍之前的博客内容,否则很难顺利实现这里的代码。

library(aroma.affymetrix)#调入库
verbose<-Arguments$getVerbose(-8,timestamp=TRUE)#设置返回信息时间格式
chipType<-"HuEx-1_0-st-v2"#设置芯片类型
cdf<-AffymetrixCdfFile$byChipType(chipType,tags="core,A20071112,EP")#设置芯片CDF文件
cs<-AffymetrixCelSet$byName("projectname_qiuworld",cdf=cdf)#设置CEL数据位置
bc <- RmaBackgroundCorrection(cs,tag="core")#设置RMA背噪参数
csBC <- process(bc,verbose=verbose)#DO RMA
qn <- QuantileNormalization(csBC, typesToUpdate="pm")#四分位标准化参数,其实是标准RMA的一个步骤
csN <- process(qn,verbose=verbose)#DO normalization
plmTr <- ExonRmaPlm(csN, mergeGroups=TRUE)#设置使用PLM模型获取基因水平表达数据
fit(plmTr,verbose=verbose)#获取基因水平表达数据
plmEx <- ExonRmaPlm(csN, mergeGroups=FALSE)#设置使用PLM模型获取外显子水平表达数据
fit(plmEx,verbose=verbose)#获取外显子水平表达数据
rs <-calculateResiduals(plmTr,verbose=verbose)#计算残基,用于绘制质量控制的NUSE图和RLE图
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)#NUSE图
plotRle(qamTr)#RLE图
cesTr <- getChipEffectSet(plmTr)#准备读入基因水平表达值
trFit <- extractDataFrame(cesTr, units=NULL, addNames=TRUE)#读入数据
cesEx <- getChipEffectSet(plmEx)#准备读入外显子表达值
exFit <- extractDataFrame(cesEx, units=NULL, addNames=TRUE)#读入数据
d.gene <- log2(trFit[,6:ncol(trFit)])#对数据做log2转换
row.names(d.gene)<-as.character(trFit[,1])#设置probeset_id为行号
d.exon <- log2(exFit[,6:ncol(trFit)])#对数据做log2转换
row.names(d.exon)<-as.character(exFit[,2])#设置probeset_id为行号

使用Affymetrix Power Tools分析获取基因水平以及外显子水平的信号

apt-probeset-summarize –a rma-sketch –p HuEx-1_0-st-v2.r2.pgf -c HuEx-1_0-st-v2.r2.clf -m HuEx-1_0-st-v2.r2.core.mps –qc-probesets HuEx-1_0-st-v2.r2.qcc -o OUT_GENE *.CEL#获取基因水平表达值
apt-probeset-summarize -a rma-sketch –p HuEx-1_0-st-v2.r2.pgf -c HuEx-1_0-st-v2.r2.clf -s HuEx-1_0-st-v2.r2.core.ps –qc-probesets HuEx-1_0-st-v2.r2.qcc -o OUT_EXON *.CEL#获取外显子水平表达值

得到表达值后读入R。

d.exon <- read.table("./OUT_EXON/rma-sketch_core_exon.summary.txt",sep="\t", header=T, row.names=1)
d.gene <- read.table("./OUT_GENE/rma-sketch.summary_core_gene.txt",sep="\t", header=T, row.names=1)

使用JETTA分析获取基因水平以及外显子水平的信号
需要注意的是JETTA将使用完全不同的方法来处理,所以它不会进入使用bioconductor去除假阳性信号,而是直接使用它自身的MADS和ASA方法来确定alternative splicing events.

library(jetta)
LIB = "/path/to/library/HuEx-1_0-st-v2/"
CELs = "/path/to/CELs/"
CLF = paste(LIB, "HuEx-1_0-st-v2.r2.clf",sep="")
PGF = paste(LIB, "HuEx-1_0-st-v2.r2.pgf",sep="")
PSA = paste(LIB, "HuEx-1_0-st-v2.na31.hg19.probeset.csv",sep="")
param.tr = jetta.do.expression(CELs,CLF,PGF,ProbesetAnnotationFile=PSA,SummarizationMethod="median_polish",SummarizationTarget='tc')
tcexpr = jetta.read.expression(param.tr)

Pages: 1 2 3 4

10 thoughts on “bioconductor系列教程之三–Affymetrix 外显子(exon)分析综合

  1. Reply Joe 6月 11,2012 4:03 上午

    您好!之前您写的primeview.db帮助非常大,非常感谢!

    于此同时,我发现 在运用nsFilter时 似乎不能载入primeview的环境,不知如何解决,特向您请教。谢谢!!

  2. Reply keke 8月 19,2012 12:17 上午

    请问一个粗浅的问题:Affymetrix 外显子(exon)芯片与传统的3‘芯片相比,有没有太多的优势?为什么现今仍有很多研究者使用传统的3‘芯片?为什么传统的3‘芯片没有被淘汰?

    • Reply admin 8月 19,2012 6:48 上午

      这个问题很难回答啊。我感觉原因有以下几点:
      1,设计不同。3‘IVT芯片(基因芯片)的对象是基因3’末端(覆盖大约3‘端600bp),每个基因只有1至2组探针组(probeset)。所以它的样品一般是以polyA为引物产生的cDNA。而exon array(外显子芯片)覆盖整个基因的所有外显子,其样品是mRNA,而后是以不同引物(而不是单一的polyA)逆转录产生的cDNA。所以基因芯片的研究对象只能是基因表达水平,而外显子芯片则可能是transcripts水平。要注意的是,基因芯片面对同一基因不同结尾时,会区别对待。
      2,实验步骤不同。正如前面提到的,它们主要的区别在于逆转录的引物。
      3,探针组成不同。基因芯片针的探针组由11PM和11MM探针组成,而外显子芯片则是4PM。平均每个基因,基因芯片有16个探针,而外显子芯片会有147个探针。虽然探针的总数相对于基因芯片来说提高了,但是其精确程度却很依赖相互验证,这在实验水平来说很难保证。这就意味着对其的分析会更难。
      4,成本不同。这是最主要人们在研究基因表达水平时倾向使用前者的原因。

  3. Reply keke 8月 19,2012 6:57 上午

    收益非浅,衷心感谢热心回复

  4. Reply Yuan Xun 1月 7,2015 3:36 上午

    您好!请问一个问题,我在使用oligo时,输入“affyRaw<-read.celfiles(celFiles)”,为什么总是提示:Warning message:package ‘pd.hugene.1.0.st.v1’ was built under R version 3.2.0 ?

  5. Reply David 7月 19,2015 9:54 下午

    您好,老师。首先,谢谢您的博文,内容非常好!
    请问在“确定alternative splicing events”这一步中“使用SI/LIMMA”的代码里首行里的 annot.map 之前未定义过,在我运行时会提示错误,这应该怎么改正呢?

    • Reply admin 7月 21,2015 10:26 下午

      现在对于SI/LIMMA中的laternative splicing events的处理,不再使用SI了,请自己查看limma::diffSplice函数,对你会有帮助的。

  6. Reply CCUICUCPR 3月 27,2017 12:45 下午

    想请问老师,那个自动化处理流程的密码是什么啊?

Leave a Reply

  

  

  

%d 博主赞过: