This website contains supplementary material to the paper:

A. Ramírez-Mena, E. Andrés-León, M. Jesus Alvarez-Cubero, A. Anguita-Ruiz, L.J. Martinez-Gonzalez, J. Alcala-Fdez. Explainable artificial intelligence to predict and identify prostate cancer tissue by gene expression. Computer Methods and Programs in Biomedicine 240 (2023) 107719. doi: 10.1016/j.cmpb.2023.107719

Summary:

• Differential expression analysis
• Datasets
• Final geneset overview
• Machine learning methods considered
• Explainability
• References

Differential expression analysis

After downloading FASTQ files from the TCGA-PRAD project, and obtaining the raw counts for every sample, differential expression analysis was carried out as explained in the paper. Detailed data generated during this processed is available below for download:

Description	#genes	Download
DEGs for TCGA-PRAD (full population, satisfying |logFC|>1 and FDR <0.05)	1991
DEGs for TCGA-PRAD (paired samples, satisfying |logFC|>1 and FDR <0.05)	1332
PRAD-DEGs: intersection of DEGs for TCGA-PRAD from paired and unpaired population (satisfying |logFC|>1 and FDR <0.05)	1065
55TOP-PRAD-DEGs	55
PRAD-DEGs-PROST-int-CANCER	114
47-PCa-Genes	47

Datasets

In general, all PC studies of significant size have an imbalance between cases and controls, which is mainly caused by the fact that control biopsies have to be taken from non-tumorous areas of patients with low tumor grade, or are negative biopsies taken from patients with high PSA at the time who may later develop the disease. Healthy individuals never have prostate biopsies taken for research purposes, as this would be ethically unacceptable because it is an invasive technique with side effects. Therefore, the strategy used to select the different populations is state-of-the-art when working with prostate cancer biopsies.

In addition, the TCGA-PRAD population allows raw data to be downloaded after permission has been granted, which was key during the training stage of our algorithm. Other databases, such as GSE183019, which was used as a validation cohort, contains a balanced number of cases and controls, but is smaller, includes only pre-processed data, and lacks ethnic diversity.

The population origin of the data is also a relevant factor to consider, as we wanted to include a wide variety of populations in the study, but no available dataset contained a balanced ethnicity proportion. This is another reason why we adopted the TCGA PRAD repository, which, although unbalanced, provides us with both healthy and prostate cancer-affected tissue samples (some of which are paired), ethnic diversity, and raw RNASEQ data from all samples, which are good conditions for applying our analysis pipeline prior to algorithm training. In fact, despite the complexity we had to manage due to the balanced training dataset, we were able to verify that our algorithm performs similarly in populations with different ancestries, demonstrating its ability to generalize its behavior.

The following table shows the datasets considered in this study. Along with the dataset name, the following fields are shown: the purpose of the dataset for this work TRAINING/VALIDATION; the number of samples in each dataset, distinguising between tumoral samples (T) and controls (NT); observations and a link to the data source.

Dataset	Purpose	Samples	Observations	Download
TCGA-PRAD	Model and parameters training	498T vs 52NT (52 paired)	Access controlled
GSE22260	Validation	20T vs 10NT (10 paired)
GTEX	Validation	245NT
GSE183019	Validation	84T vs 84NT
GSE114740	Validation	10T vs 10NT (paired)

One of the strengths of this system is its ability to classify samples into healthy/tumour samples independently of the genetic ancestry of the individuals in the study. This is an important problem to those systems trained with purely genomic data.

To demonstrate this, non-tumour samples from the GEO datasets: GSE114740, GSE22260 and 10 TCGA samples were clustered, against samples from the 1000 genomes consortium project (1KGP) [1]. In this consortium we found a total of 2504 samples sequenced at a minimum depth of 30x and classified according to their genetic ancestry. No clear information about the population is provided in the GEO datasets. For the TCGA samples, according with the clinical information, 8 white (not Hispanic or Latino) and two black or African American were selected.

To make this clustering more precise, 15,020 polymorphisms known as AIM (ancestry informative markers), have been used, in order to allow a classification into populations [2]. In detail, for each of the non-tumoral samples from the three datasets, fastq files were aligned using BWA [3]. Picard [4] was used to remove duplicate sequences, as well as to sort the reads and index the resulting bam file. These were processed with GATK (Genome Analysis Tool Kit) [5] to identify variants. In the case of the 1000 Genomes data, the variant files provided by the consortium were used. After filtering all the SNPs not related with the AIMs mentioned above with intersectBed (from bedtools) [6], plink was used to cluster the samples based on their genetic ancestry [7].

As shown in Fig. S1, we have 5 populations from 1KGP named AFR (Africans), EUR (Europeans), AMR (Admixed American: Colombians, Peruvians, Mexicans and Puerto Ricans), EAS (East Asians) and SAS (South Asians). From the GSE114740 study, we find that the 10 patients cluster very closely with East Asians, as would be expected given that the samples come from an hospital in Beijing. In the case of the GSE22260 study, the 10 healthy tissue samples cluster close to Europeans as expected given these are samples from North Americans and thus, of European ancestry, with two exception, one sample close to east Asians and another next to AMR/SAS. Similarly, the TCGA samples are from North Americans and therefore also cluster close to Europeans with the exception of the two black or African American samples that cluster close to the African individuals.

