Study design

We report association modeling, biomarker discovery, biochemical enrichment analysis and topological network visualization of plasma metabolomic, fecal bacterial metagenomic and clinical data from 50 ME/CFS patients and 50 healthy controls. 562 plasma metabolites were assessed using targeted and untargeted mass spectrometry platforms. Figure 1 shows the pipeline for metabolomic data processing and the statistical methods used.

Study population and plasma collection

Subjects included 50 cases and 50 controls from the Chronic Fatigue Initiative (CFI) cohort41 recruited at four sites across the US who met the 1994 CDC Fukuda15 and/or Canadian consensus criteria for ME/CFS16. Controls were frequency-matched to cases on age, sex, race/ethnicity, geographic/clinical site and season of sampling42. All ME/CFS subjects (n = 50) completed standardized screening and assessment instruments including medical history and symptom rating scales, had a physical examination and provided blood samples. Controls (n = 50) of the CFI cohort study42 had been found to be free of: self-reported ME/CFS or ME/CFS symptoms or other conditions deemed by the recruiting physician to be non-representative of a healthy control population including substance abuse in the prior year and any history of self-reported psychiatric illness; antibiotics in the prior three months; immunomodulatory medications in the prior year; and clinically significant findings on physical exam or screening laboratory tests. All participants provided informed written consent in accordance with protocols approved by the Institutional Review Board at Columbia University Medical Center. All participants consented to phlebotomy, collection of stool samples and completion of clinical questionnaires.

Blood samples were collected into BD VacutainerTM Cell Preparation Tubes (CPT with sodium citrate anticoagulant) between June and October 2014, and centrifuged to pellet red blood cells. The plasma was shipped to Columbia University at 4 °C. After aliquoting, samples were stored at −80 °C until thawed for metabolomic analyses. All samples were analyzed within 2 years of collection.

Clinical assessments and medical history

Clinical symptoms and baseline health status were assessed on the day of physical examination and biological sample collection from both cases and control subjects using the following surveys: the Short Form 36 Health Survey (SF-36), the Multidimensional Fatigue Inventory (MFI), DePaul Symptom Questionnaire (DSQ) and Pittsburgh Sleep Quality Index (PSQI43). Table S5 describes the individual surveys with the specific questions listed. The SF-36 includes the following subject-reported evaluations about current health status: physical and social functioning, physical and emotional limitations, vitality, pain, general health perceptions and mental health change44. The MFI comprises a 20-item self-reported questionnaire focused on general, physical and mental fatigue, activity and motivation45. Cognitive function was tested based on the DSQ questionnaire data46 and was scored using a standard cognitive disturbance definition as well as a modified definition based on a subset of questionnaire variables. Sleeping disturbances linked to ME/CFS were tested and scored based on DSQ and PSQI questionnaire items. Each instrument was transformed into a 0–100 scale to facilitate combination and comparison, wherein a score of 100 is equivalent to maximum disability or severity and a score of zero is equivalent to no disability or disturbance.

IBS co-morbidity was based on self-reported diagnosis of IBS on the medical history form. IBS was diagnosed in 24 of the 50 ME/CFS patients (48%). One control subject of out 50 reported a diagnosis of IBS (2%). Three ME/CFS cases and 1 control were newly diagnosed with IBS.

Metabolomics

