ML Classification of Cancer Types Using High Dimensional Gene Expression Microarray Data

1. Introduction

Cancer is a disease which can start almost anywhere in the human body, in which some of the body’s trillion cells grow uncontrollably and spread to other parts of the body. There are over 200 types of cancer such as colon, liver, ovarian and breast, and so on [1,2]. In 2023, 1,958,310 new cancer cases and 609,820 cancer deaths were projected in the United States. [3]. This prompts a clear understanding of the underlying mechanism and characteristics of this potentially fatal disease alongside identifying the most significant genes responsible for it.

Microarray data analysis has been a popular approach for diagnosing cancer, and DNA microarray is a technology used to monitor large numbers of various gene expressions at the same time [4,5]. Gene expression analysis can assure medical experts whether a patient suffers from cancer within a relatively shorter time than traditional methods. Analysis of gene expression requires the identification of informative genes [6] and whereas [7] demonstrates that gene expression classification or cancer classification is the process of identifying informative genes that can be used to predict new sample classes.

Gene selection, however, is generally regarded as much more problematic in multi-class situations (where there are three or more classes to be differentiated) [8,9]. The microarray genes expression data constitutes many highly correlated genes for just a small sample size and have high levels of noise in them. The small number of cancer samples compared with the number of features can degrade the performance of the classifier and increase the risk of over-fitting. Machine Learning (ML) techniques are used as an aim to model the progression rate and treatment of cancer patients. Moreover, ML-based classifiers are widely used in classification of cancer sub-types Various supervised and unsupervised ML techniques which have been adopted to identify the most significant genes [10,11]. Studies [10,11,12,13] show that these techniques suffer from overfitting and multi-collinearity problems due to noise, large number of genes, and small sample size. The unsupervised learning algorithms such as the hierarchical clustering [14], [15], K-means clustering [16] etc., have been used to identify genes which are responsible for cancer. These techniques did not identify the most significant genes resulting in low classification accuracies. Hence it may be beneficial to use feature selection methods which can address the challenges arising from high data dimensionality and small sample size.

Over the years, several studies have been carried out based on feature selection. Principal Components Analysis (PCA) is an exploratory multivariate statistical technique for simplifying complex data sets [17,18,19]. It has been used in a wide range of biomedical problems, including the analysis of microarray data in search of outlier genes [20] as well as the analysis of other types of expression data [21,22].

Degroeve, De Baets, Van de Peer and Rouz´e (2002) created a balanced training data set by randomly selecting 1000 positive instances and 1000 negative, and also created a test data with 281 positive and 7505 negative instances and another test data set with 281 positive and 7643 negative instances; they used SVM classifier, a naive Bayes classifier and a traditional method for feature selection for predicting splice site and obtained improved performance. Precision obtained for these datasets ranged in 93-98% range, but the recall and F1-measures ranged in 25-49% range [23]. Peng, Li and Liu (2006) compared various methods of gene selection over four microarray gene expression datasets and showed that the hybrid method works well on the four datasets [24].

Sharma and Paliwal (2008) used Gradient LDA method for three small microarray gene expression datasets: acute leukemia, small round blue-cell tumor (SRBCT) and lung adenocarcinoma and have obtained higher accuracies than some competing methods [25]. Bar-Joseph, Gitter and Simon (2012) provided a discussion of how time-series gene expression data is used for identification of activated genes in biological processes, and also describe how basic patterns lead to gene expression programs [26]. Dwivedi (2018) used the method of Artificial Neural Network (ANN) for classification of acute cases of lymphoblastic leukemia and myeloid leukemia and reported over 98% overall classification accuracy [27]. Sun et al. (2019) used the genome deep learning method to analyze 6083 samples from the Whole Exon Sequencing mutations with 12 types of cancer and 1991 non-cancerous samples from the 1000 Genome Project and obtained overall classification accuracies ranging in 70% - 97% [28]. A survey of feature selection literature for gene expression microarray data analysis based on a total of 132 research articles [29] was conducted by Alhenawi, Al-Sayyed, Hudaib and Mirjalili (2022). Khatun et al. (2023) developed an ensemble rank-based feature selection method (EFSM) and a weighted average voting scheme to overcome the problems posed by high dimensionality of microarray gene expression data [30]. They obtained overall classification accuracies of 100% (leukemia), 95% (colon cancer), and 94.3% for the 11-tumor dataset. Osama, Shaban and Ali (2023) have provided a review of ML methods for cancer classification of microarray gene expression data; data pre-processing and feature selection methods including filter, wrapper, embedded, ensemble, and hybrid algorithms [31].

Kabir et al. (2023) compared two different dimension reduction techniques—PCA, and autoencoders for the selection of features in a prostate cancer classification analysis. Two machine learning methods—neural networks and SVM—were further used for classification. The study showed that the classifiers performed better on the reduced dataset [32]. Another study Adiwijaya et al. (2018) utilized Principal Component Analysis (PCA) dimension reduction method that includes the calculation of variance proportion for eigenvector selection followed by the classification methods, a Support Vector Machine (SVM) and Levenberg-Marquardt Backpropagation (LMBP) algorithm. Based on the tests performed, the classification method using LMBP was more stable than SVM [7].

Based on previous research, the general scheme in the process of classification of microarray data for the detection of proposed cancer can be conducted via preprocessing the data and dimensionality reduction followed by gene classification. In this study, the step for dimensionality reduction was performed using a Principal Component Analysis (PCA) followed by ML classification technique RF and LDA.

2. Materials and Methods

In this article, we have used the Linear Discriminant Analysis (LDA) classifier [33] and the random forest (RF) classifier [34] on an 801 rows x 20531 columns (genes) dataset of patients with five cancer types: BRCA, KIRC, COAD, LUAD and PRAD; the dataset has no missing values. Variables in this dataset are RNA-Seq gene expression levels measured by illumina HiSeq platform. The variables are dummy named gene XX. This dataset (gene expression cancer RNA-Seq) was downloaded from the UCI Machine Learning Repository [35]. The statistical software package R (2023) was used for all data analyses and visualizations [36].

The LDA and RF methods were first attempted on the raw 801x20531 dataset, but the R code produced stack overflow error due to very large number of columns. We then computed Principal Component (PC) scores [37] of the data and performed the 5-level classification on an increasing number of PC's and obtained excellent classification results using just the first two components PC1 and PC2.

BRCA (Breast Cancer gene) genes produce proteins that help repair damaged DNA and are referred to as tumor suppressor genes since certain changes in these genes can cause cancer [38]. People born with a certain variant of BRCA tend to develop cancer at early ages. Chang, Dalpatadu, Phanord and Singh [39] fitted a Bayesian Logistic Regression model for prediction of breast cancer using the Wisconsin Diagnosis Breast Cancer (WDBC) data set which was downloaded from the UCI Machine Learning Repository; precision, recall and F1-measures of 0.93, 0.89, and 0.91 were reported for the training data, and 0.87, 0.91, 0.89 for the test data, respectively.HER2 protein accelerates breast cancer cell growth and HER2 positive patients when treated with medicines which attack the HER2 protein. Gene expression patterns of HER2 are quite complex and pose a challenge to pathologists. Cordova et al. (2023) developed a new interpretable ML method in immunohistochemistry for accurate HER2 classification and obtained high precision (0.97) and high accuracy (0.89) using immunohistochemistry (IHC) and fluorescence in situ hybridization (FISH) data [40].

Kidney renal cell carcinoma (KIRC) is the most prevalent type of kidney cancer, with survival rate of less than 5 years and 338,000 estimated number of new cases each year [41]. ICD profile of KIRC. Wang et al. (2023) correlated the immunogenic cell death (ICD) of KIRK with the heterogeneity and therapeutic complexity which is useful for developing optimal immunotherapy strategy for KIRC patients [42].

