二,使用bioconductor去除假阳性信号
在这里我们需要使用之前获得的基因水平和外显子水平的表达值,依照《下》中所引述的affymetrix白皮书中的几条标准来去除假阳性。
基本筛选(必要):1,基因水平去除只在基中一组当中表达的数据。因为如果其中一组原始信号强度根本无法与背景区别的话,那么经过信号标准化之后,强行将它们的信号拉到可以比较的水平,于是就会出现以噪音为信号的错误。2,去除外显子水平中至少在一组中无法从背景中分离的数据。原因同前。因为标准化可能放大噪音。
而为了做到这一步,我们需要得到detection above background (DABG)的值。在xps中也有直接得到xps值的方法,但这里还是只讲述使用APT获取dabg值,因为xps只是在R中实现相同的算法而已。
apt-probeset-summarize -a dabg -p HuEx-1_0-st-v2.r2.pgf -c HuEx-1_0-stv2.r2.clf -b HuEx-1_0-st-v2.r2.antigenomic.bgp -o ./OUT_DABG *.CEL |
我们将在OUT_DABG中得到名为dabg.summary.txt的文件。
之前我们得到的表达值无论是基因水平还是外显子水平,都是使用的RMA算法,并且将它放在d.gene以及d.exon当中。如果是使用的PLIER方法的话,我们需要对其值进行log2转换
d.exon <- log2(d.exon+16)#加上16是为了防止出现负值 d.gene <- log2(d.gene+16)
以下会分几个步骤依次以core probesets, undetected probesets以及corss-hybridizing probesets标准来去除假阳性。在这里我们需要假设我们拥有10个CEL文件,其中前5个为对照,后5个为实验组。
ctrl=c(1:5) expr=c(6:10) half=3 numberOfSample=10 #########filter for core probesets########### dabg <- read.table("./OUT_DABG/dabg.summary.txt", sep="\t", header=T,row.names=1)#读取DABG值 dabg.core <- dabg[match(row.names(d.exon), row.names(dabg)),]#只保留和exon结果对应的部分 #########filter for undetected probesets########## count.det <- function(x){length(which(x<0.05))}#定义一个函数,用于统计DABG结果中P<0.05的个数,因为我们只保留其个数大于50%的样品数的值 group1.det <- apply(dabg.core[,ctrl], 1, count.det)#对对照组计算 group2.det <- apply(dabg.core[,expr], 1, count.det)#对实验组计算 x <- sort(union(which(group1.det>=half), which(group2.det>=half)))#保留对照组以及实验组当中任何一组中认为其结果是高于背景的。 d.exon.fil <- d.exon[x,]#对结果正式过滤 #########filter for cross-hybridizing probesets########### annot <- read.table("HuEx-1_0-st-v2.na31.hg19.probeset.csv", sep=",", header=T)#读取注释文件 annot.core <- annot[match(row.names(d.exon.fil), annot[,1]),]#只保留和exon结果对应的部分 keep <- which(annot.core$crosshyb_type==1)#只保留cross-hybridizing值等于1的 ids <- annot.core[keep,1]#获取需要保留的probeset_id d.exon.fil2 <- d.exon.fil[match(ids, row.names(d.exon.fil)),]#保留相应的exon结果并保存 write.table(d.exon.fil2, "./OUT_EXON/rma-sketch.summary_core_exon_filtered.txt", sep="\t", quote=F,row.names=T) #########filter for genes undetected in both groups########### keep <- intersect(row.names(dabg.core), annot[,1])#只保留有注释的部分 dabg.core2 <- dabg.core[match(keep, row.names(dabg.core)),] dabg.core2[,numberOfSample+1] <- annot[match(row.names(dabg.core2), annot$probeset_id), "transcript_cluster_id"]#附加基因水平ID gene.ids <- levels(as.factor(dabg.core2[,ncol(dabg.core2)]))#获取基因水平ID count.exon.det <- function(x){length(which(x<0.05))/length(x)}#定义函数,统计DABG值超半数 gene.detection <- matrix(nrow=length(unique(dabg.core2[,ncol(dabg.core2)])), ncol=numberOfSample)#预先设计表格来统计一个样品一个transcript_cluster内有是否有半数可以被检测到 row.names (gene.detection) <- gene.ids#设置表格行号为基因水平ID colnames(gene.detection) <- colnames(d.gene)#设置表格列号为样品名 for(i in 1:numberOfSample)#进行统计 { gene.detection[,i] <- tapply(dabg.core2[,i], dabg.core2[,ncol(dabg.core2)],count.exon.det) } count.gene.det <- function(x){length(which(x>=0.5))}#定义一个函数,用于统计同一组样品中有多少超半数是可以被检测到的 genes.det.group1 <- apply(gene.detection[,ctrl], 1, count.gene.det)#对对照组计算 genes.det.group2 <- apply(gene.detection[,expr], 1, count.gene.det)#对实验组计算 keep.genes <- which((genes.det.group1>=half)&(genes.det.group2>=half))#只保留两组都被检测到的情况,否则会认为是基因水平表达差异,而非外显子水平的差异。 keep.gene.ids <- row.names (gene.detection)[keep.genes]#读取需要保留的基因ID y <- match(keep.gene.ids, row.names(d.gene)) y <- y[-which(is.na(y)=="TRUE")]#去除空值 d.gene.fil <- d.gene[y,]#保留需要保留的基因,并保存 write.table(d.gene.fil, "./OUT_GENE/rma-sketch.summary_core_gene_filtered.txt", sep="\t", quote=F,row.names=T)
您好!之前您写的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值如何操作啊?谢谢。
这个我帮不到你了。我已经都忘记怎么操作的了。