接到学习任务，要求搞清楚两个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)) |
从图中我们可以看出其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(Mipgj,φg)
其中Mi是库的大小，就是样品i中reads的总数；φg就是基因< i>g的分散程度；pgj是i样品所属的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) |
接下来的步骤就是获取差异表达基因了。
> 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") |
您好！又翻墙过来请教您了。
我最近用的一款新品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来做的。
谢谢！
@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增加很多，结果不是很好。如果你非要保留那个点，只能补实验了。
您好！
按照您的方法，情况大为好转，但是现在还是出现问题，
显示：错误于.qc.file.setQCArray(l) : Array name should be a single string.
此外: 警告信息：
In readLines(fl) :
读’D:/R-2.15.0/library/simpleaffy/extdata/primeviewcdf.qcdef’时最后一行未遂
@Joe, 关于这个qcdef文件，你可以直接将hgu133plus之类的qcdef改名即可。因为它们使用了共同的qc define.但是我也遇到了一些问题，导致qc无法正常运行。可能这个问题还需要问原作者如何解决了。实在是对此暂时无能为力了。
博主您好，请问一下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端吗
谢谢
alpha1 和alpha2用于定义absent/present call。 只需要3‘端是因为RNA降解是从5’至3‘方向。