Application of Transcriptome–Based Gene Set Featurization for Machine Learning Model to Predict the Origin of Metastatic Cancer

1. Methods

1.1. Data Source

The data sources for training included transcriptome data for primary cancer collected from the OCP and TCGA databases. OCP data were created via micro-array with normalization [20], whereas the TCGA data, downloaded from Firehose (gdac.broadinstitute.org), were created via RNA-Seq and normalized to RNA-Seq by expectation-maximization (RSEM). Total 27,941 samples for 17 cancer tissue type were selected for training (Table 1). For external validation, 103 formalin-fixed paraffin-embedded (FFPE) tissue samples from patients collected by Kangbuk Samsung Medical Center (KBSMC) from 2018 to 2022. Patient samples were collected after informed consent was obtained, and the study followed the guidelines of the Declaration of Helsinki, and received approval from the Institutional Review Board (KBSMC 2022-11-018). Additionally, public datasets comprising 107 samples across 7 cohorts were obtained from the Gene Expression Omnibus (GEO) to further validate our findings (Table 2).

1.1.1. RNA Sequencing and Gene Expression Profiling

mRNA was extracted from FFPE tissue samples from patients using an RNeasy FFPE Kit (Qiagen, Hilden, Germany), according to the manufacturer’s instruction. In summary, FFPE tissue sections were deparaffinized by treatment with deparaffinization solution and lysed by proteinase K digestion followed by heat treatment. Next, the supernatants treated with DNase were added to Buffer RBC and ethanol to adjust binding conditions for RNA. The samples were applied to the RNeasy MinElute spin column, where the total RNA bound to the membrane and contaminants were efficiently washed away. RNA is then eluted in RNase-free water. The RNA concentration was determined using a NanoDrop (Thermo Fisher Scientific, MA, USA). Subsequently, RNA-seq was performed by Macrogen (Seoul, Korea) on the Illumina (CA, USA) RNA-Seq platform for paired-end sequencing, employing the SureSelectXT RNA Direct Reagent Kit (Agilent, CA, USA). The raw FASTQ RNA-Seq data were trimmed using the Trimmomatic-0.39-1 tool [21]. Alignment and quantification were performed using STAR 2.7.8a and RSEM 1.3.3 with the GRCh38.105 genome reference [22,23].

1.2. Featurization and Feature Selection

Figure 2A shows that the average expression levels of each gene between the TCGA and OCP groups for breast cancer are different in both range and distribution pattern. To integrate the different characteristics of these data, we converted the gene-wise information of the transcriptome to the gene set dimension. We created the gene set enrichment scores for all samples, and statistical values were extracted using the Kolmogorov-Smirnov test with 8,300 gene sets from the "Hallmarker gene sets", "C2 curated gene sets", "C6 oncogenic signature gene sets", and "C8 cell type signature gene sets" obtained from MsigDB (www.gsea-msigdb.org, v2023.2) [24,25,26,27,28,29,30,31,32,33]. Positive or negative signs were assigned based on the directionality of the expression difference.

After computing the 8,300 gene set scores, we investigated which of these could provide a comprehensive representation of the OCP and TCGA data. In feature selection, the receiver operating characteristic (ROC) based feature selection approach can be used effectively as a tool to evaluate individual features [34,35]. In particular, for the binary-class problem, the single feature classifier constructed from feature $f_i$ can establish an appropriate threshold $\theta$. If $x\ge \theta$, x is classified as the TCGA class. And $x<\theta$, x is classified as the OCP class. If one feature $f_i$ has an area under the ROC curve (AUC) value for a single-feature classifier farther from 0.5 than another feature $f_j$, we can say that $f_i$ is more discriminative than $f_j$ in the two classes [34]. We calculated AUC for each gene set feature to discriminate between the OCP cohort and TCGA cohort. The Figure 2B shows a histogram of the AUC values calculated for each gene set feature. An AUC of 0.5 means that the classifier has no discriminative capacity at all, which means that OCP and TCGA can be used comprehensively when using features around 0.5. Therefore, feature selection was performed at intervals of ± 0.05 around a baseline of AUC 0.5, and these selected feature groups were used for modeling.

We found that simply converting to gene set scores significantly normalized the ranges and distribution patterns of the OCP and TCGA cohorts. The Figure 2C-1) shows all 8300 gene sets used, and shows that TCGA and OCP data are more highly correlated than at the gene level. Moreover, filtering gene sets that show differences between groups by AUC range can achieve a greater correlation (Figure 2C-2). To examine the distribution through T-SNE clustering for major cancer types such as lung, stomach, and breast cancer, compatibility between TCGA and OCP can be observed (Figure 2D). The dimension reduction for features was performed using T-SNE from Scikit-learn to visualize the distribution of samples and groups in a 2-dimensional space, further aiding in the understanding and optimization of the model. The Figure 2D-1 is when all 8300 gene sets are used, and D-2 is when 1249 gene sets from the AUC range of 0.45-0.55 are used. In Figure 2D-2, using only 1249 features selected by AUC, we can see that the features are more focused on the characteristics of the cancer type than the data type.

