bioconductor系列教程之二简单的Affymetrix 外显子(exon)分析流程下(综合应用)

3.筛选数据

在进行数据筛选之前,我们必须弄清楚一个理想的处显子替换剪切的数据会是什么样子的:

  1. 基因表达水平相对于背景来说较高,有很好的信噪区分度。
  2. 基因表达水平在所有样品中表现稳定。
  3. 同一基因的大部分外显子信号在组间差异比较小。
  4. 最理想的情况是一个基因中人一个外显子替换剪切,这样可以有效减小分析的复杂度。
  5. 对于外显子替换剪切在组内差异很小,最好是几乎没有,而组间差异明显。

那么什么样子的数据会导致不理想的外显子替换剪切分析结果呢?我们知道,对于依靠显色或者荧光信号强弱来区分翻译或者转录水平之种实验,一个最基本的假设就是信号强度与表达水平成线性关系。但是现实中多数的现实数据都会呈现一个大致的S型关系,也就是说,当信号强度处于的高低两端时还使用线性关系来分析并不能代表其真实的表达水平,而只有信号强度处于合适的量程之内才比较可信。我们如何来区别一个信号的是否在量程之内呢?有人说,这个简单嘛,规定好上限和下限就好了。但是,对于杂交这种实验来讲,如果探针不同,其杂交能力就不同。我们使用摄影的一个概念:宽容度。宽容度是指,感光材料在摄影过程中按正比关系记录景物亮度反差的暴光量范围。比如说一个感光材料可以分辨出十一级(对数级)的光强度差别,那么我们就认为它的宽容度是11。不同的杂交能力就意味着不同的宽容度,也就是对于信号强度上限及下限不同。而实验组间或者组内针对同一探针组,如果其信号分布范围正好部分高出了上限,或者低出了下限,那就是会引起信息损失,会带来错误。

几种可能给出假阳性结果的干扰示例
几种可能给出假阳性结果的干扰示例

为了让我们最终的分析结果可靠,我们要尽量地去除假阳性结果,所以很有必要对经过预处理的结果进一步进行筛选。筛选分成三步

  • 基本筛选(必要):1,基因水平去除只在基中一组当中表达的数据。因为如果其中一组原始信号强度根本无法与背景区别的话,那么经过信号标准化之后,强行将它们的信号拉到可以比较的水平,于是就会出现以噪音为信号的错误(上图左)。2,去除外显子水平中至少在一组中无法从背景中分离的数据。原因同前。因为标准化可能放大噪音(上图右)。
  • 标准筛选(推荐):1,去除外显子水平具有高杂交亲和力的探针组数据。原因同前,因为标准化可能放大噪音(上图中)。2,规定基因水平的背景。3,去除exon/gene信号比过大的数据,这其实就是本类筛选中1的另一种检测标准。4,去除exon/gene信号比过小的数据,这意味着exon信号可能未遵守线性方程,这其实是上类筛选中2的另一种检测标准。
  • 高级筛选(可选):1,去除基因水平表达差异过大的数据,其实是基本筛选中1的另一种检测。2,对外显子的探针组局限至核心探针组。3,只搜索基因水平信号较强的数据。4,增加在预测中感兴趣的探针组的权重。5,去除基因水平生物噪音过小的探针组,也就是杂交信号分步过窄的探针组。6,境加已知的替换剪切的关注度(或者说权重)。

图8中显示的这几组错误,都是因为无法区分信噪而导致的。所以可以考虑使用一些诸如小波变换或者傅立叶转换之类的手段进行预处理以求信噪分离。这个做成的话估计可以写文章发表。而高级筛选中对于权重的设计也可以考虑设定成标准化的流程。

下面是具体怎么做:

基本筛选(必要):2,去除外显子水平中至少在一组中无法从背景中分离的数据

> dabg<-read.table("pm-only.dabg.summary.txt",header=T,sep="\t",dec=".")
> dabg[1:3,]
  probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
1     2590411                      0.04083                        0.13424
2     2609210                      0.52343                        0.67600
3     2616957                      0.00861                        0.00159
  X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
1                      0.02966                      0.03579
2                      0.22479                      0.06348
3                      0.00213                      0.00295
  X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