Final geneset overview

Below is an brief description of our final geneset 47-PCa-Genes, including their Gene_ID, a brief summary and their location coordinates. The full table, including an extensive description of every gene can be downloaded by clicking the following link:

Table S1. Brief description of the 47-PCa-Genes geneset, including Gene_ID, summary and location.

Machine learning methods considered

A brief description of every ML method considered in this work is shown below, including a short explanation of the parameteres trained for each one:

• Random forest [8]: This is a supervised machine learning algorithm which can be used both for regression and classification. It works by building a forest of trees with different samples and attributes and its output is based on the majority vote of those trees (for classification) or their average result (for regression).
  The "mtry" parameter for this method has been tuned by running a stratified 5-fold cross-validation.This parameter controls the number of variables randomly sampled as candidates at each tree split.
• XGBoost [9]: This is a powerful machine learning algorithm that combines weak predictive models, typically decision trees, to create a stronger and more accurate model. It uses a gradient boosting framework, iteratively improving the model by minimizing the residual errors. XGBoost is known for its scalability, speed, and performance in various tasks such as classification, regression, and ranking. It incorporates regularization techniques to prevent overfitting and provides interpretability features, like feature importance, to gain insights from the data. The following parameters were tuned:
  • nrounds: The number of boosting rounds or decision trees to build.
  • eta: The learning rate, controlling the contribution of each tree.
  • max_depth: The maximum depth of each decision tree in the ensemble.
  • subsample: The fraction of observations to randomly sample for each tree.
  • colsample_bytree: The fraction of features to randomly sample for each tree.
  • min_child_weight: The minimum sum of instance weight needed in a child for further partition.
  • gamma: The minimum loss reduction required to make a further partition on a leaf node.
• KNN [10]: This is a supervised machine learning algorithm which can be used both for regression and classification. It works by computing the distance between a query and all the samples present in the data, selecting then the k nearest instances to the query and averaging their labels (for regression) or voting for the most frequent label (if it is a classification problem).
  The "k" parameter for this method has been tuned by running a stratified 5-fold cross-validation.
• rpart [11]: This is a supervised machine learning algorithm which can be used both for regression and classification. The rpart library, available for R, implements recursive partitioning and is used to build classification and regression trees (CART). According to their creators, it is an implementation of most of the functionality of the 1984 book by Breiman, Friedman, Olhsen and Stone. A cart model's output is a decision tree, where each node represents a prediction for the outcome variable and each fork is a split based on a predictor variable.
  The "complexity parameter" (CP)  argument for this method has been tuned by running a stratified 5-fold cross-validation. The role of this parameter is to prune any split that does not improve the fit by, at least CP.

As described in the article, a stratified 5-fold CV with 5 repeats was used for each method and data balancing strategy during the training stage. Below, the results for every iteration for each method are displayed (only the best class/balancing strategy for each method is shown). Given the confusion matrix shown in Fig. S2 for a binary classification problem with two classes (positive and negative), where True Positives (TP; an example in which the model correctly predicted the positive class), False Positives (FP; an example in which the model mistakenly predicted the positive class), False Negatives (FN; an example in which the model mistakenly predicted the negative class) and True Negatives (TN; an example in which the model correctly predicted the negative class) are represented, we describe the metrics considered in this study.

• G-mean: This metric provides a combination of sensitivity and specificity by computing their geometric mean.
• Sensitivity: Measures what proportion of actual positives are correctly predicted as positive. Sensitivity = (True Positive)/(True Positive + False Negative)
• Specificity: This measure is defined as the proportion of actual negatives, which got predicted as the negative class. Specificity = (True Negative)/(True Negative + False Positive)
• AUC: Calculates the area under a Receiver Operator Characteristic(ROC) curve for a classifier. A ROC plots True Positive Rate (TPR) agains False Positive Rate (FPR) at various threshold values. AUC measures the area under this curve.
• F1: F-score, also called the F1-score, is a measure of a model’s accuracy on a dataset. It is used to evaluate binary classification systems, which classify examples into ‘positive’ or ‘negative’. The F-score is a way of combining the precision and recall of the model, and it is defined as the harmonic mean of the model’s precision and recall.

Note that RF and XGBoost are two ensemble learning methods that aim to improve prediction accuracy by combining the results of several simpler decision tree models. Both models can handle high-dimensional data, are robust to outliers and noise, and use regularization techniques to control model complexity and prevent overfitting. However, both models provided statistically similar results in the experimental study. For the sake of clarity for the reader, we decided to include only the results obtained by RF in our work because it obtained slightly better average results in 3 of the 5 quality measures (G-mean, specificity, and F1) used in our analysis.

The following table shows the mean across the 25 test sets for each balancing strategy in the test sets during the experimental analysis.

Below are several tables with the results for each iteration and repeat during the training stage on the discovery population (TCGA-PRAD) for each algorithm considered with the class balancing strategy that worked best for each one.

