一,分析获取基因水平以及外显子水平的信号
使用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)
您好!之前您写的primeview.db帮助非常大,非常感谢!
于此同时,我发现 在运用nsFilter时 似乎不能载入primeview的环境,不知如何解决,特向您请教。谢谢!!
请问一个粗浅的问题:Affymetrix 外显子(exon)芯片与传统的3‘芯片相比,有没有太多的优势?为什么现今仍有很多研究者使用传统的3‘芯片?为什么传统的3‘芯片没有被淘汰?
这个问题很难回答啊。我感觉原因有以下几点:
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,成本不同。这是最主要人们在研究基因表达水平时倾向使用前者的原因。
收益非浅,衷心感谢热心回复
您好!请问一个问题,我在使用oligo时,输入“affyRaw<-read.celfiles(celFiles)”,为什么总是提示:Warning message:package ‘pd.hugene.1.0.st.v1’ was built under R version 3.2.0 ?
这是版本冲突提示。有些时候可能会导致错误的结果。
您好,老师。首先,谢谢您的博文,内容非常好!
请问在“确定alternative splicing events”这一步中“使用SI/LIMMA”的代码里首行里的 annot.map 之前未定义过,在我运行时会提示错误,这应该怎么改正呢?
现在对于SI/LIMMA中的laternative splicing events的处理,不再使用SI了,请自己查看limma::diffSplice函数,对你会有帮助的。
想请问老师,那个自动化处理流程的密码是什么啊?
已发送,请查收。
请问这一步使用APT获取dabg值如何操作啊?谢谢。
这个我帮不到你了。我已经都忘记怎么操作的了。