1                        0.07204                        0.11710
2                        0.43592                        0.18317
3                        0.00193                        0.00239
> pma<-dabg
> nrow(dabg)
[1] 1411399
> for(j in 2:7){
+ index<-which(dabg[,j]<0.05)
+ pma[index,j]<-1
+ pma[-index,j]<-0
+ }
> exon.p<-array(dim=c(nrow(pma),2))
> exon.p[,1]<-rowSums(data.matrix(pma[,c(2,4,5)]))
> exon.p[,2]<-rowSums(data.matrix(pma[,c(3,6,7)]))
> pma[1:5,]
  probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
1     2590411                            1                              0
2     2609210                            0                              0
3     2616957                            1                              1
4     2391254                            0                              0
5     3499280                            0                              0
  X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
1                            1                            1
2                            0                            0
3                            1                            1
4                            0                            0
5                            0                            0
  X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
1                              0                              0
2                              0                              0
3                              1                              1
4                              1                              0
5                              0                              0
> exon.p[1:5,]
     [,1] [,2]
[1,]    3    0
[2,]    0    0
[3,]    3    3
[4,]    0    1
[5,]    0    0
> for(i in c(1,2)){
+ index<-which(exon.p[,i]>1)
+ exon.p[index,i]<-1
+ exon.p[-index,i]<-0
+ }
> exon.p[1:5,]
     [,1] [,2]
[1,]    1    0
[2,]    0    0
[3,]    1    1
[4,]    0    0
[5,]    0    0
> index.exon.p<-which(exon.p[,1]&exon.p[,2])
> index.exon.p[1:10]
 [1]  3  8  9 12 13 16 19 22 24 30
> exons.main.plier.sum<-read.table("plier-gcbg-sketch.summary.txt",header=T,sep="\t",dec=".")
> exons.main.plier.sum.filtered<-exons.main.plier.sum[index.exon.p,]
> nrow(exons.main.plier.sum.filtered)
[1] 507859
> nrow(exons.main.plier.sum)
[1] 1411399

基本筛选(必要):1,基因水平去除只在基中一组当中表达的数据。

> gene.level<-read.table("quant-norm.pm-gcbg.iter-plier.summary.txt",header=T,sep="\t",dec=".")
> colnames(gene.level)
[1] "probeset_id"                    "X0306_YH148_EX_exMCF7_r1.CEL"  
[3] "X0306_YH149_EX_exMCF10A_r1.CEL" "X0406_YH156_EX_exMCF7_r2.CEL"  
[5] "X0406_YH157_EX_exMCF7_r3.CEL"   "X0406_YH158_EX_exMCF10A_r2.CEL"
[7] "X0406_YH159_EX_exMCF10A_r3.CEL"
> nrow(gene.level)
[1] 22011
> nrow(exon.p)
[1] 1411399
> nrow(pma)
[1] 1411399
> rownames(exon.p)<-pma[,1]
> exon.p[1:5,]
        [,1] [,2]
2590411    1    0
2609210    0    0
2616957    1    1
2391254    0    0
3499280    0    0
> core.mps<-read.table("HuEx-1_0-st-v2.r2.dt1.hg18.core.mps",header=T,sep="\t",dec=".")
> gene.level.mps<-merge(gene.level,core.mps,by.x="probeset_id",by.y="probeset_id")
> gene.level.mps[1:3,]
  probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
1     2315554                     376.4905                       445.0123
2     2315633                     411.0840                       514.0779
3     2315674                     234.6524                       157.0173
  X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
1                     337.7216                     334.4438
2                     326.9009                     362.2654
3                     205.3857                     186.6247
  X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
1                       547.5689                       495.2287
2                       314.4923                       304.3485
3                       128.8736                       131.8188
  transcript_cluster_id
1               2315554
2               2315633
3               2315674
                                                                                                                                     probeset_list
1                                                         2315588 2315589 2315591 2315594 2315595 2315596 2315598 2315602 2315603 2315604 2315607 
2                                                                                 2315638 2315639 2315640 2315641 2315642 2315643 2315644 2315645 
3 2315690 2315692 2315697 2315707 2315709 2315713 2315714 2315716 2315718 2315719 2315720 2315725 2315726 2315727 2315728 2315729 2315730 2315731 
  probe_count
