Integrating multi-platform gene expression data and machine learning assisted biomarker discovery in colorectal cancer

2.1 Inclusion criteria

We integrate transcriptomic datasets from RNA-Seq and microarray technologies. The search term “Colorectal cancer” was used to retrieve relevant studies from the Gene Expression Omnibus (GEO) database. The dataset was selected using the inclusion criteria listed below:

• a)

  The datasets must originate from Homo sapiens.

• b)

  The datasets must specifically exclude cell line-based investigations and only include total RNA extracted from Colorectal tissue and nearby healthy regions.

• c)

  The two distinct dataset groups must adhere to a standard sample-collecting procedure.

• d)

  Samples having mutations, induced gene expression, therapeutic interventions, medication treatments, or gene knockdowns shouldn’t be included in the datasets.

By following the above-mentioned selection criteria, we choose five microarray datasets—GSE25070, GSE37182, GSE41328, GSE74602, and GSE8671. The data platform for GSE25070, GSE37182, and GSE74602 datasets were Illumina HumanRef-8 v3.0 expression bead chip, while GSE41328 and GSE8671 studies involved Affymetrix Human Genome U133 Plus 2.0 Array. Similarly, two studies for RNA-Seq data selected for the study include, GSE142279 and GSE50760. A detailed representation of individual datasets, along with the sample numbers, is represented in Table 1.

2.2 Data processing

We used GEO query (Davis and Meltzer 2007) R package from Bioconductor to acquire the different microarray technologies and platforms from the GEO repository. Data processing makes an important impact since it turns the raw data matrix into a clean dataset, which enhances the quality of information. Standard transcriptomic data preprocessing steps include imputation of missing data, filtering, transformations, and sample-based normalization. The accuracy, precision, and robustness of analyses can be greatly increased by using the right tools, techniques, and datasets. We utilized the arrayQualityMetrics R package that makes use of six quality tests to eliminate low-quality samples and outliers from microarray datasets. The duplicate samples from the RNA-seq datasets were removed. The Robust multi-array average (RMA) algorithm normalized the highly variable data due to the different platforms used in our study (Bolstad et al., 2003). It normalizes the microarray data by quantiles, adjusts the background, and summarizes the data. The rma function from R package “affy” (Gautier et al., 2004) and “oligo” (Carvalho et al., 2013) was used for pre-processing the data generated from Affymetrix microarrays while data from Illumina platforms were processed using “lumi” (Du et al., 2008) package. The annotate R package (Gentleman et al., 2013), was used for gene annotations and conversion of EntrezIDs to official gene symbols. Alternatively, the biomaRt R package was used to annotate the Ensemble IDs and convert them into gene symbols once the expression raw count matrix file from RNA seq data was obtained from its source repository (Smedley et al., 2015; Howe et al., 2021). The cqn package corrected the GC content bias and normalized the data (Hansen et al., 2012), while the NOISeq package was used to compute the gene expression values (Tarazona et al., 2015).

2.3 Gene expression integration

The average and median statistical summaries can summarize all the genes with the same identifier, which is essential for a consistent data analysis. Moreover, log transformation was applied to datasets that have not been scaled to log2 depth. The transcripts were pooled by taking the average of each series separately.

The depth of each microarray platform was examined with consistent representation of multiple datasets based on platform IDs. Considering the images from microarray scanners are digitized at 16 bit-depth with known saturation limits, the RNA-seq-derived expression values can cover a wider numeric range. An unsupervised monotonic linear rescaling to match the dynamic ranges of each dataset was applied. This method preserved rank information and fold-change ordering, but only as a unit change, and was done before quantile normalization and ComBat batch effects correction. This approach aligns the intensity levels across platforms and provides point estimates for microarray intensities aligned to RNA-seq distributions. A dynamic-range harmonization (bit-depth–aware scaling) was applied depending on the maximum value of each series. This step rescales each dataset to a common [0,1] range. It does not change any value, affect ranks within a dataset, or use outcome labels. A similar approach has been used before while integrating heterogeneous gene expression data (Castillo et al., 2017; Galvez et al., 2019). This led to a bit depth between platforms: 16-bit for Human Genome U133 Plus 2.0 Array (GSE41328 and GSE8671), 16-bit for the HumanAll platform (GSE25070, GSE37182, and GSE74602), 20-bit for the HiSeq 2000 (GSE50760), and 24-bit for the HiSeq X Ten (GSE142279). Each platform selected a specific bit depth to standardize their sample data. Data integrity was maintained by selecting only genes that appeared in all datasets for the integration process. Batch effect correction was performed through the ComBat method from the sva R package after filtering common genes to create a uniform distribution of samples across all batches. The inter-array normalization of expression values across samples used the limma R package normalizeBetweenArrays function to ensure sample comparison and avoid classification errors in machine learning later. The process ensures every sample matches the others while quantile normalization applies the same empirical distribution to all samples. The integrated dataset was formed by merging n samples from N classes with p common genes.

