bioconductor系列教程之二简单的Affymetrix 外显子(exon)分析流程(上) 9

搜索替换剪切(Searching for Alternative Splicing)

搜索方法其实有很多,这会在以后的教程当中介绍。这里就只讲讲splicing index方法的实现。

我们假设一个基因在样品A当中表达量是样品B当中的两倍,那么参与该基因的所有外显子(exon)都会在样品A当中侦测量是样品B中的两倍。如下图中所示样品1当中基因水平的表达量差异是12倍左右,如果针对其中一组探针组,它的外显子表达强度(exon intensity, NI)就是探针组的强度除以基因水平的表达强度。得到NI值之后,我们就可以使用NI的对数比值来求得Splicing Index(SI)值了。

splicing index算法

splicing index算法

这个算法简单吧。但其实它又不简单。我们首先要问的问题就是怎么才能求得这个基因水平的表达强度?Affymetrix公司提供了一种算法,称为iterPLIER,其原理是先用PLIER算一遍所有核心探针组外显子的表达强度,然后对同一基因内的这些外显子做相关性分析,最接近的22组探针组拿出来再用PLIER算法算一遍,然后再将其中最接近的11组拿出来再用PLIER算法算一遍,最后再平均一下,得到该基因的表达强度。

这个办法很不错,但是R当中并没有实现这一算法的现成的程序。要不就直接从Affymetrix公司下载他们的应用软件。但那已经超出了我希望的范围。那我们有没有其它的办法呢?其实,基因水平的表达强度这个问题转换一下,就是外显子强度的平均水平。有文章表明,其实对于大多数基因来讲,探针组的选择对最后基因表达定量影响不大,但对于个别基因,就影响很大。我们考虑到有替换剪切的存在,那么假设同一基因假设对于同一基因其表达了的显子的强度应该一致。所以我们可以简单地对外显子进行一下分类处理,如果同一基因下的外显子表达量很接近,那么我们就以它们的平均值来代表其基因水平的表达强度。为了实现它,我们对核心探针组规定如下算法:

如果同一基因的外显子探针,我们对它的表达强度进行K-means clustering分析,只分两组,一组是在该基因中表达了的,另一组是在该基因中没有表达的。把其中高的那个平均值直接当成基因表达定量。

显然这种算法相当的粗糙,因为实际上即使在相同的样品内也存在替换剪切这一过程,而并不是简单的要么表达要么不表达这样子。这里只是示例,我就不严格了吧。这个算法显然是费时费力的。所以我们需要缩小一下范围,我们只对我们感趣的基因做一下这种分析,也就是我们之前选择出来的limma.sigs.filtered中外显子涉及的基因进行一下这类分析。

> limma.sigs.gene.id.unique <- unique(limma.gep.map$gene_id)
> length(limma.sigs.gene.id.unique)
[1] 2568
> x.pli.eset<-exprs(x.pli)
> len<-length(limma.sigs.gene.id.unique)
> sigs.gene.level.exprs<-array(dim=c(len,12))
> colnames(sigs.gene.level.exprs)<-colnames(x.pli.eset)
> rownames(sigs.gene.level.exprs)<-limma.sigs.gene.id.unique
> for(i in 1:len){
+ probesets<-exonic(gene.to.probeset(as.character(limma.sigs.gene.id.unique[i])))
+ j=length(probesets)
+ if(j>2) { j=2 } else { if(j>0) { j=1 } }
+ if(j>0){
+ fit<-kmeans(x.pli.eset[rownames(x.pli.eset) %in% probesets,],j)
+ sigs.gene.level.exprs[as.character(limma.sigs.gene.id.unique[i]),]<-as.numeric(fit$centers[j,])
+ cat(i,"\t")
+ }
+ }

好了,我们拿到了基因水平的表达定量,我们再对每个组做个平均。同时我们准备好外显子的平均值。

