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)