edgeR案例学习:deepSAGE分析 6

接到学习任务,要求搞清楚两个R包,edgeR和REDseq。这一篇就主要讲述edgeR中的案例之一,对应edgeR包中说明文件的案例12

要了解更多关于SAGE(Serial Analysis of Gene Expression),请阅读DeepSAGE—digital transcriptomics with high sensitivity, simple experimental protocol and multiplexing of samples

安装及调用edgeR

> source("http://bioconductor.org/biocLite.R")
BiocInstaller version 1.4.3, ?biocLite for help
> biocLite("edgeR")
BioC_mirror: http://bioconductor.org
Using R version 2.15, BiocInstaller version 1.4.3.
Installing package(s) 'edgeR'
trying URL 'http://www.bioconductor.org/packages/2.10/bioc/bin/macosx/leopard/contrib/2.15/edgeR_2.6.0.tgz'
Content type 'application/x-gzip' length 1558515 bytes (1.5 Mb)
opened URL
==================================================
downloaded 1.5 Mb
 
 
The downloaded binary packages are in
	/var/folders/Dj/Dj+bWjS7HxiNJ0kYFeKdTE+++TI/-Tmp-//RtmpcteSIZ/downloaded_packages
> library(edgeR)
Loading required package: limma

我们需要做准备一个名为targets.txt的文本文件,用于对样品的说明。接下来的工作就是读取数据。

> library(GEOquery)
Loading required package: Biobase
Loading required package: BiocGenerics
Welcome to Bioconductor
 
    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.
 
Setting options('download.file.method.GEOquery'='auto')
> gset<-getGEO("GSE10782",GSEMatrix = F)
File stored at: 
/var/folders/Dj/Dj+bWjS7HxiNJ0kYFeKdTE+++TI/-Tmp-//RtmpBoLGWD/GSE10782.soft.gz
Parsing....
Found 9 entities...
GPL9185 (1 of 9 entities)
GSM272105 (2 of 9 entities)
GSM272106 (3 of 9 entities)
GSM272318 (4 of 9 entities)
GSM272319 (5 of 9 entities)
GSM272320 (6 of 9 entities)
GSM272321 (7 of 9 entities)
GSM272322 (8 of 9 entities)
GSM272323 (9 of 9 entities)
There were 50 or more warnings (use warnings() to see the first 50)
> slotNames(gset)
[1] "header" "gsms"   "gpls"  
> x<-slot(gset,"gsms")
> setwd("~/Documents/shRNA/edgeR")
> for(i in 1:length(x)){
+	name<-Accession(x[[i]]);
+	data<-Table(dataTable(x[[i]]));
+	write.table(data,paste(name,".txt",sep=""),quote=F,sep="\t",row.names=F)
+ }
> targets<-read.delim("targets.txt",stringsAsFactors=F)
> targets
          files group                          description
1 GSM272105.txt  DCLK transgenic (Dclk1) mouse hippocampus
2 GSM272106.txt    WT          wild-type mouse hippocampus
3 GSM272318.txt  DCLK transgenic (Dclk1) mouse hippocampus
4 GSM272319.txt    WT          wild-type mouse hippocampus
5 GSM272320.txt  DCLK transgenic (Dclk1) mouse hippocampus
6 GSM272321.txt    WT          wild-type mouse hippocampus
7 GSM272322.txt  DCLK transgenic (Dclk1) mouse hippocampus
8 GSM272323.txt    WT          wild-type mouse hippocampus
> d<-readDGE(targets)
> colnames(d)
[1] "1" "2" "3" "4" "5" "6" "7" "8"
> colnames(d)<-c("DCLK1","WT1","DCLK2","WT2","DCLK3","WT3","DCLK4","WT4")


我们现在已经构建了一个DEGList——d,它里面有targets文件中提供的一共八个样品的数据。这时,我们可以接着就分析其中有差异的部分了。然而,d当中包含了80万个特异的短序,这对于我们来说,虽然完整,但其中过多的都是只出现了几次的短序,这会大大增加最终的结果中的false positive的比例。所以我们需要对数据进行过滤。过滤的办法就是将出现次数少于5次/Million reads的短序都过滤掉。edgeR提供的cpm函数可以很方便地做到这一点。

