EsoDetect: computational validation and algorithm development of a novel diagnostic and prognostic tool for dysplasia in Barrett’s esophagus

Dataset search

An exhaustive search for public datasets containing gene expression data related to BE, including normal esophageal epithelium, NDBE, BE with different degrees of dysplasia (LGD and HGD) and EAC was performed in the following databases: PubMed (https://pubmed.ncbi.nlm.nih.gov/), Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/), Sequence Read Archive (SRA, https://www.ncbi.nlm.nih.gov/sra), and European Genome-Phenome Archive (https://ega-archive.org/). For the diagnostic application, the aim was to distinguish between NDBE and LGD BE. For the prognostic application, non-progressed BE (nonP-BE) and progressed BE (P-BE) data was studied. P-BE was defined as a BE adjacent to EAC. A summary of the methodology used is represented in Fig. 1 and described in detail below.

Data pre-processing

In this study, raw RNA-seq data from projects GSE193946, GSE58963, and E-MTAB-4054 were obtained from the Sequence Read Archive (SRA) and the European Genome-Phenome Archive (EGA). We processed the data using a Docker environment equipped with Kallisto version 0.46.1 (docker image: jlnetosci/kallisto:v0.46.1), which facilitated the pseudo-alignment of the reads against the Homo_sapiens.GRCh38.cdna.all.release-107 reference transcriptome from Ensembl. Post-alignment, the transcript abundance estimates generated by Kallisto were imported into the R programming environment using the tximport package. This allowed transcript-level data to be transformed into gene-level counts, which were subsequently analyzed for differential expression. The combined data was filtered for low-expressed genes using the filterByExpr function in EdgeR (Robinson, McCarthy & Smyth, 2009), resulting in a dataset of 20,608 genes for downstream analysis. Samples were then normalized using the Trimmed Mean of M-values (TMM) normalization method and differential expression analysis was performed using EdgeR. For downstream analysis, including feature selection and classifier training, log-transformed CPM normalized values were used, which were subsequently corrected for batch effects using the ComBat function from the sva package (Leek et al., 2012).

In this study, microarray data was sourced from the Gene Expression Omnibus (GEO) database using the GEOquery package available in the R software. The data included accessions GSE1420, GSE36223, GSE13083, GSE37200, GSE34619, GSE26886, GSE39491, GSE100843, and an additional dataset from Watts et al. (2007). Data was loaded and normalized using both the affy and oligo packages in R, depending on the array platform. The CEL files were read and processed using the frma function for robust multi-array average (RMA) normalization. Probe-level data was annotated and collapsed to gene-level data using Bioconductor annotation packages hugene10sttranscriptcluster.db, hgu133a.db, hgu133plus2.db, and hgu133a2.db along with the WGCNA package. Finally, the resulting gene expression data was merged into a single dataset for downstream analysis, with additional annotations indicating BE progression status. For prognostic application, 25 genes selected in previous work to distinguish nonP-BE from P-BE (Cardoso et al., 2016), were used in this study—ACTN1, C1S, CCN1 (alias CYR61), CDH1, CEBPB, CEBPD, COL4A1, CTSB, DKK3, DUSP1, IER3, JUN, LAMC1, PLPP3, RBPMS, SNAI1, SNAI2, SPARC, TNS1, TRMT112, TP53, TWIST1, VWF, WWTR1 (alias TAZ) and ZEB1. Box plots representing normalized expression values were generated using the ggplot2 (v3.4.0) and ggsignif (v0.6.4) R packages. Statistical analysis was performed using one-way ANOVA, followed by a post hoc Tukey’s ‘Honest Significant Difference’ test, both from the R stats package (v4.1.1). When ANOVA assumptions were not met, a Kruskal–Wallis Rank Sum Test (R stats package v4.1.1) was performed, followed by a post hoc Dunn’s Kruskal-Wallis Multiple Comparisons test (FSA R package v0.9.3). The significance threshold was set at p-value < 0.05.

Threshold selection and determination of individual predictive power

For the distinction between NDBE and LGD (diagnostic) or nonP-BE and P-BE (prognostic) a Thresholding function was applied to the expression levels of each selected gene to determine an expression threshold. Performance metrics such as accuracy, recall (or sensitivity), precision (or positive predictive value—PPV), specificity, negative predictive value (NPV), and false positive rate (FPR) were calculated for each threshold, considering the known class of the samples. For the diagnostic application, other feature selection methods (Lasso, mutual information (MI) criteria, recursive feature elimination (RFE), SelectKBest) were also applied to narrow down the most informative features that appeared at least twice in one of the methods. The threshold that yielded the highest F1-score was selected. Based on this metric, genes were ranked and the top 16 (diagnostic) and top 13 (prognostic) were considered for downstream analysis. Two additional genes—TP53 and CDH1—were also included in the downstream analysis of both prognostic and diagnostic gene sets.

Algorithmic analysis and evaluation of performance metrics

Gene expression values were used for algorithm training. Several classes of classifiers, with shown applicability to microarray and RNAseq data (Jabeen, Ahmad & Raza, 2018; Peixoto et al., 2023; Pirooznia et al., 2008), such as logistic regression (LR), naive Bayes (NB), K-nearest neighbours (KNN), and support vector machines (SVM) (with linear and radial basis function kernels), were implemented with default hyperparameters in Python programming language (v3.10.0), using the scikit-learn package (v1.0.1). We specifically chose LR (linear, interpretable), NB (probabilistic, fast), KNN (distance-based, non-parametric), SVM (linear margin-based), SVM (RBF kernel, nonlinear margin-based) models as they cover probabilistic, distance-based, linear and non-linear margin-based approaches which are widely used because of good performance on small, imbalanced gene-expression cohorts in the microarray/RNA-seq literature. A leave-one-out cross-validation procedure was used to evaluate the diagnostic or prognostic value of all possible combinations of genes (from n = 2 up to all selected diagnostic or prognostic genes). This involved leaving out one sample at a time for validation while using the remaining samples to create a balanced training set. The Synthetic Minority Oversampling Technique (SMOTE) was employed from the imbalanced-learn (v0.8.1) package. For LR, KNN, and SVM, features were standardized (scaled and centered) using scikit-learn’s standard scaler module by subtracting the mean and scaling to the unit variance. Performance metrics such as accuracy, precision (PPV), recall (sensitivity), NPV, and precision and specificity were calculated and recorded for each full iteration of the validation strategy. The top-performing algorithms were chosen by maximizing performance metrics (accuracy, specificity, precision, recall, NPV, and F1-score, Table 1). The most frequent models, with the highest F1-score, were chosen to further select the best classifiers for both diagnostic and prognostic applications. The most frequently occurring genes (frequency ≥ 50 %) within the selected classifiers were chosen as features. Subsequently, the performance metrics were calculated using a decremental number of features, and the median value and standard deviation of each group of decremental subsets of genes were computed.

In-vitro gene expression analysis

Reverse transcription—quantitative real-time polymerase chain reaction (RT-qPCR)

For 1-Step RT-qPCR, reactions were performed in triplicate, using the TaqPath 1-step RT-qPCR Master Mix (#A15300, Thermo Fisher Scientific, Waltham, MA, USA) with a final reaction volume of 10 µL. Each reaction containing one µL of template, 0.25 µM of probe and 0.5 µM of each primer (Primers and probes used for RT-qPCR are listed in Table S2). Data acquisition and analysis were conducted using the QuantStudio Design & Analysis Software v1.5.1 software, using the cycling program: UNG incubation at 25 °C—2 min, Reverse Transcription at 50 °C—15 min, followed by Polymerase activation at 95 °C—2 min and 40 cycles of Amplification at 95 °C—3 s and 58 °C—30 s. To normalize gene expression levels, the geometric mean of the reference genes (PGK1, ELF1, and RPL13A) was subtracted from cycle threshold (Cq) of the target genes.

Diagnosis and prognosis dataset selection

For the development of the diagnostic application, 13 RNAseq-based datasets were identified, of which only three had publicly available clinical data—GSE58963 (MacCarthy et al., 2014), E-MTAB-4054 (Maag et al., 2017), GSE193946 (Zhang et al., 2022a; Zhang et al., 2022b)—and were therefore included in the present study. The BE data contained in each dataset is represented in Table 2. In total, data from 61 samples—comprising 21 NDBE, 40 LGD BE and 27 HGD—were included in the study.

For the prognostic application, 16 microarray datasets were identified, but only those generated on an Affymetrix platform were included in the downstream analysis to facilitate data merging. A total of nine microarray datasets were analyzed, including three previously analyzed by Cardoso et al. (2016)—GSE1420 (Kimchi et al., 2005), (Watts et al., 2007), and GSE13083 (Stairs et al., 2008) and six new ones, namely GSE36223 (Ostrowski et al., 2007), GSE37200 (Silvers et al., 2010), GSE34619 (Di Pietro et al., 2012), GSE26886 (Wang, Ma & Kemmner, 2013), GSE39491 (Hyland et al., 2014), and GSE100843 (Cummings et al., 2017). In total, data from 200 samples—representing 135 nonP-BE and 65 P-BE—were included in the study as shown in Table 2.

Identification of differentially expressed genes in a diagnostic and prognostic setting

In this study, we aimed to identify diagnostic biomarkers that can distinguish between ND-BE and LGD-BE. For this purpose, we utilized three RNAseq datasets (as listed in Table 2). Low-expression genes were excluded from each dataset, resulting in the inclusion of 20,608 genes in our analysis. After normalization, we conducted differential expression analysis between LGDBE and NDBE (Fig. S1A), and HGDBE and NDBE (Fig. S1B) using EdgeR’s quasi-likelihood approach. This approach accounted for disease staging and batch effects from the three different datasets as factors in the model (Table S3, Fig. S1C).

Following the differential expression analysis, we identified 30 biomarkers through a systematic selection process. First, we selected differentially expressed genes (DEGs) with an absolute log fold change (logFC) of ≥ 1 between LGDBE and NDBE, with a false discovery rate (FDR) of <0.05. From these DEGs, we filtered for genes that showed the same direction of expression change in the HGDBE vs. NDBE comparison (FDR < 0.05), resulting in 14 genes (Fig. S1A). Second, we identified DEGs in the HGDBE vs. NDBE comparison with an absolute logFC of ≥2 (FDR < 0.05). Among these genes, we selected those that also exhibited the same direction of expression change in the LGDBE vs. NDBE comparison (considering p-value < 0.05 for significance), resulting in 16 genes (Fig. S1B). This two-step filtering strategy ensured that the selected biomarkers not only had significant differential expression but also consistent expression patterns across different stages of disease progression. Given the strong batch effect observed (see Fig. S1), there was a risk of losing biologically relevant genes in the LGDBE vs. NDBE comparison due to this variation. To mitigate this problem, we also performed separate analyses of the EMTAB_4054 (Table S4) and GSE58963 (Table S5) datasets. We employed the glmRobust pipeline to independently identify differentially expressed genes between the LGDBE and NDBE groups within each dataset. From these separate analyses, we identified an additional 13 genes with an absolute logFC greater than 1 and an FDR < 0.05. These genes were consistently found in both datasets and exhibited the same direction of expression change (Fig. S2). Moreover, these genes showed consistent directional changes in the previous HGDBE vs. NDBE comparison. Thus, they were also included in the biomarker list (Table S6). Given their established role in the biology of BE and EAC, we also included two additional genes—TP53 and CDH1—in the downstream analysis, resulting in a total of 45 candidate genes for distinguishing between NDBE and LGD.

For the prognostic set of biomarkers, we re-analyzed 25 genes that we had previously identified to have prognostic value (Cardoso et al., 2016), namely ACTN1, C1S, CCN1 (alias CYR61), CDH1, CEBPB, CEBPD, COL4A1, CTSB, DKK3, DUSP1, IER3, JUN, LAMC1, PLPP3, RBPMS, SNAI1, SNAI2, SPARC, TNS1, TP53, TRMT112, TWIST1, VWF, WWTR1 (alias TAZ) and ZEB1. For validation purposes, we added six independent datasets to the three datasets we originally analyzed. We observed significant differential gene expression (adj. p-value < 0.05) between P-BE and nonP-BE categories for most of the genes of interest, except for CDH1, DKK3, SNAI2, and WWTR1.

Application of a thresholding function for the selecting genes with the highest predictive value

To each selected gene, we applied a Thresholding function, to determine a gene expression threshold for distinguishing different levels of gene expression between groups of samples with distinct diagnosis (NDBE vs. LGD-BE) or with distinct prognosis (nonP-BE vs. P-BE). We defined the best individual threshold of gene expression for each selected gene, 45 for diagnosis and 25 for prognosis, reflecting the individual predictive value of each gene. Genes were then ranked by the harmonic mean of recall and precision (F1-score) to ensure accurate selection. This procedure identified the top 15 genes for predicting the malignant progression of BE lesions with F1-score above 0.67 (Table S7). From the top 45 diagnostic genes, including CDH1 and TP53, genes with higher expression values (log2CPM above 1) were filtered. To further refine a list of candidates, we used several feature selection methods: Lasso, MI criteria, RFE, and SelectKBest. Additionally, feature correlation analysis was conducted to identify and eliminate highly correlated features (Pearson’s correlation coefficient > 0.9). Hence, for diagnostic purposes, we further narrowed down the selection to genes that were chosen at least twice in one of the feature selection methods and F1-score above 0.7, which identified the top 16 genes for diagnosing dysplasia in the context of BE (Table S8).

Building on our identification of the top diagnostic and top prognostic genes using the Thresholding function and various feature selection methods, the ROC curve results (Fig. 2) further validate their predictive power using a logistic regression classifier. In both diagnostic (Fig. 2A) and prognostic (Fig. 2B) contexts, most of the genes have area under the curve (AUC) values above 0.50 (random chance line), with some reaching individual values of 0.90, demonstrating a substantial predictive value of the selected genes in both diagnostic and prognostic applications.

SVM algorithms were the best for both diagnostic and prognostic applications

The diagnostic and prognostic gene groups were utilized to train the most effective diagnostic and prognostic algorithms. Various classifiers—logistic regression (LR), naive Bayes (NB), K-nearest neighbors (KNN), linear support vector machine (LSVM), and radial basis function support vector machine (RBF SVM)—were examined using increasing combinations of genes, ranging from n = 2 up to the total number, for diagnostic and prognostic applications.

The algorithms were ranked based on their performance metrics for each application (see Table 1). However, no algorithms optimized all performance metrics for both applications. Nevertheless, the LSVM algorithms emerged as the best for diagnostic purposes, maximizing the F1-score and accuracy (refer to Fig. 3A and Table 1).

For the prognostic application, a similar trend was observed, where the RBF SVM type performed best according to the F1-score and accuracy metrics (refer to Fig. 3B and Table 1).

The study found that among the selected types of algorithms, those with an F1-score above 0.96 included 1,115 LSVM for diagnostic and 5,794 RBF SVM for prognostic. The analysis identified the most frequent genes (over 50 %) across the best-performing algorithm class (Tables S9 and S10). Ultimately, ten genes were selected for identifying LGD BE using a LSVM algorithm: IGHV3-43, SLC38A4, PLLP, CELA3A, IGHV4-31, TMPRSS5, TP53, NR4A1, ATF3, IFI27. For identifying P-BE, ten genes were selected using an RBF SVM algorithm: SNAI1, C1S, DUSP1, CEBPB, COL4A1, ZEB1, CEBPD, CCN1, LAMC1 and TWIST1.

The performance of each selected algorithm (LSVM for diagnosis and RBF SVM for prognosis) was evaluated using the most frequent genes (10 for diagnosis and 10 for prognosis) as features. To test different random states while avoiding algorithm bias, 100 runs were performed for each algorithm with the same features. Table 3 presents the mean values and respective standard deviations (SD) for each performance metric. All performance metrics were above 0.90, except for specificity for the LSVM diagnostic algorithm. The low standard deviations (below 0.05) indicated an increase in the predictive value of each algorithm when the selected genes were combined.

Finally, the performance of the two algorithms was evaluated by gradually decreasing the number of selected genes (Fig. 4). The diagnostic algorithm showed a decrease in performance after the removal of just one gene (Fig. 4A). In contrast, the prognostic algorithm showed noticeable changes only after the removal of four genes (Fig. 4B).

In-vitro validation of key diagnostic and prognostic biomarkers

We conducted a validation study of the panel of biomarkers to distinguish between different stages of BE progression. Specifically, we performed RT-qPCR analysis to compare the expression levels of these biomarkers in differente cell lines: metaplasia (BAR-T and BAR-T10), dysplasia (CP-B, CP-C and CP-D), and EAC (OE33, KYAE-1 and ESO26). Each biomarker was tested with three technical replicates in each cell line.

For dysplasia diagnosis, we analyzed the expression of biomarkers in both metaplasia and dysplasia cell lines (Fig. 5). In evaluating EAC prognosis, we compared the expression levels between metaplasia and EAC cell lines (Fig. 6). Normalized expression values against reference genes (PGK1, ELF1 and RPL13A) highlighted significant differences in key markers. For instance, biomarkers such as IFI27 and ATF3 differentiated metaplasia from dysplasia with statistically significant p-values (p = 0.009 and p = 0.003, respectively), revealing their potential utility in dysplasia diagnosis. Similarly, CEBPB, SNAI1 and CCN1 (alias CYR61) genes showed significant expression changes between metaplasia and EAC (p-values of 0,031, 0.022 and 0.038, respectively). This supports their relevance for EAC prediction. The observed differential expression patterns suggest that these biomarkers serve as valuable molecular tools for early detection of dysplasia and the risk of progression to EAC, facilitating timely clinical intervention.

Interestingly, some of the top-performing genes, namely IGHV3-43, IGHV4-3 1, IGHV3-53, and PGC, showed no detectable expression in cell lines. Since these genes ranked high according to the diagnostic algorithm, we hypothesized that their expression originates from immune cells, typically absent in cell lines. To investigate this further, we specifically tested these genes in tissue samples from BE patients with and without dysplasia. Contrary to cell lines, the expression of these genes was detectable in patient samples, supporting the notion that immune cells- may play a critical role in BE progression (Fig. 5).