Discovery and predictive modeling of urine microbiome, metabolite and cytokine biomarkers in hospitalized patients with community acquired pneumonia

Pneumonia is the leading cause of infectious related death costing 12 billion dollars annually in the United States alone. Despite improvements in clinical care, total mortality remains around 4%, with inpatient mortality reaching 5–10%. For unknown reasons, mortality risk remains high even after hospital discharge and there is a need to identify those patients most at risk. Also of importance, clinical symptoms alone do not distinguish viral from bacterial infection which may delay appropriate treatment and may contribute to short-term and long-term mortality. Biomarkers have the potential to provide point of care diagnosis, identify high-risk patients, and increase our understanding of the biology of disease. However, there have been mixed results on the diagnostic performance of many of the analytes tested to date. Urine represents a largely untapped source for biomarker discovery and is highly accessible. To test this hypothesis, we collected urine from hospitalized patients with community-acquired pneumonia (CAP) and performed a comprehensive screen for urinary tract microbiota signatures, metabolite, and cytokine profiles. CAP patients were diagnosed with influenza or bacterial (Streptococcus pneumoniae and Staphylococcus aureus) etiologies and compared with healthy volunteers. Microbiome signatures showed marked shifts in taxonomic levels in patients with bacterial etiology versus influenza and CAP versus normal. Predictive modeling of 291 microbial and metabolite values achieved a + 90% accuracy with LASSO in predicting specific pneumonia etiology. This study demonstrates that urine from patients hospitalized with pneumonia may serve as a reliable and accessible sample to evaluate biomarkers that may diagnose etiology and predict clinical outcomes.