> cpm.d <- cpm(d)
> d<-d[rowSums(cpm.d>1)>=4,]
> dim(d)
[1] 44881     8

现在数据已经完全准备好了。
我们还是需要看一点数据的质量,看看是否同一组内的样品可以聚类,做一个简单的PCA分析。

> plotMDS(d,xlim=c(-2,1))

MDS plot

MDS plot


从图中我们可以看出其wildtype基本上堆在一起,而DCLK的基本上在一起,所以我们可以认为样品数据组间区分度良好。

然后估计分散程度。至于为什么要估计分散程度,在Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation一文当中有详细的阐述。要义就是对于每个样品中的每个基因的表达情况被检测到都是有一个分散程度的。如果我们不区别对待基因的分散程度的话,就会增加false positive的结果。而在edgeR: a Bioconductor package for differential expression analysis of digital gene expression data一文中有一个比较简约的解释。对于任一样品i中的任一个基因g,它的分布假设符合负二项式分布(negative binomial (NB) distribution)。

Ygi~NB(Mipgjg)

其中Mi是库的大小,就是样品i中reads的总数;φg就是基因< i>g的分散程度;pgji样品所属的j组中基因g 的相对丰度。对于NB分布,期待值µgi=Mipgj,而variance为µgi(1+µgiφg),这两个值都是可以计算的,于是我们可以得到我们感兴趣的pgj值。这样就可以比较不同组之间的表达差异了。当然,在φg趋近于0的时候,也就是样品无穷大,而分散度却很小的时候,其分布就符合泊森分布(Poisson distribution)了。

> d<-estimateCommonDisp(d)
> names(d)
[1] "samples"           "counts"            "common.dispersion" "pseudo.alt"        "conc"              "common.lib.size"  
> d$common.dispersion
[1] 0.1514558
> sqrt(d$common.dispersion)
[1] 0.3891732
> d <- estimateTagwiseDisp(d)
> names(d)
[1] "samples"            "counts"             "common.dispersion"  "pseudo.alt"         "conc"               "common.lib.size"    "prior.n"            "tagwise.dispersion"
> plotBCV(d)
> abline(h=sqrt(d$common.dispersion), col="firebrick", lwd=3)

biological coefficient of bariation

biological coefficient of bariation of abundance


接下来的步骤就是获取差异表达基因了。

> et <- exactTest(d,pair=c("WT","DCLK"))
> topTags(et,n=5)
Comparison of groups:  DCLK-WT 
                      logFC   logCPM       PValue          FDR
TCTGTACGCAGTCAGGC  9.400251 5.362400 9.727206e-25 4.365667e-20
CATAAGTCACAGAGTCG  9.825456 3.477457 3.361606e-18 7.543612e-14
GCTAATAAATGGCAGAT  3.194647 5.819106 7.615461e-16 1.139298e-11
ATACTGACATTTCGTAT -4.321802 4.330351 5.854056e-15 6.568397e-11
CCAAGAATCTGGTCGTA  3.908138 3.539008 5.137517e-14 4.611538e-10
> detags<-rownames(topTags(et,n=5)$table)
> cpm(d)[detags,order(d$samples$group)]
                      DCLK1      DCLK2      DCLK3      DCLK4       WT1        WT2       WT3       WT4
TCTGTACGCAGTCAGGC  59.58111  31.540402 178.807638  50.677933  0.000000  0.2810365  0.000000  0.000000
CATAAGTCACAGAGTCG  24.94959  24.045653  23.570098  10.749864  0.000000  0.0000000  0.000000  0.000000
GCTAATAAATGGCAGAT 144.11180 100.242268  53.642291 109.034340 12.791475  8.9931672  3.390888 12.093130
ATACTGACATTTCGTAT   1.86191   1.561406   3.251048   1.535695 32.120815 64.0763160 13.563552 33.096987
CCAAGAATCTGGTCGTA  26.06673  20.610560  19.099907  19.964034  0.852765  1.4051824  0.000000  2.227682
> summary(results<-decideTestsDGE(et,p=0.085))
   [,1] 
