A reproducible approach to high-throughput biological data acquisition and integration

Introduction

The challenges of “big data” handling in modern biology have been highlighted repeatedly in recent publications (Lynch, 2008; Marx, 2013; Staff, 2011). The need for reproducible analyses of individual and integrated genome-scale data has resulted both in significant progress (Chan et al., 2008; Waldron, Coller & Huttenhower, 2012; Wirapati et al., 2008) and in substantial setbacks, such as the widely-publicized failed colon and ovarian cancer clinical trials (Dressman et al., 2007; Garman et al., 2008; Potti et al., 2006; Samuel Reich, 2011). Reproducibility in experimental biology is enforced by a combination of best practices and technical aids built up over decades of research: standardized protocols, reagents, recording practices, and equipment that ease the organizational burden necessary for high-quality research. Although parallel best practices have begun to be widely available for the equivalent processes in bioinformatic and data science investigations (Charaniya, Hu & Karypis, 2008; Noble, 2009; Wilson et al., 2014), they are not yet widely supported by standardized methods or computational environments enforcing reproducibility, which we develop here.

Reproducibly executing any single high-throughput biological study can be challenging, while reproducibly integrating multiple large-scale data sources presents an even greater hurdle. Meta-analysis of multiple genome-scale studies of the same biological system has proven to be a fruitful route for successful translational research, with over 500 such meta-analyses of human tumor gene expression alone now published (Chan et al., 2008; Wirapati et al., 2008). However, retrieving, standardizing, and curating such data is nontrivial, time-consuming, and error-prone. Integration of heterogeneous data can be an even more onerous process. For example, GEO (Barrett et al., 2011) and the BioGRID (Stark et al., 2011) each have over 2,000 citations as individual data sources, but fewer than 100 publications cite both (Warde-Farley et al., 2010). This lack of automatable systems for large-scale biological data handling and retrieval means that reliably incorporating such data into an analytical project is precluded without extensive computational expertise, and even in the best case, it is prone to error and challenging to reproduce (Gentleman, 2005; Ioannidis et al., 2009). This is particularly true when incorporating data from more than one source or type, or when linking data from public repositories to additional datasets from one or more individual labs. An example application, such as the identification of disease biomarkers, may involve establishing candidate markers by synthesis of published gene expression data, linking promising candidate markers to proteomic experiments, establishing putative biological mechanisms through databases of regulatory and protein interactions, and finally validating only the best candidates manually in the lab (Lu et al., 2005; Rifai, Gillette & Carr, 2006; Van ’t Veer et al., 2002). The advantages of such holistic approaches are obvious, but greatly expand the complexity of analysis.

Data retrieval systems have recently been developed that facilitate this process for consistent, large-scale access to individual ’omics repositories. The Gene Expression Omnibus (GEO) (Barrett et al., 2011; Barrett et al., 2013) and ArrayExpress (Rustici et al., 2013), for example, represent two of the largest functional genomics resources for transcriptional data, and both feature elegant web-based query systems (e.g., the Gene Expression Atlas (Kapushesky et al., 2012)) and programmatic access (e.g., GEOQuery (Davis & Meltzer, 2007)). Similarly, molecular interaction databases such as IntAct (Kerrien et al., 2012; Orchard et al., 2014) and the Biological General Repository for Interaction Datasets (BioGRID) (Chatr-Aryamontri et al., 2013; Stark et al., 2011) typically provide access to selected subsets of protein–protein interaction networks by way of taxonomic, technological, or temporal metadata. Each of these interfaces is database-specific, although they provide their data downloadable in the PSI-MI standard format; in general, expression data retrieved programmatically in R must be manually converted and combined with interaction data. Although some systems such as DAS (Dowell et al., 2001) provide a uniform interface for multiple ’omic data types, accessing data from most public repositories presents a major challenge in the compatibility of the resulting data formats, mapping between different gene identifier conventions, and the curation of inconsistently-annotated experiments. Manual parsing and conversion of experimental attributes and values is tedious and error-prone (Baggerly & Coombes, 2009), and the large number of normalization and preprocessing methods in use make comparison of different existing studies difficult or even impossible (Sims, 2009). This suggests the need for a unified approach capable of automatically retrieving and reproducibly processing data from multiple different ’omics repositories.