1          43
2          32
3          69
> nrow(gene.level)
[1] 22011
> nrow(gene.level.mps)
[1] 22011
> gene.p<-array(dim=c(nrow(gene.level.mps),2))
> for(j in 1:nrow(gene.level.mps)){
+ probeset_list<-unlist(strsplit(as.character(gene.level.mps[j,"probeset_list"]),"\\ "))
+ if(sum(as.numeric(exon.p[probeset_list,1]))>(length(probeset_list)/2)){gene.p[j,1]<-1}else{gene.p[j,1]<-0}
+ if(sum(as.numeric(exon.p[probeset_list,2]))>(length(probeset_list)/2)){gene.p[j,2]<-1}else{gene.p[j,2]<-0}
+ }

标准筛选(推荐):1,去除外显子水平具有高杂交亲和力的探针组数据。如果探针会和多个位点杂交的话,也就是说其特异性较低的话,那么其结果并不可靠。所以我们需要去除这类的探针。而affymatrix公司也对他们的芯片做了评估,对于每个探针,都对其杂交特异性区分为三类,1)高特异性,只和特异性的位点杂交,大约超过80%的探针都属于这一类;2)低特异性,倾向于多个位点杂交;3)混合型,介于前两者之间。

          1 = unique    - All probes in the probe set perfectly match
                          only one sequence in the putatively
                          transcribed array design content. The vast
                          majority (>80%) of probe sets are unique.  

          2 = similar   - All the probes in the probe set perfectly
                          match more than one sequence in the
                          putatively transcribed array design content.

          3 = mixed     - The probes in the probe set either perfectly match or
                          partially match more than one sequence in
                          the putatively transcribed array design content. 

这些信息可以从affymatrix的注释文件当中得到,每个注释文件都会有crosshyb_type列。其实它就是cross-hybridization potential of the probe set。他的值就是上面所示的数字。

> annotation<-read.table("HuEx-1_0-st-v2.r2.dt1.hg18.csv",header=T,sep=",",dec=".")
> dim(annotation)
[1] 1404693      36
> annotation[1:5,]
  probeset_id psr_id exon_cluster_id transcript_cluster_id seqname strand start
1     2315101      1               1               2315100    chr1      +  1788
2     2315102      2               2               2315100    chr1      +  2520
3     2315103      3               3               2315100    chr1      +  2927
4     2315104      4               4               2315100    chr1      +  3229
5     2315105      5               4               2315100    chr1      +  4013
  stop probe_count crosshyb_type number_independent_probes    level bounded
1 2030           4             2                         4 extended       0
2 2555           4             2                         1 extended       0
3 2958           4             3                         1     full       1
4 3751           4             2                         4 extended       0
5 4231           4             2                         2 extended       0
  NoBoundedEvidence fl mrna est vegaGene vegaPseudoGene ensGene sgpGene
1                 0  0    1   2        0              0       0       0
2                 0  0    1   2        0              0       0       0
3                 0  0    0   0        0              0       0       0
4                 0  0    1   0        0              0       0       0
5                 0  0    1   0        0              0       0       0
  exoniphy twinscan geneid genscan genscanSubopt mouse_fl mouse_mrna rat_fl
1        0        0      0       0             0        0          0      0
2        0        0      0       1             1        0          0      0
3        0        0      0       0             1        0          0      0
4        0        0      0       0             0        0          0      0
5        0        0      0       0             0        0          0      0
  rat_mrna microRNAregistry rnaGene mitomap has_cds is_normgene_exon
1        0                0       0       0       0                0
2        0                0       0       0       0                0
3        0                0       0       0       0                0
4        0                0       0       0       0                0
5        0                0       0       0       0                0
  is_normgene_intron
1                  0
2                  0
3                  0
4                  0
5                  0
> crosshyb_type<-annotation[,c(1,10)]
> crosshyb_type[1:5,]
  probeset_id crosshyb_type
1     2315101             2
2     2315102             2
3     2315103             3
4     2315104             2
5     2315105             2
> crosshyb.filtered<-crosshyb_type[crosshyb_type[,2]==1,]
> crosshyb.filtered[1:5,]
   probeset_id crosshyb_type
