A Systems Biology Analysis Unfolds the Molecular Pathways and Networks of Two Proteobacteria in Spaceflight and Simulated Microgravity Conditions

2.1. Gene set enrichment analysis of expression data

Gene set enrichment analysis (GSEA) was performed with GenePattern (Reich et al., 2006), a multitool platform for genomic analysis. The GSEA module in GenePattern was run with global gene expression data and corresponding gene sets. The gene sets (ontology) for both Pseudomonas and Salmonella were obtained from the KEGG database. The expression data sets (from ArrayExpress, see the result section for details) were converted to an appropriate file format as required in GenePattern. To generate the ranked list of genes, we applied a “signal-to-noise” metric for ranking, which may be represented by (M _microg − M _1g )/(SD _microg + SD _1g ), where M _microg and M _1g represented the mean expression values of individual genes among experiments in microgravity and 1g conditions, respectively. The SD represented the corresponding standard deviation. In one case (Pseudomonas in RWV E-GEOD-16970), we used a “difference in class” metric (difference in mean values of gene expression) due to a lower number of samples (<3). For enrichment analysis, we applied running sum statistics in order to use a weighted Kolmogorov–Smirnov-like statistic scoring scheme (as described in Subramanian et al., 2005). In brief, the score was evaluated by walking down the ranked gene list (RL), and if genes in a gene set (GS) matched (hit) in the ranked list, the score would increase. It would decrease in the case of a nonmatch (miss). \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ P_ { { \rm { hit } } } } \left( { { \rm { GS } } , i } \right) = \mathop \sum \limits_ { \scriptstyle { g_j } \, \in\, { \rm { GS } } \hfill \atop \scriptstyle j \le i \hfill } { { \frac { { { \left\vert { { r_j } } \right\vert } ^p } } { { N_R } } } } \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{ P_ { { \rm { miss } } } } \left( { { \rm { GS } } , i } \right) = \mathop \sum \limits_ { \scriptstyle { g_j } \, \notin\, { \rm { GS } } \hfill \atop \scriptstyle j \le i \hfill } { \frac { 1 } { { N - { N_H } } } } \end{align*} \end{document}

where g_j represents the j ^th gene in RL = {g_i , … , g_N }, r(g_j ) = r_j p = 1 and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}{N_R} = \mathop \sum \limits_{j \, \in\, { \rm{GS}}} {{{ \left\vert {{r_j}} \right\vert }^p}} \end{align*} \end{document}

The enrichment score was calculated as the highest deviation of P _hit − P _miss from zero. To estimate the p value through null hypothesis testing, genes of gene sets were randomly permuted 1000 times, keeping the total number of genes in each set constant, and GSEA was performed for each permuted gene set. Normalized enrichment score (NES) values from this permuted set were compared with the original one to estimate the p values for each NES. The corresponding false discovery rate (FDR) q values were estimated from the distribution of the p values. We ran the same GSEA on each data set three times to check the consensus among various runs.

2.2. Consensus non-negative matrix factorization—clustering of gene expression data

The consensus non-negative matrix factorization (CNMF) algorithm was run on an expression data set of leading edge genes in the NMFC module of GenePattern. Clustering was based on the expression levels on the 1g and microgravity (simulated or space) conditions. The details of CNMF can be found elsewhere (Brunet et al., 2004). Briefly, taking the gene expression data set as a positive matrix L of size M × N, with a number of clusters C, the non-negative matrix factorization (NMF) computes iteratively for matrices F and G such that L = FG with F × C and C × G sizes. Starting from a randomly initialized matrix of size F and G, the matrices are iteratively updated to minimize a divergent functional as described in Lee and Seung (1999). We set the number of iterations per clustering as 2000 times. A consensus algorithm is run on top of the NMF to evaluate the mathematical stability for the number of clusters (C). For a given number of NMF clusters, the original data are perturbed and simulated by applying an in-built resampling tool in GenePattern. This assesses a consensus among multiple runs of NMF clustering on those simulated (perturbed) data. The resultant consensus matrix (summary of the consensus) is evaluated for a cophenetic coefficient, which determines the pairwise distance from the original data set. The highest number of mathematically stable clusters is determined from the change of the cophenetic correlation coefficient as a function of the number of clusters. The point after which the coefficient value falls sharply is estimated as the maximum possible number of mathematically stable clusters.

2.3. Mapping the genes in a protein network database

The genes that belong to each stable cluster were mapped on the STRING protein association network database to extract protein-protein association (PPA) maps. The mapping was performed by simply uploading the list of genes from each cluster on the STRING database and assigned the associated organism for those genes. This resulted in PPA networks. On each network, the interaction enrichment score and associated KEGG pathways with p and q value were determined from the in-built functional tools.

2.4. Orthology analysis

The orthology analysis of genes from 15 common altered pathways between Pseudomonas aeruginosa PAO1 and Salmonella enterica serovar Typhimurium was performed by using STRING. Taking Pseudomonas as a reference, the orthologus genes in Salmonella were searched, and the degree of protein sequence conservation between these two species was reported as a heat map.

3.1. GSEA reveals cellular and metabolic pathways altered in microgravity

The GSEA algorithm (Subramanian et al., 2005) uses a running sum statistics on the global gene expression data based on known functional sets of genes (ontology) to determine whether particular cellular or metabolic pathways are more likely to be involved with microgravity in a statistically significant way. The gene expression data for P. aeruginosa and S. typhimurium under microgravity (Wilson et al., 2007, 2008; Crabbe et al., 2010, 2011) were extracted from ArrayExpress (EMBL-EBI). The data accession numbers of these studies are E-GEOD-22684 (Crabbe et al., 2011), E-GEOD-16970 (Crabbe et al., 2010), and E-GEOD-8573 (Wilson et al., 2007, 2008). Crabbe et al. (2011) examined the effect of spaceflight and RWV on Pseudomonas. Crabbe et al. (2010) examined the effect of RWV and RPM on Pseudomonas. Similarly, Wilson et al. (2007) and Wilson et al. (2008) examined the effect of spaceflight conditions on Salmonella in rich and minimal media, respectively. To our knowledge, those are the only global gene expression data in microgravity for both species, which is available in the public domain. We then classified all genes of known cellular and metabolic functions of P. aeruginosa and S. typhimurium into 102 and 103 gene sets, respectively, by using the KEGG database. Next, we performed the GSEA, which determines whether a gene set is enriched on the top or bottom of a ranked list of all genes in a gene expression data set. Genes were ranked according to difference in the mean expression between microgravity and normal gravity divided by the summation of the standard deviation of that particular gene among replicates in each of the two experimental conditions (Fig. 2A). The histogram of the ranked genes with distribution of the score is shown in Fig. 2B. Next, enrichment scores for each gene set were calculated with running sum statistics, where the score increases when genes in the gene set are encountered in the ranked list (Fig. 2C). NES values were calculated by normalization of the enrichment scores based on the number of genes in each gene set. Unlike differential gene expression analysis, the p value and FDR q value are not assigned for an individual gene in the ranked gene list. The p value and FDR q value of an entire gene set are calculated in GSEA. A significance test was performed by estimating p values for these NES values by using a null hypothesis approach (Subramanian et al., 2005), and the FDR q was calculated by comparing the tails of NES distributions for each gene set. The resultant gene sets with positive and negative NES are defined as upregulated and downregulated, respectively.

OPEN IN VIEWER

FIG. 2.

GSEA of P. aeruginosa under spaceflight conditions. (A) Heat map of genome-wide expression for each biological replicate in space and on Earth. Genes are ordered according to a ranked list of genes that were classified based on the difference in expression between space and Earth gravity conditions (difference in the mean gene expression value between space and Earth is divided by the standard deviation in gene expression among the replicates in space and on Earth). (B) Gene correlation profile of the ranked gene list with space and Earth. (C) Representation of statistically significant pathways along with their enrichment scores and the position of the pathway genes in the ranked list. The heat map in the bottom of each pathway represents the gene correlation profile as shown in (B). (Color graphics available at www.liebertonline.com/ast)

Figure 2C shows the GSEA for two upregulated (ABC transporter and valine, leucine, isoleucine degradation) and two downregulated (ribosome and pyrimidine metabolism) pathways in P. aeruginosa under spaceflight conditions. For each set of experiments, a large number of cellular and metabolic pathways (gene sets) associated with microgravity conditions were altered (upregulated and downregulated) with statistically significant p and q values (p < 0.05, q < 0.25). The high NES values (absolute) were generally associated with high statistical significance (Supplementary Fig. S1; Supplementary Data are available online at www.liebertonline.com/ast), suggesting the robustness of the analysis. The numbers of altered pathways for each data set were shown in Fig. 3. Those pathways were tabulated in Table 1. The detailed NES, p value, and q value for each of those pathways in each experiment were tabulated in Supplementary Table S1.

OPEN IN VIEWER

FIG. 3.

Number of altered cellular and metabolic pathways identified by GSEA on different studies. The numbers 1, 2, 3, and 4 with each experiment represent data sets analyzed from the references Crabbe et al. (2010), Crabbe et al. (2011), Wilson et al. (2007), and Wilson et al. (2008), respectively. (Color graphics available at www.liebertonline.com/ast)

3.2. Comparative analysis suggests involvement of new and overlapping pathways in microgravity

For P. aeruginosa under simulated microgravity, a previous study (Crabbe et al., 2010) suggested seven altered functional categories, namely, AlgU and Hfq regulon, quorum sensing, motility, stress, anaerobic metabolism, and ATP synthesis. Similarly, in spaceflight conditions (Crabbe et al., 2011), the suggested categories are Hfq regulon, anaerobic metabolism, virulence factor, ribosome synthesis, and dehydrogenation of succinate. Two studies on S. typhimurium under space conditions (Wilson et al., 2007, 2008) showed involvement of Hfq regulons, ribosome, iron utilization, biofilm, virulence, motility, and several genes without any assigned pathways. Our analysis showed a higher number of altered pathways, which include a large number of cellular and metabolic pathways that have not been reported previously. The new pathways found in this study are tabulated in italics in Table 1. We observed that in the work of Crabbe et al. (2011) five KEGG pathways were altered. In our analysis, with the same data set, we found 37 altered pathways, which included all five pathways indicated by Crabbe et al. (2011). It was interesting to find that, in all four cases, metabolic pathways outnumbered the cellular processes.

However, GSEA runs only for the genes that are assigned to functional classes. Therefore, unclassified proteins (around 50% of the total genes including hypothetical proteins) in KEGG are not enriched in our analysis. Though the Hfq gene was found to be downregulated in the RNA degradation pathway for both species in our results, there is no functional class assigned for Hfq regulon genes in KEGG. When we added two gene sets—Hfq regulon genes upregulated (480 genes) and Hfq regulon genes downregulated (334 genes) for Pseudomonas (Sonnleitner et al., 2006), which were used by Crabbe et al. (2010) and Crabbe et al. (2011) to identify the Hfq regulon genes in microgravity—GSEA enriched the downregulated and upregulated Hfq regulons in space with p < 0.001 and q ≤ 0.1. However, these two gene sets are broad in nature and are not rigorously curated or specific, and were therefore excluded from further analysis. For Salmonella in space, the Hfq regulons were divided in the multiple functional categories (Wilson et al., 2007). KEGG does not classify any of those functions except ribosomal proteins, and our analysis showed downregulation of ribosomal proteins as was the case in the work of Wilson et al. (2007).