Previous related efforts for reproducible research include (1) general workflow systems for broad bioinformatic analysis that can be tuned for data retrieval, (2) tools enabling programmatic access to public databases, and (3) curated alternatives to the largest public experimental databases. General workflow systems, such as Kepler (Barseghian et al., 2010) and Taverna (Hull et al., 2006), facilitate the construction of reproducible programmatic workflows. In many cases they include programmatic access to biological databases, but subsequent data and metadata format harmonization must be implemented separately by the user. Tools that allow programmatic access to public databases, such the specialized interfaces provided by GEOquery (Davis & Meltzer, 2007) and Ensembl BioMarts (Kinsella et al., 2011), are among the most common. GEOquery is specific to automated queries and downloads from one particular database, and while queries can be written across multiple BioMart databases, even the latter tool includes no facilities for normalization, integration, or meta-analysis of the resulting independent query responses. Finally, curated alternative data management and analysis platforms, such as EZArray (Zhu & Xu, 2008), BxArrays™, and ONCOMINE (Rhodes et al., 2004), provide web interfaces in which a user may perform only pre-selected analyses of the provided gene expression datasets. Another open data management and access platform, inSilicoDB (Taminau et al., 2011), allows programmatic download of gene expression data curated by maintainers or contributors from the community, although users may only retrieve one experiment at a time. None of these tools facilitate bulk download of ’omics data from multiple sources for efficient processing on local computing resources, and most enable only per-repository data acquisition. A logical next step for large scale, multi’omic data acquisition for meta-analysis and heterogeneous data integration is again a batch-driven data acquisition system spanning multiple sources.

Motivated by this need to access and integrate diverse sources of experimental results, we developed an approach for reproducible processing of high-throughput and heterogeneous data. The system, implemented as the ARepA (Automated Repository Acquisition) software environment, enabled us to rapidly and systematically retrieve, integrate, and analyze large ‘omics datasets from multiple sources in a consistent and robust fashion. Here, we first undertook a study of the NFκB signaling pathway in the context of human prostate cancer by integrating ten tissue-specific human expression datasets across six platforms. This not only recovered but extended the currently known molecular network. Next, due to burgeoning recent interest in microbiome studies and the considerable knowledge gaps regarding host processes are dependent upon the gut microbiota, we investigated transcripts that were differentially expressed in the germ-free mouse colon, identifying and integrating six gene expression datasets. PPAR-alpha, RXR-alpha, and several other immune- and carbohydrate metabolism-related genes were detected as differentially expressed in this gene expression meta-analysis. Finally, we integrated protein–protein interactions, gene expression, and regulatory networks for three very diverse model bacteria (Bacillus subtilis, which is Gram-positive; Escherichia coli, which is Gram-negative, and Pseudomonas aeruginosa, which is Gram-negative and pathogenic), correctly identifying differentially-conserved peptide secretion system pathways.

Each of these three studies required rapid, reproducible, and diverse data collection and integration, and thus demonstrates the utility of ARepA for any lab wishing to combine its own data with public data repositories. We provide an implementation of ARepA that can access ’omics data from multiple data sources in a consistent and robust fashion. Furthermore, it can be configured to process data only from certain sources, of specified types, or from particular organisms of interest, depending on the user’s interest and requirements. Currently, it supports seven repositories: the Gene Expression Omnibus (GEO) (Barrett et al., 2011; Barrett et al., 2013), IntAct (Kerrien et al., 2012; Orchard et al., 2014), BioGRID (Stark et al., 2011), RegulonDB (Gama-Castro et al., 2011), STRING (Szklarczyk et al., 2011), Bacteriome (Su et al., 2008), and MPIDB (Goll et al., 2008) (Table 1), from which data are retrieved, processed and standardized for any organism of interest. At present, this produces more than 2 M interactions from ∼600 K physical interactions (Bacteriome, BioGRID, Intact, and MPIDB), ∼4 K regulatory interactions (RegulonDB) and 1.5 M functional associations (STRING), as well as more than 2.7 K expression data sets covering more than 800 K samples (GEO). When retrieving data from six model organism, ARepA currently produces >630 K interactions and >15 K gene expression data sets (>530 K samples), accompanied by corresponding metadata and derived molecular networks (Table 2 and Fig. 1).

