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

确定差异表达的探针组

Exonmap提供了几个很简单而又不错的函数:pc(), fc()和tt()。pc()可以比较任意两组之间的表达差异(基于对数的标准化的值),fc(),tt()可以提取pc()计算后的表达倍数变化(fold-changes)和p值。
一般的,如果我们会考虑设置fc()的绝对值大于1(这意味着2倍的表达差异),tt()的值小于10-4

> ab.pc.rma <- pc(x.rma,"group",c("a","b"))
> cd.pc.rma <- pc(x.rma,"group",c("c","d"))
> ab.fc.pc.rma <- fc(ab.pc.rma)
> cd.fc.pc.rma <- fc(cd.pc.rma)
> ab.tt.pc.rma <- tt(ab.pc.rma)
> cd.tt.pc.rma <- tt(cd.pc.rma)
> summary(ab.fc.pc.rma)
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-10.430000  -0.326100  -0.005890  -0.001951   0.319100  10.740000 
> summary(ab.tt.pc.rma)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.1077  0.3607  0.4010  0.6656  1.0000 
> summary(cd.fc.pc.rma)
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-11.290000  -0.343400   0.014070  -0.001101   0.360300  10.460000 
> summary(cd.tt.pc.rma)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.1292  0.3799  0.4149  0.6773  1.0000 
> keep <- (abs(ab.fc.pc.rma)>1)&(ab.tt.pc.rma<1e-4)
> ab.sigs<- featureNames(x.rma)[keep]
> length(ab.sigs)
[1] 13262
> keep <-(abs(cd.fc.pc.rma)>1)&(cd.tt.pc.rma<1e-4)
> cd.sigs<- featureNames(x.rma)[keep]
> length(cd.sigs)
[1] 6501

当然,我们也可以使用limma来计算差异表达的探针组。和之前的教程不同的是,这里计算的是探针组,而不是基于基因的。但实际上limma并不区分这些。

> library(limma)
> design<-model.matrix(~ -1+factor(c(rep(c(1:4),each=3))))
> colnames(design) <- c("A","B","C","D")> design
   A B C D
1  1 0 0 0
2  1 0 0 0
3  1 0 0 0
4  0 1 0 0
5  0 1 0 0
6  0 1 0 0
7  0 0 1 0
8  0 0 1 0
9  0 0 1 0
10 0 0 0 1
11 0 0 0 1
12 0 0 0 1
attr(,"assign")
[1] 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$`factor(c(rep(c(1:4), each = 3)))`
[1] "contr.treatment"

> cont.matrix<-makeContrasts(AvsB=A-B,CvsD=C-D,levels=design)
> fit<-lmFit(x.pli,design)
> fit1<-contrasts.fit(fit,cont.matrix)
> fit2<-eBayes(fit1)
> result<- decideTests(fit2,lfc=1)
> summary(result)
      AvsB    CvsD
-1   23222   22045
0  1364822 1369460
1    23145   19684
> keep <- result@".Data" !=0
> summary(keep)
    AvsB            CvsD        
 Mode :logical   Mode :logical  
 FALSE:1364822   FALSE:1369460  
 TRUE :46367     TRUE :41729    
 NA's :0         NA's :0  

> limma.ab.sigs <- rownames(result[keep[,'AvsB'],])
> length(limma.ab.sigs)
[1] 46367
> limma.cd.sigs <- rownames(result[keep[,'CvsD'],])
> length(limma.cd.sigs)
[1] 41729

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 博主赞过: