Bacterial transcriptome reorganization in thermal adaptive evolution
© Ying et al. 2015
Received: 18 June 2015
Accepted: 3 October 2015
Published: 16 October 2015
Evolution optimizes a living system at both the genome and transcriptome levels. Few studies have investigated transcriptome evolution, whereas many studies have explored genome evolution in experimentally evolved cells. However, a comprehensive understanding of evolutionary mechanisms requires knowledge of how evolution shapes gene expression. Here, we analyzed Escherichia coli strains acquired during long-term thermal adaptive evolution.
Evolved and ancestor Escherichia coli cells were exponentially grown under normal and high temperatures for subsequent transcriptome analysis. We found that both the ancestor and evolved cells had comparable magnitudes of transcriptional change in response to heat shock, although the evolutionary progression of their expression patterns during exponential growth was different at either normal or high temperatures. We also identified inverse transcriptional changes that were mediated by differences in growth temperatures and genotypes, as well as negative epistasis between genotype—and heat shock-induced transcriptional changes. Principal component analysis revealed that transcriptome evolution neither approached the responsive state at the high temperature nor returned to the steady state at the regular temperature. We propose that the molecular mechanisms of thermal adaptive evolution involve the optimization of steady-state transcriptomes at high temperatures without disturbing the heat shock response.
Our results suggest that transcriptome evolution works to maintain steady-state gene expression during constrained differentiation at various evolutionary stages, while also maintaining responsiveness to environmental stimuli and transcriptome homeostasis.
Evolutionary experimentation is a powerful tool for the exploration of how living organisms adapt to environmental change and is commonly applied in studies of evolution, typically focusing on changes in genome sequences. The results of these studies have provided experimental evidence to substantiate or revise numerous theories of evolution [1, 2]; the evolutionary mechanisms were generally explained by changes in DNA sequences and/or some particular genes [3–6].
Nevertheless, adaptive evolution also works to optimize cellular gene expression; thus, transcriptome evolution is a relevant area of study. In addition to studies analyzing the transcriptome response to stimuli such as heat shock and nutritional stress [7–11], recent studies report the effects of highly abstract phenomena on global transcriptome parameters. For example, studies showing gene transcriptional changes that were correlated with growth rates, performed in both yeast [12–15] and bacteria [16–18], as well as studies of the trade-off relationship for gene expression responsible for growth fitness and the stress response [18, 19]. Our previous studies identified transcriptome-wide growth-induced transcriptome reorganization  and demonstrated that a similar global coordination was also found for stochastic adaptation independent of regulatory mechanisms .
Global optimization of gene expression is critical for living cells to achieve novel adaptive states in new environments and must occur not only during transient adaptive responses but also during long-term adaptive evolution . The pioneering transcriptome studies using experimental evolution provided detailed transcriptional changes for genes that are specifically involved in evolution [21, 22]. Studies of parallel evolution in different environments usually found transcriptome diversity and a correlation to genetic background [23–25]. Although a number of studies have reported the general features of transcriptome evolution [6, 23–27], it is unclear whether transcriptome evolution requires alteration of the genomic response to acquire an adaptive state.
As a pilot investigation, we compared transcriptomes of the heat shock response and thermal adaptation. We previously performed experiments in thermal adaptive evolution with a laboratory Escherichia coli strain. We clarified the mechanism of genome evolution, i.e., positive selection vs. neutral evolution ; however, how transcriptomes reorganize during the evolution of thermal adaptation remained unknown. The heat shock response is a critical mechanism that living cells use to tolerate high temperatures/thermal stress [28, 29]. In E. coli, this mechanism is largely dependent on feedback regulation from the sigma factor 32 gene (rpoH) and heat shock proteins such as GroEL/ES and DnaK [29, 30]. It is intriguing to know whether Escherichia coli cells alter this responsive machinery to reach adaptive states, or maintain heat shock responsivity using other expression patterns during evolution.
To address this question, we performed a microarray analysis with the representative ancestor (Anc) Escherichia coli cells, and three evolved cell populations (41B, 43B and 45 L), which were acquired at varied evolutionary periods and were well adapted to growth temperatures of 41.2, 43.2 and 44.8 °C, respectively. To distinguish the thermal stress response states and the states of adaptation to high temperatures, we obtained the transcriptomes of both the exponential growth phase at the different growth temperatures (designated as the steady states) and the heat shock response (designated as the responsive states). Multilevel analyses were conducted to investigate the evolution of thermal adaptive transcriptomes. Overall, our results found that long-term evolutionary adaptation to high temperatures was different for the transient response to thermal stress. We summarized these results by proposing a potential mechanism that illustrates transcriptome evolution accompanied by genome evolution.
Cell culture and growth
The Escherichia coli cells were cultured in 5 mL of minimal medium M63, supplemented with 2 mM leucine and 25 μg/mL kanamycin, at corresponding temperatures. The culture conditions were exactly the same as previously reported . Cells were counted using flow cytometers, either the FC500 (Beckman) or the FACSCanto II (Becton-Dickson). Cell concentrations were calculated as the ratio of gated particles representing the number of Escherichia coli cells carrying the reporter gene gfp (green fluorescent protein) and beads of known concentrations, as previously described [18, 31]. The growth rates were calculated using the initial and final cell concentrations and the culturing time, as described elsewhere [18, 20, 31]. The initial and final cell densities were maintained at approximately 1 × 104 and 2 × 108 cells/ml, respectively. The culturing times were varied from 14 to 23 h. Cells within the exponential growth phase were collected for microarray analysis.
Heat shock experiments
The condition for the heat shock experiment was a 5-min incubation following a temperature upshift from 36.9 °C to 44.8 °C, as previously described . Exponentially growing cells at approximately 2 × 108 cells/mL were rapidly transferred to an adjacent water bath incubator (Personal-11EX; Taitec) set at 44.8 °C. Following a 5-min incubation at 44.8 °C, the cell culture was immediately poured into a cold phenol-ethanol solution to prepare the samples for microarray analysis. Each heat shock experiment was performed separately to enable precise control of the timing of the heat shock to ensure that the mRNA levels accurately reflected the stress response.
Sample preparation and microarray
Three independent cell cultures were used for the microarray analysis to acquire the mean expression under each culture condition. RNA sample preparation, microarray analysis with the Affymetrix GeneChip system, and data extraction were all based on the finite hybridization (FH) model [32, 33] and were performed as previously described [18, 20]. A single-nucleotide tilling array (high-density DNA microarray) that covers the entire Escherichia coli genome  was used. Thus, single-nucleotide substitution bias within ORF (open reading frame) regions could be disregarded in our evaluation of gene expression level. The results of 32 arrays for four bacterial strains in steady and responsive states were used for the analysis, which included triplicates of exponential growth at regular (36.9 °C) and evolution-stimulating (41.2, 43.2 or 44.8 °C) temperatures and duplicates of heat shock conditions.
Data normalization and annotations
Expression data sets from biological replicates were shown as mean values. The raw data sets were subjected to global normalization, resulting in a common median value of zero (logarithmic value) in all data sets, as previously described . The data sets of a total of 4383 genes were used in the analysis. Both the normalized expression data sets and the raw CEL files were deposited into the NCBI Gene Expression Omnibus database under the GEO Series accession number GSE61749. The full datasets of gene names, gene categories , and gene regulations (TF) were downloaded from GenoBase, Japan (https://dbarchive.biosciencedbc.jp/jp/genobase-v6/desc.html) and RegulonDB v8.0  (http://regulondb.ccg.unam.mx), as described previously [11, 18].
Binomial tests were performed to evaluate the statistical significance of extracted gene groups. These statistical analyses were carried out using free software packages available from the Broad Institute (http://www.broadinstitute.org). All statistical tests and computational analyses, except for the gene set enrichment analysis (GSEA), were performed using R . Gene sets enrichment analysis (GSEA)  was performed as previously described [11, 18]. The gene categories and gene regulatory links (TF) comprising more than 15 genes were used for the enrichment analysis and bimodal tests. Principal component analysis (PCA) [39, 40], which classifies expression patterns according to gene expression level variance between the culturing conditions, was performed as described previously [18, 20].
Results and discussion
Overview of transcriptome evolution for thermal adaptation
To identify changes in gene expression for thermal adaptation, four Escherichia coli strains were acquired during experimental evolution that we performed previously . We used an ancestor (Anc) population and evolved cell populations 41B, 43B and 45 L, which were adapted to the temperatures 41.2, 43.2 and 44.8 °C, respectively (Fig. 1a). Their growth rates were examined at both the regular (_r) and evolutionary high (_e) temperatures. The results showed that the evolved cells demonstrated highly recovered growth fitness at their corresponding evolutionary temperatures (Fig. 1b), consistent with our previous report . These growing cells were collected for microarray analysis, such that growth rates corresponded to transcriptomes at the exponential growing phases analyzed in this study (Additional file 1: Figure S1). In addition, we acquired the transcriptomes of cells undergoing a heat shock response (_hs, a transient responsive state induced by a temperature increase from 36.9 to 44.8 °C) to identify whether thermal adaptation changed responsivity to a temperature increase (Additional file 1: Figure S1).
An overview dendrogram based on hierarchical clustering analysis clearly shows that global gene expression patterns could be divided into two main clusters, the responsive and the steady states, which corresponded to the transcriptomes during heat shock response and the exponential growing phases, respectively (Fig. 1c). These results suggest that the transcriptome reorganization that occurs during thermal adaptive evolution differs from the reorganization that occurs during a heat shock response. In addition, the transcriptomes of cells in the early and the late evolutionary processes were roughly separated, within the steady-state cluster. This tendency suggests that genetic changes (i.e., genome mutations) played a role in the transcriptional changes during the evolutionary process at 43.2 °C (from 41B to 43B). According to the genome sequence analysis performed previously , a comparable number of genome mutations responsible for transcriptional regulation were newly accumulated in each strain (summarized in Additional file 1: Table S1). Overall, these results suggest that global transcriptome optimization might underlie transcriptome evolution.
Equivalent magnitudes of responsive changes
The enrichment analysis consistently found similar patterns of transcriptional change in most gene categories among all strains, although the genes in the categories of partial information and carrier showed significant differences between Anc and 45 L (Fig. 2b, left). A similarity in the patterns of transcriptional change was also observed for transcriptional regulation. Despite the small number of regulatory networks identified, such as the gene networks regulated by iscR and fhlA, most gene regulation patterns maintained identical directional changes among all strains (Fig. 2b, right). We note that such comparable changes in gene expression were also found in the transcriptional networks controlled by the global regulators harboring nonsynonymous substitutions or deletion mutations (e.g., lrp and oxyR; Additional file 1: Table S1). These results indicate that thermal adaptive evolution did not disturb the mechanism of heat shock response and that the evolved competence at high temperatures might have been acquired using other pathways.
Differentiated transcriptomes at steady states
Gene set enrichment analysis (GSEA) showed that transcriptome reorganization for thermal adaptation did not randomly occur but followed directional changes (Fig. 3c). The changes in gene expression compared in Fig. 3a represent the transcriptome differences derived from genotypes (i.e., the genomes harboring all the mutations fixed during the thermal adaptive evolution) and similar differences in Fig. 3b are derived from the growth temperatures. Gene function enrichment analysis was performed on the two types of differences, the genotype- and the growth temperature-mediated changes. Roughly, the changes in gene expression mediated by the genotypes were inversely correlated to changes mediated by the growth temperatures (Fig. 3c). The resultant patterns could be manually divided into four clusters (C1–C4). C1 showed significant inverse correlations; C3 and C4 showed gradual directional changes. Such directional changes were also found in regulatory networks (Additional file 1: Figure S3). The inverse transcriptional changes reflected a homeostasis of the transcriptome, which is supported by the finding that the expression of genes responding to a temperature increase (i.e., heat shock genes) returned to normal levels in all evolved cells, regardless of the high temperatures (Additional file 1: Figure S2B).
Thermal adaptive transcriptomes between the regular and the heat shock states
Functional enrichment analysis was performed for the genes positively and negatively loaded on the four main PCs (top 5 % each, 439 genes in each PC). The enrichment gene analysis (p < 0.001) failed to identify any essential functions (Fig. 4b, Additional file 1: Figure S4). The result suggests that transcriptome adaptation is not simply achieved by a central function of gene activity but is largely mediated by genes with diverse functions. The enrichment of transcriptional networks analysis resulted in largely varied transcription factor regulatory links (TF, p < 0.001) (Fig. 4b). The results showed that only a few regulatory networks were enriched in the most weighted components, PC1 and PC2, indicating that transcriptional networks participated as hubs in the transcriptome, and considerably contributed to thermal adaptation. In comparison, more regulatory links were enriched in the lower weighted components, PC3 and PC4. Differentially expressed regulatory networks were most abundantly identified in PC4, which represented variation in genotypes. This is consistent with the fact that the mutated regulators were largely identified in PC4 and, furthermore, demonstrated that genome mutations (lrp, rpoD and oxyR) did contribute to transcriptional changes at regulatory levels. In addition, of the mutated genes, most were nonsynonymous or InDel (insertion or deletion) mutations, found at considerable frequencies in PC1 and PC4 (Additional file 1: Figure S6). These findings indicate that genome evolution shaped the transcriptomes at PC1, which correlates with growth direction, and at PC4, which correlates with genotype.
Negative epistasis in transcriptome reorganization
Despite the weak correlations between genotype and heat shock, a negative epistasis in transcriptome reorganization was clearly detected (Fig. 5b). An approximate 20–30 % reduction in transcriptional changes was estimated to be the result of the simultaneous occurrence of genome mutations and heat shock, compared to their additive occurrence. This result is highly consistent with a previous finding that reported a ~30 % reduction in transcriptional changes . The negative epistasis was most significant for genes that were highly weighted in the main PCs, e.g., a ~50 % reduction in overlapping genes from PCs1–4 (Additional file 1: Table S2). This negative epistasis explains the fact that changes in gene expression were largely inversely related to genotype and temperature (Fig. 3c). Taken together, the changes in gene expression attributed to genetic alteration tended to cancel the changes attributed to the environmental perturbation. Such annulment might enable the cells to maintain the cellular homeostasis of the transcriptome. Although the underlying mechanism was not clearly resolved, our results suggest that adaptive evolution might involve a process of repressing gene expression changes triggered by multiple factors.
A potential path of transcriptome evolution for thermal adaptive
Thermal adaptation in Escherichia coli cells differentiated steady-state transcriptomes under both normal and evolution-stimulating high temperatures; nevertheless, the heat shock response was maintained with a high sensitivity to thermal stress, as found in the ancestor population. Thermal adaptive evolution directed the transcriptomes to novel steady states different from regular steady states or responsive states. These findings indicate that long-term evolution does not alter existing response machineries but rather adjusts gene expression homeostasis to adaptive steady states.
Availability of supporting data
The raw data is available at the GEO database under:
We thank Junko Asada for technical assistance. This work was partially supported by JST KAKENHI No. 23231061 (to TY) and in part by a research grant from the Institute for Fermentation, Osaka, Japan (to BWY).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Elena SF, Lenski RE. Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation. Nat Rev Genet. 2003;4(6):457–69.View ArticlePubMedGoogle Scholar
- Barrick JE, Lenski RE. Genome dynamics during experimental evolution. Nat Rev Genet. 2013;14(12):827–39.PubMed CentralView ArticlePubMedGoogle Scholar
- Philippe N, Crozat E, Lenski RE, Schneider D. Evolution of global regulatory networks during a long-term experiment with Escherichia coli. Bioessays. 2007;29(9):846–60.View ArticlePubMedGoogle Scholar
- Barrick JE, Yu DS, Yoon SH, Jeong H, Oh TK, Schneider D, et al. Genome evolution and adaptation in a long-term experiment with Escherichia coli. Nature. 2009;461(7268):1243–7.View ArticlePubMedGoogle Scholar
- Kishimoto T, Iijima L, Tatsumi M, Ono N, Oyake A, Hashimoto T, et al. Transition from positive to neutral in mutation fixation along with continuing rising fitness in thermal adaptive evolution. PLoS Genet. 2010;6(10):e1001164.PubMed CentralView ArticlePubMedGoogle Scholar
- Wang L, Spira B, Zhou Z, Feng L, Maharjan RP, Li X, et al. Divergence involving global regulatory gene mutations in an Escherichia coli population evolving under phosphate limitation. Genome Biol Evol. 2010;2:478–87.PubMed CentralView ArticlePubMedGoogle Scholar
- Patten CL, Kirchhof MG, Schertzberg MR, Morton RA, Schellhorn HE. Microarray analysis of RpoS-mediated gene expression in Escherichia coli K-12. Mol Genet Genomics. 2004;272(5):580–91.View ArticlePubMedGoogle Scholar
- Durfee T, Hansen AM, Zhi H, Blattner FR, Jin DJ. Transcription profiling of the stringent response in Escherichia coli. J Bacteriol. 2008;190(3):1084–96.PubMed CentralView ArticlePubMedGoogle Scholar
- Brockmann-Gretza O, Kalinowski J. Global gene expression during stringent response in Corynebacterium glutamicum in presence and absence of the rel gene encoding (p)ppGpp synthase. BMC Genomics. 2006;7:230.PubMed CentralView ArticlePubMedGoogle Scholar
- Jozefczuk S, Klie S, Catchpole G, Szymanski J, Cuadros-Inostroza A, Steinhauser D, et al. Metabolomic and transcriptomic stress response of Escherichia coli. Mol Syst Biol. 2010;6:364.PubMed CentralView ArticlePubMedGoogle Scholar
- Ying BW, Seno S, Kaneko F, Matsuda H, Yomo T. Multilevel comparative analysis of the contributions of genome reduction and heat shock to the Escherichia coli transcriptome. BMC Genomics. 2013;14:25.PubMed CentralView ArticlePubMedGoogle Scholar
- Regenberg B, Grotkjaer T, Winther O, Fausboll A, Akesson M, Bro C, et al. Growth-rate regulated genes have profound impact on interpretation of transcriptome profiling in Saccharomyces cerevisiae. Genome Biol. 2006;7(11):R107.PubMed CentralView ArticlePubMedGoogle Scholar
- Castrillo JI, Zeef LA, Hoyle DC, Zhang N, Hayes A, Gardner DC, et al. Growth control of the eukaryote cell: a systems biology study in yeast. J Biol. 2007;6(2):4.PubMed CentralView ArticlePubMedGoogle Scholar
- Brauer MJ, Huttenhower C, Airoldi EM, Rosenstein R, Matese JC, Gresham D, et al. Coordination of growth rate, cell cycle, stress response, and metabolic activity in yeast. Mol Biol Cell. 2008;19(1):352–67.PubMed CentralView ArticlePubMedGoogle Scholar
- Fazio A, Jewett MC, Daran-Lapujade P, Mustacchi R, Usaite R, Pronk JT, et al. Transcription factor control of growth rate dependent genes in Saccharomyces cerevisiae: a three factor design. BMC Genomics. 2008;9:341.PubMed CentralView ArticlePubMedGoogle Scholar
- Nahku R, Valgepea K, Lahtvee PJ, Erm S, Abner K, Adamberg K, et al. Specific growth rate dependent transcriptome profiling of Escherichia coli K12 MG1655 in accelerostat cultures. J Biotechnol. 2010;145(1):60–5.View ArticlePubMedGoogle Scholar
- You C, Okano H, Hui S, Zhang Z, Kim M, Gunderson CW, et al. Coordination of bacterial proteome with metabolism by cyclic AMP signalling. Nature. 2013;500(7462):301–6.PubMed CentralView ArticlePubMedGoogle Scholar
- Matsumoto Y, Murakami Y, Tsuru S, Ying BW, Yomo T. Growth rate-coordinated transcriptome reorganization in bacteria. BMC Genomics. 2013;14:808.PubMed CentralView ArticlePubMedGoogle Scholar
- Lopez-Maury L, Marguerat S, Bahler J. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nat Rev Genet. 2008;9(8):583–93.View ArticlePubMedGoogle Scholar
- Murakami Y, Matsumoto Y, Tsuru S, Ying BW, Yomo T. Global coordination in adaptation to gene rewiring. Nucleic Acids Res. 2015;43(2):1304–16.PubMed CentralView ArticlePubMedGoogle Scholar
- Cooper TF, Rozen DE, Lenski RE. Parallel changes in gene expression after 20,000 generations of evolution in Escherichia coli. Proc Natl Acad Sci U S A. 2003;100(3):1072–7.PubMed CentralView ArticlePubMedGoogle Scholar
- Riehle MM, Bennett AF, Lenski RE, Long AD. Evolutionary changes in heat-inducible gene expression in lines of Escherichia coli adapted to high temperature. Physiol Genomics. 2003;14(1):47–58.View ArticlePubMedGoogle Scholar
- Fong SS, Joyce AR, Palsson BO. Parallel adaptive evolution cultures of Escherichia coli lead to convergent growth phenotypes with different gene expression states. Genome Res. 2005;15(10):1365–72.PubMed CentralView ArticlePubMedGoogle Scholar
- Le Gac M, Cooper TF, Cruveiller S, Medigue C, Schneider D. Evolutionary history and genetic parallelism affect correlated responses to evolution. Mol Ecol. 2013;22(12):3292–303.View ArticlePubMedGoogle Scholar
- Puentes-Tellez PE, Kovacs AT, Kuipers OP, van Elsas JD. Comparative genomics and transcriptomics analysis of experimentally evolved Escherichia coli MC1000 in complex environments. Environ Microbiol. 2014;16(3):856–70.View ArticlePubMedGoogle Scholar
- Vijayendran C, Barsch A, Friehs K, Niehaus K, Becker A, Flaschel E. Perceiving molecular evolution processes in Escherichia coli by comprehensive metabolite and gene expression profiling. Genome Biol. 2008;9(4):R72.PubMed CentralView ArticlePubMedGoogle Scholar
- Pelosi L, Kuhn L, Guetta D, Garin J, Geiselmann J, Lenski RE, et al. Parallel changes in global protein profiles during long-term experimental evolution in Escherichia coli. Genetics. 2006;173(4):1851–69.PubMed CentralView ArticlePubMedGoogle Scholar
- Guisbert E, Yura T, Rhodius VA, Gross CA. Convergence of molecular, modeling, and systems approaches for an understanding of the Escherichia coli heat shock response. Microbiol Mol Biol Rev. 2008;72(3):545–54.PubMed CentralView ArticlePubMedGoogle Scholar
- Yura T, Nakahigashi K. Regulation of the heat-shock response. Curr Opin Microbiol. 1999;2(2):153–8.View ArticlePubMedGoogle Scholar
- Hartl FU, Martin J. Molecular chaperones in cellular protein folding. Curr Opin Struct Biol. 1995;5(1):92–102.View ArticlePubMedGoogle Scholar
- Ying BW, Tsuru S, Seno S, Matsuda H, Yomo T. Gene expression scaled by distance to the genome replication site. Mol Biosyst. 2014;10(3):375–9.View ArticlePubMedGoogle Scholar
- Ono N, Suzuki S, Furusawa C, Agata T, Kashiwagi A, Shimizu H, et al. An improved physico-chemical model of hybridization on high-density oligonucleotide microarrays. Bioinformatics. 2008;24(10):1278–85.PubMed CentralView ArticlePubMedGoogle Scholar
- Ono N, Suzuki S, Furusawa C, Shimizu H, Yomo T. Development of a physical model-based algorithm for the detection of single-nucleotide substitutions by using tiling microarrays. PLoS One. 2013;8(1):e54571.PubMed CentralView ArticlePubMedGoogle Scholar
- Suzuki S, Ono N, Furusawa C, Kashiwagi A, Yomo T. Experimental optimization of probe length to increase the sequence specificity of high-density oligonucleotide microarrays. BMC Genomics. 2007;8:373.PubMed CentralView ArticlePubMedGoogle Scholar
- Riley M, Abe T, Arnaud MB, Berlyn MK, Blattner FR, Chaudhuri RR, et al. Escherichia coli K-12: a cooperatively developed annotation snapshot--2005. Nucleic Acids Res. 2006;34(1):1–9.PubMed CentralView ArticlePubMedGoogle Scholar
- Salgado H, Peralta-Gil M, Gama-Castro S, Santos-Zavaleta A, Muniz-Rascado L, Garcia-Sotelo JS, et al. RegulonDB v8.0: omics data sets, evolutionary conservation, regulatory phrases, cross-validated gold standards and more. Nucleic Acids Res. 2013;41(Database issue):D203–213.PubMed CentralView ArticlePubMedGoogle Scholar
- Ihaka R, Gentleman R. R: A language for data analysis and graphics. J Comput Graph Stat. 1996;5(3):299–314.Google Scholar
- Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.PubMed CentralView ArticlePubMedGoogle Scholar
- Venables WN, Ripley BD. Modern Applied Statistics with S. New York: Springer; 2002Google Scholar
- Mardia KV, Kent JT, Bibby JM. Multivariate Analysis. London: Academic Press; 1979.Google Scholar