Using STAMP to identify SEED subsystems which are differentially abundant between CandidatusAccumulibacter phosphatis sequences obtained from a pair of enhanced biological phosphorus removal (EBPR) sludge metagenomes(data originally described in Parks and Beiko, 2010).

STAMP (Statistical Analysis of Metagenomic Profiles) is a software package for analyzing metagenomic profiles (e.g., a taxonomic profile indicating the [...]]]> STAMP

Using STAMP to identify SEED subsystems which are differentially abundant between CandidatusAccumulibacter phosphatis sequences obtained from a pair of enhanced biological phosphorus removal (EBPR) sludge metagenomes(data originally described in Parks and Beiko, 2010).

STAMP (Statistical Analysis of Metagenomic Profiles) is a software package for analyzing metagenomic profiles (e.g., a taxonomic profile indicating the number of marker genes assigned to different taxonomic units or a functional profile indicating the number of sequences assigned to different biological subsystems or pathways) that promotes ‘best practices’ in choosing appropriate statistical techniques and reporting results. It encourages the use of effect sizes and confidence intervals in assessing biological importance. A user friendly, graphical interface permits easy exploration of statistical results and generation of publication quality plots for inferring the biological relevance of features in a metagenomic profile. STAMP is open source, extensible via a plugin framework, and available for all major platforms.

Announcements

• Mar. 6, 2014: STAMP v2.0.2 released. Adds support for BIOM profiles including those generated by PICRUSt andQIIME.
• Feb. 9, 2014: STAMP v2.0.1 released. Adds heatmap plot. Fixes bug in multi-group plots causing DP values to be scaled by an 100.
• Sept. 2, 2013: OS X user may be experiencing problems running STAMP due to a bug in the latest version of matplotlib. Click here for details.
• May 24, 2012: STAMP v2.0.0 (release candidate 6) released. Improved handling of degenerate cases and unclassified taxa in STAMP profiles.
• Previous announcements

Mailing list

• Join our mailing list to keep informed about STAMP developments.

Documentation

• Quick installation instructions (Microsoft Windows, Apple’s Mac OS X, Linux, Command-line interface)
• User’s Guide
• FAQs
• Version history

Downloads

Please uninstall previous versions of STAMP before installing a new release.

• STAMP v2.0.2 setup package for Microsoft Windows (~35MB)
• OS X and Linux users can follow the instructions to install from source
• STAMP GitHub Repository
• Previous versions

Examples

• Image gallery
• Example datasets

Citing STAMP

If you use STAMP in your research, please cite:

Parks, D.H. and Beiko, R.G. (2010). Identifying biologically relevant differences between metagenomic communities. Bioinformatics, 26, 715-721. (Abstract)

Contact Information

STAMP is in active development and we are interested in discussing all potential applications of this software. We encourage you to send us suggestions for new features. Suggestions, comments, and bug reports can be sent to Rob Beiko (beiko [at] cs.dal.ca). If reporting a bug, please provide as much information as possible and a simplified version of the data set which causes the bug. This will allow us to quickly resolve the issue.

Funding

The development and deployment of STAMP has been supported by several organizations:

• Genome Atlantic
• The Dalhousie Centre for Comparative Genomics and Evolutionary Bioinformatics, and the Tula Foundation
• Killam Trust
• The Natural Sciences and Engineering Research Council of Canada
• The Dalhousie Faculty of Computer Science

This page discusses how to load GEO SOFT format microarray data from the Gene Expression Omnibus database (GEO) (hosted by the NCBI) into R/BioConductor. SOFT stands for Simple Omnibus Format in Text. There are actually four types of GEO SOFT file available:

GEO Platform (GPL) These files describe a particular type of microarray. They [...]]]> http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/r/geo/

This page discusses how to load GEO SOFT format microarray data from the Gene Expression Omnibus database (GEO) (hosted by the NCBI) into R/BioConductor. SOFT stands for Simple Omnibus Format in Text. There are actually four types of GEO SOFT file available:

GEO Platform (GPL)
These files describe a particular type of microarray. They are annotation files.

GEO Sample (GSM)
Files that contain all the data from the use of a single chip. For each gene there will be multiple scores including the main one, held in the VALUE column.

GEO Series (GSE)
Lists of GSM files that together form a single experiment.

GEO Dataset (GDS)
These are curated files that hold a summarised combination of a GSE file and its GSM files. They contain normalised expression levels for each gene from each sample (i.e. just the VALUE field from the GSM file).

As long as you just need the expression level then a GDS file will suffice. If you need to dig deeper into how the expression levels were calculated, you’ll need to get all the GSM files instead (which are listed in the GDS or GSE file).

To me, it was natural to ask: How can I turn a GEO DataSet (GDS file) into an R/BioConductor expression set object (exprSet)? (answer) And while we’re at it, how to load the GEO Platform annotation (GPL file) too? (answer)

In the MOAC Module 5 assignment, the approach taken was to sanitize the data by hand, allowing it to be loaded into R with a simple call to the read.table command. Its a good idea to look at the raw files to understand what you are dealing with, but surely there is a more elegant way…

It turns out there are several existing GEO parsers, but one stands out above all others: Sean Davis’GEOquery (released roughly December 2005).

Installing GEOquery

Assuming you are running a recent version of BioConductor (1.8 or later) you should be able to install it from within R as follows:

> source("http://www.bioconductor.org/biocLite.R") > biocLite("GEOquery") Running bioCLite version 0.1.6 with R version 2.3.1 ...

For those of you on an older version of BioConductor, you will have to download and install it by hand fromhere.

If you are using Windows, download GEOquery_1.6.0.zip (or similar) and save it. Then from within the R program, use the menu option “Packages”, “Install package(s) from local zip files…” and select the ZIP file.

On Linux, download GEOquery_1.6.0.tar.gz (or similar) and use sudo R CMD INSTALL GEOquery_1.6.0.tar.gz at the command prompt.

Loading a GDS file with GEOquery

Here is a quick introduction to how to load a GDS file, and turn it into an expression set object:

library(Biobase) library(GEOquery) #Download GDS file, put it in the current directory, and load it: gds858 <- getGEO('GDS858', destdir=".") #Or, open an existing GDS file (even if its compressed): gds858 <- getGEO(filename='GDS858.soft.gz')

I’m using GDS858 as input. The SOFT file is available in compressed form here GDS858.soft.gz, but GEOquery takes care of finding this file for you and unzipping it automatically.

Loading this file from the hard disk takes about two minutes on my laptop.

There are two main things the GDS object gives us, meta data (from the file header) and a table of expression data. These are extracted using the Meta and Table functions. First lets have a look at the metadata:

> Meta(gds858)$channel_count [1] "1" > Meta(gds858)$description [1] "Comparison of lung epithelial Calu-3 cells infected ..." > Meta(gds858)$feature_count [1] "22283" > Meta(gds858)$platform [1] "GPL96" > Meta(gds858)$sample_count [1] "19" > Meta(gds858)$sample_organism [1] "Homo sapiens" > Meta(gds858)$sample_type [1] "cDNA" > Meta(gds858)$title [1] "Mucoid and motile Pseudomonas aeruginosa infected lung epithelial cell comparison" > Meta(gds858)$type [1] "gene expression array-based"

Useful stuff, and now the expression data table:

> colnames(Table(gds858)) [1] "ID_REF" "IDENTIFIER" "GSM14498" "GSM14499" "GSM14500" [6] "GSM14501" "GSM14513" "GSM14514" "GSM14515" "GSM14516" [11] "GSM14506" "GSM14507" "GSM14508" "GSM14502" "GSM14503" [16] "GSM14504" "GSM14505" "GSM14509" "GSM14510" "GSM14511" [21] "GSM14512" > Table(gds858)[1:10,1:6] ID_REF IDENTIFIER GSM14498 GSM14499 GSM14500 GSM14501 1 1007_s_at U48705 3736.9 3811.0 3699.6 3897.6 2 1053_at M87338 343.0 500.3 288.3 341.3 3 117_at X51757 120.9 34.3 145.8 110.5 4 121_at X69699 1523.8 1281.1 1281.9 1493.4 5 1255_g_at L36861 51.6 15.9 45.9 8.1 6 1294_at L13852 253.2 164.8 200.0 205.2 7 1316_at X55005 199.6 250.7 290.3 218.6 8 1320_at X79510 81.7 13.4 13.9 88.7 9 1405_i_at M21121 18.9 5.6 11.0 9.5 10 1431_at J02843 99.7 74.5 72.6 114.8

