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

二,使用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)

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

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

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

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

    1. 这个问题很难回答啊。我感觉原因有以下几点:
      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. 您好!请问一个问题,我在使用oligo时,输入“affyRaw<-read.celfiles(celFiles)”,为什么总是提示:Warning message:package ‘pd.hugene.1.0.st.v1’ was built under R version 3.2.0 ?

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

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

发表评论

电子邮件地址不会被公开。 必填项已用*标注