Material and Methods

Case studies

Human prostate cancer data retrieval and molecular network integration

We screened ARepA’s GEO metadata database for human prostate cancer and prostate tissue conditions matching “prostate cancer” in human gene expression profile studies with at least 6 samples each, identifying ten relevant datasets from six different platforms (Affymetrix, Agilent, and CNIO, Table S1). Identifying, processing, and standardizing these datasets in ARepA, accompanied by computing co-expression networks, required <1 h running time in total (on a 2.9 GHz Intel Core i7 16 GB machine). These data were then meta-analyzed to a single prostate cancer specific functional network using unsupervised data integration averaging across normalized co-expression values (z-scores) (Huttenhower et al., 2009) (Text S1). Next, we used this integrated network to predict genes highly functionally related to NFκB gene family (NFκB1, NFκB2, RELA, RELB, REL) in prostate cancer by extracting a high-confidence subgraph based on the probabilistic graph search algorithm as described in Huttenhower et al. (2009); Myers et al. (2005). This starts from a user-defined set of query genes and identifies k additional neighbours in the network that are connected with high confidence to the original query genes (Huttenhower et al., 2008). Here, we defined a query gene set of twelve genes from the NFκB signaling pathway in BioCarta (Table S4) and a neighbourhood size of k = 10, resulting in a high-confidence NFκB signaling network containing 22 genes in total as illustrated in Fig. 2.

Murine expression data retrieval and meta-analysis

We conducted a differential expression meta-analysis of genes and pathways up- or down-regulated in the germ-free (without microbes) murine gut (Fig. S1 and Text S1). We used six case-control datasets containing intestinal tissue and germ-free versus wild type comparisons, identified through ARepA’s metadata screen by matching “germ-free,” “wild type,” and “intestinal tract” in mouse gene expression profiles with at least four samples each (Table S2). This again required <1 h total running time (on a 2.9 GHz Intel Core i7 16 GB machine). To integrate the resulting ARepA R data files, we first computed log fold changes and confidence intervals (2.5% and 97.5%) for all genes between germ-free gut and wild type gut within each dataset using the R/limma package (Smith, 2005), resulting in ∼3,600 differentially expressed genes. Next, we performed a meta-analysis using the R/metafor package (Viechtbauer, 2010) on the six datasets by applying a random-effects model on the fold changes with default options, fitting the model with the restricted maximum-likelihood estimator (REML). The false discovery rate was controlled by the Benjamini–Hochberg method (Benjamini & Hochberg, 1995). We finally tested all genes and their resulting meta-p-values for gene set enrichment (Subramanian et al., 2005) in KEGG and BioCarta pathways. This resulted in two pathways that were significantly enriched for upregulated genes under germ-free conditions (using 1,000 permutations, Fig. 3), while 15 pathways were enriched for downregulated genes (Table S3).

Data integration for bacterial model organisms

We retrieved all gene expression data (GEO) and gene interaction data (IntAct, MPIDB, RegulonDB, Bacteriome and STRING, Table 2A) from ARepA for Bacillus subtilis, Escherichia coli and Pseudomonas aeruginosa. These were processed as above and standardized onto a shared gene identifier (UniRef90). Within each bacterium, these data were integrated into a species-specific functional network using an unsupervised data integration averaging across normalized co-expression values (z-scores) (Huttenhower et al., 2009) and extracting a high-confidence subgraph based on the probabilistic graph search algorithm (Huttenhower et al., 2009; Myers et al., 2005) (Text S1). This query (Huttenhower et al., 2008) started from a user-defined set of genes of the sec and tat genes and genes from the Type I, Type II, Type III, Type V, and Type VI secretion systems (Table S4) and used a neighbourhood size of k = 5 for each subgraph extraction. As Gram-positive bacteria use only the general secretory (sec) pathway and the tat pathway, we used two tat genes and nine sec genes as the query gene set in B. subtilis. In contrast, E. coli, a Gram-negative bacterium,  employs an additional secretion system (Type II secretion system) and has one type I gene (TolC); thus, we used four tat genes, twelve sec genes, one type I gene, and twelve type II genes as a query gene set in E. coli. P. aeruginosa, another Gram-negative bacterium, additionally utilizes the type I, type III , type V, and type VI secretion systems, so we defined our query gene set of three tat genes, nine sec genes, two type I genes, eleven type II genes, twelve type III genes, one type V gene, and eight type VI genes for retrieving a species-specific high-confidence network in P. aeruginosa. We recovered networks of 16 genes for B. subtilis (Fig. 4A), 34 genes for E. coli (Fig. 4B), and 49 genes for P. aeruginosa (Fig. 4C), which showed clustering of secretion types conserved across species.

Results

Motivated by the need to reproducibly access and unify diverse data resources, we developed a novel approach for the retrieval, normalization, and meta-analysis of high-throughput biological data. In this study, we applied our system, ARepA (Automated Repository Acquisition), to three integrative “big data” analyses: an integration of gene networks in human prostate cancer to identify interactors with the NfκB pathway, a differential expression meta-analysis of genes and pathways up- or down-regulated in the germ-free murine colon, and an integrative analysis of heterogeneous data from three model microbes to identify differentially conserved peptide secretion system pathways.

ARepA is implemented as a publicly available software package, currently providing reproducible and automated data and metadata processing for the Gene Expression Omnibus (GEO) (Barrett et al., 2011), IntAct (Kerrien et al., 2012), BioGRID (Stark et al., 2011), RegulonDB (Gama-Castro et al., 2011), STRING (Szklarczyk et al., 2011), Bacteriome (Su et al., 2008), and MPIDB (Goll et al., 2008). Investigators can retrieve normalized data and metadata in standardized formats from these repositories for any organism of interest. For six model organisms studied here in detail (human, mouse, yeast, E. coli, B. subtilis, and P. aeruginosa), it currently yields data comprising more than 630 K interactions: 250 K physical interactions (Bacteriome, BioGRID, IntAct, and MPIDB), and 380 K functional associations (STRING), in addition to approximately 15 K expression data sets covering more than 530 K samples (GEO) with standardized formats and gene identifiers. Metadata and network data are also automatically produced for each dataset and can be summarized in a portable R data file (Fig. 1).

Standardized and reproducible high-throughput data retrieval and integration

As evidenced by these studies, ARepA provides an analysis platform that is modular and hierarchical. The system includes access to multiple repositories; within each, hierarchical data retrieval first queries all experiments from a particular source, then all datasets within each experiment, and finally all metadata for each dataset, in a manner similar to the design of the ISA standard (Rocca-Serra et al., 2010). Data from each repository can be retrieved independently if desired, as can data from individual datasets, and a user can customize multi-layered data processing pipelines in a consistent, templated manner if desired. Each repository retrieves multiple experiments (typically organized by taxon ID and publication ID), and each experiment can retrieve multiple platforms or datasets. Information is shared vertically, and data passed from parent to child analyses can include taxonomic information, desired dataset types, and repository-specific configuration information, such as type of normalization functions and machine learning algorithms. An ARepA user can generally request that all modules beneath a particular repository be executed, or that only a single submodule (e.g., one dataset) is processed.

As an example of dedicated module implementation, the GEO module first downloads raw gene expression data (SOFT or SeriesMatrix files) as well as platform annotations from the Gene Expression Omnibus (Fig. 5). The number of platforms used by each dataset is determined, and data for each platform is built recursively. The expression data is parsed, normalized and imputed (Huttenhower et al., 2008), and gene IDs are standardized using GeneMapper. Finally, metadata is created for each dataset, combining automatically-derived information from the SOFT files with manually curated tables if available. The final outputs are a standardized and consistent gene expression matrix encoded as tab-delimited text, corresponding metadata, a normalized co-expression network, and an R data file containing an expression set. Other dedicated modules, such as those for interaction data, operate similarly by first downloading and parsing the requested data and directly generating tab-delimited text networks for each dataset (typically separated by taxon ID and publication ID). Standardization of interaction data and metadata are performed as described for the GEO module.