17     2315127             1
22     2315128             1
24     2315131             1
26     2315133             1
29     2315146             1
> exons.main.plier.sum.filtered[1:5,]
   probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
3      2616957                     66.08712                       91.02172
8      3842613                     63.75740                       30.61370
9      3871433                    142.23688                      112.46411
12     3452435                     37.20707                       44.19509
13     2580386                      9.52387                       19.23181
   X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
3                      75.91781                     76.48327
8                      36.41760                     44.08072
9                     121.33124                    125.38457
12                     12.31370                     38.58102
13                      5.35404                     16.48731
   X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
3                       105.39617                       91.21448
8                        47.89569                       73.11166
9                       174.84539                      204.31540
12                       42.15767                       17.42073
13                        5.49444                       17.81327
> exons.main.plier.sum.filtered.1<-merge(exons.main.plier.sum.filtered,crosshyb.filtered,by.x="probeset_id",by.y="probeset_id")
> exons.main.plier.sum.filtered.1[1:5,]
  probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
1     2315133                     24.91217                       38.93058
2     2315165                    419.12063                      479.39810
3     2315201                     98.01441                      178.25163
4     2315356                     28.38951                       44.84634
5     2315358                     10.72641                       13.53168
  X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
1                     25.37772                     69.32379
2                    312.38263                    299.15070
3                     53.59437                     80.99036
4                     58.10497                     54.67421
5                     21.59089                     20.03721
  X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL crosshyb_type
1                       46.70432                       46.66366             1
2                      243.39897                      307.17719             1
3                       78.63091                       50.00538             1
4                       62.39389                       44.07864             1
5                       27.99601                       19.34531             1
> exons.main.plier.sum.filtered.1$crosshyb_type<-NULL
> exons.main.plier.sum.filtered.1[1:5,]
  probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
1     2315133                     24.91217                       38.93058
2     2315165                    419.12063                      479.39810
3     2315201                     98.01441                      178.25163
4     2315356                     28.38951                       44.84634
5     2315358                     10.72641                       13.53168
  X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
1                     25.37772                     69.32379
2                    312.38263                    299.15070
3                     53.59437                     80.99036
4                     58.10497                     54.67421
5                     21.59089                     20.03721
  X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
1                       46.70432                       46.66366
2                      243.39897                      307.17719
3                       78.63091                       50.00538
4                       62.39389                       44.07864
5                       27.99601                       19.34531
> nrow(exons.main.plier.sum.filtered.1)
[1] 374888
> nrow(exons.main.plier.sum.filtered)
[1] 507859

标准筛选(推荐):2,规定基因水平的背景。这个过程上分为三步,1)设定背景噪音域值,大于域值的设定为P;2)如果一组样品中超过半数的值为P,标记该组为P;3)如果两组都为P则保留,否则过滤掉。

> summary(gene.level.filtered)
  probeset_id      X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
 Min.   :2315674   Min.   :3.300e-04            Min.   :3.000e-04             
 1st Qu.:2790626   1st Qu.:5.220e+01            1st Qu.:4.658e+01             
 Median :3309112   Median :2.137e+02            Median :2.067e+02             
 Mean   :3233009   Mean   :6.433e+02            Mean   :6.880e+02             
 3rd Qu.:3707383   3rd Qu.:6.941e+02            3rd Qu.:7.829e+02             
 Max.   :7385683   Max.   :1.530e+04            Max.   :1.534e+04             
 X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
 Min.   :2.700e-04            Min.   :3.200e-04           
 1st Qu.:5.029e+01            1st Qu.:5.195e+01           
 Median :2.171e+02            Median :2.198e+02           
 Mean   :6.111e+02            Mean   :6.219e+02           
 3rd Qu.:6.892e+02            3rd Qu.:7.056e+02           
 Max.   :1.610e+04            Max.   :1.568e+04           
 X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
 Min.   :4.000e-04              Min.   :1.800e-04             
 1st Qu.:4.602e+01              1st Qu.:4.652e+01             
 Median :1.993e+02              Median :1.986e+02             
 Mean   :5.913e+02              Mean   :5.941e+02             
 3rd Qu.:6.770e+02              3rd Qu.:6.836e+02             
 Max.   :1.475e+04              Max.   :1.379e+04 