-1   750
0  42601
1   1530
> detags<-rownames(d)[as.logical(results)]
> plotSmear(et,de.tags=detags)
> abline(h=c(-2,2),col="dodgerblue")
DGE data can be visualized as ‘MA’ plots (log ratio versus abundance)

DGE data can be visualized as ‘MA’ plots (log ratio versus abundance)

6 thoughts on “edgeR案例学习:deepSAGE分析

  1. Reply Joe 5月 2,2012 7:31 上午

    您好!又翻墙过来请教您了。
    我最近用的一款新品Affy的PrimeView,估计是比较新,bioconductor中很多包都没有。其中affy中qc()时总是显示,“错误于setQCEnvironment(cdfn) : Could not find array definition file ‘ primeviewcdf.qcdef ‘. Simpleaffy does not know the QC parameters for this array type.See the package vignette for details about how to specify QC parameters manually.”
    我又尝试用makecdfenv去做个包,但是似乎不成功,这个问题迟迟没有得到解决。请教您如何解决?
    另外,在后续结果进行配对分析时,如果有的样本 某一个时间点缺失了,如何去处理呢?之前我是用duplicatecorrelation来做的。

    谢谢!

    • Reply admin 5月 2,2012 11:41 上午

      @Joe, 先在你的R中获取simpleaffy的目录,
      system.file(“extdata”,package=”simpleaffy”)
      然后在其中新建一个名为primeviewcdf.qcdef的文件,内容为
      array primeviewcdf
      alpha1 0.05
      alpha2 0.065
      spk bioB AFFX-r2-Ec-bioB-3_at
      spk bioC AFFX-r2-Ec-bioC-3_at
      spk bioD AFFX-r2-Ec-bioD-3_at
      spk creX AFFX-r2-P1-cre-3_at
      ratio actin3/actin5 AFFX-HSAC07/X00351_3_at AFFX-HSAC07/X00351_5_at
      ratio actin3/actinM AFFX-HSAC07/X00351_3_at AFFX-HSAC07/X00351_M_at
      ratio gapdh3/gapdh5 AFFX-HUMGAPDH/M33197_3_at AFFX-HUMGAPDH/M33197_5_at
      ratio gapdh3/gapdhM AFFX-HUMGAPDH/M33197_3_at AFFX-HUMGAPDH/M33197_M_at

      关于第二个问题,你可能需要去除那个时间点了。之前我做过这种非平衡的设计,但是最终的结果是p.value增加很多,结果不是很好。如果你非要保留那个点,只能补实验了。

  2. Reply Joe 5月 3,2012 9:16 下午

    您好!
    按照您的方法,情况大为好转,但是现在还是出现问题,
    显示:错误于.qc.file.setQCArray(l) : Array name should be a single string.
    此外: 警告信息:
    In readLines(fl) :
    读’D:/R-2.15.0/library/simpleaffy/extdata/primeviewcdf.qcdef’时最后一行未遂

    • Reply admin 5月 4,2012 7:52 上午

      @Joe, 关于这个qcdef文件,你可以直接将hgu133plus之类的qcdef改名即可。因为它们使用了共同的qc define.但是我也遇到了一些问题,导致qc无法正常运行。可能这个问题还需要问原作者如何解决了。实在是对此暂时无能为力了。

  3. Reply ZiSu 4月 29,2014 6:26 上午

    博主您好,请问一下primeviewcdf.qcdef文件中的
    alpha1 0.05
    alpha2 0.065
    代表什么,还有就是
    spk bioB AFFX-r2-Ec-bioB-3_at
    spk bioC AFFX-r2-Ec-bioC-3_at
    spk bioD AFFX-r2-Ec-bioD-3_at
    spk creX AFFX-r2-P1-cre-3_at
    难道算bioB,bioC,bioD,cre的浓度只需要3端吗
    谢谢

Leave a Reply

  

  

  

%d 博主赞过: