A distinct glycerophospholipid metabolism signature of acute graft versus host disease with predictive value

Graft versus host disease (GvHD), caused by alloreactive donor T cells, is a major factor that limits the success of allogeneic hematopoietic cell transplantation (alloHSCT) (1). An estimated 40%–60% of alloHSCT recipients are affected by GvHD (2, 3). Currently, the diagnosis of acute GvHD (aGvHD) is mainly based on clinical symptoms pertaining to 1 or more target organ systems (4, 5). Glucocorticoids and immunosuppressive agents are used as the first-line treatment for GvHD. However, generalized immunosuppression is liable to lead to opportunistic infections. Therefore, there has been an increased interest in identifying biomarkers of aGvHD and understanding the pathogenesis of aGvHD so as to guide treatment decision-making and direct more intensive clinical surveillance of high-risk patients (6–14).

Extensive research has been dedicated to identifying blood biomarkers for aGvHD to improve diagnosis and facilitate personalized treatment. Biomarkers for aGvHD have been studied thoroughly at the facets of cytokines (15–17), immune cells (18–21), single nucleotide polymorphisms (22), miRNAs (23–27), and chemokines (12, 28–31). However, no single biomarker predicts outcomes across different transplant practices or patient characteristics. Up to now, the best combination was ST2 and REG3α. Studies have shown that ST2 and REG3α levels at the initiation of GvHD therapy and during the first month after transplantation may help improve risk stratification for treatment-resistant GvHD (11, 12, 30). An increasing body of evidence suggests a strong connection between aGvHD and metabolism (32, 33). Cell metabolism determines the fate and function of T cells, which play a key role in the pathogenesis of aGvHD. Metabolomics is an important component of system biology that allows us to measure metabolic alterations that reflect the genetic, epigenetic, and environmental determinants of cellular physiology (34, 35). The stearic acid/palmitic acid (SA/PA) ratio detected on day 7 after HSCT was shown to be a useful biomarker of aGvHD and relapse (36). However, the influence of gene changes on metabolism and the gene-metabolite networks involved in aGvHD are not well characterized. Integration of metabolomics and transcriptomics may yield further insights into the pathogenesis of aGvHD. Indeed, it may help uncover the complex regulatory networks involving genes and metabolic pathways in aGvHD than either approach alone.

In the present study, we performed global metabolic analysis in alloHSCT recipients with or without aGvHD to identify critical metabolites that may help predict aGvHD. We performed transcriptomic profiling in a cohort of patients (training set) to identify significantly altered genes between patients with and without aGvHD. Integrative analysis of metabolomics and transcriptomic data revealed that glycerophospholipid (GPL) metabolism may play an important role in the pathogenesis of aGvHD. These relevant metabolic and transcriptomic alterations were further verified in a separate cohort of patients (validation set). Subsequently, we combined the metabolic and transcriptomic data to identify significantly disturbed pathways at both the metabolic and transcriptional levels. Lastly, we developed a model including 5 highly connected GPL metabolites to predict the probability of aGvHD.

Metabolic alterations in aGvHD patients. The overall study design is shown in Figure 1, and the clinical characteristics of patients are summarized in Table 1. There were no significant differences with respect to age, sex, diagnosis, disease status, conditioning regimen, or type of transplantation between the 2 sample sets or between patients with and without aGvHD in each set. To explore the metabolic profile of aGvHD, we used a metabolomic approach based on liquid chromatography−mass spectrometry (LC-MS) to identify differentially expressed metabolites in the aGvHD and control groups (without aGvHD) in the training set. The metabolomic profile of the aGvHD group at +15 days after transplantation was distinctly different from that of the control group, as evidenced by robust orthogonal partial least squares discriminant analysis (OPLS-DA) models established in the training (R²X_cum = 0.523; R²Y_cum = 0.745; and Q²_cum = 0.634) and validation sets (R²X_cum = 0.519; R²Y_cum = 0.684; and Q²_cum = 0.576). The scores plot demonstrated a clear separation between cases with and without aGvHD, and no overfitting of the OPLS-DA models were found (Supplemental Figure 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.129494DS1).

Figure 1

Experimental flow chart.

Table 1

Characteristics of patients