Urine cytokines. We interrogated 34 cytokines in the urine of healthy controls, patients infected by influenza, S. pneumo, or S. aureus. Out of 34 cytokines tested, we detected 17 cytokines in the urine and 11 of those cytokines showed differential presentation among the four groups of participants (Table 1). Pneumonia caused by S. aureus differed from the healthy controls for all 11 cytokines and in general the level of cytokines detected in influenza patients was consistently lower than the levels in patients with bacterial pneumonia and only differed from healthy volunteers for IFNγ, IL-6, IL-18, eotaxin, IP-10 and MCP-1. Four of the cytokines demonstrated a significant difference between the 3 types of pathogens; IFNγ (P = 0.005), IL-18 (P = 0.0052), MCP-1 (P = 0.0029) and SDF-1 (P = 0.0451). The remaining cytokines that were present in the urine but did not differ from each other or healthy volunteers were IL-1β, IL-1α, IL-1RA, IL-22, IL-27 and IL-8. We did not detect IL-12p70, IL-13, IL-2, IL-5, GM-CSF, IL-10, IL-17A, IL-21, IL-23, IL-9, IFNα, IL-31, IL-7, TNFβ, MIP-1β, RANTES or TNFα in the urine in any of the groups.
Urine microbiome. Total DNA was extracted from urine samples for quantitative PCR of 16S copy numbers, which demonstrated that healthy volunteers exhibited higher bacterial DNA copy numbers compared with bacterial or viral pneumonia patient urine samples (Fig. 1A); potentially influenced by the initial administration of antibiotics that all cause pneumonia patients receive upon admission. Similarly, compared with healthy  www.nature.com/scientificreports/ ing between healthy and influenza samples, while patients infected with bacterial pathogens, S. aureus and S. pneumoniae, clustered even more distinctly from healthy volunteer samples (RDA significance: Variance = 12.74, F = 1.13, P = 0.05) (Fig. 2B). A network analysis of detected taxa also demonstrated more tightly clustering of taxa found in healthy volunteer associated (Fig. 2C). Specific taxa associated with healthy volunteers included the genus Lactobacillus within Firmicutes, while pneumonia patients demonstrated elevated levels of the class Gammaproteobacteria, family Enterobacteriaceae, and the genera Clostridium and Streptococcus (Figs. 2D; S1B, S2A,B). Hierarchical taxonomic composition for each patient group are summarized in Fig. S1. To identify specific taxanomic differences between groups, we further employed linear discriminate analysis of effect size (LEfSe) between experimental groups. Initially, we compared all 30 patients with CAP to healthy volunteers ( Fig. S3A). At the phylum level, Proteobacteria was identified as significant in pneumonia samples while Synergistetes was identified in healthy controls based on LDA scores. At the genus level, Clostridium and Sutterella were identified in case with CAP while Lactobacillus, Prevotella, Magasphaera, Dorea, Vibrio, and Coprococcus were the most significantly enriched. Cladogram projection of these differences demonstrated the CAP case samples clustered primarily within the phyla Proteobacteria, while the healthy volunteers were more taxonomically distributed (Fig. S3B). The analysis was repeated after regrouping samples based on viral vs bacterial pathogen, which showed changes at the order level, where Bifidobacteriales was abundant in healthy controls, Enterobacteriales was abundant in bacterial pneumonia samples, and Sphingomonadales was abundant in influenza samples (Fig. S3C). A final regrouping determined comparisons of healthy vs S. aureus and S. pneumoniae (Fig. S3D) and healthy vs influenza (Fig. S3E), where healthy samples consistently displayed greater levels of the family Rikenneliaceae and the order Bifidobacterium. Interestingly, the level of detectable Streptococcus was most elevated in S. pneumoniae case samples, while Syntrophobacter and Delftia were most elevated in patients with S. aureus (Fig. S2B).
Urine metabolites. Metabolites were extracted from 50 µl of urine and subjected to ultra-high-performance liquid chromatography coupled high resolution mass spectrometry (UPLC-HRMS) 87 known metabolites were detected in the forty urine samples and identified using known masses (± 5 ppm) and retention times (Δ ≤ 1.5 min). Creatinine is considered the best internal standard to correct for urine volume variations as its rate of elimination is independent of urine flow and urine volume and creatinine concentration are inversely proportional 35,36 . Thus, to ensure that observations were directly comparable peak intensity was normalized to creatinine. Then these data were compared to unnormalized data to make sure there was no masking of biologically relevant changes by normalization (DNS). As is the convention in metabolomics we first used unsupervised multivariate statistical analysis to determine the dataset structure and relationships between groups. www.nature.com/scientificreports/ To evaluate the group trends, sample uniformity and identify potential outliers, multivariant principal component analysis (PCA). The variation were explained by F1 and F2 with a cumulative percent variability of 78.56% spread among the patient groups (i.e. Healthy 53.5 and 6.1, IAV 13.8 and 70.0, S. aureus 8.6 and 22.4, and S. pneumoniae 30.0 and 1.2 percent per F1 and F2 respectively). Adding a third component marginally increased the cumulative percent variability to 85.96%. The two component PCA analysis shows good separation between CAP patients and the healthy group (Fig. 3, circles). Likewise, the high-risk classes (IV-V) and low risk (I-III) centroids showed clear separation (Fig. 3, squares). We then used unsupervised clustering of both metabolites and individuals, they were clustered independently using k-means clustering followed by ascendant hierarchical clustering based on Euclidian distances. The data matrix's was rearranged according to the corresponding clustering with spatial relationship proportional to similarity among patient samples or metabolites (Fig. 4). These clusters were also represented via a dendrogram displayed vertically for metabolites and another horizontally for patients. We find the healthy volunteers centered and groups nicely together (red) as did the IAV (Blue) while the bacterial pneumonia samples were interspersed together (gray and orange) (Fig. 4). Consistent with the PCA analysis the high-risk groups tended to be close together on the far left or right (Fig. 4, brown bars).
Next, we employed a simple one-way ANOVA with Tukey's honestly significant difference test (Tukey's HSD) with Benjamini-Hochberg post hoc correction (XLSTAT) to identify 6 metabolites with significant differences among patient groups ( Table 2). Adenosine 5′-phosphosulfate (APS) was the most significant metabolite with differential concentration based on pneumonia, it was significantly higher in the urine of healthy volunteers. In humans, all APS is converted to 3′-phosphoadenosine 5′-phosphosulfate (PAPS) for the sulfonation of glycosaminoglycans, proteins, peptides, lipids, bile acids, xenobiotics and steroids [37][38][39] . Guanidoacetic acid was also significantly higher in healthy volunteers and is a precursor to creatine, metabolite in the Urea cycle as well as metabolism www.nature.com/scientificreports/ of amino groups of several amino acids including glycine, serine, threonine, arginine and proline. 2,3-dihydroxybenzoate is a conjugate base of 2,3-dihydroxybenzoic acid that is increased after consumption of nutrients (e.g. cranberry juice) or aspirin and is also a biomarker of OH radicals 40,41 . Succinate was significantly decreased in patients with pneumonia in our studies as well as two other metabolite profiling studies of pneumonia from human pleural fluid and mouse urine infected with S. pneumoniae 18,42 . We found citrate and succinate, metabolites related to the citric acid cycle, to be significantly reduced in all three groups with CAP (Table 2). Reduced citrate levels have previously been reported in plasma of patients with pneumonia and in mouse urine 18,27 . Likewise, Adamko et al. observed decreases in both citrate and succinate in urine from children with bacterial and viral respiratory infections 18,27,43 . Conversely, we found uridine to be significantly increased in the urine of all patients with pneumonia. This is in keeping with previous reports of uridine transiently increasing in the lung and BAL fluid of mice with viral pneumonia from influenza 42 . It is worth noting these are highly abundant analytes whose values are relative to the peak intensity of creatinine in each sample (creatinine mean across samples was 2.753e + 009). Taken together these metabolites represent likely candidates for including among the signature biomarkers of pneumonia. www.nature.com/scientificreports/  www.nature.com/scientificreports/ We applied a supervised four component partial least squares discriminant analysis (PLS-DA) to distinguish between patient groups and identify differentially expressed variables. The correlation map of the first two components reveals a clear separation of the healthy individuals and group (solid and open grey circles respectively) from the pneumonia patients (Fig. 5A). The index values of the Variable Importance in Projection (VIP) from the PLS-DA were then used to identify 9 metabolites with VIP scores over one (Fig. 5B). However, the overall fit of this model was not robust (i.e. low Q 2 values), indicating the quality of the fit varies a lot depending on the metabolite. Likewise, the R 2 values were around 0.3 suggesting the components generated by the PLS regression did not summarize either the X or Y variables well. Thus, we revised this analysis by first centering and reducing the explanatory variables before starting the PLS-DA calculations (PLS-DA VCR ). The quality of the PLS-DA VCR was improved (i.e. Q 2 cumulative 0.083-0.378). While the Q 2 value is positive, thus has predictive relevance, it remains somewhat low suggesting the quality of the fit of this model varies a lot depending on the metabolite. The PLS-DA VCR also improved the regression's ability to summarize both the X and Y variables (i.e. R 2 Y 0.552 and R 2 X 0.483) resulting in better separation of pathogen groups (Fig. 5C). However, this produced a large number of metabolites, 35% of those identified, with VIP scores above 1 (Fig. 5D).
predictive modeling. There were 291 variables including two demographics such as gender and age and 185 OTUs detected in urine samples, 17 cytokines, and 87 metabolites of 40 subjects. Note that we included OTUs that were observed for at least two subjects. First, we implemented multi-class classification with 5-folds cross validation to distinguish between four subject categories using a total of 291 predictors. However, none of the machine learning model provided desirable accuracy (all < 47.5%). While none of the 6 metabolites identified by Tukey HSD failed the Dunnett's test, this method overlooked 12 metabolites that passed (Supp. Table 2; Fig. S4A). Given the post hoc correction method was for false discovery rate, it is not surprising that this expression analysis resulted in no Type I error but high levels of Type II Errors (Supp. Table 3). The initial PLS-DA www.nature.com/scientificreports/ identified nine metabolites with VIP > 1, of these only citrate and taurine showed significant differences (Supp. Table 2 and Fig. S4A). Further, the PLS-DA analysis produced the most, 16, Type II errors (Supp. Table 3). The PLS-DA VCR analysis identified 31 metabolites with VIP > 1 (Supp. Table 2). Nineteen of the metabolites identified with PLS-DA VCR analysis passed the individual analysis thereby improving the Type II errors when compared to the PLS-DA. However, it misidentified 13 metabolites resulting in the largest number of Type I errors of any of the models (Supp. Table 3). Lasso Model 1 identified seven potential biomarkers that distinguished healthy from CAP patients, all of which passed the Dunnett's test (Supp. Table 2). Model 1 also produced the least errors with more positive identifications (Supp. Table 3). Thus Lasso model 1 did not require reducing explanatory variables as was done in the PLS-DA VCR analysis, that resulted in the greatest level of Type I errors, while producing the least Type I or II errors. It is important to note that in the first iteration, our model pulled out several predictions that are significantly altered and represent abundant metabolite markers in the urine. We then implemented LASSO logistic regression with fivefold cross-validation and a total of 13 OTUs, 2 cytokines, and 13 metabolites were found to be discriminating between different subject categories as listed in Supp. Table 4.
Model 1 identified one OTUs and three metabolites to distinguish heathy subjects from S. aureus, S. pneumoniae, and Influenza (Model 1 providing AUC with 95% CI of 0.98; 0.94-1.00). Model 2 identified six OTUs, two cytokines, and four metabolites to distinguish Influenza from S. aureus and S. pneumoniae (Model 2 providing AUC of 1.00). Model 3 identified six OTUs, one cytokine, and six metabolites to distinguish S. aureus from S. pneumoniae (Model 3 providing AUC of 1.00). The confusion matrices with performance indicators for each model is presented in Table 3.
In recursive implementation of LASSO regression in three steps, we identified a total of 28 OTUs, cytokines, and metabolites to classify subjects into their actual categories. However, model 1 assumes the subject is not healthy and model 3 assumes the subjects not healthy nor influenza. Therefore, to develop a model that can be implemented on subjects without any assumption on their status, by using these selected 28 predictors, which is significantly smaller than the original dataset with 291 predictors, we readdressed the multi-class classification problem. We implemented various machine learning algorithms and found that Ensemble Method (Ensemble Method: Subspace, Learner Type: Discriminant, Number of Learners: 30, Subspace Dimension: 13) provided the highest overall classification accuracy of 85.0% (Table 4). Note that the parameter setting of the final model was fixed across folds and there was no parameter optimization implemented.  www.nature.com/scientificreports/ Table 4 shows that most of the classification error is due to misclassification between two bacterial groups. When we merge two bacterial groups (Table 5), we found that one can distinguish between healthy, influenza and bacterial categories with very high accuracy of 92.5% and with perfect positive predictive values for healthy and influenza subjects.

Discussion
The search for accurate predictors of infection and disease remains an important frontier in the era of omics and personalized medicine. In the setting of CAP, which is the leading cause of infectious related death in the United States, various methods have been employed using serum and other clinical samples to confirm and determine severity of respiratory infection, including the quantification of cytokines and eicosanoids. However, these approaches have achieved limited sucess and are not widely implemented for patient care. For the clinician, even the basic differential diagnosis between viral vs bacterial pathogens remains difficult in CAP cases due to the overlapping symptomatic presentations. Considering the dynamic nexus of host immunological, metabolic, and microbial networks, we moved beyond the search for single or limited numbers of biomarkers by instead comprehensively profiling urine cytokines, the microbiome, and the metabolome in samples collected from newly admitted pneumonia patients with either influenza or bacterial pathogens compared with healthy volunteers. Urine was chosen as a non-invasive sample since it is readily obtainable in the in-patient and out-patient setting.
There have been few studies to determine the utility of measuring cytokine markers in urine as potential biomarkers of infectious disease states. Out of the 34 cytokines we measured in urine, 11 were significantly different between groups. Patients with bacterial pneumonia exhibited the greatest elevation and number of cytokines in urine that differed from healthy controls; S. aureus pneumonia patients differed in all 11 cytokines with IL-4, IL-15, Gro-α, MIP-1α and SDF-1 uniquely elevated in those patients. Whereas patients with influenza exhibited the lowest levels of cytokines detectable in urine that differed compared to healthy controls. In several studies IL-6 levels in the serum have been found to correlate with either bacterial infections or disease severity in patients with CAP [44][45][46][47] . Our results also detected increased IL-6 in the urine which was elevated in all pneumonias compared to healthy controls. However, based on the predictive modeling, IL-15 and IL-18 in combination with the metabolome and microbiome data may be more useful for distinguishing bacterial vs viral pneumonias.
Analysis of the urine microbiome demonstrated a complex community under both healthy and infectious states, including the presence of slow growing Lactobacillus and Corynebacterium, consistent with recent reports on urine microbiome. While urinary tract bacterial populations were historically overlooked outside of the context of infection, next generation sequencing techniques enabled culture independent insights into these communities. Recent findings demonstrate major shifts in the urine microbiome community under diseases of the urinary tract, including 48 urolithiasis and certain cancers, however investigation of the urine microbiome as signatures of disease at distant body sites has not been employed. We observed elevated diversity and evenness in urine samples from pneumonia patients compared with healthy volunteers, including consistently elevated Enterobacteriaceae and Clostridium. Furthermore, community composition data suggested greater dissimilarity in the urine microbiome in patients with bacterial pneumonia than in those with influenza ( Fig. 2A-C).
Similar to the microbiome, the metabalome of CAP patients clustered from healthy volunteers (Figs. 3, 4), suggesting a divergence in metabolite profiles under an infectious state. Drivers of this differential clustering included loss of numerous metabolites, including citrate/isocitrate, succinate, guanidoacetic acid, N-acetylglutamine, among others, compared with healthy volunteers. No metabolites were specific to influenza infection alone compared with the other groups. However, methyladenosine, uridine and 2-dehydro-d-gluconate were elevated under bacterial infection compared with influenza and healthy controls. These divergent profiles could be useful in determining pathogen kingdom.
Since the microbiome and metabolite profiles are influenced by the environment, nutrition, age, and lifestyle of the host, in addition to genetics, these concatenated profiles provide a unique snapshot of individual health 19,[22][23][24][25][26][27][28][29][30][31] Indeed, while these complex profiles can be examined independently, changes in the collective abundance patterns of metabolites and microbes may indicate deeper homeostatic disturbances, which may be reflected through changes in interleukin signaling. The membership of the bodies microbial communities have dynamic interconnected relationships with one another and the host that change under states of disease and stress. Therefore, the microbiome and metabolome complement to serve as personalized readout of individual health. The ability to detect rapid measurable changes in these profiles in response to challenges, such as infection, would be a novel systems biology approach to personalized medicine. For instance, the components of the metabolome and microbiome are physically or stoichiometrically co-related, leaving precise abundance patterns that may accurately reflect discreet pathophysiological states 23,34 . Utilization of meta-biomarkers, such as the www.nature.com/scientificreports/ microbiome and metabolome, would represent a distinct shift away from the majority of clinical biomarkers currently in use, even in the era of genomics, transcriptomics, and proteomics 23,34 . After combining all urine meta-biomarkers, totaling 385 data points, we performed machine learning models by implementing LASSO logistic regression in three unique models. Model 1 aimed to distinguish healthy subjects from pneumonia; Model 2 to distinguish between bacterial (S. aureus and S. pneumoniae) pathogens from Influenza; and Model 3 to distinguish between S. aureus and S. pneumoniae. For each predictive model, we implemented a fivefold cross validation process to avoid overfitting. Specifically, the data were split into five distinct folds where 4 folds were used for model testing and the remaining for validation. By repeating this process five times by changing the test fold, we identified a total of 28 predictors, including two cytokines, 13 microbial taxa, and 13 metabolites that provided a predictive power of 92.5% in distinguishing patient groups.
There are several limitations to our study, including the total sample size of 40 individuals. While we were able to detect consistent changes in our meta-biomarkers, larger studies with greater numbers and groups that included other pathogens may improve the resolution of our predictive models in determining unique signatures of pneumonia or other respiratory diseases. Moreover, despite cross-validation yielded high predictive accuracy, there is a need to both validation of the data on a larger cohort and on a more diverse external cohort to be able to claim broader generalizability. Another limitation was related to a characteristic of the clinical standard of care, where all pneumonia patients in this study were placed on antibiotics upon admission to the hospital. Future studies may attempt to compare urine samples collected from individual before and after the implementation of antibiotics. On the other hand, the inclusion of patients with influenza acted as a unique control group for the bacterial groups, since all patients were placed on antibiotics. We observed large perturbations in the metabiomarkers in bacterial groups compared with the influenza group, suggesting that the changes were indeed driven by the pathogen and not a general response to infection. A final limitation that should be noted was the imbalance of gender between experimental groups, where the healthy volunteers were 90% female while the pneumonia groups were 30-50% female, and future work should place emphasis on larger and balanced gender composition between groups.
Here we describe a comprehensive profile of urine meta-biomarkers, including the microbiome, metabolome, and cytokines in pneumonia patients. Using these biomarkers, we achieved high success in predicting pneumonia pathogens. Depending on the infectious pathogen identified in each patient, distinctly different immune profiles were observed in cytokine profiles, and even larger shifts were observed in the metabolite and microbiome profile, especially in response to bacterial infections. This study provides a proof of concept that urine samples, which are easily accessible in outpatient and inpatient settings, could provide additional diagnostic insights to patient infectious status and future risk factor for complication.