RF (undersampling strategy)
 

Table S2. Detailed results for each iteration in the training stage (5-fold CV with 5 repeats) for RF (undersampling strategy) over the discovery population.

KNN (hybrid strategy)
 

Table S3. Detailed results for each iteration in the training stage (5-fold CV with 5 repeats) for KNN (hybrid strategy) over the discovery population.

rpart (undersampling strategy)
 

Table S4. Detailed results for each iteration in the training stage (5-fold CV with 5 repeats) for rpart (undersampling strategy) over the discovery population.

Although RF is statistically proven in our article to be performing better than the other methods, in Fig. S3 we compare our final classifier with other classifiers built in this study across external populations. Only algorithms achieving at more than 0.4 in both sensibility and specificity are represented for simplicity. This graph also proves that our classifier based on RF also works better when it comes to generalizing the results in the validation populations.

Fig. S3. Classifiers performance across external databases, the discovery population is also represented in a different background. Labels show sensibility,specificity values for each classifier and database.

Explainability (SHAP)

Understanding why a model makes a certain prediction can be as crucial as the prediction's accuracy in many applications. However, the highest accuracy for large modern datasets is often achieved by complex models that even experts struggle to interpret, such as ensemble or deep learning models, creating a tension between accuracy and interpretability. In response, various methods have recently been proposed to help users interpret the predictions of complex models, but it is often unclear how these methods are related and when one method is preferable over another. To address this problem, SHAP (SHapley Additive exPlanations) assigns each feature an importance value for a particular prediction. SHAP, is a game theoretic approach to explain the output of any machine learning model. It connects optimal credit allocation with local explanations using the classic Shapley values from game theory and their related extensions. Shapley values are approximating using Kernel SHAP, which uses a weighting kernel for the approximation, and DeepSHAP, which uses DeepLift to approximate them.

The table below shows the 47-PCa-Genes geneset ordered according to their importance in our final classifier for the discovery population (TCGA-PRAD) as well as the validation populations.

Table S5. Gene list ordered by their SHAP importance in both the discovery and validation databases employed in this study.

The following table shows, for our discovery population (TCGA-PRAD) the mean, mininum, maximum and median contributions of each gene, ranked according to their importance. These values have been calculated for PCa and non-PCa affected tissue.This table can also be downloaded below.

Table S6. Summary of gene contributions to the classifier output for PCa and Non-PCa samples. Mean, minimum, maximum and median values are shown for each sample group in the discovery population.

Different SHAP attribution graphs are shown in Fig. S4 for the complete set of genes in the classifier for the TCGA-PRAD database. Every dot represents a gene for a specific sample and its color is related to its expression level (red=high; blue=low).

Fig. S4a. Shap attribution for every gene in the full set of samples included in the TCGA-PRAD database.	Fig. S4b. Shap attribution for every gene included in the T samples for the TCGA-PRAD database.	Fig. S4c. Shap attribution for every gene included in the NT samples for the TCGA-PRAD database.

Two different boxplots are shown in Fig. S5. On the left, the boxplot for Shap attribution is displayed for both T and NT samples separately. Similarly, on the right boxplot for gene expression is represented for T and NT samples (gene expression values have been previously normalized).

Fig. S5a. Shap attribution boxplot for TCGA.	Fig. S5b. Gene expression boxplot for TCGA.

References

1. C. Genomes Project et al., "A global reference for human genetic variation," Nature, vol 526, pp. 68-74, Sep. 2015.
2. F.D. Carmona et al., "A large-scale genetic analysis reveals a strong contribution of the HLA class II region to giant cell arteritis susceptibility," Am J Hum Genet, vol 96, pp. 565-580, Apr. 2015.
3. H. Li and R. Durbin, "Fast and accurate short read alignment with Burrows-Wheeler transform," Bioinformatics, vol 25, pp. 1754-1760, Jul. 2009.
4. B. Institute, Picard Tools, http://broadinstitute.github.io/picard/, 2018.
5. A. McKenna et al., "The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data," Genome Res, vol 20, pp. 1297-1303, Sep. 2010.
6. A.R. Quinlan and I.M. Hall, "BEDTools: a flexible suite of utilities for comparing genomic features," Bioinformatics, vol 26, pp. 841-842, Mar. 2010.
7. S. Purcell et al., "PLINK: a tool set for whole-genome association and population-based linkage analyses," Am J Hum Genet, vol 81, pp. 559-575, Sep. 2007.
8. T. K. Ho, "Random decision forests," in Proceedings of 3rd international conference on document analysis and recognition, vol. 1. pp. 278–282, 1995.
9. Chen T, Guestrin C. XGBoost: A Scalable Tree Boosting System. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining [Internet]. p. 785–94, 2016.
10. N. S. Altman, "An introduction to kernel and nearest-neighbor nonparametric regression," The American Statistician, vol. 46, no. 3, pp. 175–185, 1992.
11. T. Therneau and B. Atkinson, "rpart: Recursive Partitioning and Regression Trees," r package version 4.1-15, [Online]. Available: https://CRAN.R-project.org/package=rpart, 2019.