Plasma-based lipidomics reveals potential diagnostic biomarkers for esophageal squamous cell carcinoma: a retrospective study

Chemicals and reagents

High performance liquid chromatography (HPLC) grade methanol, acetonitrile (ACN), isopropanol (IPA), methyl-tert-butyl ether (MTBE), and ammonium acetate were purchased from Merck (Darmstadt, Germany). Distilled water was purchased from Wahaha Group Co., Ltd. (Hangzhou, China). Formic acid (HPLC) was purchased from Roe Scientific Inc. (Delaware, USA).

Participants and sample collection

This retrospective study analyzed plasma samples obtained from Zhejiang Cancer Hospital (Hangzhou, China) between December 2010 and December 2012. Plasma samples were collected from 40 pathologically diagnosed ESCC patients and 31 healthy controls (HCs). Table 1 presents the basic clinical information of the participants. Blood was collected from individuals who had fasted overnight and transferred into vials pre-treated with the anticoagulant reagent (ethylenediaminetetraacetic acid disodium potassium salt). Plasma was obtained by centrifuging the blood at 2400xg for 8 min. The samples were then stored at −80 °C until analysis. All procedures involving human participants were conducted in accordance with the ethical standards set by the Ethics Committee of Zhejiang Cancer Hospital (IRB- 2019-66), following the principles of the 1964 Helsinki Declaration and its subsequent amendments or comparable ethical standards. Furthermore, since the samples used in our study came from the biobank of Zhejiang Cancer Hospital, the Medical Ethics Committee of Zhejiang Cancer Hospital waived the need for informed consent.

Lipid extraction from plasma

Plasma lipid extraction was conducted referring to previous publication (Yang et al., 2020) with some modifications. Briefly, 300 µL of chilled methanol was added to 40 µL of plasma, followed by mixing for 1 min. Then, 1 mL of MTBE was added and the mixture was shaken at a frequency of 60 Hz for 1 h at room temperature. After that, 250 µL of water was added, and incubated on ice for 10 min. An aliquot of 400 µL of the upper organic phase was transferred to a new tube after separation by centrifugation at 16,200xg at 4 °C for 15 min, and subsequently was dried using a concentrator. The dried residue was reconstituted by mixing with 80 µL ACN-IPA-water (65:30:5, v/v/v) followed by centrifugation at 16,200xg, at 4 °C for 15 min, and an aliquot of 60 µL of supernatant was transferred to sample vial. Finally, 5 µL of supernatant was loaded to liquid chromatography-mass spectrometry (LC-MS) analysis.

The pooled quality control (QC) plasma samples were generated by combining equal aliquots of plasma from each individual sample, which were then dispensed into 40 µL volumes. The extraction process employed for these pooled samples was identical to that used for the individual sample pretreatment.

LC-MS analysis

LC-MS analysis was conducted according to previous research (Yang et al., 2022). In brief, ultimate 3000 UHPLC system coupled with Q exactive orbitrap mass spectrometer (both from Thermo Fisher Scientific, Waltham, MA, USA) was used for lipidomics analysis. Chromatographic separation was carried out on an Acquity UPLC BEH C18 column (2.1 mm × 100 mm, 1.8 µm, Waters, Milford, MA, USA). Solvent A consisted of a mixture of ACN/water (3:2, v/v) containing 0.1% (v/v) formic acid and 10 mM ammonium acetate, while solvent B was composed of IPA/ACN (9:1, v/v) with the same additives. The flow rate was 0.3 mL/min, and column temperature was set at 50 °C. The elution condition was set at 0.0−1.5 min, 32% B; 1.5 min-15.5 min, 32%–85% B; 15.5–15.6 min, 85%–97% B; 15.6–18.0 min, 97% B; 18.0–18.1 min, 97%–32% B; 18.1–20.0 min, 32% B. The settings for the mass spectrometer included a capillary voltage of 3.0 kV and a capillary temperature of 300 °C. The sheath gas flow rate was set to 50 Arb. The auxiliary gas had a flow rate and temperature of 15 Arb and 320 °C, respectively. The scan range was set at m/z 100–1200. The full scan MS had a resolution of 70,000 and an AGC target of 3 × 10⁶. The data-dependent MS/MS had a resolution of 17,500 and an AGC target of 1 × 10⁵. The normalized collision energy was set to 30, 40, and 50 eV, respectively.

The analytical procedure employed a full scan mode to collect data from all samples in the batch. For qualitative quality control (QC) samples, data-dependent acquisition (DDA) mode was utilized. To ensure consistent performance and accuracy during the analysis, QC samples were interspersed within the sample injection sequence. The sequence commenced with three consecutive QC samples, followed by the inclusion of one QC sample every ten samples. The sequence concluded with another three consecutive QC samples. This approach helped to correct for mass spectrometry signal fluctuations and maintain reliable data quality throughout the analysis.

Metabolomics data analysis