2.4 Feature selection

Feature selection is a crucial step in genomic datasets, where it streamlines the classification process by eliminating superfluous features. Feature selection reduces the workload of the classifier and, as a result, increases the classification accuracy by identifying the most important characteristics for high-dimensional datasets. We used LASSO regression to decrease the number of genes or features that helped in further analyzing the data. The training dataset with 10-fold CV was subjected to the LASSO model using the R package glmnet version 4.1-3. The L1 penalty term’s inherent property allows LASSO to precisely reduce some of the estimated coefficients to zero. Consequently, characteristics with non-zero regression coefficients can be automatically chosen as informative (Tibshirani 1996). The data were randomly divided into ten sets; nine were used for training, and the remaining set was used for testing. Initially, using the cv.glmnet package, 10-fold CV was used to compute the penalty regularization parameter lambda (Ahmed et al., 2022). The final model was constructed using glmnet, with the lambda value corresponding to the maximum area under the curve (AUC, lambda. min) selected for regularization. The expression data of genes with non-zero coefficients were used to create the final LASSO model.

2.5 GO-KEGG pathway enrichment analysis

Kyoto encyclopedia of genes and genomes (KEGG) and Gene ontology (GO) pathway enrichment analysis were performed to investigate the biological pathways and functions of the genes obtained in the previous step. Functional analysis was conducted by database for annotation, visualization, and integrated discovery (DAVID), with the KEGG pathways and GO terms significantly enriched, explained by P < 0.05 after an adjusted p-value cut-off (Kanehisa and Goto 2000).

2.6 Classification through traditional machine learning classifiers

We employed five classification models to evaluate the performance of Lasso regularized regression, namely random forest (RF) (Rigatti 2017), support vector machines (SVM) (Suthaharan and Suthaharan 2016), eXtreme Gradient Boosting (XGBoost) (Chen et al., 2015), naive bayes (NB) (Webb et al., 2010), and k-nearest-neighbors (KNN) (Guo et al.,). We applied these ML models using the R environment with caret (Kuhn 2011), xgboost (Chen et al., 2015), pROC (Robin et al., 2021), e1071 (Meyer et al., 2019), and random forest (RColorBrewer and Liaw, 2018) packages. To ensure reproducibility, we fixed random seeds to 123. A stratified 80:20 split of the integrated dataset into training and test sets, with class proportions preserved in each split. The test set was held out and used only once for final evaluation. A 10-fold cross-validation (K-fold CV, K = 10) was applied to the training set to tune the hyperparameters for each classification approach. Hyperparameters were chosen based on the mean AUROC (area under the receiver operating characteristic curve) from the inner 10-fold CV. The 1-SE rule was used to select the simplest setup that stayed within one standard error of the best option. Specifically, the parameters that were modified include: σ (Laplace) and usekernel (Gaussian Naïve Bayes) in Naïve Bayes, the booster, objective, nrounds, max_depth, eta, gamma, colsample_bytree, min_child_weight, and subsample in XGBoost, number of trees for RF, k for k-NN, and γ for support vector machines (SVM) and random forest (RF). The objective was to select the best values for each procedure to improve model performance on unseen data. The classification models were developed using the selected genes and the training dataset and subsequently evaluated using the test dataset. We made pairwise comparisons of model performance on the testing dataset using DeLong’s test for AUROC and McNemar’s test for paired accuracy. We used the fixed thresholds from the training data to calculate F1 performance via bootstrap 95% confidence intervals (B=2000). Holm adjustment was used in the case of more than one pairwise comparison.

2.7 Deep learning models

We used two deep-learning models, i) 1-dimensional convolutional neural network (1D-CNN), ii) multilayer perceptron (MLP), to train our integrated data and extract important features. For deep learning models, we used features obtained from LASSO regression. The deep neural network models were trained using the Adam optimizer. The learning rate was selected via a small grid using a training-only internal validation split. To minimize overfitting and improve generalization, L2 weight decay, dropout, and early stopping on validation loss parameters were employed. The binary classification uses a softmax activation function in its output layer. The testing dataset was not changed until the final evaluation was done. The selected features were then used for classification by deep learning models. We load the reticulate library (v1.41) and configure R to use a Python installation. The deep neural network TensorFlow (v2.16) (Developers 2022) was used in combination with the Keras (v2.1.5) R package (Arnold 2017) to design a neural network.

