Barcode of patients with DM in the TCGA-HNSC, TCGA-BLCA, and TCGA-PAAD projects were retrieved from Genomic Data Commons (GDC) Data Portal, and gene expression data, clinical data, and histopathology images of these patients were downloaded. The few number of available cases of DM is a common challenge in most studies on metastasis. BLCA, and PAAD were selected based on the number of cases of DM in these cancer types that are available on The Cancer Genome Atlas (TCGA), and the number of cases with available corresponding genomic, image, and clinical data. A few cancer types were dropped during this selection process due to the absence of some variables in their records. For example, Skin Cell Carcinoma (SKCM) was eliminated due to missing “number_of_lymph_nodes_positive_by_HE” variable which is present in BLCA, PAAD, and HNSC datasets. HNSC was added based on the interest of one of the authors. See Supplementary Material S1 for file ID and barcode of samples. BLCA (N = 80) records had 59 male, and 21 female patients. Average age in PAAD (N = 58; Male = 32; Female = 26), and HNSC(N = 24; Male = 21; Female = 3) records is 63, and 61 respectively. The race of patients in the BLCA cohort are (White = 66; Black or African American = 6; Asian = 5), and (White = 47; Black or African American = 3; Asian = 5) for PAAD, while HNSC had (White = 19; Black or African American = 4; Asian = 1). The T, and N staging, and Number of lymph nodes positive by hematoxylin and eosin (HE) staining in each cancer type is shown in Table 1.

We designed a study that utilized machine learning techniques to simultaneously identify biomarkers of DM in each of BLCA, PAAD, and HNSC, and investigate if genes involved in DM are similar across the three cancer types. Each cancer type dataset was pre-processed separately. Equal numbers of complementary gene expression data of samples without DM (BLCA = 80; PAAD = 58; HNSC = 24) were downloaded for each of the cancer types, and we ensured that there were no overlapping samples between the groups with DM, and those without DM. After extracting Transcript Per Million (TPM)) normalized values of protein-coding genes from the individual records, all samples in a group were merged to create a table, and NaN values were replaced with zero. Other initial preprocessing steps include log to base ten transformation of data to adjust for the wide and non-linear values, removal of genes with a value of zero in greater than 80% of the samples, and selection of only top 10000 genes with highest variance across all samples.

To identify important genes that are involved in DM, we utilized the Random Forest (RF) algorithm with an optimization technique described in (Mori Y. et al., 2021). RF is a supervised ML algorithm that creates multiple decision trees from bootstrapped samples data and randomly selected features to output a result, that is, based on a voting system. It works well with high dimensional data, and is a commonly utilized technique in gene expression data analysis. The optimization technique we employed consists of running the RF algorithm over (N = 1000) iterations to classify samples as DM positive or DM negative, and selecting the top (K = 100) most important genes at each of the N instances. This is followed by ranking the important genes overall based on how many times each gene appears (i.e., frequency) in K over the N iterations. The algorithm was used to classify samples in each of the groups based on the presence or absence of DM (Figure 1). We selected the top 50 overall highest ranked genes for BLCA, PAAD, and HNSC. Also, the three cancer types datasets were combined to derive the “ALL3” group, and the above steps were repeated to select the top 50 highest ranked genes in this group as well (Figure 1).

To assess the strength of the selected biomarkers for prediction of DM, and to investigate if genes involved in DM are similar across the three cancer types, first, we used only the selected genes in each of the groups as features to train multiple ML models for the task of DM prediction. Support Vector Machine (SVM), K-Nearest Neighbor (KNN) and RF models were trained on datasets from each of the four (BLCA, PAAD, HNSC ALL3) groups (Figure 1), and the metrics were evaluated. Next, we looked for overlaps between different combinations of the selected genes in the four groups. Furthermore, we trained new ML (SVM, RF, and KNN) models on BLCA, PAAD, and HNSC datasets, but with genes selected from the combined dataset of the three cancer types (i.e., ALL3 group) (Figure 1). Lastly, we created a list of union of selected genes from each of the three cancer types (i.e., BLCA + PAAD + HNSC), and these were used as variables to predict presence or absence of DM in each of the three cancer dataset.