Materials and methods
Sample processing. Urine samples were collected using sterile technique and were aliquoted separately for cytokines, microbiome analysis, and metabolites 15,16,36 . For multiplex cytokine assays and microbiome analysis, the urine samples were centrifuged at 10,000 × g for 10 min and then analyzed as described below. For Metabolite analysis, samples (50 μl for urine) were extracted with 1.3 ml of extraction solvent (40:40:20 HPLC grade methanol, acetonitrile, water with 0.1% formic acid) pre-chilled to 4 °C in a cold room and incubated for 20 min at − 20 °C. The samples were centrifuged for 5 min (16.1 rcf) at 4 °C. The supernatants were transferred to new 1.5 ml centrifuge tubes and pellets were resuspended with 200 µl of extraction solvent. Extraction was allowed to proceed for 20 min at − 20 °C and all supernatants collected in glass vials. Vials containing the collected supernatant were dried under a stream of N2 until all the extraction solvent had been evaporated. Residue was resuspended in 300 µl of sterile water and transferred to 300 µl autosampler vials. Samples were immediately placed in autosampler trays for mass spectrometric analysis.

Multiplex for cytokines and statistical analysis. The levels of a panel of inflammatory mediators in
urine samples were measured using a 34-plex ProcartaPlex Multiplex Immunoassay according to manufacturer's instructions (Invitrogen, Carlsbad CA, USA). Cytokine standards were prepared to determine the concentration of cytokines in the samples. The samples were run on a Millipore Magpix instrument and analyzed with xPO-NENT 4.2 software. For data analysis, a five-parameter logistic curve fitting method was applied to the standards and the sample concentrations extrapolated from the standard curve. The results were normalized to creatinine as measured by Creatinine Detection Kit (Thermofisher, Waltham, MA).
Kruskal-Wallis test was used to compare median level of cytokine among the four groups of samples. The median cytokine levels between urine samples from health volunteers and those from influenza, S. pneumo or S. aureus were tested via Wilcoxon sum rank test. The p values were not adjusted for multiple comparisons. All analyses were performed using R-3.4.0 (https ://www.R-proje ct.org/).
16S rRNA-based PCR, ilumina library preparation, and data analysis. To  www.nature.com/scientificreports/ GGG AGG CAG CAG T-3′; Reverse: 5′-GGA CTA CCA GGG TAT CTA ATC CTG TT-3′) and an in-house standard to generate a standard curve. To assess bacterial community structure, primers specific for 16S rRNA V4-V5 region (Forward: 338F: 5′-GTG CCA GCMGCC GCG GTAA-3′ and Reverse: 806R: 5′-GGA CTA CHVGGG TWT CTAAT-3′) that contained Illumina 3′ adapter sequences, as well as a 12-bp barcode, were used. Sequences were generated by an Illumina MiSeq DNA platform at Argonne National Laboratory and analyzed by the program Quantitative Insights Into Microbial Ecology (QIIME) 50 . Operational Taxonomic Units (OTUs) were picked at 97% sequence identity using open reference OTU picking against the GreenGenes database. OTUs were quality filtered based on default parameters set in the open-reference OTU command in QIIME and sequences were rarified to an average sampling depth of 7,084 reads per sample. Representative sequences were aligned via PyNAST and taxonomy was assigned using the RDP Classifier. Processed data were then imported into Calypso 8.84 for further analysis and data visualization 51 . Alpha diversity was assessed using observed Shannon index and Eveness. Network analyses were generated with Speaman's correlations. Positive correlations were FDRadjusted at P < 0.05 and presented as network edges. OTUs generated in QIIME were finally analyzed using linear discriminant analysis (LDA) effect size (LEfSe) where non-parametric factorial Kruskal-Wallis sum-rank testing (P < 0.05) identified significantly abundant taxa followed by unpaired Wilcoxon rank-sum test to determine LDA scores > 2 was considered significant. Dendrograms of LEfSe display taxonomic distribution of significant taxa 52 53 . The samples were run with a spray voltage was 3 kV. The nitrogen sheath gas was set to a flow rate of 10 psi with a capillary temperature of 320 °C. AGC (acquisition gain control) target was set to 3e6. The samples were analyzed with a resolution of 140,000 and a scan window of 85-800 m/z for from 0 to 9 min and 110-1,000 m/z from 9 to 25 min. Solvent A consisted of 97:3 water:methanol, 10 mM tributylamine, and 15 mM acetic acid. Solvent B was methanol. The gradient from 0 to 5 min is 0% B, from 5 to 13 min is 20% B, from 13 to 15.5 min is 55% B, from 15.5 to 19 min is 95% B, and from 19 to 25 min is 0% B with a flow rate of 200 µl/min 53 . Files generated by Xcalibur (RAW) were converted to the open-source mzML format 54 via the open-source msconvert software as part of the ProteoWizard package 55 . Maven (mzroll) software, Princeton University 56, 57 was used to automatically correct the total ion chromatograms based on the retention times for each sample. Metabolites were manually identified and integrated using known masses (± 5 ppm mass tolerance) and retention times (Δ ≤ 1.5 min). Unknown peaks were automatically selected via Maven's automated peak detection algorithms.
Multivariate statistical analysis for the MS/MS data was performed using XLSTAT software (Addinsoft, New York, NY) interfaced with Excel (Microsoft Corporation, Redmond, WA). The average coefficient of variation (C.V.) was 0.395 (± 0.211). Thus an inclusion criterion for technical replicates were applied based on C.V. ≤ 0.5 resulting in 11 exclusion (i.e. 11 technical replicates in duplicate and 29 in triplicate) resulting in C.V. average of 0.288 (± 0.114). To ensure that observations were directly comparable and to account for the biofluid concentration peak intensity was normalized to creatinine (these data were compared to unnormalized data to make sure there was no masking of biologically relevant changes by normalization). To evaluate the group trends, sample uniformity and identify potential outliers unsupervised multivariant principal component analysis (PCA) was employed. The variation were explained by F1 and F2 with a cumulative percent variability of 78.558 and a marginal increase to 85.958% with the addition of F3. These data were then independently k-means clustered followed by ascendant hierarchical clustering based on Euclidian distances. The data matrix's was rearranged according to the corresponding clustering with similarity proportional to a closer spatial relationship for patient sample columns and metabolite rows. 21 metabolites with less than 0.25 standard deviation were eliminated to simplify the graph. These clusters were also represented via a dendrogram displayed vertically for metabolites and another horizontally for patients. The data values of the permuted matrix were replaced by corresponding color intensities based on interquartile range with color scale of red to green through black resulting in a heat map. Patient identifiers and risk categories were replaced by color bars. XLSTAT expression analysis was then used to determine metabolite significance between groups using one-way ANOVA with Benjamini-Hochberg post hoc correction and found to have with significant differences using Tukey's honest significance test (Tukey HSD) for multiple comparisons. Partial least squares discriminant analysis (PLS-DA) was then applied to separate patient groups and identify metabolites with corresponding variable importance in the projection (VIP) values above 1. The four component PLS-DA was then rerun with the variables centered and reduced prior to analysis to improve the model quality and identify the corresponding VIP. Each identified metabolite was then analyzed individually using ANOVA then the means of pneumonia samples were then compared to healthy controls (Dunnett's multiple comparisons test) prior to and post outlier removal performed with PRISM software (Graphpad, San Diego, CA www.nature.com/scientificreports/ decompositions were normalized in a way that sum of all detected OTUs are equal to 1 for each subject. We then combined identified OTUs, cytokines and metabolites as predictors of four subject categories. Our first approach is to implement multi-class classification using various machine learning algorithms such as random forest, ensemble trees, support vector machines, k-nearest neighborhood. However, small sample size (n = 40), multiple output categories (four subject groups) and expected larger number of predictors (OTUs and metabolites) made predictive modeling very challenging. Therefore, as an alternative approach, we implemented recursive binary classification in three steps and obtained three different models at each step. First, Model 1 is to distinguish healthy subjects from the other three disease categories, Model 2 to distinguish between bacterial (S. Aureus and S. Pneumoniae) disease from Influenza, and Model 3 to distinguish between S. Aureus and S. Pneumoniae. Considering the small sample size and large number of predictors, for each of the three models, we first implement Least Absolute Shrinkage and Selection Operator (LASSO) 58 logistic regression 59 . LASSO is statistical method retraining strong features of both subset selection and ridge regression. It implements ordinary least squares subject to sum of absolute values of the regression coefficients being less than a predetermined constant value. Logistic regression LASSO is an extension of LASSO for an output variable with binomial distribution. By implementing LASSO, some of the regression coefficients are shrink to take a valued of zero therefore only variables with non-zero regression coefficients remain in the model. By taking advantage of LASSO, we will first identify OTUs and metabolites that are the most effective in distinguishing between our subject categories in Model 1, Model 2, and Model 3. Next, we combined all selected OTUs and metabolites in readdress multi-class classification problem using the machine learning approaches mentioned.
For each predictive model, we implemented a stratified 5-folds cross validation process to avoid overfitting. To ensure unbiasedness of our cross-validation strategy, we split the entire cohort into five distinct each including a same number of subjects from each category. Next, a model built on using four folds of data and tested on the remaining fold. By repeating this process five times by changing the test fold, we obtain a predict class labels for each subject using a model that is trained on other subjects. We did not implement cross-validation with the goal of hyper-parameter tuning and optimization, instead, we used cross-validation (1) to evaluate the variability (or stability) of the predictive models from one subset to another (2) to evaluate the model performance on an unseen dataset. Therefore, we did not transfer parameter setting from one fold to another, instead, we fixed model parameters across each folds. We compared different machine learning algorithms based on model performance statistics such as specificity, sensitivity, and positive and negative predictive values. ethics statement. The usage of human samples was approved and performed in accordance with the regulations and guidelines set by the Univeristy of Louisville Insitutional Review Boards and the Human Subjects Protection Program. Samples were obtained from the University of Louisville Infectious Disease Biorepository (IRB # 13.0001) and de-identified metadata were used for analysis under the Biomarkers study (IRB # 17.0601). All patients provided written informed consent for sample biorepository storage and subsequent use in research studies.

Data availability
The datasets generated during microbiome analysis in this study are available as raw sequence reads at NCBI Sequence Read Archive: https ://submi t.ncbi.nlm.nih.gov/subs/sra/SUB76 20442 /overv iew. Datasets generated during metabolite analysis in this study are available in the Metabolights Database at https ://www.ebi.ac.uk/ metab oligh ts/MTBLS 1722. www.nature.com/scientificreports/