2.8 1-Dimensional convolutional neural network

A convolutional neural network (CNN) is a type of deep learning model that automatically learns spatial levels of features through convolutional layers. A standard CNN architecture consists of convolutional and pooling layers, as well as regularization mechanisms such as batch normalization and dropout, to enhance generalization. First, the gene expression data were normalized, and then replaced into the three-dimensional tensor organization, ready to be processed serially by a 1D-CNN model. The CNN network consists of two convolutional layers, each with 64 and 128 filters, respectively. Both layers utilize ReLU activation and a kernel size of 3. The resulting tensor is then flattened, passed through a dense layer (64 units with the ReLU activation function) (LeCun et al., 2015). The model was trained in 100 epochs with a 0.0005 learning rate before early stopping prevented overfitting. Multiple evaluation metrics, including accuracy alongside sensitivity and specificity, F1-score, and area under the receiver operating characteristic curve (AUC), were used to assess the model’s performance.

2.9 Multilayer perceptron

MLP deep neural network contains three dense layers with 128, 64, and 32 hidden units that use ReLU activation functions. The model was trained at 100 epochs with 128 batches per epoch and 20% validation data (Popescu et al., 2009). The first convolutional layer weights are converted into absolute values to determine the feature importance. Each feature receives an importance score through the total of absolute weights across all filters. The features were ranked from highest to lowest score, with the top-ranked feature receiving the highest score and vice versa. The 25 leading features appeared in bar plots and were stored for further analysis.

2.10 Model validation using TCGA GDC dataset

We performed validation of the models using transcriptomic RNA-Seq data of CRC samples from the Cancer genome atlas (TCGA) Genomic data commons (GDC) portal (https://portal.gdc.cancer.gov/). The TCGAbiolinks R package (Colaprico et al., 2016) was used to download the data, and the Ensemble IDs were mapped to HGNC gene symbols via biomaRt package (Durinck et al., 2009). The validation on TCGA was performed on the original, unaltered data without balancing the data. Sample types were identified using short letter codes and labeled as Tumor or Normal. We applied a variance stabilizing transformation (VST) (Kelmansky et al., 2013) with DESeq2 (Love et al., 2014) to normalize the count data.

2.11 Common features across all models

To analyze intersections of the top 20 Features extracted via five ML and 2 Deep Learning models, we utilized the recently developed UpSet analysis (Lex et al., 2014) as visualization of intersections of multiple sets is not possible by classic Euler or Venn diagrams. The UpSet could analyze 2¹², i.e., 4096 intersections.

2.12 Survival analysis

The gene expression profiling interactive analysis 2 (GEPIA2) web server was used to assess the prognostic significance of the significant genes retained after the final stage of data filtering (Tang et al., 2019). The platform’s default settings were implemented to see the prognostic significance of the genes across the TCGA Colon adenocarcinoma (COAD) dataset. Significant features retrieved using artificial intelligence (AI) models and had prognostic value (Log-rank test P-value < 0.05) for either disease-free survival (DFS) or overall survival (OS) were chosen to identify COAD prognostic genes. Additionally, GEPIA2 compared the expression levels were compared between tumor and normal tissue samples and assessed the dysregulation of COAD prognostic genes.

3.1 Feature selection

The curse of dimensionality problem was addressed by the gene selection process from the expression data, which includes just n samples out of thousands of p genes. We reduced the dimension by removing non-essential genes and retaining the essential genes by utilizing an effective feature selection method. In our analysis, we utilized LASSO regression for our feature selection step. Prior to the start, our integrated normalized dataset contains 7880 features or genes and 416 samples, with 213 Normal samples and 203 Tumor samples. LASSO regression screening identified a total of 42 gene characteristics of CRC (Fig. 2a). The best Lambda for our model was 0.005, and the minimum square error (MSE) value was 0.234 (Fig. 2b) (Table S1).

Close

Fig. 2.

LASSO model fitting and selection of optimal regularization parameter. (a) coefficient profiles plotted as a function of the logarithm of the regularization parameter (log λ). each curve represents the trajectory of a predictor’s coefficient across varying λ values. The first vertical dashed line corresponds to the λ value that minimizes the cross-validation error, while the second represents the largest λ within one standard error of the minimum, offering a more regularized and parsimonious model. (b) ten-fold cross-validation curve showing binomial deviance (y-axis) versus log(λ) (x-axis). red dots represent the mean cross-validation error, with error bars denoting ±1 standard error. The two vertical dashed lines correspond to the same λ values as in (a), guiding model sparsity and selection.

Export to PPT

3.2 GO and KEGG enrichment analysis of core genes

The biological functions and associated pathways of genes identified via LASSO regression were determined via GO and KEGG enrichment analyses. The analyses were categorized into three parts, which include Biological processes (BP), cellular components (CC), and molecular functions (MF). The top five enriched terms from each category were summarized in Table 2. In the BP category, core genes were primarily associated with chemotaxis, lipid transport, and general transport processes. For CC, the most significant terms were related to the secretion of proteins into the extracellular environment. In the MF category, the top enriched terms included heparin binding, cytokine activity, and monooxygenase activity.

Along with the GO keywords, a total of four significant KEGG pathways (p-value < 0.05) were identified. The most notable pathways were cytokine-cytokine receptor interaction, viral protein interaction with cytokines and cytokine receptors, ABC transporters, and folate transport and metabolism.

3.3 Classification results through traditional machine learning models

The features identified through the Lasso model were given as data input to various traditional ML models to assess the performance of these genes when new samples were presented based on classification accuracy. We used five different classifiers RF, SVM, XGBoost, Naïve Bayes, and kNN. The optimum hyperparameters used in our models are: mtry = 1 for RF; SVM-Kernel = radial, cost = 0.50 and sigma = 0.03 with 135 support vectors; nrounds = 500, max_depth = 3, eta = 0.05, gamma = 0, colsample_bytree = 0.6, min_child_weight = 1 and subsample = 0.5 parameters for XGBoost; laplace = 0, use kernel = TRUE and adjust =1 for Naïve Bayes and k = 11 for kNN. The ML models developed were evaluated based on accuracy, precision, recall, and AUC. The results showed that Random Forest and XGBoost with the highest diagnostic accuracy (98.78%) followed by kNN (97.56%), SVM (96.56%), and Naïve Bayes (95.36%). Table 3 displays the validation (10-fold cross-validation on the training dataset) and the test results obtained from the use of these five algorithms.

The AUC evaluates the performance of a binary classification model. The receiver operating characteristic - area under curve (ROC-AUC) curve for all the models comes in the range of >90%, representing excellent models (Fig. 3). Each model extracts important features/genes and ranks them according to their importance, which is represented and shown graphically in Fig. 4. For SVM and k-NN, we used the recursive feature elimination (RFE) method to select important features. RFE is a wrapper-based feature selection method that begins by considering all features. It repeatedly fits the model, ranks the features according to their importance, and then prunes away the least significant ones. This process continues until a preset number of features remain or until model performance no longer improves. The RF model used Mean Decrease in Gini (Gini Importance) method, which measures how much a feature contributes to reducing impurity (Gini Index) across all decision trees in the RF model. The XG-Boost used Gain-based feature importance method that measures how much a feature improves the purity of splits in the decision trees. The Naïve bayes used mean differences in conditional probabilities method that assumes conditional independence of features where each feature contributes independently to classification based on how different its distributions (mean values) are across classes.

Close

Fig. 3.

Receiver operating characteristic (ROC) curves of five classifiers; SVM, XGBoost, Naïve Bayes, k-NN, and RF, evaluated using 10-fold cross-validation. The AUC values for all models exceed 0.97, indicating excellent discriminatory power. Random Forest achieved the highest AUC (0.999), followed closely by Naïve Bayes (0.988), XGBoost (0.989), k-NN (0.987), and SVM (0.978).

Export to PPT

Close

Fig. 4.

Important features extracted from five different machine learning models. (a) The top important features are extracted via the random forest model using the ini importance method. (b) The most important features are extracted via an svm model using the recursive feature elimination method. (c) The most important features extracted via Naïve Bayes model using mean differences in conditional probabilities. (d) The top important features are extracted via XGBoost model using the gain-based feature importance method. (e) The top significant features extracted via kNN model using the recursive feature elimination method. (f) The top important features extracted via 1D-CNN using the absolute weights of the first convolutional layer. The importance score for each gene (feature) was computed as the sum of the absolute values of its weights across all convolutional filters. (g) The most significant genes extracted via MLP.

Export to PPT

The RF model extracts 42 genes with the most important features. The SVM model extracts 15 genes with the most important features. The XGBoost model selected 42 genes with the most significant features. The Naïve Bayes model obtains 42 features as the most valuable features. Finally, the k-NN model acquires 33 genes as the most significant features.

Comparisons across five classifiers showed small pairwise differences and no statistically significant differences (all adjusted DeLong p≥0.99 and all adjusted McNemar p≥0.99) after Holm correction. The largest absolute AUROC difference was approximately 0.010, and the largest absolute accuracy and F1 differences were approximately 0.024. Point estimates suggested Naïve Bayes and random forest had slightly higher AUROCs than the SVM, and XGBoost and k-NN had slightly lower AUROCs than the SVM, but these differences were statistically insignificant (Table S2).

3.4 Classification results through deep learning models

We also employed two deep-learning neural network models (MLP and 1D-CNN). A regression-based prediction model called the MLP, along with a 1D-CNN, was put forth to classify gene expression data. The ID-CNN and MLP models, when trained on LASSO regression feature-extracted data, give an accuracy of 98.81% and 97.62% respectively (Fig. 5). 1D-CNN shows 98.79% F1 score, 97.67% Sensitivity, 100% Specificity, and 99.21% AUC. Similarly, the results of MLP show 97.67% F1 score, 97.67% Sensitivity, 97.56% Specificity, and 99.32% AUC (Table 3).

Close

Fig. 5.

Training and validation performance of (a) 1D-CNN and (b) MLP models. Both models show decreasing training loss and stable validation loss, indicating effective learning and convergence. Accuracy curves reveal strong generalization. The use of LASSO regularization promotes model sparsity and may enhance feature selection without compromising performance.

Export to PPT

3.5 Validation of models

The diagnostic validity of the seven Models (5 traditional machine learning and 2 deep learning) constructed was validated using an external dataset retrieved from TCGA database. The TCGA-colorectal adenocarcinoma (COAD) dataset contains 481 Tumor samples and 43 Normal samples. The DESeq2 variance stabilizing transformation was applied to the raw counts to achieve normalization of the gene expression data. All the models achieved an accuracy of greater than 90% showing the model’s accuracy and reliability (Table S2).

3.6 Upset analysis to identify overlapping genes

We performed the UpSet analysis to visualize commonly regulated genes extracted by five different ML models and two deep learning models. The top features were selected from each model. The results from UpSet analysis revealed that seven genes, Carbonic Anhydrase 7 (CA7), ATP binding cassette subfamily A member 8 (ABCA8), somatostatin (SST), myomesin 1 (MYOM1), CC motif chemokine ligand 23 (CCL23), procollagen C-endopeptidase enhancer 2 (PCOLCE2), and CXC motif chemokine ligand 10 (CXCL10), as overlapping features present across all models (Fig. 6).

Close

Fig. 6.

UpSet plot illustrates the overlap of selected features among five classification algorithms. The plot shows the intersection sizes of features selected by different classifiers. Vertical bars represent the number of shared features across different combinations of algorithms, while the matrix below indicates which algorithms contribute to each intersection. The horizontal bars on the left indicate the total number of features selected by each algorithm.

Export to PPT

3.7 Prognostic analysis of overlapping features

We compared the expression level of these hub genes in the COAD dataset using GEPIA2 database, comparing 275 tumor and 349 normal tissues. The results showed that the expression level of all essential genes was significantly downregulated for PCOLCE2, CA7, ABCA8, SST, MYOM1, and CCL23 (Figs. 7a-f) and significantly upregulated for CXCL10 (Fig. 7g)

Close

Fig. 7.

The box plot for hub genes expression level of PCOLCE2, CA7, ABCA8, SST, MYOM1, CCL23, and CXCL10 (a-g) in tumors indicated by red block and control samples showed by grey block. The red block shows the expression level of genes in tumor samples, which is significantly differentiated from normal samples. *P < 0.05. (h) Prognostic analysis of PCOLCE2 in colon adenocarcinoma using the GEPIA2 portal for overall survival. The red line denotes high gene expression samples, whereas the blue line signifies low gene expression samples. The research utilized the log-rank test to establish statistical significance while setting the threshold at P < 0. 0085. The X-axis shows time intervals while the Y-axis records patient survival percentages.

Export to PPT

Close

Fig. 8.

The schematic framework of the proposed workflow.

Export to PPT

We also perform survival analysis to analyze the role of significant features and the progression of colon adenocarcinoma. The key genes identified were analyzed for overall survival and disease-free survival to find parameters that were linked with prognostic relevance. The survival analysis showed that PCOLCE2 demonstrated a significant association with disease-free survival, with a log-rank p-value of 0.0085 (Fig. 7h). This suggests that PCOLCE2 expression levels are significantly correlated with disease recurrence in prostate cancer patients. The hazard ratio (HR) of 1.9 indicates that patients with high PCOLCE2 expression have nearly double the risk of disease progression compared to those with low expression. The survival curve for the high-expression group remains consistently below that of the low-expression group, further supporting the observation that elevated PCOLCE2 expression is associated with worse prognosis and a higher likelihood of disease recurrence.