Three untargeted metabolomic assays and 1 targeted assay for 562 metabolites from over 20 biochemical pathways were performed with gas chromatography time-of-flight (GCTOF) and liquid chromatography–tandem mass spectrometry (LC-MS/MS) instruments by the West Coast Metabolomics Center at University of California, Davis, USA. Sample preparation and metabolomic analyses were described previously47,48,49. Mass spectrometry metabolomic data have been deposited in the Metabolomics Workbench (http://www.metabolomicsworkbench.org) with project identifier PR000576.

Statistical analyses

Mass spectrometry data were initially filtered for internal standard metabolites (quality control spike-in metabolites) and unannotated/unknown metabolites. Annotated metabolites were further filtered out if more than 50% of samples (50 individuals out of 100) showed non-detectable or missing profiles. For each metabolite, we created a dummy variable for whether the values were missing, and ran chi-squared tests between this dummy variable and the ME/CFS status variable. We found no significant difference in prevalence of missing metabolite values between ME/CFS and controls. For each of the remaining metabolites, missing values were imputed using 50% of its smallest available value. Normalization was performed by dividing each metabolite profile by the total sum of metabolite intensities per sample50. Base-10 log transformation was applied to limit outlier effect. Keeping features on a positive domain, all data points were multiplied by a factor of 1.0e + 09 before log transformation. Data for each metabolite were then scaled by the control standard deviation.

The nonparametric Mann-Whitney U test and adjusted univariate logistic regression were applied to identify potential predictors differentiating ME/CFS from controls. For the binary outcome of IBS subgroups vs. controls, logistic regression models were adjusted for all frequency matched variables and BMI. For the binary outcome of ME/CFS vs. controls, additional adjustment for IBS was included. The metabolites with a p-value below 0.05 in both U test and adjusted univariate logistic regression were then used to develop a representative set of variables by random forests17. The top 10 random forests ranked metabolites were fitted as predictors in the predictive multivariate logistic regression models51,52. Table S6 shows the descriptive value of chemical compounds differed in ME/CFS, subgroups and ME/CFS female subjects only. 5-MT was excluded from random forests analysis as it is confounded by the use of antidepressants in 50% of cases. In-sample receiver operating characteristic (ROC) curves were plotted and area under the curve (AUC) was measured to assess the goodness-of-fit of the models. To examine the predictive performance of the multivariate logistic regression models, random resampling cross-validation was performed with 1000 iterations. Data were randomly split into a training set (80%) and a test set (20%) within each iteration.

Correlations between bacteria, metabolites and disease scores were examined using Spearman correlations.

Data were analyzed and visualized with Matlab (R2013a, The Mathworks Inc., MA), R (version 3.3.3) and SPSS (version 24, IBM, NY). All p-values were 2-tailed.

MetaMapp network mapping and enrichment analysis

MetaMapp networks were created using the method described in the ref.18. Cytoscape software version 3.4.0 was used to visualize the networks with overlaid statistical results. Networks were visualized using the organic layout algorithm. Node size was mapped to fold-change differences in the case vs. control groups. Node color was mapped to the direction of the between-group differences. Structurally identified metabolites were grouped into 21 chemical groups. Metabolites were manually classified because many compounds, including complex lipids, are poorly covered in pathway maps from major databases such as KEGG, Reactome or MetCyc. Metabolite groups were used as input for an enrichment analysis using Fisher exact test in R. P-values for the enrichment were corrected for multiple hypothesis testing using the Benjamini-Hochberg false discovery rate (FDR) method53 controlling the FDR at 0.05 level.

Topological data analyses

Metagenomic and metabolomic data including plasma metabolite levels, bacterial composition and inferred metabolic pathways, plasma immune profiles and health symptom severity scores were integrated for TDA using the AYASDI platform (Ayasdi, Menlo Park, California). AYASDI represents high-dimensional, complex biological data sets as a structured 3-dimensional network54. Each node in the network comprises 1 or more subject(s) who share variables in multiple dimensions. Lines connect network nodes that contain shared data points. Unlike traditional network models wherein each node reflects only a single sample, the size of a node in the topological network was proportional to the number of variables with a similar profile. We built a network comprised of 100 samples and 1900 variables (562 variables representing plasma metabolites, 574 representing fecal bacterial relative abundances at different taxonomic levels, 586 variables representing predicted bacterial metabolic pathways, 61 variables reflecting immune molecules, 81 variables representing different ME/CFS fatigue and other symptom score/health questionnaire items and information on co-morbidities; and 36 demographic variables). All variables were weighted equally. Jaccard and variance-normalized Euclidean distance methods were used as the distance metrics; a range of filter lenses (neighborhood lens 1 and 2, MDS coordinate 1 and 2, Metric PCA 1 and 2) was used to identify networks. A metric represents a notion of similarity (or distance) between rows in the data. The Jaccard metric measures the dissimilarity of asymmetric information on non-binary variables. Each row must be a collection of objects (no order necessary) and this metric computes the Jaccard score (1 - intersection/union) between the two sets. The Jaccard metric treats the columns as delimiters. The rows are treated as a list of objects in the set. The variance-normalized Euclidean metric is a variant on the Euclidean that accounts for the scale choices. For each variable, this metric finds its mean and standard deviation, and rescales the value of the coordinate by subtracting the coordinate mean and dividing by the corresponding standard deviation. A lens is a filter that converts the dataset into a vector. The neighborhood lenses generate an embedding of high-dimensional data into two dimensions by embedding a k-nearest neighbors graph of the data. A k-nearest neighbors graph is generated by connecting each point to its nearest neighbors. The MDS coordinate 1/2 and Metric PCA 1/2 lenses compute a variant of PCA coordinate lenses for data that does not use the Euclidean metric. The AYASDI maps the data into a Euclidean space using the rows of the distance matrix as the coordinates and then performs PCA.

Standard statistical methods were applied to define the primary variables of these networks.

Data and materials availability

Data has been uploaded to the Metabolomics Workbench (http://www.metabolomicsworkbench.org) with project identifier PR000576.