Next, to discern common pathways related to microgravity, if there is some association, we compared the number of overlapping pathways among different data sets (Fig. 4). It was observed in our analysis that the number of overlapping pathways between P. aeruginosa and S. typhimurium in space was much higher than the number between P. aeruginosa in space and simulated microgravity (RWV and RPM) conditions. Ribosome, RNA degradation, protein export, flagellar assembly, methane metabolism, toluene degradation, oxidative phosphorylation, TCA cycle, glycolysis, purine metabolism, and pyrimidine metabolism were among the common altered pathways between P. aeruginosa and S. typhimurium in space. We tested the degree of protein sequence conservation in pathway-specific orthologous genes for all the 15 common altered pathways between the two species (Supplementary Fig. S2). The results show a high degree of sequence conservation. The effect of spaceflight on both the species, therefore, may have a common cellular origin.

OPEN IN VIEWER

FIG. 4.

Venn diagrams among the experiments represent the common altered pathways in microgravity. (A) P. aeruginosa in space, RWV, and RPM. Valine, leucine, and isoleucine degradation pathway was common among space, RWV, and RPM. Flagellar assembly, protein export, ribosome, RNA degradation, TCA cycle, and oxidative phosphorylation were increased in RWV and RPM but were decreased in spaceflight conditions. (B) P. aeruginosa in space, S. typhimurium in space, rich media, and S. typhimurium in space, minimal media. Ribosome, RNA degradation, protein export, flagellar assembly, methane metabolism, toluene degradation, oxidative phosphorylation, TCA cycle, glycolysis, purine metabolism, and pyrimidine metabolism were among the common altered pathways. The numbers 1, 2, 3, and 4 with each experiment represent data sets analyzed from the references Crabbe et al. (2011), Crabbe et al. (2010), Wilson et al. (2007), and Wilson et al. (2008), respectively. (Color graphics available at www.liebertonline.com/ast)

We found that the valine, leucine, and isoleucine degradation pathway was common among space, RWV, and RPM for P. aeruginosa (Fig. 4A). In contrast, several of the pathways, which were upregulated in simulated microgravity for P. aeruginosa, were downregulated in space for both the species (Table 1). The gene expression of flagellar assembly, protein export, ribosome, RNA degradation, TCA cycle, and oxidative phosphorylation was increased in RWV and RPM, but was decreased in spaceflight conditions.

3.3. CNMF clusters the leading edge genes based on the gene expression patterns

Next, we performed a leading edge subset analysis to identify the significant genes (leading edge genes) that contributed the most for a pathway to be statistically significant. The mapped genes before the peak value of enrichment scores in Fig. 2C are the leading edge genes. Those genes for all experiments were tabulated in Supplementary Table S2. To evaluate whether a set of patterns in leading edge gene expression exists, we performed an unsupervised CNMF clustering method (Lee and Seung, 1999; Brunet et al., 2004) for the experiments, where the combined number of leading edge genes among altered pathways was more than 100. The parts-based local representation is a unique property of NMF that allows one to extract functional and more meaningful subsets from a matrix data set in comparison to other methods such as principal component analysis (PCA) and vector quantization (VQ) (Lee and Seung, 1999). NMF has also outperformed all conventional clustering methods including K-means and singular value decomposition (SVD) in extracting functional relationships in yeast gene expression data (Kim and Tidor, 2003). A consensus clustering algorithm together with NMF (CNMF) was successfully applied to classify cancer samples from global gene expression data, where other clustering methods like hierarchical clustering, self-organizing maps, and principle component analysis have failed (Brunet et al., 2004). This has been used successfully in effectively interpreting several gene expressions and other large-scale biological data (Devarajan, 2008). CNMF was applied separately for upregulated and downregulated pathways. We identified the maximum number of mathematically stable clusters in each case by plotting the cophenetic coefficient as a function of the number of clusters (k). Figure 5A shows such a correlation derived from upregulated leading edge genes of P. aeruginosa in spaceflight conditions. The cophenetic coefficient determines the mathematical stability of the clusters. Its value declined sharply after a certain number of clusters, which is indicated by an arrow (Fig. 5A). This number suggests the maximum number of mathematically stable clusters, and we can then hypothesize the possible number of patterns in gene expression. The cophenetic coefficient correlation plots for other cases are shown in Supplementary Fig. S3. Figure 5B shows the maximum number of stable clusters for all the studies. We have shown that the CNMF was able to distinguish the expression patterns between ribosomal genes and various metabolic genes by classifying them into different clusters in each and every case (Fig. 6, Supplementary Figs. S5 and S6). This validated that the unsupervised CNMF clustering was able to differentiate the gene expression patterns in a biologically meaningful way.