Data analysis was performed according to existing articles (Yang et al., 2022) with some modifications. Briefly, the research utilized ProteoWizard’s msconvert tool (https://proteowizard.sourceforge.io/download.html) to convert RAW format data into mzXML format data. The R package xcms was then employed for detecting and extracting ion features, which included tasks such as peak picking and retention time correction. To correct signal shifts, the R package statTarget utilized QC-based random forest signal correction (QC-RFSC). Subsequently, ion feature filtration was conducted. In this step, variables were retained if they had a non-zero value in at least 80% of samples within any single group. However, variables in QC samples with a relative standard deviation (RSD) greater than 30% were excluded. Imputation was performed using the K-nearest neighbors (KNN) algorithm. Prior to chemometrics analysis, the detected ions in each sample belonging to the same class were normalized by setting the sum of their peak areas to 100,000. This rigorous approach ensured that only the most consistent and reliable features were retained for further analysis, thereby significantly improving the overall robustness and reliability of the lipidomics data.

Train set samples were used to discover the differentially expressed plasma lipids between ESCC and HC groups. Firstly, unsupervised principal component analysis (PCA) was employed to visualize the overall separation trend of all samples based on the ion features. Subsequently, a supervised partial least squares discriminant analysis (PLS-DA) was utilized to assess the classification ability of these ion features, yielding a variable importance in projection (VIP) value for each ion feature. Furthermore, the statistical significance of the ion features between the ESCC and HC groups was evaluated using a two-tailed Student’s t-test, and the Benjaminii-Hochberg false discovery rate (FDR) was also calculated. The criteria for defining differentially expressed ion features were as follows: VIP >1.0, adjusted p-value (FDR) <0.05, and fold change (FC) >1.50 or <0.667.

The annotation of lipids was as previously described (Tsugawa et al., 2020). Briefly, the RAW format data was first converted to the abf format using Abf converter (https://www.reifycs.com/abfconverter/). Then, MS-DIAL software ver.4.9 (http://prime.psc.riken.jp/compms/index.html) was used to perform feature detection on all ions with the following parameters: the tolerance of MS1 and MS2 were set 0.01Da and 0.025Da; identification score cut off was set 80%; in positive mode, [M + H]⁺, [M + NH4]⁺, [M + Na]⁺ and [M + H-H2O]⁺ were selected as the adduct types; in negative mode, [M - H]⁻ was selected as the adduct type; the retention time was tolerance set to 0.05 min and MS1 tolerance was set to 0.015 Da in all ions feature alignment option.

The differential lipids between ESCC and HC groups in both train set and test set were illustrated using heatmap.

Diagnostic significance of plasma lipid

In order to assess the diagnostic value of plasma lipids, receiver operating characteristic (ROC) curve analysis was performed for each differentially expressed lipid in the train set using the R package pROC. This analysis allowed calculation of AUC values. The top nine lipids with the highest AUC values in the ROC curve were selected as the most diagnostic plasma lipids.

To determine the optimal cutoff value in the train set, the threshold was set at the point where the Youden Index (Sensitivity + Specificity − 1) was maximized. This approach aims to find the threshold that maximizes the difference between the true positive rate (sensitivity) and the false positive rate (1 − specificity), striking a balance between sensitivity and specificity.

To evaluate the diagnostic value of the plasma lipids, the cutoff value determined in the train set was used to predict the classification of samples in the test set. The performance of each lipid was then assessed using a confusion matrix, which allowed calculation of diagnostic metrics such as sensitivity, specificity, and accuracy. By following this procedure, the study aimed to identify the most informative plasma lipids for diagnostic purposes and evaluate their performance in classifying samples in both the train and test sets.

In addition to single lipid-based diagnostic models, multiple lipid-based ML models were investigated, including partial least squares (PLS) and random forest (RF) from the caret package, and support vector machine (SVM) from the e1071 package. ROC curves were plotted for each model using the train set, and their prediction performance was evaluated in the test set using confusion matrix calculations.

Relationship between plasma lipid levels and clinical features

To establish the relationship between differentially expressed plasma lipid levels and clinical features such as sex, age, drinking history, smoking history, lymph node metastasis (LNM), and TNM stage, statistical tests were conducted. The Kruskal-Wallis test or Wilcoxon test was employed to compare the lipid levels among different groups in all samples. The p-value of less than 0.05 was considered statistically significant, indicating a significant association between the lipid levels and the clinical features being examined.

Basic characteristics of the participants

This study enrolled 40 patients pathologically diagnosed ESCC and 31 healthy controls (HCs). Age and sex were matched between the two groups without any statistically significant difference, as shown in Table 1. Among the ESCC patients, nine were in stage I, 15 were in stage II, 16 were in stage III while no patients were in stage IV. The samples were randomly divided into a training set (ESCC = 28, HC = 20) and a test set (ESCC = 12, HC = 11).

Metabolic shift between ESCC and HC

Metabolic features were analyzed in train set and depicted in Fig. 1. Following peak picking, retention time alignment, grouping, and signal shift correction, a total of 41,028 ion features were obtained, including 22,591 positive ions and 18,437 negative ions. Based on these metabolic features, PCA and PLS-DA score plots were generated to investigate the differences between ESCC and HC groups. The results exhibited a significant separation between the ESCC and HC groups, indicating notable dysregulation in the plasma lipid profile of ESCC patients (Figs. 1A and 1B). The reliability of the PLS-DA result was further validated using a permutation test (n = 20) (Fig. 1C). Furthermore, based on the criteria of differential lipid (FDR < 0.05; VIP > 1.0; FC > 1.50 or FC < 0.667), a total of 5,899 differential metabolic features were shown in the volcano plot, including 2,654 up-regulated features and 3,245 down-regulated features (Fig. 1D).

Differentially expressed plasma lipids in ESCC

A total of 99 differentially expressed plasma lipids were identified and found to exhibit significant differences between individuals from ESCC and HC. Among these lipids, 15 were upregulated (FC >1.5) and 84 were downregulated (FC <0.667) in the plasma of ESCC patients compared to HC. The detailed information of these lipids is summarized in Table 2. Heatmap displaying the expression patterns of the differentially expressed lipids in the train set (Fig. 2A) revealed evident differences between ESCC and HC. Among these lipids, approximately three-fourths exhibited a downward trend in ESCC patients compared to HC. Similar expression patterns were observed in the test set (Fig. S1). The proportion of lipid classifications was presented in the pie chart (Fig. 2B), including 21 fatty acids (FAs), 22 glycerolipids (GLs), 37 glycerophospholipids (GPs), 18 sphingolipids (SPs), and one sphingomyelin (SM). GPs represented the largest proportion of differential lipids at 37%, with 34 downregulated lipids and three upregulated lipids. Among the upregulated differential lipids (Fig. 2C), there were eight FAs, three GPs and SPs, and one GL. Among the downregulated differential lipids (Fig. 2D), there were 13 FAs, 21 GLs, 34 GPs, 15 SPs, and one SM.

Diagnostic performance of the differential plasma lipids

ROC curve analysis revealed that plasma lipids show promise as diagnostic biomarkers for ESCC. The top nine lipids with the highest AUC values were FA 15:4, FA 27:1, FA 28:7, FA 28:0, FA 36:0, FA 39:0, FA 42:0, FA 44:0, and DG 37:7. Eight of these nine lipids achieved an AUC value of 1.00, indicating excellent diagnostic accuracy. FA 15:4 and DG 37:7 showed an up-regulated trend in ESCC samples, while the remaining lipids showed a down-regulated trend (Fig. 3).

In testing step, the confusion matrix charts (Fig. 4) were constructed for a more intuitive representation of the diagnostic performance of the top nine lipids. All the nine lipids depicted in the chart exhibited prediction accuracy exceeding 85%, with six exceeding 95%. And seven of these nine lipids achieving 100% diagnostic efficiency in detecting tumors. The top 15 lipids, ranked by their prediction accuracy in the test set, were selected. Table 3 presents their AUC values in the train set, as well as prediction accuracy, sensitivity, specificity, precision, and recall values in the test set. Among these lipids, 13 belonged to FAs, and 11 lipids (FA 27:1, FA 28:7, FA 28:0, FA 39:0, FA 42:0, FA 44:0, FA 22:2, FA 36:0, FA 25:0, FA 19:1, FA 27:0) achieved prediction accuracy exceeding 0.90. The results for the remaining differential lipids can be found in the Table S1.

Multiple lipid-based models were constructed using the training set data. The ROC curves in Figs. 5A, 5B and 5C showed high AUC values of 0.990, 0.990, and 0.980, indicating excellent prediction performance. These results were consistent with the performance of individual differential lipids. The accuracy of the models was further validated in the test set, achieving accuracies of 95.7% (Figs. 5D, 5E and 5F). These findings demonstrate the effectiveness of the multiple lipid-based models as accurate diagnostic tools for ESCC.

Relationship between plasma lipid levels and clinical features

The results revealed that 14 differentially expressed plasma lipids were associated with age (Fig. S2), 14 ones with sex (Fig. S3), 18 ones with smoking history (Fig. S4), and only one with drinking history (Fig. S5).

Compared the ESCC patients with and without lymph node metastasis, there were 11 differential lipids showed significant difference (P-value <0.05), including FA 14:0, FA 15:4, FA 27:1, FA 28:7, FA 27:0, FA 42:0, DG 37:7, PC O-16:1_22:6, PC O-17:1_21:4, SM d18:1/23:0, PI 18:0_20:3. Among them, six were up regulation and five were down regulation in ESCC patients with lymph node metastasis. Boxplots illustrating the expression levels of these 11 lipids were shown in Fig. 6. Out of the 99 differential lipids, 10 lipids—FA 14:0, FA 15:4, FA 28:7, FA 27:0, FA 28:0, DG 37:7, PE O-17:1_22:6, PE O-18:3_22:6, PC O-16:1_22:6, and SM d18:1/23:0—exhibited significant differences among the three stage groups (P-value <0.05). Boxplots depicting the expression levels of these 10 lipids are shown in Fig. 7. Among them, FA 15:4 displayed the most prominent upregulated trend in stage III compared to stages I and II, with the P-value of 7.27e−07. Only one of these 10 lipids, PE O-18:3_22:6, exhibited statistically significant differences in pairwise comparisons across stages I, II, and III, with an increase in expression levels in stage II but declining in stage III patients.