Finally, data retrieved in this environment can be integrated or meta-analyzed either heterogeneously (e.g., our bacterial study) or homogeneously (e.g., the human prostate cancer network integration and murine gut differential expression studies) by procedures such as fixed or random effects meta-analysis or network integration. The former can be directly applied to identify which genes are consistently up- or down-regulated in a given organism under selected experimental conditions. This environment facilitates the identification of relevant experiments, collection of associated datasets, and application of standardized R/Bioconductor meta-analysis environments such as metaphor (Viechtbauer, 2010) to gene expression data from various laboratories and platforms. Network integration studies, in contrast, can be performed for data from any source of interest (e.g., co-expression, protein–protein interactions, regulatory interactions, posttranslational modifications, etc.) by, effectively, meta-analyzing gene interactions (edges) rather genes (nodes). The simplest approach averages across normalized interaction scores (as illustrated in the human prostate cancer network integration and bacterial study) and can be extended using more sophisticated probabilistic (Franceschini et al., 2013; Lee et al., 2004; Myers & Troyanskaya, 2007; Zuberi et al., 2013) or discriminative (Chikina et al., 2009; Guan et al., 2012) models. Both gene and interaction level integrations can be performed from within ARepA directly from the standardized metadata and formatted data.

Discussion

ARepA was developed to provide an automatic, consistent, reproducible, and high-throughput method of accessing, retrieving, standardizing, processing, and integrating diverse, heterogeneous data resources, suitable for large integrative analyses and meta-analyses in translational and bioinformatic research studies. To illustrate ARepA’s utility in such studies, we have described here three biological studies demonstrating the utility of ARepA for integrative analyses. When focusing on human data, ARepA’s metadata database allowed us to identify and standardize ten human prostate cancer gene expression datasets from GEO, which were subsequently meta-analyzed, while a subsequent co-expression network analysis recovered new candidate genes with potential roles in the NFκB signaling pathway in prostate cancer. A similar example in mouse integrated six gene expression datasets selected by querying ARepA for metadata indicating germ-free and intestinal tissue conditions. Finally, multiple data types from three model microbes were integrated to assess functional differences in peptide secretion systems.

In this study, we correctly recovered and extended the molecular network of the NFκB signaling pathway in human prostate cancer by using this environment to identify ten tissue-specific human expression datasets and subsequently meta-analyze them across six different array platforms. The total computing time needed to identify, normalize, and integrate these datasets was approximately 1 h. In contrast, performing a comparable process using Taverna (Hull et al., 2006) required 8 total main steps: (1) downloading gene expression data from GEO, (2) loading and accessing gene expression data (maxd software: maxdLoad2 and maxdBrowse (Hancock et al., 2005)), (3) attaching (curated) metadata (maxdLoad2), (4) normalizing and pre-processing the raw data (Taverna/Biomoby/runPreprocessing and AffyArrayNormalization (Kawas, Senger & Wilkinson, 2006) or RShell (Wassink et al., 2009)), (5) identifying and downloading dataset specific mapping files, (6) mapping from probeset to gene identifier (Taverna/Biomoby/GeneExpressionAnnotation, BeanShell scripting language), (7) generating co-expression networks (Huttenhower et al., 2008), and (8) building an integrative expression set object. Running a similar process using GEOquery (Davis & Meltzer, 2007) required seven total main steps: (1) downloading expression data from GEO (R/GEOquery/getGEO), (2) normalizing and pre-processing the raw data (R/affy/RMA (Gautier et al., 2004; Irizarry et al., 2003)), (3) identifying and downloading dataset specific mapping files, (4) mapping from probeset to gene identifier for non-standard microarray platforms, (5) attaching (curated) metadata, (6) generating co-expression networks (Huttenhower et al., 2008), and (7) building an integrated expression set object. Using ARepA requires none of these steps, as it is specifically constructed for acquisition, standardization, query, and integration tasks, resulting in co-expression/interaction networks and expression sets without the need for additional intervention.