> threshold<-5
> gene.level.filtered.p<-array(dim=c(nrow(pma),6))
> for(j in 2:7){
+ index<-which(gene.level.filtered[,j]>threshold)
+ gene.level.filtered.p[index,j-1]<-1
+ gene.level.filtered.p[-index,j-1]<-0
+ }
> gene.level.filtered.pp<-array(dim=c(nrow(gene.level.filtered),2))
> gene.level.filtered.pp[,1]<-rowSums(data.matrix(gene.level.filtered.p[,c(1,3,4)]))
> gene.level.filtered.pp[,2]<-rowSums(data.matrix(gene.level.filtered.p[,c(2,5,6)]))
> gene.level.filtered.pp[1:5,]
     [,1] [,2]
[1,]    3    3
[2,]    3    3
[3,]    3    3
[4,]    2    1
[5,]    3    3
> for(i in c(1,2)){
+ index<-which(gene.level.filtered.pp[,i]>1)
+ gene.level.filtered.pp[index,i]<-1
+ gene.level.filtered.pp[-index,i]<-0
+ }
> gene.level.filtered.pp[1:5,]
     [,1] [,2]
[1,]    1    1
[2,]    1    1
[3,]    1    1
[4,]    1    0
[5,]    1    1
> index.gene.level.filtered.pp<-which(gene.level.filtered.pp[,1]&gene.level.filtered.pp[,2])
> gene.level.filtered.1<-gene.level.filtered[index.gene.level.filtered.pp,]
> gene.level.filtered.1[1:5,]
   probeset_id X0306_YH148_EX_exMCF7_r1.CEL X0306_YH149_EX_exMCF10A_r1.CEL
2      3063807                    111.45444                       87.23597
4      2949118                    175.25352                      256.74699
6      4014076                     22.78219                      101.42396
9      3719161                    673.76710                      925.63541
10     2916345                   2592.87403                      544.81952
   X0406_YH156_EX_exMCF7_r2.CEL X0406_YH157_EX_exMCF7_r3.CEL
2                     136.65633                     79.67514
4                     277.17128                    247.13024
6                      18.26672                     17.13998
9                     716.63370                    647.56955
10                   2388.37326                   2402.37036
   X0406_YH158_EX_exMCF10A_r2.CEL X0406_YH159_EX_exMCF10A_r3.CEL
2                         71.7941                       81.91057
4                        262.8664                      242.69974
6                        128.8538                      139.77980
9                        885.6453                      944.66380
10                       515.2965                      466.86786
> nrow(gene.level.filtered.1)
[1] 12544
> nrow(gene.level.filtered)
[1] 13129
> write.table(gene.level.filtered.1,"CC.genes-core.i-plier.sum.txt",sep="\t",row.names=F,quote=F)
> write.table(exons.main.plier.sum.filtered.1,"CC.exons-main.plier.sum.txt",sep="\t",row.names=F,quote=F)
cel_files	group_id
X0306_YH148_EX_exMCF7_r1.CEL	brain
X0306_YH149_EX_exMCF10A_r1.CEL	heart
X0406_YH156_EX_exMCF7_r2.CEL	brain
X0406_YH157_EX_exMCF7_r3.CEL	brain
X0406_YH158_EX_exMCF10A_r2.CEL	heart
X0406_YH159_EX_exMCF10A_r3.CEL	heart

使用MiDAS来判断alternative splicing。

../bin/apt-midas --cel-files cel_ids.txt -g CC.genes-core.i-plier.sum.txt -e CC.exons-main.plier.sum.txt -m HuEx-1_0-st-v2.r2.dt1.hg18.full.mps -o out -pvalues -fstats -normalized -stabilize 8.0

将结果读入R来排个序。

> midas <-read.table("midas.pvalues.txt", sep="\t", header=T)
> head(midas)
> o<-order(midas$pvalue)
> midas.ordered<- midas[o,]

至此,整个筛选工作结束。

发表评论

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