OPEN IN VIEWER

FIG. 5.

CNMF clustering on leading edge genes. (A) The cophenetic coefficient as a function of the number of clusters, where the CNMF algorithm was applied on the expression profile of upregulated leading edge genes of P. aeruginosa in spaceflight. The arrow indicates the maximum number of mathematically stable clusters possible in this particular example. (B) Maximum number of mathematically stable clusters for other experiments. Clusters with upregulated and downregulated genes are shown in black and red bars, respectively. P. aeruginosa in space, P. aeruginosa in RWV, S. typhimurium in space (rich media), and S. typhimurium in space (minimal media) represent the data analyzed from the references Crabbe et al. (2011), Crabbe et al. (2010), Wilson et al. (2007), and Wilson et al. (2008), respectively. (Color graphics available at www.liebertonline.com/ast)

OPEN IN VIEWER

FIG. 6.

Protein-protein association networks of P. aeruginosa from leading edge genes of an individual CNMF cluster. Each node represents a protein, and the line connecting the nodes (edge) represents the functional association. The relative thickness of each line signifies the confidence level of such association within a network. The functional annotations of different parts of the network are marked on the figure. (A) Cluster 1 from upregulated pathway genes in space condition. (B) Cluster 1 from downregulated pathway genes in space condition. (C) Cluster 2 from downregulated pathway genes in space condition.

3.4. Mapping of the patterned gene expression reveals the protein networks affected by microgravity

To understand the plausible regulatory and functional connections among genes, pathways, and microgravity, we mapped the genes of distinct CNMF cluster(s) with the STRING protein network database (Franceschini et al., 2013). Figure 6 shows the PPA networks for one upregulated and two downregulated clusters of P. aeruginosa in spaceflight conditions. The network maps for all other studies are shown in Supplementary Figs. S4–S7. It was observed that the clusters showed either highly interactive protein-protein networks (p < 10⁻⁸) or almost no interactions (Table 2 and Figs. 6 and S4–S7). The interaction networks were further enriched with KEGG pathways by using a built-in tool within STRING. The built-in tool requires at least nine interactions to assign any function. Five or higher numbers of interactive proteins and the enrichment p value < 0.001 were taken for further analysis. The details of the PPAs of all clusters with underlying functions are tabulated in Table 2. We observed the network connections between benzoate degradation and butanoate metabolism, between ABC transporters and two-component systems, between purine and pyrimidine metabolism and flagellar assembly, and between flagellar assembly and protein exports in various networks in P. aeruginosa (Tables 1 and 2 and Fig. 6). We have shown network connections between glycolysis and TCA cycle in S. typhimurium and between oxidative phosphorylation and flagellar assembly pathways for both P. aeruginosa and S. typhimurium. Interestingly, all those connections have been observed in other studies (Saier and Ramseier, 1996; Tobisch et al., 1999; Kapatral et al., 2004; Macnab, 2004; Maurer et al., 2005; Metzler-Zebeli et al., 2010; Dintner et al., 2011; Rooks et al., 2014) with several bacterial species not related to microgravity.