Based on the criteria of variables important to projection (VIP) value > 1 and P < 0.05 (2-tailed Student’s t test), 38 metabolites were identified as potential biomarkers (Supplemental Table 1). These metabolites pertained to GPL metabolism, fatty acid metabolism, amino acid metabolism, citrate cycle, sphingolipid metabolism, purine metabolism, and bile acid metabolism.

Integration of metabolomics and transcriptomics revealed altered GPL metabolism pathway during pathogenesis of aGvHD. To discover the essential metabolites, data pertaining to 38 potential biomarkers in the training and validation sets were used to build correlation networks (Figure 2A). Five highly connected metabolites (lysophosphatidylcholine [LysoPC] [18:1], C16:0 PAF, LysoPC [P-16:0], LysoPC [18:2], and LysoPC [22:6]) with the brightest colors appeared in the center of the networks in both sets. Interestingly, these metabolites were all GPLs, which indicated that GPL metabolism may be notably changed in aGvHD. Figure 2B illustrates the altered pathways based on metabolic profiling data. We identified GPL metabolism as the top hit because of its lower P value and greater pathway impact; this suggested that GPL metabolism is one of the most significant pathways.

Figure 2

Integration of metabolomics and transcriptomics identified a unique set of GPL metabolism pathways in aGvHD. (A) Correlation analysis using profiling data from 38 metabolites in the training (n = 57) and validation sets (n = 50). (B) Overview of the pathway analysis based on metabolites alteration in +15 day samples after transplantation. (C) Pathway enrichment analysis of transcriptomic profiles associated with aGvHD.

To trace the upstream variations of metabolome, we conducted transcriptional profiling of the same specimens. Of the 57,773 genes profiled, 1,148 were significantly altered in patients with aGvHD (t test, FDR corrected P < 0.05, |fold change| > 1.5). Figure 2C shows the top 30 significant pathways based on transcriptomic pathway analysis. After metabolic pathway analysis and transcriptomic pathway analysis, we found 2 differentially expressed pathways at both the metabolomic and mRNA expression levels in patients with aGvHD: GPL metabolism and inositol phosphate metabolism (Figure 2). Furthermore, transcriptomic pathway analysis (Figure 2C) showed that the GPL metabolism pathway was one of the most significant pathways with Bonferroni correction P value (Q value) < 0.001; this indicated that it is a critical pathway related to the development of aGvHD. We identified 13 significantly altered genes associated with GPL metabolism. Moreover, gene-related pathways known to play a key role in aGvHD physiopathogenesis were also differentially expressed between patients with and without aGvHD; these included Th1 and Th2 differentiation, Th17 differentiation, and MAPK signaling pathways.

Glycerophospholipolytic metabolites and enzymes were significantly altered during aGvHD. Next, we investigated the significant metabolic alterations in GPL metabolism. Compared with patients without aGvHD, the following 17 metabolites involved in the GPL metabolism were identified in this study: LysoPC (16:0), LysoPC (20:3), LysoPC (17:0), LysoPC (18:2), LysoPC (22:6), LysoPC (18:1), LysoPC (P-16:0), LysoPE (0:0/20:1), LysoPE (0:0/22:4), PC (16:0/22:6), PC (13:0/0:0), PC (14:0/18:1), PC (15:0/20:5), Phosphatidylethanolamine (PE) (20:4/20:3), PE (16:0/18:2), Diacylglycerol (21:0/22:4/0:0) (DG [21:0/22:4/0:0]), and platelet activating factor C16:0 (C16:0 PAF) (Supplemental Table 1). All of these were differently expressed in aGvHD plasma in both the training and validation sets (P < 0.05). Compared with patients without aGvHD, those with aGvHD exhibited significant accumulation of PCs and PEs, along with reduction of LysoPCs, LysoPEs, DG, and C16:0 PAF. Of note, the levels of metabolites were completely or partially restored to normal levels with the remission of aGvHD in both the training and validation sets (Supplemental Table 2).

