Metabolite-based machine learning model predicts success of chemoimmunotherapy for advanced lung squamous cell cancer

We retrospectively screened 3124 patients treated at Shanghai Chest Hospital from July 2020 to July 2021, and finally included 79 patients who met the criteria and ended follow-up on April 31, 2024. All patients provided written informed consent, which was approved by the Ethics Committee and Institutional Review Board of Shanghai Chest Hospital (Reference number: LS1808). The inclusion criteria were: (1) histologically or cytologically confirmed squamous cell carcinoma of the lung; (2) stage IIIB to IV according to the ninth version of TNM, which was not suitable for radical surgery; (3) having measurable lesions; and (4) receiving first-line PD-1 (programmed cell death 1) inhibitor combination chemotherapy. (5) Eastern Cooperative Oncology Group (ECOG) performance status (PS) score between 0 and 1. Exclusion criteria were: (1) presence of driver gene mutations; (2) surgical treatment after immunotherapy; (3) active infection, such as HIV, hepatitis B, hepatitis C, etc.; (4) incomplete clinical data or did not complete the necessary systemic examination; (5) patients with other primary active malignancies; (6) patients with a history of autoimmune diseases, severe cardiopulmonary dysfunction or other serious complications.

We collected clinical information on patients' age, sex, smoking history, TNM stage, PD-L1 expression, ECOG PS, etc. via the hospital 's electronic system of medical records. Patients were treated with PD-1 inhibitors combined with chemotherapeutic agents administered every 3-4 weeks as a cycle until disease progression or serious adverse reactions or death. Immunologic medications included pembrolizumab, tislelizumab, cariselizumab, and sintilimab at a single dose of 200 mg. According to the patient 's body surface area and tolerance, the specific chemotherapy drug administration regimen was generally platinum-based doublet, platinum drugs included carboplatin, cisplatin, lobaplatin or nedaplatin, drugs used in combination with platinum drugs included paclitaxel drugs, gemcitabine or docetaxel; for patients not suitable for platinum drugs, gemcitabine combined with vinorelbine or gemcitabine combined with docetaxel is given. Patients' condition will be assessed before each cycle by chest computed tomography (CT), abdominal ultrasound, bone scan, brain enhanced magnetic resonance imaging (MRI) or positron emission tomography-computed tomography (PET-CT), and progression will be judged by at least one professional radiologist and clinical medicine.

Our serum samples were obtained with prior approval from the Shanghai Chest Hospital Clinical Biobank and could be investigated with these samples. Serum samples were previously stored in a -80°C freezer, and we obtained serum samples and performed untargeted metabolomic analysis within 4 weeks before patients received immune combination therapy by consulting patients' blood collection time records.

Serum samples stored in a -80°C freezer were preprocessed for analysis, first thawed on ice and vortexed, and 50 μL of sample was mixed with 300 μL of extraction solution (acetonitrile: methanol = 1:4, volume ratio) in a 2 mL microfuge tube. Vortex again, centrifuge, collect 200 μL of supernatant, place the collected supernatant at -20°C for 30 mins, then centrifuge to obtain 180 μL of supernatant for liquid chromatography-mass spectrometry (LC-MS) analysis (19).

Pretreated serum samples were analyzed by LC-MS in the next step, where the specific setting of the chromatographic part was the column Waters ACQUITY UPLC BEH C18, 1.8 µm particle size, 2.1 mm inner diameter x 100 mm length; the column temperature was 40°C, the flow rate was 0.4 mL/min, and the volume of sample injected into the column was 2 μL. Gradient elution was performed in a solvent system containing 0.1% formic acid in water (mobile phase A) and 0.1% formic acid in acetonitrile (mobile phase B) in steps. That is, 5% mobile phase B started at 0 min; within 11 min, a linear gradient to 90% mobile phase B; 90% mobile phase B was held for 1 min; return to 5% mobile phase B within 0.1 min and hold for 1.9 min to quickly return to starting conditions and prepare for the next injection.

The mass spectrometer was TripleTOF 6600, the data acquisition software was Analyst TF 1.7.1 (Sciex, Concord, ON, Canada), the mass spectrometry acquisition mode was information-dependent acquisition (IDA) mode, and appropriate parameters were set to ensure high sensitivity and resolution of LC-MS analysis. In addition, we prepared a pooled human serum sample from all study participants as a consistent quality control (QC) sample. Then, a blank solvent was used as a negative control. These QC samples were evenly distributed throughout the analysis batches. Specifically, they were inserted every 10 samples to monitor the stability and consistency of the analytical process. For the pooled human serum, we monitored the relative standard deviation (RSD) of the peak areas for major metabolites. An RSD threshold of 30% was set, and any metabolite exceeding this threshold was excluded from further analysis. For the blank solvent, we checked for any carryover effects by ensuring no significant peaks were observed that could indicate contamination from previous samples. It can be used to assess the quality of the data and also to monitor within-run and between-run variability during the analysis. Sample handling and repeatability of the analytical method can be assessed by comparing the results of QC samples analyzed at different time points.

Raw data files generated by the mass spectrometer such as mzML or raw format were converted to mzXML format using ProteoWizard software. Peak extraction was performed using XCMS software and extracted peaks were aligned to correct for minor differences in retention times due to variations in experimental conditions. Peak areas were corrected using the "Support Vector Regression" (SVR) method, and peaks with greater than 50% missing in each group of samples were filtered. Metabolite identification information is then obtained through various libraries such as self-built libraries and public databases as well as using metDNA methods.

The primary endpoint of this study is overall survival (OS), defined as the time from the start of the first treatment until the patients' death from any cause. The secondary endpoint was progression-free survival (PFS), defined as the time from the start of the first treatment to the patients' tumor progression (in any way) or death (in any way). We noticed that several papers mentioned 24-month survival rates and used this as a basis to assess the efficacy of lung cancer immunotherapy (20, 21), so combining the actual survival of 79 patients, we divided patients into Response (R = 41) and Non-Response (NR = 38) groups using 24-month cutoff values. Serum was collected from patients within four weeks prior to chemoimmunotherapy for untargeted metabolomics analysis. Combined with variable impact projection (VIP) and P values, differential metabolites were initially selected, and then the differential metabolites finally included in the model were further selected by LASSO and RF methods (22). Patients were randomly divided into training and validation sets at a ratio of 7:3, modeled using multiple machine learning methods such as RF, SVM, and LR, and then receiver operating characteristic (ROC) curves were plotted to assess model performance. The best model tested by different machine learning methods was selected, then the risk score was calculated, and the total population was divided into low-risk and high-risk groups using the median as the cutoff. Comparison of prognostic differences between high-risk and low-risk groups further assessed model performance.

Baseline characteristics of patients were compared using Chi-square test or Fisher 's exact test. OS and PFS were calculated using Kaplan-Meier method and log-rank test. Univariate and multivariate Cox regression were used to assess the impact of different factors on OS and PFS. ROC curves were utilized to evaluate the predictive capability of each factor. Specifically, the model's classification ability was demonstrated by plotting the relationship between the true positive rate (sensitivity) and the false positive rate (1-specificity). The area under the curve (AUC) was subsequently calculated, with a higher AUC value indicating a stronger predictive ability of the model. Specifically, we followed the procedure outlined below: we generated 1000 bootstrap samples from our dataset. For each bootstrap sample, we calculated the AUC. We sorted the AUC values obtained from the bootstrap samples and determined the 2.5^th and 97.5^th percentiles to obtain the 95% confidence interval (CI). Let AUC denote the area under the ROC curve, and let CI_lower and CI_upper denote the lower and upper bounds of the 95% CI, respectively. The results can be expressed as: 95% CI=[CI_lower,CI_upper].

Orthogonal partial least squares discriminant analysis (OPLS-DA) was used to demonstrate differences between groups, which was accomplished using the MetaboAnalystR package in R software. The volcano plot primarily serves to illustrate the relative content disparities of metabolites between two sets of samples, alongside the statistical significance of these differences. Student's t-test was employed to analyze variables with a normal distribution and equal variance between the two groups, whereas the Mann-Whitney test was utilized for those with a non-normal distribution. Significant differences in selecting differential metabolites were determined by VIP > 1 and p-value < 0.05. The annotation of metabolic pathways, encompassing the differential metabolites, was executed with the aid of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. LASSO regression was conducted using the glmnet, foreign, and tidyr packages in R, while the RF method was executed via the varSelRF package in R. The R package used by SVM is mainly e1071 package, and the R package used by LR is mainly glm package.

All data analyses were conducted utilizing R version 4.4.1 and SPSS version 26.0, and Adobe Illustrator 2022 was used for picture drawing integration. A p-value less than 0.05 in a two-tailed test was deemed statistically significant.

The overall flow of our study was shown in Figure 1. There were 79 patients in our study who received combination therapy as their first line of treatment. Their baseline clinical characteristics such as age, sex, smoking history, ECOG PS score, TNM stage, number of metastatic organs and PL-L1 expression were shown in Table 1. As can be seen from the Table 1, our selected population was concentrated at ≥ 65 years, male, with a history of smoking, ECOG PS 0, TNM stage IV, and PD-L1 expression ≥ 1%. In Table 2, after we divided 79 patients into NR (OS < 24 months) and R (OS ≥ 24 months) groups according to OS, we then compared the clinical characteristics of the two groups and performed a chi-square test and found that these baseline characteristics were evenly distributed in both groups without statistical difference (P > 0.05).

In order to initially screen for differential metabolites between the two groups, we first verified the rationality of grouping, and the OPLS-DA model showed that the two groups could be significantly separated (Figure 2A), demonstrating that patients with different survival periods had different metabolic profiles. Based on the VIP obtained from the OPLS-DA model (biological replicates ≥ 3), combined with the P-value of univariate analysis (biological replicates ≥ 2), we could initially screen 117 differential metabolites between the two groups (VIP > 1 and P < 0.05). Among them, as shown in volcano plot (Figure 2B), 70 differential metabolites were up-regulated and 47 differential metabolites were down-regulated. In addition to performing statistics on VIP values and P-values for differential metabolites, we also calculated FC values for differential metabolites and plotted radar plots for the top 10 metabolites with the largest difference, i.e., the largest absolute log₂FC value. In Figure 2C, these 10 metabolites were Moupinamide, 15(R)-17-Phenyl trinor prostaglandin F2alpha isopropyl ester, Guanabenz acetate, 5,6-dihydroxy-2-(4-hydroxyphenyl)-7-methoxy-8-[3,4,5-trihydroxy-6-(hydroxymethyl)oxan-2-yl]-4H-procen-4-one, 5-Hydroxypseudobaptigenin, Lenticin, Dexmedetamide, Esprocarb, Asn-Tyr, and Asn-Val-Phe-Lys. In addition, we also compared metabolic pathways, and KEGG metabolic pathways were significantly differentially enriched in metabolites between NR and R groups. The closer the P-value is to 0, the more significant the enrichment (Figure 2D). The first five pathways with the smallest size from small to large were Chemical carcinogenesis-receptor activation, Breast cancer, Progesterone-mediated oocyte maturation, Oocyte meiosis, and Cortisol synthesis and secretion.

In our study, the Figure 3A revealed that the median OS was 13 months (95% CI, 10.0-16.0) in the NR group, which was significantly shorter than that in the R group (P < 0.0001). In addition, median PFS was 6.0 months (95% CI, 5.0-7.0) versus 19.0 months (95% CI, 12.7-25.3) in NR group and R group, respectively (P < 0.0001). Therefore, OS was also significantly shorter in the NR group than in the R group and the data was shown in Figure 3B. To identify which clinical factors influence patient survival, we performed univariate Cox regression analysis and found no features to be significantly associated with survival (Figure 3C).

Because the univariate Cox analysis was not significant, we did not include clinical characteristics in the model. Our preliminary statistical screening with p < 0.05, VIP > 1 resulted in a large number of metabolites, 117. To further screen more reliable biomarkers, we first selected 46 metabolites using LASSO (Figures 4A, B) and then identified metabolites using RF, and their cross-metabolites were candidate metabolites, and a total of 8 differential metabolites were identified as the best biomarkers. To avoid over-fitting and false positives, 500 feature selections and 10 cross-validations were performed. These 8 differential metabolites were Mevalonate, 1,2−Ethanedithiol, 4alpha−Hydroxymethyl−4beta−methyl−5alpha−cholesta−8,24−dien−3beta−ol, Arg−Asp−Leu−Tyr−Ser, Donhexocin, Glu−Cys−Ala, TG(10:0/14:0/a−15:0)[rac], D−3−Hydroxykynurenine, respectively, and their predictive performance was verified by AUC of ROC curves. As can be seen from Figure 4C, the predictive performance of individual differential metabolites is not sufficiently satisfactory (AUC < 0.75). However, the model formed by their combination had good predictive performance, with AUC values reaching 0.915, (95% CI, 0.857-0.972), suggesting that this model may be used as a predictive panel to predict the efficacy of immunotherapy combined with chemotherapy in patients with advanced lung squamous cell carcinoma (Figure 4D).

Subsequently, patients (n = 79) were randomly divided into training set (n = 55) and validation set (n = 24) in a 7:3 ratio, and various machine learning models were employed for datasets consisting of these eight metabolites, including RF model, SVM model, and LR model. The comparison of the relative content of these eight metabolites between the NR and R groups was generally consistent in the training and validation sets, but most of them were not significantly different (Figure 5), still suggesting that single metabolites were not able to be biomarkers. However, the RF, SVM, LR models constructed by the panel of eight metabolites showed good predictive performance, with AUCs of 0.973, 0.938, and 0.934 in the training set and 0.944, 0.897, and 0.900 in the validation set (Figure 6), with the RF model having the best predictive performance. The accuracy, precision, recall and F1 values of RF model in the training set were 0.98, and the accuracy, precision, recall and F1 values in the validation set were 0.96, 1.00, 0.96 and 0.98, respectively.

In Figure 7A, to further investigate the relationship between relative combined metabolite panel content and patient survival outcomes, we obtained risk scores by constructing a risk model and divided patients into low-risk (n = 40) and high-risk (n = 39) groups using the median of risk scores as the cutoff value. All baseline factors and risk factors were included in the multivariate Cox analysis (Figure 7B), and age > 65 years was significantly associated with shorter overall survival (HR, 2.990; 95% CI, 1.440-6.380; P = 0.005), whereas the low-risk group had significantly longer survival than the high-risk group (HR, 0.280; 95% CI, 0.120-0.640; P = 0.003). Multivariate Cox regression analysis indicated that risk score and age were independent prognostic factors for OS. Figures 7C, D present OS and PFS comparisons between the high-risk and low-risk groups. The median OS in the high-risk group was 18.0 months (95% CI, 15.0-21.0), which was significantly shorter than that observed in the low-risk group (P = 0.00023). In the high-risk group, median PFS was 7.0 months (95% CI, 5.9-8.1), which was significantly lower than PFS in the low-risk group, which was 15 months (95% CI, 10.0-20.0), and the difference was statistically significant (P < 0.0001). Time-dependent ROC curves indicated AUCs of 0.79, 0.90, and 0.73 at the 1-year, 2-year, and 3-year intervals, respectively, for the risk factors in forecasting OS (Figure 6E). These results indicated that our model could potentially predict the survival outcomes of patients.