> chip_MCF7 = colnames(x.pli.eset)[1:3]
> chip_MCF10A = colnames(x.pli.eset)[4:6]
> chip_100MCF7 = colnames(x.pli.eset)[7:9]
> chip_100MCF10A = colnames(x.pli.eset)[10:12]
> exon.xpr.MCF7<-rowMeans(x.pli.eset[rownames(x.pli.eset) %in% limma.sigs.filtered, chip_MCF7])
> length(exon.xpr.MCF7)
[1] 16569
> exon.xpr.MCF10A<-rowMeans(x.pli.eset[rownames(x.pli.eset) %in% limma.sigs.filtered, chip_MCF10A])
> exon.xpr.100MCF7<-rowMeans(x.pli.eset[rownames(x.pli.eset) %in% limma.sigs.filtered, chip_100MCF7])
> exon.xpr.100MCF10A<-rowMeans(x.pli.eset[rownames(x.pli.eset) %in% limma.sigs.filtered, chip_100MCF10A])
> gene.xpr.MCF7<-rowMeans(sigs.gene.level.exprs[,chip_MCF7])
> gene.xpr.MCF10A<-rowMeans(sigs.gene.level.exprs[,chip_MCF10A])
> gene.xpr.100MCF7<-rowMeans(sigs.gene.level.exprs[,chip_100MCF7])
> gene.xpr.100MCF10A<-rowMeans(sigs.gene.level.exprs[,chip_100MCF10A])

很明显外显子比对应的基因数应该远小于外显子的数目。我们现在需要做的事情就是对每个外显子找到其对应的基因,而后计算其NI值及SI值。

> ni.MCF7<-exon.xpr.MCF7[limma.sigs.gene.s[,"IN1"]]/gene.xpr.MCF7[limma.sigs.gene.s[,"gene_id"]]
> ni.MCF7[1:5]  2316685   2316686   2316728   2318582   2318583 
0.8534844 0.7478648 0.6979099 1.1653441 1.0928776 
> ni.MCF10A<-exon.xpr.MCF10A[limma.sigs.gene.s[,"IN1"]]/gene.xpr.MCF10A[limma.sigs.gene.s[,"gene_id"]]
> ni.100MCF7<-exon.xpr.100MCF7[limma.sigs.gene.s[,"IN1"]]/gene.xpr.100MCF7[limma.sigs.gene.s[,"gene_id"]]
> ni.100MCF10A<-exon.xpr.100MCF10A[limma.sigs.gene.s[,"IN1"]]/gene.xpr.100MCF10A[limma.sigs.gene.s[,"gene_id"]]
> si.ab<-log2(ni.MCF7)-log2(ni.MCF10A)
> si.cd<-log2(ni.100MCF7)-log2(ni.100MCF10A)
> all.si.ab<-cbind(limma.sigs.gene.s,as.matrix(si.ab))
> colnames(all.si.ab)
[1] "IN1"     "gene_id" ""       
> colnames(all.si.ab)[3]<-"SI"
> all.si.cd<-cbind(limma.sigs.gene.s,as.matrix(si.cd))
> colnames(all.si.cd)[3]<-"SI"

Pages: 1 2 3 4 5 6 7 8

9 thoughts on “bioconductor系列教程之二简单的Affymetrix 外显子(exon)分析流程(上)

  1. Pingback: bioconductor系列教程之三–Affymetrix 外显子(exon)分析综合 | 糗世界

  2. Reply s223554 2月 21,2013 3:54 上午

    现在exonmap包不能用了啊,read.exon()有替代吗?

  3. Reply admin 2月 21,2013 6:44 上午

    你可以参考本博的其他教程使用xps,oligo或者aroma,jetta来分析。

  4. Reply wxy 10月 10,2013 5:42 上午

    xmapcore_homo_sapiens_60已下载在本机如何将其添加入mysql服务器。我的系统是win7,mysql也安装在本机上。谢谢。

  5. Reply yylin1984 2月 8,2015 8:12 下午

    您好,当安装”xmapcore”包的时候,提示说不适用于当前版本的R,不知道有没有解决办法?

  6. Reply 7月 24,2015 2:21 上午

    欧老师,又来麻烦您了。这篇博客中的数据下载地址不能打开了,因为我的电脑内存不够,我在按照bioconductor系列教程之二简单的Affymetrix 外显子(exon)分析流程中(质量控制)步骤操作的时候,
    +- rawData/
    | +- /
    | | +- / <– must match exactly a chip type folder under annotationData/
    | | +- CEL files
    中的cel files可以用什么文件代替呢?

Leave a Reply

  

  

  

%d 博主赞过: