Machine learning identifies a 5-serum cytokine panel for the early detection of chronic atrophy gastritis patients

2.1Ethics statement

All studies involving human participants were reviewed and approved by Institutional Ethics Committee of Nanjing Medical University (KY21069). The patients provided their written informed consent prior to enrollment.

2.2Study population inclusion and exclusion criteria

The target population, aged between 40 and 69 years, were invited to participate in the study from January to December of 2021. In China, the incidence of GC rises rapidly after the age of 40 [30]. Additionally, as this study was a community-based GC screening project, individuals over the age of 70 were excluded as they have a higher risk of complications associated with undergoing gastroscopy and biopsy. Individuals were excluded from the study for any of the following reasons: history of gastric surgery, including endoscopic mucosal resection or submucosal dissection, coagulation disorders or severe cardiovascular or cerebrovascular diseases, liver, kidney, or psychiatric disorders. In addition, patients were excluded if taking proton pump inhibitors within 2 weeks because these can influence PG Levels in the serologic tests. In addition, patients who routinely take an antiplatelet drug such as aspirin were excluded, as these drugs can cause gastrointestinal bleeding during the endoscopic biopsy procedure.

2.3Sample size calculation

We based our sample size calculation on obtaining a sensitivity and specificity of 0.7. From prior experience approximately 30–40% of patients screened at our hospital have CAG. Using this information, we calculated that for a sensitivity of 0.7 and with 35% positivity rate we would require 230 patients for this study [31].

2.4Patient characteristics

A total of 247 samples, including 157 CSG patients (the control group for this study) and 90 CAG patients, met the inclusion criteria and were included in our study. The overall male to female ratio was 103:144. The median age was 58. The clinical characteristics and cytokines levels of the two cohorts are summarized in Tables 1 and 2 respectively.

2.5Human blood sample collection and serum preparation

Blood samples were obtained immediately prior to the endoscopy procedure for 247 patients as part of a community-based screening program for gastric cancer between January and December of 2021. The blood was collected into two vacutainers (5 ml into each) containing K2 Ethylenediamine tetraacetic acid (EDTA), one sample was used for multiplex cytokine concentration detection by Meso Scale Discovery (MSD) and the other for pepsinogen detection by enzyme-linked immunosorbent assay (ELISA). The blood samples were mixed by gentle inversion and kept at room temperature (RT) for no more than 2 hours, they were then centrifuged at 3,000 × g for 10 min at 4^∘C, The upper one-third of the serum supernatant was collected. The serum samples were aliquoted and stored at -80^∘C until sample analysis.

2.6Biopsy collection and histological reporting

All participants underwent a biopsy procedure during the gastroscopy procedure. Two biopsies were taken, one from the antrum and one from the small curvature of the stomach for histology examination. We did not do the full 5-point biopsy according to the OLGIM (Operative Link on Gastric intestinal metaplasia assessment) evaluation system because a severity grading of atrophy was not required for this study.

2.7Diagnosis of the CSG and CAG by histological examination of biopsy samples

Chronic Superficial Gastritis (CSG) was characterized by the gastric mucosa having lymphocytes and plasma cell infiltration without changes in atrophy or intestinal metaplasia by pathological analysis. Gastric Atrophy (GA) refers to the reduction of intrinsic glands in the gastric cavity, the thinning of gastric mucosa, and the shallowing of gastric pits by pathological analysis. Gastric Intestinal Metaplasia (GIM) is the replacement of gastric mucosa epithelial cells by intestinal mucosa epithelial cells (goblet cells, pan cells and absorptive cells) by pathological analysis [32]. In this study, we collectively referred to GA and GIM as Chronic Atrophic Gastritis (CAG) [33].

There are two types of CAG: a gastric antrum predominant type in patients infected with H.pylori, and an autoimmune type limited to the gastric body and fundus [34]. Due to the rare occurrence of the autoimmune type, a distinct subgroup of patients with specific biomarkers gastrin G17 and pepsinogen [35], none of the autoimmune type were included in this study.

2.8Serum detection of Pepsinogen I (PGI) and Pepsinogen II (PGII)

Serum expressions of Pepsinogen I (PGI) and Pepsinogen II (PGII) were determined using Enzyme-linked immunosorbent assay (ELISA). A total of 5 ml of blood was collected from each of the 247 patients, and the serum was separated from the blood by centrifugation at 3000 × g for 10 minutes at 4^∘C. The PGI and PGII concentrations (in nanograms per liter) were measured using an ELISA kit (Biohit Co., Ltd., Finland) immediately according to the manufacturer’s instructions. Briefly, both standard and serum samples were added to the ELISA plates (100 μl/well), followed by a 2-hour incubation at 37^∘C. After that, detection solution A was added (100 μl/well) and incubated at 37^∘C for 1 hour. The plates were washed 3 times, and detection solution B was added (100 μl/well) and incubated at 37^∘C for 1 hour, followed by 5 washes. Tetramethylbenzidine (TMB) substrate was added (90 μl/well), and the plates were incubated at 37^∘C for 15–25 minutes. Finally, the reaction was terminated by adding termination solution (50 μl/well), and the optical density (OD) value was read out at 450 nm.

2.9Cytokine concentration quantification by multiplex Meso Scale Discovery (MSD) enzyme linked immunoassay

The following 40 cytokines: Tumor necrosis factor-α (TNF-α). Tumor necrosis factor-β (TNF-β), Eotaxin, Eotaxin-3, interferon-inducible protein-10 (IP-10), Stromal cell-derived factor 1α (SDF-1α), Macrophage migration inhibitory factor (MIF), Macrophage Inflammatory Protein-5 (MIP-5), Interferon gamma (IFN-γ), Macrophage Inflammatory Protein-3α (MIP-3α), MCP-4, Eotaxin-2, Interleukin-12p70 (IL-12p70),Macrophage-derived chemokine (MDC), Macrophage inflammatory protein-1β (MIP-1β), Interleukin-13(IL-13), Interleukin-15(IL-15), Interleukin 17A(IL-17A), Interleukin-1RA (IL-1RA), Interleukin-1β (IL-1β), Interleukin-2 (IL-2), Interleukin-4 (IL-4), Interleukin-5 (IL-5), Interleukin-6 (IL-6), Interleukin-8 (IL-8), Interleukin-10 (IL-10), Interleukin 12 (IL-12/IL-23p40), Interleukin-15 (IL-15), Interleukin 16 (IL-16), Interleu-kin-1α (IL-1α), Interleukin7 (IL-7), Interleukin-1RA (IL-1RA), Interleukin-2Ra (IL-2Ra), Interleukin-3 (IL-3), Interleukin-33 (IL-33), Granulocyte-macrophage colony-stimulating factor (GM-CSF), Epithelial-neutro-phil activating peptide (ENA-78), and Vascular endothelial growth factor A (VEGF-A) were measured using the Meso Scale Discovery platform (MSD) according to the manufacturer’s instructions (product insert of K15067, MSD, Rockville, MD).

MSD is an indirect binding quantitative electrochemiluminescence (ECL) method designed to detect multiple antibodies in human serum simultaneously. In the MSD platform, biological reagents (antigens or antibodies) are coated on the carbon electrode at the bottom of a micro titer plate. Bound samples are then incubated with a ECL labeled (SULFO-TAG) detection antibody. Following application of an electrical current across the carbon electrode, the SULFO-TAG reacts with the Ru(bpy)32+ reagent and Tripropylamine (TPA) catalyst in the detection buffer and resulting in the emission of light that can be quantitated at 620 nm with a CCD camera.

2.10Data pre-processing and statistics

Three cytokines, IL-2, IL-3 and IL-1α were not detected in 42.5%, 69.2% and 45.6% of patients, and therefore, were excluded in subsequent machine learning analyses. For the remaining 37 cytokines, the average missing rate is 0.52% (the median is 0). The few missing values for these cytokine concentrations were filled in using the lower limit of detection (LLOD) of the MSD assay divided by 10, reflecting the fact that the concentration of these cytokines is too low to be detected. Student t-test was used to assess the statistical significance of cytokine levels between the different groups. An overall adjusted p-value of 0.05 using Benjamini-Hochberg for multiple testing correction was considered as the significance cut-off for all statistical analyses. All statistical analyses and machine learning algorithms for this study were performed using R version 4.1.0 and Python version 3.8.8 [36]. In order to reduce overfitting in the downstream analyses from the feature selection and machine learning, we randomly divided the whole dataset into a training dataset (Dataset1) and a validation dataset (Dataset2) for five times, while maintaining the proportion of CAG and CSG for each division (Supplementary Table 1). The Dataset2 here is considered as a validation dataset since it is a dataset the algorithm did not include in the training process and was only used for validation. In detail, the Dataset1 includes 125 CSG and 72 CAG, while Dataset2 includes 32 CSG and 18 CAG.

2.11Feature selection

To reduce over-fitting and improve the generalizability of the model feature selection was employed [37]. The machine learning model then only utilizes the selected features as its input.

Feature selection was performed using the Boruta algorithm [27, 28] and other redundant or non-important features were discarded. Boruta was used to select the features that correlated with the dependent variable. Boruta feature selection was performed on each of the 5 training datasets and the common features from all five data splits were then used by the ML algorithm. The Boruta method is a wrapping algorithm that is based on the random forest (RF) method. It shuffles the original real features to construct shadow features, then joins the real features and shadow features into a feature matrix for training. Finally, the feature importance score of shadow features is used as the reference base to select the feature set truly related to the dependent variable from the real features (see Supplementary Method for more detail).

The Boruta program that we used in this work was Python package Boruta (version 0.3), the default parameters were selected for execution.

2.12Machine Learning Algorithms: LightGBM

LightGBM is a well-known gradient boosting framework that uses a tree-based learning algorithm [38, 39, 40]. It can be regarded as an improved version of Gradient Boosting Decision Tree (GBDT), which uses the negative gradient of the loss function of the current Decision Tree (DT) as an approximation of the residual and fits a new DT recursively. It is designed to be distributed and efficient by incorporating two new technologies, gradient based one side sampling (GOSS) and exclusive feature bunching (EFB) [41], both of which greatly improve the efficiency and ensure the accuracy of classification.

In this study, the features are modeled using LightGBM and the parameters are continuously optimized by means of a grid search features (see Supplementary Methods for more detail). LightGBM program (version 3.3.5) is implemented through a Python package, Scikit-learn [36].

2.13Model evaluation

Model evaluation is a crucial step in machine learning to assess the performance of a predictive model. There are various metrics that can be used to evaluate the model’s performance, and here we mainly use, area under the curve (AUC), area under the precision recall curve (PR-AUC) and accuracy. AUC is the area under the receiver operating characteristic curve (ROC curve) and is a metric that evaluates the model’s ability to distinguish between positive and negative classes. The ROC curve is created by plotting the true positive rate (TPR) against the false positive rate (FPR) at different threshold values. PR-AUC, the area under the precision-recall curve, is a metric that evaluates the model’s ability to identify positive class instances while minimizing false positives. The PR-AUC is calculated as the area under the precision-recall curve. Precision is the fraction of true positives (TP) among the samples predicted as positive (TP + false positives), and recall (or sensitivity) is the fraction of true positives among all the actual positives (TP + false negatives). Accuracy is a metric that measures the ratio of correctly predicted observations (true positive + true negative) to the total number of observations. From the design of these different metrics, AUC and PR-AUC focus on the detection of the positive samples, while the accuracy metric gives the same “importance/weight” on both the positive samples and negative samples. In this study, our aim is to accurately identify the positive samples (CAGs), therefore, we evaluated the predictive models primarily based on AUC and PR-AUC.

2.14The subgroup analysis of logistic regression model

Hypertension, diabetes, coronary heart disease, stroke, and chronic hepatitis B are common chronic diseases. We examined whether these common chronic diseases were confounding factors in this study. A subgroup analysis using a logistic regression model was performed on these five confounding factors using the R software (R version 4.2.3), and then the differential effects of the 5-cytokine panel diagnosis were compared in each subgroup dataset. By examining the interaction in the subgroup regression of the generalized linear model (p for interaction), we assessed whether these confounding factors had an impact on the 5-cytokine panel diagnostic effects.

3.1Serum cytokine measurements

For all 247 patients, we collected blood samples and extracted serum. We measured the concentrations of 40 key cytokines in the serum using MSD. The raw cytokine concentration data were cleaned, as described in the Materials and Methods section, leaving data on 37 cytokines for further analysis (Table 2).

Figure 1.

Schematic showing the details of the ML-based biomarker discovery pipeline following the serum cytokine measurements.

3.2Machine learning based biomarker discovery study design

To optimally utilize the information from our patient cohort and develop a predictive model for CAG patients, we designed the process outlined in Fig. 1. Starting with the cytokine concentration matrix composed of 247 patient samples (including 157 CSGs and 90 CAGs) and 37 cytokine features, we first performed feature analysis in order to identify the most relevant features for our downstream analyses. We calculated the statistical significance for all features among the two groups, to determine the difference of these features between CSG and CAG patients (the adjusted p-value is shown in Table 2). Only two cytokines TNF-a and Eotaxin were found to be significantly different between the groups. Then, we utilized the Boruta algorithm [27, 28] to select the features (all cytokine features plus the clinical features) with the highest capacity for discrimination between CSG and CAG patients.

Once we had identified five key features, we developed a machine learning model based on the selected features by utilizing the LGBM model [41]. We divided the whole dataset (a total of 247 samples) into training set (Dataset 1) and test set (Dataset 2), and we repeated this split randomly five times to avoid over-fitting. We used Dataset1 for training with 5-fold cross-validation and Dataset2 for validation (i.e., validation in a dataset the algorithm has never seen).

Figure 2.

Feature importance and correlation analysis. (A). Heatmap of feature importance based on the Boruta feature selection method (scores range from 1 to 35), performed five independent times through five random splits of the data. Showing the results for 37 cytokines, Gender, Age, PGI, PGII, and PGI/PGII ratio. The five cytokines with the highest ranking are displayed in bold. (B). Heatmap of absolute values of Pearson correlation coefficients of features among samples. The five selected features show low correlation with each other.

3.3Feature Selection enables us to reduce the input features for our subsequent Machine Learning modeling

Clinical features (PGI/PGII ratio, age and gender, Table 1) and cytokine features (Table 2) were first analyzed using the Boruta method [27, 28]. For this feature selection procedure, we exclusively used the training set. For every one of the five independent training sets (Supplementary Table 1), we utilized Boruta for feature selection, and calculated the feature importance, ranking their capacity to distinguish CAG vs CSG patients. In turn, we selected the features which ranked first (i.e. rank of “#1”) for each separate training (see Fig. 2A and Table 3 for details, see Supplementary Table 2 for all features ranking). Subsequently, we took the intersection of these five separate feature selection attempts, which yielded a total of 5 features: IL-10, TNF- α, Eotaxin, IP-10 and SDF-1a. These features (cytokine levels) were retained as the key features that would be used for the follow-up modeling (see Table 3, in bold). We also did feature selection using another method Lasso [42], these 5 features also rank very high by Lasso (see Supplementary Table 3).

Furthermore, we compared the differences of each feature in the two groups through statistical testing (student’s t-test), and we verified that four of these five cytokine features, TNF-α, Eotaxin, IP-10 and SDF-1a, demonstrate statistically significant differences between the two groups (P< 0.05) (Supplementary Fig. 1), although two of them, IP-10 and SDF-1a, showed no significance after correction for multiple testing (see adjusted P-value in Table 2). In addition, we analyzed the correlation between different features among patients (see Fig. 2B and Supplementary Table 4) and found that some features showed strong correlations with each other. This indicating the redundancy of some features and the benefits of the feature selection performed here. The correlation plot also shows the independent predictive power of the five features we selected here as they did not correlate with each other (Fig. 2B), again shows the validity of the feature selection process. In addition, the 5 selected cytokines showed no significant correlation with age or gender demonstrating these clinical features would not confound the results.

In order to assess whether common health conditions could confound the 5-cytokine panel diagnostic ability, we conducted a subgroup analysis using a logistic regression model on common confounding factors. These results demonstrated that hypertension, diabetes, coronary heart disease, stroke, and chronic hepatitis B all had insignificant interaction values with these cytokines, with the one exception for hypertension with IP-10 levels (p= 0.009). Therefore, the 5-cytokine diagnostic panels efficiency is largely unaffected by the above confounding factors (see Supplementary Table 5).

3.4Building an accurate predictive model for CAG patients, using LGBM

Figure 3.

The training performance of the five cytokine features derived LGBM modelling based on Training1 of division 1. (A). The mean AUC curve is shown in blue. The light gray shade shows the standard deviation stemming from the five-fold-cross validation results. (B). The PRAUC mean is shown in red. The light gray shade shows the standard deviation stemming from the five-fold-cross validation results. AUC: Area Under the Curve, PRAUC: Precision Recall Area Under the Curve.

After reducing the number of the features to those that show significant predictive capacity for CAG patients, we moved on to develop a predictive model for CAG patients. For our predictive model, considering the relatively small sample size, the light gradient-boosting machine (LGBM) classifier was selected [41]. To avoid overfitting – a typical challenge in machine learning approaches – we randomly divided the whole dataset into training set (Dataset 1) and validation set (Dataset 2) for five times, as specified in the Material and Methods section. For each division, there is a training set as well as its corresponding validation set (Supplementary Table 1). We further performed five-fold-cross-validation on the training set to avoid overfitting to the maximum degree. The LGBM model included the five selected features, namely IL-10, TNF- α, Eotaxin, IP-10 and SDF-1a, as input candidate features. The classifier with the best model performance was selected by continuously iterating the optimization parameters for the training set of each division. Out of all 5 different independent runs, the training set of division 1 (sheet of “Training1” in Supplementary Table 1) is shown in Fig. 3; the performance of the other divisions (Training2 to Training5) are provided in Supplementary Fig. 2, and demonstrate similar performance. In summary, the performance of our model was quite high, with AUC = 0.85, and PRAUC = 0.81. Furthermore, we also constructed the predictive models using other 6 machine learning classifiers, i.e., Logistic Regression (LR), Decision Tree (DT), Support vector Machine (SVM), Neural Network (NN), Naive Bayes (NB), K Nearest Neighbors (KNN), and compared their performances with the model constructed by LGBM. Overall, the LGBM model’s performance is better than other 6 machine learning models (see Supplementary Table 6 and 7).

Figure 4.

The performance of the five cytokines based LGBM model in the independent, validation dataset shows high predictive capacity. (A). Area under the curve analysis. (B). Precision recall area under the curve analysis. AUC: Area Under the Curve, PRAUC: Precision Recall Area Under the Curve.

3.5Validation

The final LGBM model was validated in its corresponding validation dataset (“Validation1” in Supplementary Table 1). The obtained AUC = 0.88, PRAUC = 0.76, Accuracy = 0.78, Sensitivity = 0.78, Specificity = 0.78 show that our model’s capacity to accurately predict CAG patients is high. The validation results are shown in Fig. 4.

In the same way as the training sets (division 1), we also evaluated the performance of each of the other 4 training division derived models on their respective validation set (Supplementary Fig. 3). We have summarized the validation performance of the different training set derived model in Table 4. It clearly shows that the models’ performance across all different data partitions stays high and shows robustness. As the main goal of our work here is to accurately identify CAG patients, the metrics, “AUC” and “PRAUC” are most relevant as they focus on the true positive – the CAG patients (as opposed to using the “accuracy” metric (TP+TN)/total samples). By inspecting the variance of the models’ performance (AUC, PRAUC and Accuracy) among the five divisions, we can clearly see the limited variance, which indicates the model is very stable and robust against changes in the training datasets. Therefore, we conclude that, the derived model is valid and can be generalized.

To further demonstrate that the five features selected by the Boruta method led to the best prediction algorithm, we compared the performance of our 5-cytokine based ML model with a ML model that includes the complete set of 37 cytokines we measured (Supplementary Fig. 4) as well as a 40-feature ML model which contains all available features (37 cytokines + 3 clinical features) (Supplementary Fig. 5). In addition, we also built a model based only on the significant cytokines, (TNF- α, Eotaxin; adjusted P value < 0.05) (Table 2), with the aim to see whether our feature selection in combination with ML outperforms this simple combination (Supplementary Fig. 6). In summary, our results demonstrate that the performance of our concise 5-feature (cytokine) predictive model is comparable with the performance of the much more complex multi-feature models, and better than the simple statistical test derived marker model, which demonstrates that the five selected features suffice to empower our predictive model to accurately distinguish CAG patients.

We next examined how our model compared to white-light gastroscopy, a slightly less invasive procedure than the gold-standard gastroscopy/biopsy procedure, as it does not require the collection of biopsies. We found that clinician diagnosis of CAG using white-light gastroscopy (Supplementary Table 1) had a lower accuracy than the 5-cytokine ML model (0.67 vs 0.78) and a dramatically lower sensitivity (0.4 vs 0.78). This sensitivity is similar to a previously published study that indicated white-light gastroscopy had a sensitivity of 0.41 [43]. In contrast white-light gastroscopy had a slightly higher specificity than the 5-cytokine ML model (0.83 vs 0.78). This demonstrates that the 5-cytokine ML model outperforms white-light gastroscopy on accuracy and sensitivity in addition to being a less invasive procedure for the patients.

Figure 5.

(A). The difference of PGI/PGII between CSG and CAG patients is significant (p-value = 0.023, using student’s t-test). (B). The performance of PGI/PGII in classifying CAG with the logistic regression: the classifier does not have strong discrimination capacity. CSG: Chronic Superficial Gastritis; CAG: Chronic Atrophic Gastritis; PGI: pepsinogen I; PGII: pepsinogen II.

3.6Evaluation of the predictive power of PGI/PGII in our dataset

As previous approaches to identify CAG patients have utilized PGI and PGII measurements, we wanted to investigate how our approach compared with these more established approaches. We found the PGI/PGII ratio to be significantly different between the two sample groups in this study using a standard student t-test (P= 0.023) (Fig. 5A).

First, we employed the criteria that has been used in a number of publications [44], (PG I ⩽ 70 ng/mL and a PG I/PG II ratio of ⩽ 3) to interrogate our full dataset, we found only 3 samples met this criterion (Sample IDs: 1907649, 1908664, 1907503, see detail in sheet “Complete_data” of the Supplementary Table 1). Only one of these was a true positive, CAG (1908664), and other two were false positives, CSG (1907649, 1907503), this gives rise to a Sensitivity = 0.01 and Accuracy = 0.63 (of note, only the PG I/PG II ratio influenced the outcome in this analysis). This performance is dramatically lower than the 5-cytokine approach that we have developed. We next sought to determine whether the utility of PGI/PGII could be improved with the aid of machine learning. To properly evaluate the discriminatory capacity of PGI/PGII on CAG detection, we constructed a model based on logistic regression to predict CAG patients, using the PGI/PGII ratio. We performed five-fold cross validation by employing the entire sample set, which leads to a model with a mean AUC of 0.59 (Fig. 5B). This result strongly suggests that a model based solely on the PGI/PGII information fails to accurately discriminate CAG patients from CSG patients. Taken together, these results demonstrate that with or without machine learning, the common approach for CAG identification utilizing PGI and PGII concentrations cannot achieve the performance of our cytokine-based model.

In addition, we tested whether the addition of the remaining clinical features (age and gender) to the PGI/PGII ratio could improve the prediction performance. We therefore trained a LGBM model using these three features and evaluated the performance in its capacity to discriminate CAG vs CSG patients. These results confirm once more that the clinical features alone do not have strong predictive capacity (see Supplementary Fig. 7), and is in line with our previous findings, which demonstrated that their feature importance did not typically rank highly (Fig. 2A).

4.1Limitations of the study

It is important to note that this was a single center study and that further validation of this 5-cytokine panel machine learning algorithm in additional patient cohorts would be critical before this approach could be adopted in the clinic. Furthermore, while the incorporation of 5 cytokines in this model rather than a single cytokine measurement does increase the specificity of this test it is possible that some other cancer types or diseases may lead to a similar cytokine profile and thus lead to false positive results.

As there were no autoimmune type CAG patients included in this study, further studies would need to be conducted to demonstrate whether this 5 cytokine panel has any utility in this patient subgroup.

Finally, while most components of the GastroPanel, the alternate serum screening method, were compared to the 5-cytokine panel model, a direct comparison with the GastroPanel was not included in this study, we thus cannot definitively state that this new system is superior to the GastroPanel.