Performing this task in ARepA allowed us to recover several genes already associated with NFκB and prostate cancer (Davis, Kucuk & Sarkar, 1999; Perkins, 2012), in addition to original findings including MEN1 as a novel putative upstream regulator (Bouwmeester et al., 2004; Hacker & Karin, 2006; Heppner et al., 2001; Wen et al., 2000) and ACBD6 as a potential downstream suppressor of NFκB-activation (Oberti et al., 2011; Pringle et al., 2012; Soupene et al., 2012; Wu et al., 2006). MEN1 was linked to cancer, prostate cancer and NFκB by two main data sources. First, MEN1 co-immunoprecipitated with the NFκB proteins p50, p52, and p65 in vitro and in vivo, revealing MEN1 as a repressor of NFκB/p65-mediated transcriptional activation and a tumor suppressor protein (Heppner et al., 2001). Second, MEN1 and NFκB/p65 were co-expressed in hepatocellular carcinoma (HCC) tissues, adjacent normal liver tissues and various cell lines (including THLE-3, HepG2, Hep3B, HuH-7 and Li-7 cells) (Gang et al., 2013). Although no mechanism was previously proposed, heterozygous Men1 mutant mice (Men1^±) carrying an inactivated Men1 allele (Bertolino et al., 2003) and mixed with C57BL6/129-Sv mouse strains developed prostate cancer as compared with wild-type (Men1^+/+) littermates in aged mice (18 to 26 months) (Seigne et al., 2010). The strong functional linkage of MEN1 with IKBkG in our data, as well as p53, supports the previously suggested role for this gene product upstream of the NFκB pathway as a tumor suppressor by inhibition of NFκB-mediated transactivation. This could be further validated by, for example, siRNA knockdown of MEN1 in LAPC4 cells to determine whether MEN1 silencing results in more rapid cellular growth in a hormone-depleted state.

To the best of our knowledge, ACBD6 has not been directly linked to cancer, although its expression promotes C. trachomatis-mediated carcinogenesis in vitro (Soupene et al., 2012). A homologous protein also carrying the acyl-CoA binding domain, ACBD3, was highly expressed in primary pigmented nodular adrenocortical disease (PPNAD) tissues (Liu, Li & Papadopoulos, 2003), lung cancer tissue in treated patients that did not respond to gefitinib (Kakiuchi et al., 2004), colonic mucosa cancer tissue upon pretreatment versus post-treatment with celecoxib (Glebov et al., 2006), and bone tissue in osteosarcoma patients (Jones et al., 2012), supporting a role as a potential oncogene (Fan et al., 2010). Additionally, another homologous protein, ACBD4, was highly up-regulated in cell lines from different tumor types (including MCF-7 breast adenocarcinoma, HeLa cervical cancer, HT1080 sarcoma, KB laryngeal carcinoma, SW480 colon carcinoma, Hep-2 oral carcinoma and D54 glioma) after inducing by the histone deacetylase inhibitor valproic acid (VPA) (Chavez-Blanco et al., 2006), while Acyl-CoA synthetase (ACS), which converts long-chain fatty acids into acyl-CoA, is overexpressed in different cancer cell lines (including glioma SF268 and U251 cells) (Mashima et al., 2009), lung cancer (Pei et al., 2013), and glioblastoma (Pei et al., 2009). Based on such links between acyl-CoA binding domains, their synthetase, and cancer, as well as our data that predicted a strong functional link between ACBD6 and NFκB via its inhibitors NFκB IA and NFκBIE (Pringle et al., 2012; Wu et al., 2006) (Fig. 2), we suggest a potential role of ACBD6 in prostate cancer as an upstream suppressor protein of NFκB activation. This could be experimentally validated by screening series of prostate cancer cells (LAPC4, LNCaP) to determine whether ACBD6 is present, and in those cell lines in which it is present, siRNA knockdown of ACBD6 could be performed do to functional assays, such as proliferation in hormone sensitive (LNCaP) and hormone resistant cell lines (LNCaP-abl).