A common cancerous tumor in the digestive track is Colon adenocarcinoma (COAD) and is commonly associated with fatty acids [43]; diagnosis of COAD is difficult as there are hardly any early symptoms. Li et al. (2017) used a genetic algorithm and the k-nearest neighbors clustering method to determine genes which can accurately classify samples as well as class subtypes for a TCGA RNA-seq dataset of 9066 cancer patients and 602 normal samples [44].

Lung adenocarcinoma (LUAD) is a common form of lung cancer which also gets detected in the middle/late stages and therefore is hard to treat [45]. Yang et al. (2022) used a dataset of gene expression profiles from 515 tumor samples and 59 normal tissues and split the dataset into two significantly different clusters; they further showed that using age, gender, pathological stages, and risk score as predictors of LUAD increased the prediction accuracy measures [46]. Liu, Lei, Zhang, and Wang (2022) used cluster analysis on enrichment scores of 12 stemness signatures to identify three LUAD subtypes, St-H, St-M and St-L for six different datasets [47].

Prostate adenocarcinoma (PRAD) is common in elderly men, and patients suffering from PRAD typically have good prognosis [48]. Khosravi et al. (2021) used Deep Learning ML models on an MRI dataset from 400 subjects with suspected prostate cancer combined with histological data and reported high accuracies [49].

We will next provide brief descriptions of the methods of data analysis and the common measures of accuracy used in multi-level classification.

2.1. Principal Components Analysis (PCA)

PCA is a dimension-reduction technique which creates new and uncorrelated linear combinations of original variables (principal components); the values of the principal components are called PC-Scores and can be used in place of the original variables for further analyses such as Multiple Linear Regression (MLR) or Discrimination and Classification. Using PC-Scores instead of original variables as predictors eliminates the problem of multicollinearity.

2.2. Linear Discriminant Analysis (LDA)

LDA is itself a dimension-reduction technique which is used for separating a dataset into 2 or more subgroups, and also for classification of new data into these subgroups. LDA is typically one of the methods used for multi-level classification problems. The LDA method involves computing separating hyperplanes for classification purposes [33]. The R-code for LDA produced stack overflow errors, and we therefore used the function prcompfast of the R-package Morpho to perform PCA of the gene expression microarray dataset at hand.

2.3. Random Forest (RF)

The RF method is a decision-tree based method that can be used for classification (categorical response) or regression (continuous response) problems. It randomly selects a subset of rows (samples) and a subset of columns (features) at a time and fits decision trees a very large number of times to predict Y and then uses a voting mechanism to predict Y values.. Random forest is known to be highly accurate [50].

2.4. Training and Test Datasets

In ML literature, it is common practice to randomly split the available dataset into Training and Test datasets and report the accuracy measures of prediction for both of these datasets. Typically higher accuracy measures are obtained for the training set than the test set. The entire raw dataset was used to compute PC-Scores by using the fast-PCA method of the R-package Morpho. A dataset of 801 rows and 25 PC-scores was created, and then this dataset of PC-scores was randomly split into an 80% training set and 20% test set. The LDA and RF methods were used on the training set of PC-Scores and the accuracy measures given below were computed for both training and test sets.

2.5. Accuracy Measures for Multi-Level Classification

All of these measures are computed from the confusion matrix which is a cross-tabulation of observed Y and predicted Y values.

The commonly used accuracy measures for multi-level classifiers computed from the full confusion matrix (CM) are:

Overall Accuracy (OA) = sum of diagonal elements of CM/sum of all elements of CM

Precision, Recall and F1 are calculated for each level j (j = 1, 2, …, 5)

Precision_j = j-th diagonal element of CM/sum of j-th column of CM

Recall_j = j-th diagonal element of CM/sum of j-th row of CM

F1_j = harmonic mean of Precision_j and Recall_j

The following accuracy measures are computed for each level by calculating the one vs all binary confusion matrices:

Area Under the Curve (AUC)

Macro- and micro-averages of AUC

Explanations of the accuracy measures and computational details are provided in [51].