We further explored the gene expressions of these lipases involved in GPL metabolism, which included GPL lipase surrogate genes lecithin-cholesterol acyltransferase (LCAT, cholesterol esterifying enzyme), lysophosphatidylcholine acyltransferase 1 (LPCAT1, LysoPC acyltransferase 1), LPCAT4 (LysoPC acyltransferase 4), HRAS-like suppressor 3 (PLA2G16, phospholipase), PLCG1 (phosphor lipase C), phosphatidate phosphatase LPIN 1 (LPIN1, phosphatidic acid phosphohydrolase enzyme), DGKE (DG kinase), and CHKA (choline kinase) (Figure 3). The fluctuations in the relative signal intensities were in accordance with the aGvHD severity in both the training and validation sets (Figure 3, B and C). The relative gene expressions of these lipases in patients with and without aGvHD were also visualized using a heatmap (Figure 3A). As shown in Figure 4A, we observed a consistent trend in gene expressions of these lipases and the related GPLs. Downregulation of lysoPCs, PAF, and DG were associated with downregulation of LCAT, PLA2G16, PLCG1, and LPIN1; in addition, upregulation of PCs and PEs were associated with upregulation of LPAT1 and LPCAT4.

Figure 3

Gene expression of glycerophospholipolytic enzymes in 2 independent sets. (A–C) Heatmap showing the changes in the expression of genes involved in GPL metabolism of aGvHD and control (without aGvHD) groups. Shades of red and green represent high or low expression (see color scale) (A). Each column represents a patient with aGvHD or without aGvHD. LCAT, LPCAT1, LPCAT4, PLA2G16, PLCG1, LPIN1, DGKE, and CHKA are differentially expressed in aGvHD as compared with control in training (B) and validation sets (C). In the training set, n = 13 for aGvHD and aGvHD^Δ samples, n = 12 for without aGvHD samples. In validation set, n = 12 for aGvHD, aGvHD^Δ, and without aGvHD samples. The comparisons were made between control vs. aGvHD and between aGvHD^Δ vs. aGvHD. Box plots represent gene expression level with relative intensity (log₂) of microarray data in the training set or relative Ct value normalized with endogenous control gene U6 using qPCR data in validation set. Bars indicate median value. P value refers to 1-way ANOVA. *P < 0.05, **P < 0.01. aGvHD^Δ, aGvHD in remission after medication.

Figure 4

Schematic overview of the metabolites and major metabolic pathways, as well as pathway-related enzyme surrogates iterations, in aGvHD patients at day +15 after alloHSCT. (A) Description of the glycerophospholipid metabolism pathway related metabolites and enzyme surrogate genes changes in aGvHD. (B) Dysregulation of metabolic pathways in aGvHD. The enzyme surrogates are shown in color: red text represents increased expression of enzyme surrogates, and blue text represents decreased expression of enzyme surrogates. The yellow circles represent no detected metabolites; green circles represent decreased metabolites; and red circles represent increased metabolites.

Development of a predictive model of GPL metabolite biomarkers and its association with aGvHD. An ideal biomarker should be detectable before the onset of aGvHD. We analyzed the metabolites in plasma samples collected from aGvHD patients before the onset of aGvHD (median 15 days [range: 4–75] before aGvHD diagnosis). A distinct GPL metabolism signature was uncovered in aGvHD plasma, as evidenced by significant alteration of 17 of 38 identified plasma metabolites and their corresponding lipases in this pathway. We hypothesized that dysregulated GPL metabolism is an early event in the activation of alloreactive T cells during aGvHD pathogenesis. We further used least absolute shrinkage and selection operator (LASSO) logistic regression model to build aGvHD risk score (GRS) in the training set and identified 5 metabolites from the 17 GPL metabolites (Figure 5 and Supplemental Figure 2). We evaluated the predictive value of these 5 metabolites for the development of aGvHD (Supplemental Figure 3, A and B). Hierarchical clustering analysis revealed significantly altered metabolite concentrations (Figure 5A). Subsequently, we calculated GRS for each patient based on their individual expression levels of the 5 metabolites, where GRS = 2.2974 × v1 + 20.8171 × v2 – 102.8055 × v3 – 0.0939 × v4 – 1.1977 × v5 (Supplemental Table 3). The receiver operating characteristic (ROC) curves based on the GRS yielded satisfactory results in the training and validation sets, respectively (Figure 5B). The metabolism signature had an AUC value of 0.922 (95% CI, 0.849–0.996) in the training set and 0.896 (95% CI, 0.811–0.981) in the validation set. The optimal cut-off point of GRS in training set (–4.74) was associated with sensitivity of 88.9% and specificity of 86.7%. Using this cut-off point, patients were divided into high and low GRS groups. Patients with high GRS had OR of 32.8 for aGvHD compared with low GRS. The 5-biomarker model successfully identified HLA-matched related or unrelated alloHSCT recipients with a high risk for aGvHD (83.02% in the high-risk group and 12.96% in the low-risk group) (P < 0.0001). The algorithm performed equally well in the training set (85.71% vs. 10.34%, P < 0.0001) and the validation set (80.00% vs. 16.00%, P = 0.00002). There were no significant differences between these 2 groups with respect to clinical characteristics (Supplemental Table 4).

Figure 5

Five metabolites exhibited a strong correlation with the development and severity of aGvHD. (A) The heatmap of 5 metabolites (C16:0 PAF, LysoPC [22:6], PE [20:4/20:3], PC [13:0/0:0], LysoPE [0:0/22:4]) selected from 17 GPL metabolites. (B) ROC curve analysis of the ability of plasma metabolites, including C16:0 PAF, LysoPC (22:6), PE (20:4/20:3), PC (13:0/0:0), and LysoPE (0:0/22:4), to predict the development of aGvHD based on both the training (n = 57) and validation sets (n = 50). The area under the ROC curve was 0.922 (95% CI, 0.849–0.966) in the training set, and 0.896 (95% CI, 0.811–0.981) in the validation set. (C) Box plots showing fluctuations in relative signal intensities of 5 metabolites. Five metabolites within the panel were restored to normal levels at varying degrees after aGvHD in remission after medication. The comparisons were made between control vs. aGvHD, and aGvHD^Δ vs. aGvHD, respectively. Bars indicate median value. P value refers to 1-way ANOVA. *P < 0.05, **P < 0.01. aGvHD^Δ, aGvHD in remission after medication.

Further, we observed that the levels of 38 metabolites were largely restored to normal with amelioration of aGvHD clinical scores in both the training and validation sets (Supplemental Table 2 and Figure 5C). In particular, the levels of C16:0 PAF, LysoPC (22:6), PC (13:0/0:0), and LysoPE (0:0/22:4) in the panel were completely restored to normal after aGvHD remission. The other metabolites were also restored to varying degrees in both sets.

Subsequently, we investigated whether the metabolites exhibit aGvHD target organ specificity. LysoPE (0:0/22:4) and PC (13:0/0:0) showed a positive correlation with liver and skin, respectively. PC (13:0/0:0) showed a negative correlation with gastrointestinal aGvHD (Supplemental Figure 3, C–E).

To develop a clinically applicable method that could help predict an individual’s risk of aGvHD, we developed a nomogram based on the panel of these 5 metabolites of all patients in the 2 sets (Figure 6A). Calibration plots of the nomograms for predicting the probability of aGvHD showed that the nomograms did well as compared with an ideal model. The results indicated the reliability of the model based on GRS (Figure 6B).

Figure 6

Nomogram based on training set for predicting risk/probability of aGvHD. (A) Nomogram based on training set (n = 57) for predicting the risk/probability of aGvHD. (B) The calibration plot comparing predicted outcomes with actual outcome. Dashed line is the 45° reference line, indicating an ideal nomogram. The dotted line is the apparent accuracy of the nomogram without correcting for overfit. The solid line is the bootstrap-corrected performance of the nomogram with a scatter estimate for future accuracy. Calibration plots showed that the nomograms did well compared with an ideal model.

Lastly, we assessed whether the GRS can be used as a prognostic tool in alloHSCT. One hundred seven alloHSCT recipients were divided into high GRS and low GRS groups. Kaplan-Meier survival analysis suggested that patients in the high GRS group had significantly poorer OS as compared with those in the low GRS group (P = 0.0106, Supplemental Figure 4). On multivariate Cox regression analysis with GRS and clinical characteristics, GRS was an independent unfavorable prognostic factor for OS (P = 0.014, hazard ratio = 2.696; Supplemental Table 5).

In this study, we employed metabolic analysis, along with transcriptome analysis, to identify critical biomarkers that distinguish alloHSCT recipients with aGvHD from alloHSCT recipients without aGvHD in 2 separate cohorts. Pathway analysis of 38 significantly altered metabolites and 1,148 differentially expressed gene surrogates revealed distinctly altered GPL metabolism. Subsequently, we developed a GRS based on 5 highly connected GPLs to diagnose aGvHD (AUC 0.922 in training set and 0.896 in validation set). In addition, high GRS was correlated with poor overall survival (OS). Moreover, GPL enzyme surrogate genes also exhibited an ability to predict the development of aGvHD (Supplemental Table 6). Collectively, our data suggest that the plasma GPL signature is an independent predictive factor for aGvHD. GRS on day 15 after alloHSCT may facilitate risk stratification of HLA-matched alloHSCT recipients for aGvHD.

Recent studies in the aGvHD model revealed that alloreactive T cells may employ fatty acid or glucose to fuel their activation (32, 33), which might be related to the altered GPL metabolism observed in our results. Lipids, such as GLPs, cholesterol, and glycolipids, are the most abundant molecular components of cell membranes. With the increased need of membrane building blocks, pathogenic T cells were shown to depend on de novo fatty acid synthesis and the underlying glycolytic-lipogenic metabolic pathway for their development (37, 38). Consistent with these findings, we observed increased GPL and fatty acid metabolism at +15 days after alloHSCT in aGvHD patients prior to the onset of clinical symptoms (Figure 4B), which indicated the proliferation of alloreactive T cells. Once activation, T cells need to transit from quiescence to fast expansion. In the meantime, lipid metabolism of T cells switched from energy generation through fatty acid oxidation to fatty acid biosynthesis for membranes and signaling molecules (38–40). A previous study also reported that GPL metabolism and fatty acid metabolism were accelerated in response to the need for rapid proliferation of T cells (39), which is in accordance with our results. During activation of T cells, beta-oxidation of fatty acids is decreased while other metabolic pathways increase. As shown in Supplemental Table 1, free fatty acids (the ingredient for synthesis of GPLs) were reduced in aGvHD plasma, which indicated increased consumption of fatty acids for GPL synthesis by activated cells. Interestingly, it was reported that lipid metabolism reversed to resting state conditions and proliferation of T cells was inhibited immediately after removing the stimulus in proliferating T cells (41), which demonstrated that lipid metabolism was closely related to T cell proliferation. Our results show that the acute graft versus host reaction is already in progress by day 15 and that the increased biomarker levels may precede clinical symptoms by days or weeks. The median day for aGvHD onset was 31 days in the training set and 29 days in the validation set (Table 1). The fluctuations in relative signal intensities were in accordance with the alterations of metabolites. Because of the scarcity of T cells in peripheral blood or BM on +15 days after transplantation, we were unable to conduct a metabolomic analysis for T cells in this study.

Active proliferation of T cells entails synthesis of several cellular constituents such as membranes. As a major structural lipid component of cellular membranes, enhanced GPL biosynthesis in activated T cells has been reported in the setting of inflammatory and autoimmune diseases (42, 43); these findings are consistent with our data (Figure 4 and Supplemental Table 1). Increased plasma levels of PEs and PCs in aGvHD plasma were observed in patients with aGvHD (Figure 4B), which may suggest increased proliferative activity in alloreactive T cells. LysoPCs and LysoPEs are intermediates in the biosynthesis pathway of PCs and PEs. Reduced LysoPCs and LysoPEs and increased activity of lipases LPCAT1 and LPCAT4 in aGvHD plasma suggested increased synthesis of PC and PE (Figure 4). Downregulated expressions of lipases LCAT, PLA2G16, PLCG1, and LPIN1 eventually led to reduced levels of LysoPCs, LysoPEs, PAF, and DG. Overall, GPL levels showed a positive association with gene expressions of these lipases in plasma samples of patients with aGvHD, which suggested that dysregulation of the GPL network appeared in the very early stage of aGvHD (Figure 4). These alterations of GPL metabolism directly reflect the host histocompatibility antigens to which donor T cells respond within days of graft infusion and were significant predictors of aGvHD in univariate and multivariate analysis.

In summary, based on metabolic and transcriptomic investigation combined with pathway analysis of differentially expressed metabolites and genes, we discovered distinctive altered GPL metabolism in the early stage of aGvHD pathogenesis. Gene expressions of GPL-related lipases were significantly altered in aGvHD samples, leading to dysregulated GPLs. The GRS generated from 5 GPL metabolites had the ability to predict the development of aGvHD with an AUC of 0.922 in the training set and 0.896 in the validation set. In addition, our data suggest that high GRS on day 15 after alloHSCT was a unique feature of patients who had a high risk for aGvHD before the overt onset of aGvHD (median = 15.0 days before diagnosis). There were several limitations in our study. First, we enrolled HLA-matched alloHSCT patients to maximize the homogeneity of the study. Thus, the GRS could only be applied in HLA-matched alloHSCT settings. For other kinds of alloHSCT, further identification and verification are needed. Second, small sample size was a major limitation in our study, which made it difficult to comprehensively analyze the prediction ability of GRS in populations from different backgrounds. However, the GRS showed powerful prediction ability for aGvHD with high AUC and discrimination ability with high OR. A prospective, multicenter clinical trial is needed in the future.

To the best of our knowledge, this is the first study that identifies GRS as a potential warning signal that may help improve risk stratification for aGvHD. By integrative “Omic” analysis, we have offered both a comprehensive view of the molecular perturbations underlying the pathogenesis of aGvHD and an initial investigation of a unique metabolic and transcriptomic network that may facilitate identification of aGvHD patients at a very early stage and help institute preemptive therapy. However, owing to the small sample size in our study, the results await further validation using a larger and independent cohort.

Patients and sample collection. This study was designed to assess whether plasma-derived metabolites, measured at the 15th day after transplant, could predict aGvHD. All patients (n = 128) who underwent HLA-matched alloHSCT at the Changhai Hospital between September 2014 and June 2018 were enrolled in this study. The inclusion criteria are presented in Supplemental Methods. Twenty-one patients relapsed within 6 months after alloHSCT. The study comprised a training set and a validation set. The median duration of follow-up was 33.43 months (range, 13.60–53.77), and the last date of follow-up was December 31, 2018. The median age at alloHSCT was 36 years (range, 17–66); 43.93% patients were female. In brief, 27 allogeneic transplanted patients with aGvHD and 30 patients without aGvHD (controls) were enrolled as the training cohort, while 24 allogeneic transplanted patients with aGvHD and 26 patients without aGvHD (controls) were enrolled as the validation cohort for the metabolomics study.

Plasma samples were collected at predetermined time points for metabolomic analysis using citrate as an anticoagulant and centrifuged at 1,590 g for 15 minutes at room temperature. The peripheral blood mononucleated cells and plasma samples were collected for transcriptomic analysis at 15 days after transplantation. aGvHD and amelioration of aGvHD were diagnosed according to the generally accepted criteria (44). All patients received prophylactic therapy with cyclosporine A, mycophenolate mofetil, and short-term methotrexate (45). Clinical information, including aGvHD grade, is shown in Table 1. All plasma samples and mononucleated cells were stored at –80°C and in liquid nitrogen before analysis.

Metabolic profiling with LC-MS. Metabolic profiles of all plasma samples were obtained using reversed-phase–MS (RPLC-MS) and hydrophilic interaction chromatography–MS (HILIC-MS) platforms as described elsewhere (46). Samples were randomized prior to LC-MS analysis to decrease experimental drifts. An in-house quality control (QC) mechanism was established by pooling and mixing equal amounts of plasma samples from all enrolled subjects; these were used to control intra- and interbatch variability. Details pertaining to sample preparation and LC-MS analysis methods are provided in the Supplemental Methods.

After nonlinear alignment of the data in the time domain, automatic integration, and extraction of the peak intensities, raw data containing retention time, intensity, and m/z pairs were obtained. Then, LC-MS data acquired from RP and HILIC columns were combined as a new data set and were used to perform metabolomics profiling analysis.