Our study of host-microbiota interactions focused on identifying transcripts differentially expressed in the germ-free mouse colon. As expected based on prior knowledge of germ-free mice, there was strong differential expression of immune and carbohydrate metabolism-related genes. Meta-analysis showed that genes from the Ppar signaling pathway were highly differentially expressed under germ-free conditions. Specifically, Ppar-α, a major regulator of lipid metabolism was upregulated. Ppar signaling overall plays key roles in regulating energy metabolism, cellular proliferation and differentiation in many tissues, including the intestine, kidney, heart, and lung (Braissant et al., 1996; Nephew, Polek & Khan, 1996). Dysregulated Ppar signaling has been implicated in metabolic syndrome, for which there is an increasing appreciation for the role of gut microbiota (Aitken & Gewirtz, 2013; Tremaroli & Backhed, 2012). Prior work using Ppar-α deficient mice coupled with gnotobiotic re-derivation suggests that Ppar-α does not function in the microbiota-dependent regulation of fat storage (Backhed et al., 2004). Thus to determine its specific role in this predicted context, follow-up work could include knockout mice for these genes (inclusive of conditional mice-deficiencies within intestinal epithelial or immune cells within the gut), with or without defined microbial consortia, screened using different diets (e.g. high fat, methionine and choline deficient). The observation that the germ-free state resulted in significant perturbations to Ppar signaling pathways suggests that how microbial consortia and their metabolites affect host regulation of obesity-related inflammation, insulin responsiveness and cholesterol homeostasis should be further investigated.

In the third study, multiple data types from three model microbes (B. subtilis, E. coli, and P. aeruginosa) were integrated to correctly identify differentially conserved peptide secretion system pathways in the Gram-positive bacterium B. subtilis and the two Gram-negative bacteria E. coli and P. aeruginosa. The sec and tat pathways are highly conserved across all three bacteria and were easily retrieved during this analysis, as were type II secretion in E. coli and P. aeruginosa and the type I, III, and V, and VI pathways in P. aeruginosa (Kostakioti et al., 2005). The E. coli network was built from expression data from strain K12 (taxonomy ID 83333). While some E. coli strains (e.g., strain O157:H7) have type III secretion systems (Sharma et al., 2012), strain K12 does not. Although it has a (type I) TolC gene, it also lacks a complete HlyABCD type I secretion operon (Bakkes et al., 2010). We populated our K12 expression-based network with only the secretion systems possessed by strain K12, and correctly recovered strong sec/tat co-expression, and strong sec/tat/type II co-expression.

In P. aeruginosa, there was a strong type II/III secretion link as well as strong tat/sec coordination, while type VI secretion was only weakly connected to sec/tat. Type V secretion, which is sec-dependent (Bleves et al., 2010), was strongly linked to sec/tat. The strong coregulation of type II and III secretion in Pseudomonas is potentially due to their joint role in pathogenesis; for example, in Pseudomonas-mediated pneumonia, type II secretion is used to secrete proteases such as LasA and LasB, phospholipase H, and lipolytic enzymes which degrade host tissue; type II secretion allows Pseudomonas infections to persist in mice (Jyot et al., 2011). The toxins secreted by the type III secretion system are complimentary, but focused on immune inhibition rather than tissue destruction. For example, ExoT slows wound healing and inhibits bacterial internalization by macrophages, and ExoS induces lymphocyte apoptosis (Schulert et al., 2003). Pseudomonas that lacks type III secretion is much more vulnerable to immediate clearance by neutrophils (Sun et al., 2012). Thus, the type II and type III secretion systems in Pseudomonas work together to help it persist in the host in infection. Type VI secretion in Pseudomonas is much less well-characterized, but proteins secreted by this mechanism do not contain sec or tat signals (Pallen, Chaudhuri & Henderson, 2003). While type VI secretion has been implicated in biofilm production (Zhang et al., 2011), the operon associated with biofilm is repressed by quorum sensing (Bleves et al., 2010), so an association with sec and tat that is driven by quorum sensing for the purpose of biofilm production is unlikely.

In comparison with existing bioinformatics workflow tools like Taverna (Hull et al., 2006) and Galaxy (Giardine et al., 2005; Goecks, Nekrutenko & Taylor, 2010), ARepA is specifically tuned to download, normalize, and meta-analyze very large and diverse collections of high-throughput biological data. While users can author new modules for other analyses of the resulting integrated data, the system is tuned to facilitate this task seamlessly out of the box. Galaxy, on the other hand, focuses primarily on real-time, interactive analyses, typically of a few datasets of interest through a graphical user interface. This makes it very flexible when developing new or customized analysis workflows of particular datasets, but less easily adaptable to bulk processing of especially diverse data. Taverna, alternatively, provides an extensive platform for combining bioinformatics tools in a graphical environment, but again it does not directly facilitate bulk data download, normalization, and integration. By uniformly normalizing and providing metadata for results deposited across many high-throughput data repositories, ARepA makes it easy to integrate all relevant data types for organisms, phenotypes, diseases, or molecular targets of interest.