3. Results

PCA was run on the entire 801 rows x 20531 genes data set, and trial-and-error showed that just the first two principal components were sufficient for classification purposes. The genes with highest absolute loadings are shown in Table 9. The PC1 and PC2 scores were saved in a data file. A scatterplot of first two PC-scores for the entire dataset is shown in Figure 1. A clear separation between BRCA and KIRC cancer sub-types with some overlap between COAD, LUAD and PRAD is seen in Figure 1.

Figure 2, Figure 3, Figure 4 and Figure 5 show plots of the confusion matrices for the LDA and RF classifiers for training and test sets, respectively. Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, Table 7 and Table 8 show that all measures of multi-level accuracy are high for both training and test datasets and both LDA and RF methods.

Accuracy Measures for the LDA Classifier for Training Data

Figure 2. Confusion Matrix Plot for the LDA Classifier – Training Data.

Figure 2. Confusion Matrix Plot for the LDA Classifier – Training Data.

Table 1. Precision, Recall, F1 and AUC Measures for the LDA Classifier – Training Data.

Table 1. Precision, Recall, F1 and AUC Measures for the LDA Classifier – Training Data.

Table 2. Macro- and Micro-Averaged AUC Measures for the LDA Classifier – Training Data.

Table 2. Macro- and Micro-Averaged AUC Measures for the LDA Classifier – Training Data.

Accuracy Measures for the LDA Classifier for Test Data

Figure 3. Confusion Matrix Plot for the LDA Classifier – Test Data.

Figure 3. Confusion Matrix Plot for the LDA Classifier – Test Data.

Table 3. Confusion Matrix Plot and Accuracy Measures for the LDA Classifier – Test Data.

Table 3. Confusion Matrix Plot and Accuracy Measures for the LDA Classifier – Test Data.

Table 4. Macro and Micro averaged AUC for the LDA Classifier – Test Data.

Table 4. Macro and Micro averaged AUC for the LDA Classifier – Test Data.

Accuracy Measures for the RF Classifier for Training Data

Figure 4. Confusion Matrix Plot for the RF Classifier – Training Data.

Figure 4. Confusion Matrix Plot for the RF Classifier – Training Data.

Table 5. Precision, Recall, F1 and AUC Measures for the RF Classifier – Training Data.

Table 5. Precision, Recall, F1 and AUC Measures for the RF Classifier – Training Data.

Table 6. Macro and Micro averaged AUC for the RF Classifier – Training Data.

Table 6. Macro and Micro averaged AUC for the RF Classifier – Training Data.

Accuracy Measures for the RF Classifier for Test Data

Table 7. Precision, Recall, F1 and AUC Measures for the RF Classifier – Test Data.

Table 7. Precision, Recall, F1 and AUC Measures for the RF Classifier – Test Data.

Table 8. Macro and Micro averaged AUC for the RF Classifier – Test Data.

Table 8. Macro and Micro averaged AUC for the RF Classifier – Test Data.

In Table 9 we provide the variables (genes) with high absolute loadings on the first two PC-scores; such a table can be very useful for selection of features (genes).

Table 9. genes wit highest absolute loadings on the first two PC-scores.

Table 9. genes wit highest absolute loadings on the first two PC-scores.

4. Discussion

PCA results showed that the first 50 components (PC) cumulatively explained 71% of all variability present in the 801 × 20532 gene expression data, with the first two PC’s explaining only 26% of total variability. The first two PC’s, however, were sufficient for classification of cancer sub-types with high accuracy. This can be seen from the plot of the first two components by of cancer sub-type. LDA was able to classify each of the five cancer-subtype with high accuracies with the exception of LUAD which had a precision of 77% for the training set. The RF method was able to classify each sub-type with very high accuracy. The PCA loadings on 20532 genes were sorted in order of magnitude and genes (features) important for classification were identified. Our results are not generalizable, but the proposed classification method should be very helpful to researchers and clinicians working with gene expression microarray data of high dimensionality.