RNA isolation and transcriptome profiling. Total RNA was isolated from the peripheral blood mononucleated cells and then quantified using a Bioanalyzer (Agilent Technologies). Whole transcriptome sequencing libraries (3 pairs of aGvHD and controls, 6 libraries in total) were prepared according to the manufacturer’s instructions for the Whole Transcriptome Sample Prep Kit (Illumina) by GENEWIZ Company. After cutting off the adaptor sequence and deleting low-quality reads, all sequencing reads were mapped to the human genome using TopHat. The mapped results were processed using Cufflinks, producing a fragments per kilobase of transcript per million mapped reads (RPKM) data matrix and, with htseq-count, producing the raw hit counts data matrix. The raw hit counts gene matrices were transformed using DESeq variance stabilization algorithm and normalized using the quantile normalization function built into the Bioconductor Limma package. The normalized data matrix was then analyzed using the Limma and SAMR algorithms. The gene expression data has been deposited in the National Center for Biotechnology Information’s (NCBI) Gene Expression Omnibus (GEO) with accession number GSE131897.

GO-TermFinder was used to identify gene ontology (GO) terms that annotate a list of enriched genes with a 1.5-fold change cutoff and P ≤ 0.05. Gene function was then annotated on KAAS (KEGG Automatic Annotation Server). For pathway enrichment analysis, we used scripts in house to enrich significant differential expression gene in KEGG pathways.

Quantitative PCR. RNA from plasma samples was extracted using standard TRIzol (Invitrogen) protocol. cDNA was reverse transcribed from total RNA, and quantitative PCR (qPCR) was performed using standard TaqMan chemistry and a laser-equipped thermal cycler to detect changes in fluorescence in real time (Applied Biosystems). cDNA concentrations were calculated according to the ΔCt method, corrected for differences in PCR efficiency, and normalized to expression levels of U6. All primers for qPCR in the present study were purchased from Applied Biosystems (Supplemental Table 7).

Development of GRS. In order to construct a signature with selected GPL metabolite biomarkers for predicting aGvHD, we adopted penalized logistic regression model with the LASSO penalty to achieve shrinkage and variable selection simultaneously (47). The LASSO has been widely used for building prognostic signatures (48, 49), and the use of the LASSO logistic regression model allowed us to integrate multiple GPL metabolites into a single tool; this significantly improved the predictive accuracy as compared with accuracy using single metabolite alone. LASSO logistic regression model analysis was conducted using the “glmnet” package of the R software (version 3.4.3). The penalty parameter λ controls the amount of shrinkage. Usually, the optimal values of λ is the value that gives minimum mean cross-validated error, and a larger λ could be chosen for greater penalty at accepted cross-validated errors (50). Using the regression coefficients from LASSO logistic regression model based on the training set, we developed a GRS.

Data processing and statistical analysis. The metabolomic data were normalized using L-2-chlorophenylalanine as the internal standard. The multivariate statistical model of orthogonal partial least square discriminate analysis (OPLS-DA) was constructed with the software SIMCA-P+ (version 11.0, Umetrics). A detailed description of the identification of biomarker metabolites is provided in the Supplemental Methods. Metabolic pathway analysis was performed via MetaboAnalyst 3.5 to identify, analyze, and visualize the affected metabolic pathways. We investigated the predictive accuracy of GRS using ROC curve analysis; AUC was used to measure the accuracy. The optimal cut-off point of GRS (associated with maximal sensitivity plus specificity) was used to divide patients into 2 groups: low and high GRS. To provide the clinician with a quantitative method to predict the risk of aGvHD, we constructed a nomogram, and the performance of the nomogram was explored graphically by calibration plots. OS was calculated from the time of transplant until death from any cause. Estimated probabilities for OS were calculated by the Kaplan-Meier method, and the log-rank test was used to compare survival curves. Univariate and multivariate analysis were constructed using the Cox proportional hazard regression model. All statistical analyses were performed using the SPSS 20.0; “rms” package of R 3.4.3 was used to prepare nomogram. Two-tailed Student’s t tests were used for comparisons of 2 groups, and 1-way ANOVA was used for comparisons of multiple groups. Two-sided P values less than 0.05 were considered statistically significant. Fold changes were calculated as the average mass response (area) ratio between 2 groups.

Study approval. All procedures complied with the Helsinki Declaration standards and were approved by the Changhai Hospital IRB. The requirement for written informed consent was waived because this study used retrospective data from medical records and there were no interventions in the patients.