One of the most important application areas for an environment such as ARepA is in translational meta-analyses, which relies on large data collections from multiple sources to ensure cross-study validation (Ganzfried et al., 2013; Griffith et al., 2006). ARepA supports such efforts in two ways. First, it effectively integrates data of heterogeneous format and origin without the need for manual curation or user intervention. This is true both for data normalization and for metadata collection, the latter of which is also integrated into a single, searchable database from across multiple repositories. ARepA’s automated processing can thus be used as a first step in particularly highly targeted manual curation if desired, as has been the case in our own previous studies (Ganzfried et al., 2013; Waldron, Coller & Huttenhower, 2012). Second, ARepA provides complete provenance and reproducibility for all data retrieved, as this is similarly a priority in pre-clinical research. The data processing and normalization steps carried out by ARepA can be modified or added to by a user if desired, ensuring that exactly the same steps are carried out at each execution for each dataset. This remains practical even when handling large data collections since ARepA’s dependency-tracking architecture ensures execution only when the input data source has changed or new analyses are added.

A second potentially important environment for ARepA is as a basis for “in-house” repositories for groups supporting internal data or focusing on specific types (e.g., stem cell data). This provides a convenient way to start from public repositories by first screening metadata for areas of interest. This screen can be added to ARepA’s data retrieval process for each source repository, ensuring that only data of interest is downloaded and updated each time; the latter is important since public repositories do updated regularly and on a varying basis (e.g., GEO is updated continuously, BioGRID once a month, and STRING once to twice per year). Privately maintained data can be appended by adding a new module to ARepA that fetches data locally, allowing it to be standardized to the same format as that of the other repositories for integration. ARepA is generally designed to be customizable to individual user environments, allowing modification of data retrieval or processing scripts while still taking advantage of its dependency-tracking system. Public data are typically fetched via ftp or http; however, a user can easily acquire data using other protocols such as ssh, Aspera, or simple local copies or links. This flexibility in design allows users to integrate publically available data with privately maintained data to conduct in-silico experiments within a single consistent environment.

Since any high-throughput genomic assays are susceptible to batch effects (Leek et al., 2010), analysts performing meta-analyses should be aware of their potential impact. Batch correction methods in general either require a priori knowledge of the batches (e.g., for empirical Bayes adjustment (Johnson, Li & Rabinovic, 2007)) or are specific to the outcome variable of a particular differential expression (Leek et al., 2012), as has been previously reviewed (Chen et al., 2011; Lazar et al., 2013). Thus, batch correction is a study-specific issue that cannot be resolved solely by a high-throughput tool such as ARepA, although it provides a framework within which individual datasets can be easily handled. Although study-specific sensitivity analysis of batch correction is warranted when possible, several studies have found that traditional meta-analysis approaches are robust to batch effects and that explicit batch correction may be unhelpful for independent validation of differential expression (Re et al., 2014) or prediction models (Riester et al., 2014; Waldron et al., 2014).

Like any computational environment, this is one of a number of routes available for future expansions to ARepA. An obvious next step is the inclusion of additional repository modules, with large and well-annotated sources such as ArrayExpress (Rustici et al., 2013), Refseq (Pruitt, Tatusova & Maglott, 2005), ENCODE (Consortium, 2012), or RNA-sequencing data representing excellent candidates. This will likewise allow the integration of additional data types such as genomic sequences, regulatory elements, and high quality expression data that would broaden the set of data to be considered for subsequent meta-analyses. The addition of a graphical user interface will make ARepA more accessible for a broader user community and increase its usability, as will a unified configuration interface (e.g., a single, XML-based configuration file) to more easily change parameters. The current version of ARepA thus represents only the first step in the ongoing development of an environment for reproducible “big data” research, which we hope will be both expandable by bioinformaticians and useful for targeted biological studies.

Supplemental Information