Now, lets turn this GDS object into an expression set object (using base 2 logarithms) and have a look at it:

> eset <- GDS2eSet(gds858, do.log2=TRUE) > eset Expression Set (exprSet) with 22283 genes 19 samples phenoData object with 4 variables and 19 cases varLabels : sample : infection : genotype/variation : description > geneNames(eset)[1:10] [1] "1007_s_at" "1053_at" "117_at" "121_at" "1255_g_at" [6] "1294_at" "1316_at" "1320_at" "1405_i_at" "1431_at" > sampleNames(eset) [1] "GSM14498" "GSM14499" "GSM14500" "GSM14501" "GSM14513" [6] "GSM14514" "GSM14515" "GSM14516" "GSM14506" "GSM14507" [11] "GSM14508" "GSM14502" "GSM14503" "GSM14504" "GSM14505" [16] "GSM14509" "GSM14510" "GSM14511" "GSM14512"

GEOquery does an excellent job of extracting the phenotype data, as you can see:

> pData(eset)$infection [1] FRD1 FRD1 FRD1 FRD1 FRD440 [6] FRD440 FRD440 FRD440 FRD875 FRD875 [11] FRD875 FRD875 FRD1234 FRD1234 FRD1234 [16] uninfected uninfected uninfected uninfected Levels: FRD1 FRD1234 FRD440 FRD875 uninfected > pData(eset)$"genotype/variation" [1] control control [3] control control [5] mucoid mucoid [7] mucoid mucoid [9] motile motile [11] motile motile [13] non-mucoid, non-motile non-mucoid, non-motile [15] non-mucoid, non-motile non-mucoid, non-motile [17] non-mucoid, non-motile non-mucoid, non-motile [19] non-mucoid, non-motile Levels: control motile mucoid non-mucoid, non-motile

As with any expression set object, its easy to pull out a subset of the data:

> eset["1320_at","GSM14504"] Expression Set (exprSet) with 1 genes 1 samples phenoData object with 4 variables and 1 cases varLabels : sample : infection : genotype/variation : description > exprs(eset["1320_at","GSM14504"]) GSM14504 1320_at 6.70044

You should be able to produce a heatmap of differentially expressed genes easily enough using this page, especially as the phenotype/sub-sample information has been sorted out for you.

Loading a GPL (Annotation) file with GEOquery

In addition to loading a GDS file to get the expression levels, you can also load the associated platform annotation file. You can find this out from the GDS858 meta information:

> Meta(gds858)$platform [1] "GPL96"

So, for GDS858, the platform is GPL96, Affymetrix GeneChip Human Genome U133 Array Set HG-U133A.

Now let’s load up the GPL file and have a look at it (its a big file, about 12 MB, so this takes a while!):

library(Biobase) library(GEOquery) #Download GPL file, put it in the current directory, and load it: gpl96 <- getGEO('GPL96', destdir=".") #Or, open an existing GPL file: gpl96 <- getGEO(filename='GPL96.soft')

As with the GDS object, we can use the Meta and Table functions to extract information:

> Meta(gpl96)$title [1] "Affymetrix GeneChip Human Genome U133 Array Set HG-U133A" > colnames(Table(gpl96)) [1] "ID" "Species.Scientific.Name" [3] "Annotation.Date" "GB_LIST" [5] "SPOT_ID" "Sequence.Source" [7] "Representative.Public.ID" "Gene.Title" [9] "Gene.Symbol" "Entrez.Gene" [11] "RefSeq.Transcript.ID" "Gene.Ontology.Biological.Process" [13] "Gene.Ontology.Cellular.Component" "Gene.Ontology.Molecular.Function" 

Lets look at the first four columns, for the first ten genes:

> Table(gpl96)[1:10,1:4] ID Species.Scientific.Name Annotation.Date GB_LIST 1 1007_s_at Homo sapiens 16-Sep-05 U48705 2 1053_at Homo sapiens 16-Sep-05 M87338 3 117_at Homo sapiens 16-Sep-05 X51757 4 121_at Homo sapiens 16-Sep-05 X69699 5 1255_g_at Homo sapiens 16-Sep-05 L36861 6 1294_at Homo sapiens 16-Sep-05 L13852 7 1316_at Homo sapiens 16-Sep-05 X55005 8 1320_at Homo sapiens 16-Sep-05 X79510 9 1405_i_at Homo sapiens 16-Sep-05 M21121 10 1431_at Homo sapiens 16-Sep-05 J02843 

This shows a hand picked selection of the columns, again for the first ten genes:

> Table(gpl96)[1:10,c("ID","GB_LIST","Gene.Title","Gene.Symbol","Entrez.Gene")] ID GB_LIST Gene.Title Gene.Symbol Entrez.Gene 1 1007_s_at U48705 discoidin domain receptor family, member 1 DDR1 780 2 1053_at M87338 replication factor C (activator 1) 2, 40kDa RFC2 5982 3 117_at X51757 heat shock 70kDa protein 6 (HSP70B') HSPA6 3310 4 121_at X69699 paired box gene 8 PAX8 7849 5 1255_g_at L36861 guanylate cyclase activator 1A (retina) GUCA1A 2978 6 1294_at L13852 ubiquitin-activating enzyme E1-like UBE1L 7318 7 1316_at X55005 thyroid hormone receptor, alpha (erythroblastic...) THRA 7067 8 1320_at X79510 protein tyrosine phosphatase, non-receptor type 21 PTPN21 11099 9 1405_i_at M21121 chemokine (C-C motif) ligand 5 CCL5 6352 10 1431_at J02843 cytochrome P450, family 2, subfamily E, polypeptide 1 CYP2E1 1571 

The above all used the 12MB file GPL96.soft, but you can also get a much smaller 3MB file GPL96.annot(compressed as GPL96.annot.gz) which has slightly different information in it… see here.

Using the BioConductor hgu133a package

Instead of loading the GEO annotation file for GPL96/HG-U133A, we could use an existing annotation package from the BioConductor annotation sets, hgu133a. These libraries exist for most of the popular microarray gene chips.

First of all, we need to install the package:

> source("http://www.bioconductor.org/biocLite.R") > biocLite("hgu133a") Running bioCLite version 0.1 with R version 2.1.1 ...

Then we can load the newly installed library:

> library(hgu133a)

There is any easy way to check when this was lasted updated, and what it can translate the Affy probe names into:

> hgu133a() Quality control information for hgu133a Date built: Created: Tue May 17 13:02:12 2005 Number of probes: 22277 Probe number missmatch: None Probe missmatch: None Mappings found for probe based rda files: hgu133aACCNUM found 22277 of 22277 hgu133aCHRLOC found 20195 of 22277 hgu133aCHR found 21283 of 22277 hgu133aENZYME found 2507 of 22277 hgu133aGENENAME found 18726 of 22277 hgu133aGO found 18647 of 22277 hgu133aLOCUSID found 21747 of 22277 hgu133aMAP found 21183 of 22277 hgu133aOMIM found 15109 of 22277 hgu133aPATH found 5067 of 22277 hgu133aPMID found 21004 of 22277 hgu133aREFSEQ found 21002 of 22277 hgu133aSUMFUNC found 0 of 22277 hgu133aSYMBOL found 21303 of 22277 hgu133aUNIGENE found 21128 of 22277 Mappings found for non-probe based rda files: hgu133aCHRLENGTHS found 25 hgu133aENZYME2PROBE found 663 hgu133aGO2ALLPROBES found 5912 hgu133aGO2PROBE found 4326 hgu133aORGANISM found 1 hgu133aPATH2PROBE found 142 hgu133aPMID2PROBE found 96291

And now lets test some of those mappings on the fourth gene 121_at in the GPL file:

> Table(gpl96)[4,c("ID","GB_LIST","Gene.Title","Gene.Symbol","Entrez.Gene")] ID GB_LIST Gene.Title Gene.Symbol Entrez.Gene 4 121_at X69699 paired box gene 8 PAX8 7849

Now, what does the annotation file have to say?

> mget("121_at",hgu133aACCNUM) $"121_at" [1] "X69699" > mget("121_at",hgu133aGENENAME) $"121_at" [1] "paired box gene 8" > mget("121_at",hgu133aSYMBOL) $"121_at" [1] "PAX8" > mget("121_at",hgu133aUNIGENE) $"121_at" [1] "Hs.469728"

You will notice that there is some overlap between the information in the GEO annotation table, and thehgu133a package (which compiles its data from a range of sources). See help(hgu133a) .

You should also read this introduction, Bioconductor: Annotation Package Overview 

初读文章，我脑子里冒出的一句话是“上帝在跟我们掷骰子”，文中给出了大量的不可重复的试验，仿佛就像那些号称“具有统计学意义”（下文我再说这个所谓的“意义”）的试验结果只是若干次骰子中的一次运气好的结果而已。读完文章，我们可能不禁要问，到底是真理在缩水，还是它根本就不曾存在？下面我从四个方面来展开，分别说明人对随机性的认识、统计推断的基石、让无数英雄折腰的P值、以及可重复的统计研究。

随机变量在统计分析中占据中心地位，数学上关于随机变量的定义只是一个“干巴巴的函数”，从样本空间映射到实数集，保证从实数集上的Borel域逆回去的集合仍然在原来的sigma域中即可。随机变量的性质由其分布函数刻画。写这段话的目的不是为了吓唬你，也不是为了作八股文，而是来说明我为什么不喜欢数学的理由，对我而言，我觉得有些数学工具只是为了让自己不要太心虚，相信某时某刻某个角落有个理论在支撑你，但后果是弱化了人的感知，当然，也有很多数学工具有很强的直觉性（如果可能，我想在未来下一篇文章里面总结这些问题）。我一直认为很多人对随机性的感知是有偏差的，对概率的解释也容易掉进陷阱（参见Casella & Berger的Statistical Inference第一章，例如条件概率的三囚徒问题）。

《缩水》一文发表了很多不可重复的试验案例，我们应该吃惊吗？我的回答是，未必。举两个简单的例子：

第一个例子：很多数据分析人员都很在意所谓的“离群点”，论坛上也隔三差五有人问到这样的问题（如何判断和处理离群点），而且也有很多人的做法就是粗暴地删掉，我从来都反对这种做法。除了基于“数据是宝贵的”这样简单的想法之外，另一点原因是，离群点也许并非“异类”。离群点是否真的不容易出现？请打开R或其它统计软件，生成30个标准正态分布N(0, 1)随机数看看结果，比如R中输入rnorm(30)，这是我运行一次的结果：

30在现实中是一个比较小的样本量，你看到了什么？2.901？它接近3倍标准差的位置了。还有2.089？……如果你不知道这批数据真的是从标准正态分布中生成出来的，现在你会有什么反应？把2.9删掉？标准正态分布是一个在我们眼中很“正常”的分布，而一个不太大的样本量一次试验足以生成几个“离群点”，那么要是成千上万的试验中没能产生几项震惊世界的结果，你会怎样想？（上帝的骰子坏掉了）

另一个例子和统计学结合紧密一点，我们谈谈残差的QQ图。QQ图是用来检查数据正态性的一种统计图形，与腾讯无关，细节此处略去，大意是图中的点若呈直线状（大致分布在对角线上），那么可以说明数据的正态性比较好，因此QQ图经常被用在对回归模型残差的正态性诊断上。我的问题是，即使数据真的是正态分布，你是否真的会看见一些分布在直线上的点？若答案是否定的，那么我们就得重新审视我们对分布和随机的认识了。下图是一幅教科书式的QQ图（仍然基于30个正态分布随机数）：

“正常的”QQ图（来自R语言qqnorm(rnorm(30))）

随机性并没有这么美好，即使数据真的来自正态分布，你也有可能很容易观察到歪歪扭扭的QQ图，尤其是小样本的情况下。比如下图是50次重复抽样的正态数据QQ图，它和你想象的QQ图本来的样子差多远？