To further assess the strength of genes selected by the proposed method for predicting DM, under the same study design, we carried out an analysis to identify Differentially Expressed Genes (DEGs) between samples with DM and those without DM in each of the four groups. Instead of the TPM normalized values, raw counts of unstranded RNASeq were extracted from the TCGA RNASeq records. A sample information table was also generated. Differential Gene Expression (DGE) analysis was performed in R using the DESeq2 Bioconductor package. DESeq2 detects DEGs by normalizing raw count values of genes in the experimental groups, fitting negative binomial generalized linear models for each gene, and detecting significance by Wald test (Love et al., 2014). Threshold was set at p-adjusted value of 0.05, after initial set threshold of log fold-change of 1 and p-adjusted value of 0.05 yielded only five, one and seven DEGs in BLCA, PAAD and HNSC respectively. Thereafter, top DEGs were used as features to classify DM samples in each of the study groups (Figure 1). DEGs from this analysis and metrics of models trained on them are compared to those of genes selected via the (RF + Optimization) method.

We downloaded diagnostic pathology Whole Slide Histopathology Images (WSIs) of patients with DM in the TCGA projects of BLCA, HNSC, and PAAD from GDC Data Portal. The number of samples in each group reduced slightly after collating a list of only patients with available histopathology images, and RNASeq data, and clinical data. Number of samples with DM in BLCA, PAAD, and HNSC groups are 50, 44, and 18 respectively. Again, we downloaded random complementary WSIs (BLCA = 55; PAAD = 51; HNSA = 22) of patients without DM for each of the cancer types, and ensured no overlap between samples (Supplementary Material S1). Total number of samples in each group was split in an 80:20 ratio for models training and testing.

One WSI was preprocessed per patient. Due to the typical large size of WSIs which hovers around 100000 pixels in both horizontal and vertical axis, we isolated representative regions within each WSI for analysis. First, a maximum of 500 random non-overlapping patches of size 512 * 512 pixels were extracted from each image at 20 × magnification. These were reduced to a maximum of 200 patches after checking for a tissue area of at least 80% in each patch. To ensure selection of tumor representative regions we imitated the concept of high cellularity described in (Riasatian, 2021) which assumes that high grade tumors contain more cellular areas than normal tissues. Hence, we used a pretrained U-Net nuclei segmentation model which was trained on breast cancer histopathology images to rank the patches based on cellular content. The top 60 highest ranked patches were selected for each patient. Further random manual inspection led to the exclusion of a few more patches before Macenko normalization. At the end of initial preprocessing steps, BLCA, PAAD, HNSC, and ALL3 groups had a total of 4936, 4348, 1856, and 11154 training patches respectively, of which 25% were for validation.

A DenseNet121 model, pretrained on Imagenet data, and KimiaNet (Riasatian, 2021) were chosen for this study. This allows us to evaluate the effect of keeping weights of lower layers of a model fine-tuned on domain histopathology images (i.e., KimiaNet) on performance. DenseNet is a CNN architecture that attempts to solve the problem of varnishing gradient associated with deep neural networks by cumulatively concatenating features output of a layer within the architecture to input of subsequent upper layers, (Huang et al., 2016; 2017). In total, DenseNet121 contains 1 7 × 7 convolution, 58 3 × 3 convolution, 61 1 × 1 convolution, 4 average pooling, and 1 fully connected layers. KimiaNet is a DenseNet121 architecture Imagenet-pretrained CNN model fine-tuned on about 250,000 histopathology images. We replaced the classification layers in the models with 1 GlobalAveragePooling2D, 3 Dropout (0.1, 0.2, and 0.05), and 3 Dense (512, 32, and 1) layers, to train a new classifier head on top of the base convolutional layers, using a weakly supervised technique. The last hidden layer of the new classifier heads has 32 nodes which is subsequently used as features extractor. The modified DenseNet121 was separately trained on data from each of the four study groups (HNSC, BLCA, PAAD, ALL3) for a binary classification task of DM prediction ( Figure 2 ). Loss function was binary cross-entropy, and optimizer, Adam with a learning rate of 0.0001. Training epoch was set at 40.

For prediction on the test sets, patient level features were derived by passing all patches for each patient through the trained feature extractor models, and averaging all patch features from a slide. We trained traditional machine learning algorithm- SVM and a single layer MLP on the patient level features to classify samples as DM positive or DM negative (Figure 2). Considering this approach we have taken, HNSC group was dropped here based on small sample size.

Multimodal concatenation of patient level features derived from the CNN models to log transformed values of TPM normalized RNASeq data of selected DM associated genes from above was carried out. For BLCA, and PAAD groups, the 50 genes selected using the (RF + Optimization) technique were combined with 32 features derived from histopathology images of each patient to train classifiers, and predict presence or absence of DM (Figure 2). 150 genes, consisting of 50 each selected from BLCA, PAAD, and HNSC were combined with 32 image features to predict DM in the ALL3 group (Figure 2).

Likely due to the effect of a larger training dataset given the weakly supervised training technique that was employed, it was observed from early results that image only models trained on data from the ALL3 group generally performed better in their predictions on the test sets than other groups. Similarly, multimodal models of genomic and histopathology image were generally better than unimodal model of either data type in the ALL3 dataset, while metrics of multimodal models are generally lower than either genomic or image only models in other individual cancer dataset (i.e., BLCA, or PAAD). Therefore further experiments were carried out on the ALL3 dataset. We trained fully connected layers on KimiaNet, a DenseNet121 architecture model finetuned on about 250,000 histopathology images to extract features from the images in the ALL3 dataset. Metrics of multimodal models built from these features, and genomic data were compared to those from DenseNet121.

Four variables—Age, T stage, N stage, and Number of lymph nodes positive by HE—were selected from the clinical data, and concatenated with genomic, and histopathology features to predict presence or absence of DM in the ALL3 dataset (Figure 2). As part of the clinical data preprocessing, invalid/NaN values in the age variable were replaced by the rounded mean age of the cancer type while invalid/NaN T staging, N staging and number of positive lymph node values were replaced with 0. Min-max normalization was applied to age, and number of positive lymph node variables.

In each cancer type, samples with DM were combined with the same size of samples without DM. Also, datasets from the three cancer types were combined to derive the “ALL3” group. Following pre-processing, the optimized RF algorithm, as described in the methods section was used for binary classification of the samples (i.e., DM vs. non-DM), and simultaneous gene selection (top 50) in each of the groups (Supplementary Material S2). While RF is a popular ML based gene selection technique (Díaz-Uriarte and Alvarez de Andrés, 2006; Wenric and Ruhollah, 2018; Nivedhitha et al., 2020), to the best of our knowledge, at the time of this writing, this is the first time RF is combined with this optimization method for the task of gene selection. The top 5 genes associated with metastasis in BLCA are FKBP6, ASIC5, MAPK8IP1, F11R, and PABPC5, while CTSV, BIRC5, SERPINA7, CST2, KLHL3 are the highest ranked genes in PAAD. In HNSC, TM4SF1, C19orf18, EXPH5, FKBP2, PTPRZ1 are the highest ranked genes, and ALL3 group had CHRNG, CPT1B, CGREF1, GPR31, SPTBN5. Past publications have confirmed the activities of some of these genes in the metastasis of various cancers, either as oncogene or as tumor suppressors (Table 2). There are others whose roles in cancer metastasis are yet to be explored. A functional annotation search for the selected 50 genes on pantherdb.org revealed very similar classification patterns in all of the three cancer types. Approximately 50%–60% of the genes had no known category, and the most common functional classification of those with known categories in all the cancers are protein binding, and catalytic activity (Figure 3).

Dataset in each group was randomly split into train and test set in a 75:25 ratio, and SVM with linear kernel, KNN, and RF models were trained on the highest ranked 30, 25, 20, 15, 10, 5, 3, 2, and 1 gene from list of selected top 50 genes. Outcomes of single instance predictions, and five-fold cross-validation on the test sets were recorded (Supplementary Material S3).

To test the strength of genes selected based on the (RF + Optimization) method as predictors of DM, a DGE analysis was carried out with the DESeq2 software package to identify DEGs between the DM and non-DM samples. With a threshold adjusted p-value of 0.05, the number of DEGs in BLCA, PAAD, HNSC, and ALL3 are 229, 2142, 1100, and 658 respectively. Of these 5, 39, 19, and 13 overlap with the 50 genes selected using the ML method in each of the respective groups. Presence or absence of DM was predicted in the four groups using DEGs with the lowest 30, 25, 20, 15, 10, 5, 3, 2, and 1 adjusted p-value as variables. In almost all cases, genes selected using the ML (RF + optimization) techniques had higher evaluation metrics (Accuracy, F1-score, AUROC) than those selected via DESeq2 DGE analysis in the task of DM prediction ( Figure 4 ). Highest mean AUROC score achieved over a five-fold cross validation was 0.87, 0.92, 0.97, 0.79 in BLCA, PAAD, HNSC, and ALL3 groups respectively, compared to 0.65, 0.80, 0.92, 0.62 with DESeq2 DGE analysis, when models were trained on top 15 selected genes (Figure 4; Supplementary Material S3).