1.3. Cancer Primary Site Classification Model

To create a classification ensemble model based on machine learning, logistic regression, LightGBM, and support vector machine (SVM) were modeled using a soft voting classification approach. For this algorithm, the Python packages Scikit-learn (version 0.22.1) and LightGBM (version 3.1.1) were used. LightGBM is renowned for its high efficiency and low memory usage, which make it particularly effective for handling large datasets with high dimensionality. It also excels in terms of speed and accuracy. Additionally, SVM is advantageous for its effectiveness in high-dimensional spaces and its versatility through the use of different kernel functions, enabling it to model non-linear boundaries. The ensemble approach, which combines these powerful classifiers, leverages their individual strengths to enhance the predictive performance and achieve significant results and parameters were set via a grid search.

1.4. External Validation

For external validation in metastatic cancer, predictions were made through the model using KBSMC and GEO data, and a ranking was created based on the probability. To further validate the model in the real world, we trained the model on public data and then evaluated it on a prospectively collected sample of cancer patients from the hospital. Through this process, the top-1 prediction accuracy was calculated.

2. Results and Discussion

We demonstrate that the more integrated the data are, the better the prediction performance for datasets that from real-world clinical data, and show how to integrate the two largest sources of transcriptomic data, micro-arrays and RNA sequencing, and use them as features in machine learning model. The Figure 3A shows the average F1 scores for each model with features in the range of AUC values shown in B of Figure 2. In this case, the average F1 score is the average of the predictions across the 17 cancer primary organs. The model was trained with TCGA and then tested with cancer types from OCP data, and conversely, the model was trained with OCP data and then tested with TCGA cancer types. Overall line does not shown as medians because it has a greater number of OCP data, 19541, than TCGA, 8400 (Table 1). We found that AUC 0.5 ± 0.1 to 0.3 was associated with an average F1 score of 0.9 or higher (Figure 3A). Therefore, we selected features from this range to build our model and conduct external data validation (Figure 3B).

Models trained on TCGA data, OCP data, and a combination of both TCGA and OCP were validated using external data. The best scores for the validation sets in each model are shown in Figure 3B and Table 3. Notably, models trained on TCGA data generated through RNA sequencing exhibited significantly lower performance when validated with KBSMC data, which were also produced through RNA sequencing. However, we showed a significant improvement in performance when micro-array data from OCP were included in the training with our feature integration method. The models trained on the combined TCGA+OCP dataset demonstrated an accuracy range of approximately 0.80 to 0.86, depending on the range of feature selection (Figure 3B, Table 3). In Figure 3, C and D are the confusion matrix for the top-1 prediction of the external validation set. The AUC of the feature selection group adopted in the evaluation model ranged from 0.5 ± 0.25. Most of the answers matched the actual primary site, and for those that were incorrect, the most likely answer was a neighboring organ. For example, when the bile duct is the correct answer, the incorrect answers are for the pancreas and colorectum, and the incorrect answers for uterus are ovary and kidney. In particular, the prostate, uterus, bile duct, and pancreas, which have a relatively high incorrect answer rates, have fewer than 10 data points, making it difficult to perform sufficient validation. Our model had an average accuracy of 0.9 when calculating the accuracy for KBSMC+GEO data up to the prediction of the top-2 (Table 4).

3. Conclusion

In this study, we demonstrated the feasibility of using gene set scores for the featurization of transcriptome data generated from different platforms, facilitating the integration of such data for training purposes. This approach enabled us to achieve higher performance, underscoring the potential of integrated datasets to enhance the predictive accuracy of our models. By converting gene expression data into gene set scores, we could mitigate the variability and noise often associated with data from heterogeneous sources, thereby improving the robustness and reliability of our predictions.

Our model was validated on samples of metastatic cancer patients collected from hospitals, and while the model trained on a single of public database did not predict real-world clinical data well, our method predicted the primary site with an accuracy of better than 0.8. While the top-1 accuracy achieved in this study remained in the 0.8 range, extending the metric to include top-2 accuracy revealed a performance level of approximately 0.9(Table 4). Furthermore, our model can learn the tissue characteristics of the primary site because the answers that missed the top-1 prediction also said organs similar to the correct answer. This suggests that with more comprehensive feature selection based on the featurization approach outlined in our research, and as the volume of data across all cancer types increases, we can anticipate even greater performance. In future research, we will validate whether our model could perform as well for other non-Asian ethnicities. Cancer genomic variation varies significantly by race, and in the public database TCGA, only 672 out of 11,122 patients are Asian [19], therefore we expect our model to perform better for Americans.