To investigate if primary tumors of different cancer types share similar gene expression profiles in metastasis, first, we looked for overlap between various combinations of the selected 50 genes in each group. There was no overlap between the different combinations of BLCA, PAAD or HNC. However, in the ALL3 group, 7 genes (CGREF1, SPTBN5, TAS1R3, FAM241B, EL5, CSPG5, and MAPK8IP1), and 6 genes (CFAP45, CST2, BIRC5, MAGED4, CST6, and FAIM2) overlap with those selected in BLCA, and PAAD respectively (Figure 5; Table 3). Also, it was observed that there was at least one member of the Z NF and FKBP gene family present within the list of selected genes in BLCA, and HNSC. Jen and Wang (2016) extensively reviewed multiple studies on the role of ZNF (Zinc Finger) gene family proteins in cancer progression, and metastasis, acting both as oncogenes, and tumor suppressors. Various studies (Fong et al., 2003; Sun et al., 2021) have also implicated members of the FKBP family in cancer metastasis.

Furthermore, we trained new ML (SVM, RF, and KNN) models on BLCA, PAAD, and HNSC datasets, but with genes selected from the combined dataset of the three cancer types (i.e., ALL3 group) (Figure 1). These models performed lesser than in cases of predictions of DM using genes selected in individual cancer types. Highest mean AUROC score of 0.76, 0.82 and 0.64 was achieved in BLCA, PAAD, and HNSC datasets respectively against 0.87, 0.92, and 0.97 of models trained on selected genes derived from individual cancer types (Figure 4; Supplementary Material S3). To further confirm this specificity of genes in each cancer type, a union list of selected top 15 genes in the three cancer types - BLCA, PAAD, and HNSC was created. Results show a drop in mean AUROC scores from 0.87, 0.92, and 0.97 in BLCA, PAAD, and HNSC datasets respectively when 15 genes generated from each cancer type were used as variables to 0.83, 0.88, and 0.77 when the total of 45 genes in the list of unified genes were used as variables to predict DM (Figure 4; Supplementary Material S3).

DenseNet121 was trained on histopathology images from BLCA, PAAD, and ALL3 groups, after the HNSC group was dropped due to small sample size. One WSI was preprocessed per patient as described in the methods section. For each group, classification models—SVM, and MLP were trained on patient level features extracted by the trained CNN models, and prediction was carried out on the test set features (Figure 2). The image features were also combined with genomic features from the same patient to build multimodal (Image + Genomic) models.

Early results showed that metrics of multimodal (Image + Genomic) models built from the ALL3 dataset generally improved on image, or genomic unimodal models by 1–3 percent margin, with genomic data contributing more to the multimodal metrics by a large margin. Following a five-fold Monte Carlo cross-validation, a multimodal (Image + Genomic) SVM classifier produced a mean AUROC score of 0.77, while the corresponding image and genomic unimodal models produced a score of 0.58, and 0.74 respectively. The mean AUROC scores from MLP classifiers are (Image + Genomic = 0.81; Image = 0.59; Genomic = 0.80). In either BLCA or PAAD dataset, results of multimodal (Image + Genomic) models mostly improved on those of the image only models, however they were generally below that of genomic only models. The highest mean AUROC derived from a BLCA dataset multimodal (Image + Genomic) model is 0.62, and the corresponding values from image only, and genomic only model are 0.57, and 0.80 respectively. Similar pattern was seen in the PAAD dataset [Image + Genomic = 0.85; Image = 0.51; Genomic = 0.90) (Figure 6; Supplementary Material S4)]. These outcomes emphasize the importance of a large dataset for building a more robust image model when a weakly supervised training technique is employed.

Based on these early results, we carried out further experiments with the ALL3 dataset. To evaluate the effect of using a CNN model with lower layers pre-finetuned on histopathology domain images as features extractor, fully connected layers were trained on KimiaNet with the ALL3 dataset histopathology images, and features derived from the model were combined with genomic data to predict DM in the test set. There was generally no additional advantage observed in the results of KimiaNet over those from DenseNet121. Highest mean AUROC from the MLP multimodal model (KimiaNet Image + genomic data) is 0.77. Unimodal image and genomic models from the same dataset produced AUROC scores of 0.55, and 0.80 respectively. The same pattern was seen with the SVM classifier (Image + Genomic = 0.73; Image = 0.56; Genomic = 0.74) (Supplementary Material S4).

Lastly, four clinical variables—Age, T stage, N stage, and Number of lymph nodes positive by HE were concatenated with genomic, and histopathology features to train SVM, and MLP classifiers. This led to an improvement of the mean AUROC score of the Image + Genomic SVM model from 0.77 to a score of 0.79. Mean accuracy, and F1 scores also improved from 0.67 and 0.67 to 0.71 and 0.71 respectively. MLP classifier however produced a mean AUROC score of 0.80 against the score of 0.81 achieved with Genomic + Histopathology data (Figure 6). Combining clinical with genomic data improved the metrics of the clinical only models, however there was no marked improvement from metrics of the genomic only models (Figure 6; Tables 4, 5).