MOSBY enables multi-omic inference and spatial biomarker discovery from whole slide images

MOSBY enables multi-omic inference and spatial biomarker discovery from whole slide images | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article MOSBY enables multi-omic inference and spatial biomarker discovery from whole slide images Yasin Senbabaoglu, Vignesh Prabhakar, Aminollah Khormali, Jeff Eastham, and 6 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3938444/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 06 Aug, 2024 Read the published version in Scientific Reports → Version 1 posted 11 You are reading this latest preprint version Abstract The utility of deep neural nets has been demonstrated for mapping hematoxylin-and-eosin (H&E) stained image features to expression of individual genes. However, these models have not been employed to discover clinically relevant spatial biomarkers. Here we develop MOSBY ( M ulti- Omic translation of whole slide images for S patial B iomarker discover Y ) that leverages contrastive self-supervised pretraining to extract improved H&E whole slide images features, learns a mapping between image and bulk omic profiles (RNA, DNA, and protein), and utilizes tile-level information to discover spatial biomarkers. We validate MOSBY gene and gene set predictions with spatial transcriptomic and serially-sectioned CD8 IHC image data. We demonstrate that MOSBY-inferred colocalization features have survival-predictive power orthogonal to gene expression, and enable concordance indices highly competitive with survival-trained multimodal networks. We identify and validate 1) an ER stress-associated colocalization feature as a chemotherapy-specific risk factor in lung adenocarcinoma, and 2) the colocalization of T effector cell vs cysteine signatures as a negative prognostic factor in multiple cancer indications. The discovery of clinically relevant biologically interpretable spatial biomarkers showcases the utility of the model in unraveling novel insights in cancer biology as well as informing clinical decision-making. Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction The developments in high-throughput RNA, DNA, and protein assays as well as in digital pathology have enabled multiple different high-resolution perspectives on a given tissue. Sampled from the same tissue, these data modalities are inherently linked, due to arising from highly similar cell populations and potentially the same biological state. Successfully inferring this common biology is key for relating different data modalities, and predicting one modality from another. Deep neural nets have emerged as an effective and flexible framework to represent the common biology underlying different readouts, opening the way to ‘translate’ any data modality into another 1–5 . A high utility use case for this translation task involves digital pathology where relatively easy-to-acquire hematoxylin and eosin (H&E) stained whole slide images (WSI) can be used to infer high-throughput molecular data that is time consuming or expensive to acquire. The success of inferring the joint biology (i.e. underlying latent distribution) between WSIs and molecular data relies heavily on extracting informative features from WSIs. Contrastive self-supervised learning has made a breakthrough in computer vision by learning high quality image representations in no-annotation settings as shown in curated benchmark datasets such as ImageNet 6,7 . Contrastive learning aims to discriminate between positive and negative images where positive pairs are obtained from augmentations from the same image and negative pairs are augmentations from different images 8 . Contrastive self-supervised learning has also shown promise in identifying cancer-specific morphological features in large-scale and heterogeneous histology datasets and training feature extractor models that outperform ImageNet-pretrained networks 9–11 . However, the superior performance of self-supervised features has largely been demonstrated for classification tasks, with limited attention on regression tasks. For instance, Wang et al. 8 developed an oncology-focused feature extractor RetCCL that was trained in more than 34,000 histology images and achieved state-of-the-art performance in cancer subtype classification. Regression tasks, such as inferring gene expression levels from WSI features, have been addressed with multi-output regression networks 1,12 ; yet, self-supervised learning-based image features remain largely underexplored in this context 13 . Among molecular data modalities, the inference of bulk and single-cell transcriptomics from H&E WSI features has received considerable attention 1,2,12 ; however, an underexplored extension of this question is whether mapping the same image features to proteomic data holds greater promise due to proteins being more proximally associated with phenotype. The full power and limitations of image to molecular data ‘translations’ will become more evident as these deep learning models are expanded to learn multiple -omic data modalities such as transcriptomics, proteomics, metabolomics, and whole exome or SNP array-based DNA measurements. Here, we present the MOSBY ( M ulti- Omic translation of whole slide images for S patial B iomarker discover Y ) model that employs contrastive self-supervised learning-based (RetCCL) features from H&E images to infer high-throughput molecular data such as transcriptomics (single genes or gene signatures), proteomics, and selected DNA-based features (e.g. cancer DNA fraction, average DNA methylation). MOSBY adopts a multiple instance learning approach where it partitions gigapixel WSIs to small tiles, makes tile-level predictions as an intermediate step, and then aggregates those to obtain slide-level predictions. At the inference stage, MOSBY enables in silico spatialization of bulk omic profiles by reconstructing tiles from a given WSI. Correlation between two omic features across all tiles of a slide reveals patient-specific colocalization or spatial exclusion patterns. Expanded to the cohort, these spatial features allow identification of clinically relevant spatial biomarkers. We demonstrate the performance of MOSBY in pan-cancer TCGA data, and spatially validate results with 1) spatially resolved transcriptomic data in breast cancer, 2) serially-sectioned CD8 IHC images in bladder cancer, and 3) tumor microenvironment (TME) cell type predictions in nonsmall cell lung cancer from an orthogonal deep learning model. In contrast to a whole transcriptome approach, we focus on TME cell type marker genes for MOSBY gene models, whereby model predictions could be contrasted with validation data in ( 3 ). Signature models, on the other hand, aim to encompass a wider biology including TME cell type, metabolic, and oncogenic pathway signatures. We finally utilize signature models to showcase the derivation of biologically interpretable colocalization features consistently associated with risk and disease state in human cancers. RESULTS We developed the MOSBY model that learned a mapping from RetCCL (contrastive self-supervised learning)-based WSI features to bulk transcriptome, proteome, whole exome, SNP-array, and DNA methylome based profiles (Fig. 1a). We limited our study to 21 TCGA indications 14 with at least 200 paired RNA-seq and H&E whole slide images, resulting in a total of 12592, 10192, and 12090 images matching with transcriptome, proteome and DNA-based data respectively (breakdown by indication in Supplementary Table 1). Analyzed transcriptomic features consisted of 55 TME-related genes 15,16 (Supplementary Table 2) and 175 gene sets that covered tumor-related processes 17–32 , metabolic pathways 33 , and TME cell type or process signatures 34–37 (Supplementary Table 3). The proteome model involved 191 total and phosphoprotein antibodies from the TCGA reverse-phase protein array (RPPA) panel that focused on tumor-intrinsic and oncogenic processes 38 . Tested DNA-based features were limited to tumor purity, cancer DNA fraction, subclonal genome fraction, and average DNA methylation (all continuous measures bound to the interval [0,1]). In addition to single indication models, we also trained MOSBY in pan-tissue, pan-cancer, pan-squamous and pan-adenocarcinoma settings. Pan-tissue models consisted of LUNG (lung adenocarcinoma and squamous cell carcinoma), KIDNEY (clear cell and papillary renal cell carcinoma), BRAIN (glioblastoma and low grade glioma), and COADREAD (colon and rectal adenocarcinoma). The pan-adenocarcinoma (ADENO) model consisted of pancreatic, lung, stomach, colon, rectal, prostate, and ovarian adenocarcinomas; while the pan-SQUAMOUS model included squamous cell carcinomas of the lung, cervix, and head and neck. The ADENO and SQUAMOUS models allowed us to investigate whether histological similarities beyond tissue architecture enabled MOSBY to learn a better mapping from WSI features to multi-omic data. Similar to the HE2RNA 12 model, a 2-layer perceptron was adopted as the multi-output regression network that mapped image features to omic variables, and a maximum of 8000 image tiles per slide were used. Separate models were trained for gene, signature, protein, and DNA variables. In contrast to HE2RNA, image tiles were randomly selected from each WSI to capture an unbiased representation of the entire slide, and were all used in training without clustering. In addition, batch-normalized transcriptomic and proteomic data were used as ground truth to enable across-indication comparisons with resulting MOSBY predictions. To obtain slide-level predictions, tile-level predictions of omic features were aggregated by averaging. Contrastive self-supervised pretraining benefits prediction of omic data from H&E whole slide images The RetCCL feature extractor utilized TCGA H&E images during contrastive training, however it was not supervised by gene expression information. Hence, there are no a priori guarantees for RetCCL-based image features to predict gene expression with higher accuracy than ImageNet-based features. The success in this task hinges upon the diversity of semantic (in this case, biological) features in the representations learned by RetCCL. Thus we investigated the performance differences between MOSBY models trained with RetCCL- or ImageNet-based image features. Models were trained with 5-fold cross-validation (80/20 percent training/test set split). A random one fifth of the training set was also allocated as validation set to determine an early stopping criterion using the Spearman correlation between slide-level model prediction and ground truth omic data. Spearman correlation was preferred to mitigate the effects of outlier samples in relatively smaller cohorts. Correlation coefficients from test sets were subsequently averaged across five folds to obtain the final performance score for a given feature. RetCCL-based image features consistently led to higher cross-validated test set averages in all four omic data types compared to features extracted with ImageNet-pretrained ResNet-50 architecture (Fig. 1b-1e). The pan-cancer model (PANCAN) with RetCCL-based features achieved the highest performance for all tested data types with median cross-validated Spearman correlation of 0.722 in single genes (0.673 with ImageNet-features) (Fig. 1b), 0.693 in signatures (0.647 with ImageNet-features) (Fig. 1c), 0.549 in proteomic data (0.512 with ImageNet-features) (Fig. 1d), and 0.6 in DNA features (0.533 with ImageNet-features) (Fig. 1e). In terms of single indication models, the thyroid cancer model (THCA) achieved best performance for single gene and protein expression data sets with RetCCL-based features. Median cross-validated Spearman correlations in these models were 0.59 vs 0.524 for single gene, and 0.443 vs 0.385 for protein expression data with RetCCL and ImageNet-based features respectively (Fig. 1b, 1d). The single indication models achieving best performance for signature and DNA data were the liver and bladder cancer models respectively (LIHC and BLCA) (Fig. 1c, 1e). Across tested signatures, the LIHC model showed a median cross-validated Spearman correlation of 0.532 with RetCCL-based vs 0.45 with ImageNet-based features. For tested DNA features, the BLCA model achieved a median cross-validated Spearman correlation of 0.553 with RetCCL-based vs 0.417 with ImageNet-based features. Taken together, these results indicated that contrastive self-supervised pretraining in large-scale histology datasets benefits prediction of multiomic data from H&E-stained whole slide images. MOSBY tile level predictions are validated with spatially resolved transcriptomic data For inference tasks, MOSBY ‘full models’ were trained with data from 80% of patients, with the remaining 20% used as the validation set to determine an early stopping criterion. An MOSBY full model was trained in the TCGA breast cancer (BRCA) cohort (N = 1576 WSIs), and subsequently deployed on H&E image tiles from a publicly available spatially resolved transcriptomic breast cancer dataset 39 (referred to as ST from here on, N = 68 WSIs) that served as an independent validation cohort. Image tiles were centered around ST spot coordinates to enable one-to-one comparison between spot level ground truth data and tile level model predictions ( Methods ). For individual slides, the MOSBY signature model predictions showed the highest concordance with ground truth for ‘poor differentiation’ (Stemness_Kim_Myc, Spearman r = 0.71) and stromal features (Stroma_Estimate, Spearman r = 0.63) (Fig. 1f, Supplementary Fig. 1a). Across 68 slides, concordance was highest for a monocyte signature with a median Spearman correlation of 0.238 (Supplementary Fig. 1a), and a maximum of 0.611 (Fig. 1f). This large variation across slides was observed for all tested signatures suggesting that the quality of spatially resolved data is critical in validating MOSBY predictions. Of note, CD8 T cell infiltration and proliferation-related model predictions also showed strong concordance with ground truth, demonstrating the variety in phenotypes captured successfully by the model (Supplementary Fig. 1b). Spot level concordance for single gene expression values was overall lower than that for signatures (Supplementary Fig. 1d), potentially driven by the large degree of zero reads (i.e. dropouts) in ground truth data (Supplementary Fig. 1e). However, model predictions for genes associated with stroma, plasma cells, and epithelial features showed good performance for individual slides ( COL1A2, MZB1, EPCAM respectively) (Fig. 1g). CD68 , a myeloid marker, showed the highest median concordance in gene models (Spearman r = 0.127), which was low in magnitude but consistent with the high concordance of myeloid cells in signature models (Supplementary Fig. 1d). As MOSBY aims to discover the joint latent space between H&E WSIs and bulk omic profiles, the correlation structure in ground truth data was also utilized to validate MOSBY model predictions (Methods). This analysis revealed that the MOSBY signature-signature correlation matrix was significantly different from random for 67 out of 68 slides (empirical p-value < 0.05) (Supplementary Fig. 1c). For gene-gene correlation matrices, this value dropped to only 28 out of 68 slides (Supplementary Fig. 1f), again highlighting that computing gene signature scores is an effective strategy to deal with the dropout issue in spatially resolved transcriptomic datasets. Despite reduced overall performance in gene models, predicted gene expression levels still successfully captured specific highly correlated gene blocks derived from ground truth data, such as those among cancer-associated fibroblast (CAF) and immune cell marker genes (Supplementary Fig. 1g). An analysis across 68 slides indicated that MOSBY model predictions were particularly good at distinguishing CAF-related features from other TME population marker genes. For 63 out of 68 slides, the MOSBY gene-gene correlation matrix had a cluster significantly enriched for fibroblast marker genes, and this enrichment was higher than 89% for the large majority of the slides (Supplementary Fig. 1h). Serially sectioned CD8 IHC whole slide images validate MOSBY in silico spatialization MOSBY tile-level predictions were further validated with CD8 antibody-stained immunohistochemistry (IHC) WSIs. To this end, gene and signature models were trained with data from urothelial carcinoma immune checkpoint inhibitor (ICI) trials IMvigor210 40 and IMvigor211 41 (Methods, N = 1460 paired H&E and RNA-seq samples). H&E-CD8 double-stained slides were not available, thus serially-sectioned single-stained H&E and CD8 IHC whole slide images were computationally aligned (Methods, N = 1017 paired H&E and CD8 IHC WSIs). 454 slides that were correctly aligned based on multiple alignment metrics were subjected to a manual QC check and a stringent 95% tissue overlap threshold (Methods, Fig. 2a). The remaining 42 IHC WSIs constituted the analysis set, and were used to categorize “brownness” of pixels based on diaminobenzidine staining ( i.e. DAB mask)(Methods). Each DAB mask was split into tiles, and “brown” pixels in each tile were counted to determine tile-level CD8 IHC quantification (Methods). Pearson correlation coefficient was then computed between CD8 IHC and MOSBY-predicted tile-level values (Fig. 2a). To compare with CD8 IHC quantification, gene features CD8A, CD8B , and signature features Cibersort_CD8_T_cells and T_effector_cells were used from MOSBY model predictions. Across 42 slides, CD8B and Cibersort_CD8_T_cells showed the strongest concordance with CD8 IHC, having the highest correlation in 22 and 15 slides respectively (Fig. 2b). Overall, CD8B model predictions had the highest median correlation with CD8 IHC tile-level values (Pearson r = 0.352), with a range of [-0.006, 0.789]. The density plot for r = 0.789 correlation (N = 16859 tiles, p ≈ 0, two-tailed exact test) indicates stronger concordance for values closer to the higher and lower ends of the distribution (Fig. 2c). Visual inspection of CD8 IHC and H&E WSIs confirmed successful alignment, but also revealed minor differences between tissue sections potentially due to these being non-adjacent sections (Fig. 2d). At low magnification, visualization of DAB mask tile level quantification and CD8B model predictions indicated that the model successfully captures CD8 T cell infiltration patterns in the tumor (Fig. 2e, 2f). High magnification inspection of model predictions and DAB mask (as well as the input H&E and CD8 IHC images) confirmed that MOSBY, despite being weakly supervised on a patient level bulk RNA-seq dataset, effectively learned a CD8 T cell-specific representation in H&E images (Fig. 2g–2j). In order to gain insight into both high and low concordance cases, four representative slides were selected based on the distribution of predicted CD8B vs DAB mask correlations (near the minimum, and 25th, 50th, and 75th percentiles). Selected slides had Pearson correlation values 0.48 (N = 30032 tiles, p ≈ 0, two-tailed exact test), 0.35 (N = 22942 tiles, p ≈ 0, two-tailed exact test), 0.25 N = 11642 tiles, p ≈ 0, two-tailed exact test, and 0.011 (N = 16595 tiles, p = 0.16, two-tailed exact test) respectively (Supplementary Fig. 2a–2d). Low magnification inspection of DAB mask and MOSBY CD8B predictions indicated that model predictions captured overall CD8 T cell infiltration patterns in the first three slides despite the correlation being as low as 0.25 (Supplementary Fig. 2a–2c). For the last slide (r = 0.011), spatial patterns in DAB mask and model predictions were conspicuously discordant (Supplementary Fig. 2d). High magnification inspection of CD8 IHC and DAB mask images from this last slide revealed that DAB mask captured brown staining artifacts in regions that did not show any expected CD8-specific staining patterns (Supplementary Fig. 2e, 2f). In addition, CD8 IHC and H&E tissue sections showed misalignment in small substructures, despite an overall high score in tissue alignment (Supplementary Fig. 2e-2h). Taken together, these results showed that MOSBY successfully learned a CD8 T cell-specific representation on H&E images, and correlations as low as 0.25 may be sufficient to capture overall CD8 T cell infiltration patterns in the tumor microenvironment, suggesting the potential of this tool for patient stratification in ICI trials. Validation with pathologist annotations of cancer regions indicates epithelial gene and signature markers may be misleading in accounting for tumor-specific expression Various RNA- and protein-based epithelial markers such as EPCAM (GENE) and E-cadherin (PROTEIN) are frequently used in the literature to mark cancer cells in tissue specimens. However, epithelial markers can be expressed in both transformed and non-transformed epithelial cells. We utilized MOSBY tile level predictions in order to compare virtual spatialization of epithelial markers with pathologist tumor region annotations and assess the tumor-specificity of commonly used tumor markers. MOSBY full models were trained in TCGA lung adenocarcinoma (LUAD), and deployed on WSIs in IMpower150 43 , an independent validation cohort in nonsquamous nonsmall cell lung cancer with pathologist annotations (i.e. pen marks) demarcating cancer epithelia. Pathologist annotations were compared with MOSBY features of tumor purity (DNA-based) and epithelial markers EPCAM (GENE), E-cadherin (PROTEIN), and an epithelial gene signature by Taube et al 20 . WSI level comparison indicated that MOSBY models trained with DNA-based quantification of cancer cells (tumor purity) achieves a stronger concordance with pathologist tumor region annotations (Fig. 3a, 3b). E-cadherin closely followed DNA-based features in terms of concordance with pathologist annotations, albeit noticeably lower specificity for the cancer region in Slide 2 (Fig. 3d vs 3b). Interestingly, predicted spatialization of EPCAM transcript levels showed relatively poor concordance with cancer epithelium annotation (Fig. 3c); and the performance of the Taube et al . epithelial signature was further reduced with highest predicted tiles corresponding to non-tumor regions (Fig. 3e). Collectively, these results suggest that employing DNA-based tumor cellularity estimates is the most reliable approach for accounting for tumor-specific expression in bulk datasets; and published epithelial signatures may not be fit-for-purpose in this task. MOSBY enables investigation of intratumor heterogeneity with prediction of single gene, signature, and protein expression levels The MOSBY framework employs RNA-seq and RPPA data to predict gene (or signature) and protein expression levels respectively (Methods). The utility of the model in predicting single gene, signature, and protein expression levels and inferring intratumor heterogeneity is demonstrated in Impower150 in the context of cell proliferation features (Supplementary Fig. 3). Three representative H&E slides are shown in this figure together with pathologist annotations of cancer epithelia (Supplementary Fig. 3a), and MOSBY-predicted spatialization of cell proliferation markers Cyclin B1 (Supplementary Fig. 3b), Hallmark G2M checkpoint signature (Supplementary Fig. 3c) and MKI67 (Supplementary Fig. 3d). Even though tile-level predictions of all three cell cycle markers are highest in cancer epithelia-annotated regions, the spatialization of protein and signature markers (Cyclin B1 and Hallmark G2M checkpoint signature) demonstrate a stronger localization to cancer epithelia (Supplementary Fig. 3b, 3c, 3d). In addition, the spatialization of MKI67 expression in all three slides appears to suffer from checkerboard-like artifacts (Supplementary Fig. 3d). These results suggest that molecular data prediction from WSIs may be more successful with proteins or gene signatures, as opposed to single genes. Potential reasons for this observation are 1) proteins being closer to phenotype than genes and 2) gene signatures increasing the signal-to-noise ratio from RNA-seq data. Biologically, the predicted spatial pattern of Cyclin B1 and Hallmark G2M checkpoint signature levels show differential proliferation levels across cancer epithelia tiles (Slide 3 in Supplementary Fig. 3b, 3c), indicating MOSBY also has potential to infer intratumor heterogeneity. MOSBY predicts stroma, immune, and proliferation features with highest accuracy in TCGA Cross-validated test set performance in single-indication signature and protein models was used to determine biological features with the highest prediction accuracy, and those with 25 highest performance are highlighted in Fig. 4a, 4d. The best-performing signatures ESTIMATE Stroma and ESTIMATE Immune had median correlations 0.627 and 0.616 across all tested indications (Fig. 4a). In randomly split test sets, the Spearman correlation of ESTIMATE Stroma and ESTIMATE Immune signatures reached as high as 0.787 and 0.781 in skin cutaneous melanoma and thyroid cancer respectively (Fig. 4b). Other best-performing signatures were again largely enriched in stroma and immune features as well as general mesenchymal characteristics (such as Hallmark EMT and Taube et al mesenchymal signatures). (Fig. 4a). Of note, particular signatures known to play an important role in specific cancer indications also showed strong prediction accuracy in those pertinent settings. For instance, the Hallmark Angiogenesis signature had Spearman r = 0.746 in the liver cancer model (LIHC), and a fatty acid elongation signature had Spearman r = 0.741 in the low grade glioma setting (LGG) (Fig. 4b). After adjusting for multiple hypothesis testing with false discovery rate (adj. p < 0.05), 14 out of 21 tested indications showed significant performance (as measured by Spearman r) for more than 90% of tested signatures (total N = 175 signatures) (Fig. 4c). In contrast with our signature set that represented both tumor-intrinsic pathways and cell populations in the TME, the TCGA RPPA antibody set was heavily focused on tumor-intrinsic pathways. Therefore, the best-performing proteins had a diversity of representation from tumor-intrinsic characteristics, such as proliferation (Cyclin-B1, FOXM1), DNA repair (MSH6, PARP1), and apoptosis (cleaved Caspase7) (Fig. 4d). MSH6 and Cyclin B1 were the 2 best-performing proteins that were tested in at least 20 indications (Fig. 4d). These proteins had median Spearman r = 0.423 and 0.417 respectively across tested indications, but correlations in individual indications were as high 0.676 in PRAD for Cyclin B1, and 0.593 in LGG for MSH6. Overall, lower prediction accuracy in protein models was expected due to the lower signal-to-noise ratio in the RPPA technology compared to RNA-seq. However, in specific settings such as the PRAD model, both total and phosphoproteins showed strong prediction accuracy in test sets. Progesterone receptor (PR) and cMET_pY1235 reached Spearman r = 0.614 and 0.631 respectively in this indication (Fig. 4e). Also, the AKT/mTOR pathway showed evidence of strong prediction accuracy in the sarcoma setting where S6 and RICTOR_pT1135 antibodies showed Spearman r = 0.731 and 0.695 respectively (Fig. 4e). Despite paucity of representation in the RPPA panel, stromal and immune features were also found among the best-performing proteins such as ECM-associated Collagen-VI, Fibronectin and T cell-associated Lck (Fig. 4d). Overall, 11 out of 20 tested indications showed significant performance (adj. p < 0.1) for more than 90% of tested antibodies (total N = 191 antibodies) (Fig. 4f). In terms of single gene models, best performing features confirmed the strong prediction accuracy associated with stromal features observed above; the highest-ranking genes were marker genes of fibroblasts (e.g. LUM, COL5A1) (Fig. 4g). Top-ranking genes also included markers of macrophages (e.g. CSF1R), and T cells (e.g. CD3E). In DNA models, tumor cellularity measures (tumor purity, cancer DNA fraction) achieved better performance than subclonal genome fraction, and average DNA methylation features (Fig. 4h). Multi-indication runs may show exaggerated performance The results above demonstrated that MOSBY models trained with gene signature and protein expression data generally perform better than those trained with single genes. Thus, using signature and protein expression data, we next addressed whether single- vs multi-indication models (introduced in Fig. 1b) had more utility from a clinical perspective. We focused on performance in the unseen test set to evaluate the generalization capacity of each approach in a scenario emulating the clinical setting; i.e. predicting gene and protein expression levels for a single patient. Performance was measured, as above, with the Spearman correlation between MOSBY slide-level predictions and ground truth levels of signature or (phospho)protein expression (N = 175 for signatures, 191 for total and phosphoproteins). Multi-indication models (PANCAN, ADENO, BRAIN, KIDNEY, SQUAMOUS, LUNG) showed superior performance for both signature and protein prediction tasks, compared to single-indication counterparts (Supplementary Fig. 4a). The PANCAN model had a median cross-validated Spearman correlation of 0.693 for tested signatures (Supplementary Fig. 4a), with epithelial, proliferation, and stemness signatures showing the best prediction accuracy (cross-validated Spearman r = 0.826, 0.824, and 0.82 respectively). We focused on the Benporath_ES2 stemness signature to demonstrate the bias associated with multi-indication models. In a randomly split test set, this signature had Spearman r = 0.84 between predicted and ground truth levels in the PANCAN dataset (Supplementary Fig. 4b, left). However, the correlations in individual indications varied between 0.244–0.645 for this signature, suggesting that the PANCAN model performance may be spuriously inflated by combining multiple indications with different distributions; reminiscent of the Simpson’s paradox (e.g. OV and READ slides concentrated on the higher end, and THCA and KIRC concentrated on the lower end of the distribution; Spearman r = 0.248, 0.326, 0.41, and 0.482 for OV, READ, THCA and KIRC respectively) (Supplementary Fig. 4b, right). In cross-validated protein models, Cyclin B1, Claudin-7 and progesterone receptor (PR) showed the best prediction accuracy in the PANCAN model (Spearman r = 0.775, 0.772, and 0.749 respectively). To ask whether this performance was spuriously inflated, we again compared PANCAN and single indication models. In a randomly split test set, Cyclin B1 had Spearman correlation of 0.76 in the PANCAN model (Supplementary Fig. 4c, left), while correlations in individual indications were lower ranging from − 0.10 to 0.597. Simpson’s paradox was again evident in protein model predictions with Cyclin B1 levels in individual indications localizing to smaller parts of the pan-cancer distribution and showing lower prediction accuracy (e.g. Spearman r = − 0.103, 0.198, 0.214, and 0.263 for LIHC, THCA, READ and CESC respectively) (Supplementary Fig. 4c, right). We then investigated the potential bias in multi-indication models more systematically by comparing median test set performance (i.e. median correlation for signatures or proteins) in single vs multi-indication models (Supplementary Fig. 4e, 4g). The train/validate/test set splits in multi-indication models were implemented with indication stratification, enabling comparison with single-indication training. Results showed that neither single-indication nor multi-indication training was universally better (Supplementary Fig. 4e, 4g). However, evidence also showed that multi-indication training generally did not improve performance if a single-indication model reached a median performance of 0.45 for signatures (Supplementary Fig. 4e), and 0.35 for proteins (Supplementary Fig. 4g), with LGG as a potential exception benefiting from joint training with GBM. Pan-tissue training appeared to improve median test set accuracy in certain settings such as protein models in kidney, brain, lung, colon/rectum (Supplementary Fig. 4g), and signature models in brain (Supplementary Fig. 4e). However, at the individual signature/protein level, further heterogeneity in model performance was observed (Supplementary Fig. 4f, 4h). A subset of signatures/proteins benefited from multi-indication training whereas another subset did not. This observation prevented us from making general recommendations as to the deployment of single vs multi-indication models in a clinical setting. In circumstances where there is an interest in specific signatures/proteins, it may be suboptimal to make the model decision based on median performance. However, results in Supplementary Fig. 4e, 4g can be used as a guideline in other settings that lack a focus on specific biological features. We utilize single indication models in subsequent parts of this study in order to prevent potential bias from multi-indication training. Spatial patterns inferred from tile-level MOSBY predictions increase survival predictive power of gene signature-based models MOSBY tile-level predictions enable in silico spatialization for a tested omic feature, as well as assess spatial correlation between two tested features (positive and negative correlations indicating colocalization and spatial exclusion respectively). We define a ‘colocalization feature’ as the Pearson correlation between tile-level predictions of two omic features. As correlation coefficients were computed across all tiles on a slide, colocalization features represent slide-level as opposed to ‘local’ spatial patterns. We focused on our signature panel (N = 175) as the omic features to derive colocalization features and investigate survival associations, as the signature panel covered both tumor and non-tumor TME components. Moreover, using signatures as opposed to single genes enable the discovery of ‘biologically interpretable’ spatial biomarkers as well-designed signatures capture pathway and cell type-related gene expression with higher fidelity. For a slide, correlations from all pairwise combinations of signatures (N = 15225) (i.e. the collection of all colocalization features) are referred to as the ‘colocalization map’ from here on (Fig. 1a). A survival analysis was performed to ask whether slide-level patterns in colocalization maps harbored survival signals that could not be captured by the mere magnitude of signature expression. Three L1-regularized Cox proportional hazards regression models were fit to address this question, and the ability of the models to predict survival was computed with concordance indices (c-index) (Fig. 5a) ( Methods ). The first model only used flattened colocalization features (MOSBY predictions, N = 15225). The second Cox model only used signature expression levels (N = 175) with the goal of assessing the survival predictive power of gene expression ‘magnitude’. The third model was a joint model combining all signatures but only lasso-selected colocalization features from the first model in order to prevent colocalization features from dominating the model. Each model was run across 10 cross-validation folds to optimize the shrinkage parameter and also obtain mean and standard error estimates for the relevant c-index ( Methods ). This process was implemented separately for all tested TCGA indications, and c-index mean and standard error estimates were plotted (Fig. 5a). We observed that the joint (i.e. third) model had a higher c-index than the signature-only model in most indications. This finding revealed that slide-level spatial patterns, discovered with an inference engine such as MOSBY, have survival predictive power that could not be captured by gene expression alone. As MOSBY colocalization features are biologically interpretable, this opens the door to discovering potentially ‘actionable’ colocalization or spatial exclusion patterns that predict clinical outcomes. Moreover, the c-indices achieved in the joint model were highly competitive with or higher than those reported in the literature from end-to-end survival-trained multimodal models utilizing genomic, transcriptomic, and image datasets 44 . Of note, comparing colocalization-only and signature-only models, we also found that the total survival predictive power of slide-level spatial patterns was not as high as that of signature levels in most TCGA indications. Ovarian and rectal cancer results were an exception to this general pattern (Fig. 5a), suggesting spatial biomarkers discovered in these indications may have the greatest potential to lead to novel insights. Colocalization maps enable discovery of biologically interpretable spatial biomarkers We next investigated consistent spatial predictors of risk that were supported across multiple indications and also showed evidence of tumor specificity. In each indication, survival effects and tumor specificity of colocalization features were explored with two tests: 1) A univariate Cox regression model for survival, and 2) a Mann-Whitney test to compare tumor-normal levels of colocalization value. A potential spatial biomarker of risk was defined as a colocalization feature that was significantly associated with poor overall survival (p < 0.05) and also had elevated levels of colocalization in the tumor (adj. p < 0.05). Given these criteria, four colocalization features had evidence in four different TCGA indications to be a spatial biomarker of risk (Fig. 5b). Of these four, the colocalization between an ER stress signature (XBP1s targets ER17 22 ) and a neurotransmitter signature was associated with poor survival and malignant state in colon adenocarcinoma, lung adenocarcinoma, liver hepatocellular and ovarian cancers (Fig. 5c, 5d). In an independent non-squamous lung cancer study involving both immune checkpoint blockade (atezolizumab) and chemotherapy arms (IMpower110), this colocalization feature was also found to be associated with poor survival in the chemotherapy, but not immunotherapy arm, suggesting higher relevance as a resistance factor in chemotherapy (Fig. 5e, 5f). Of note, the ER17 and neurotransmitter signatures were not individually found to be associated with risk in any of the four mentioned indications (Supplementary Fig. 5a, 5b). Visual inspection of WSIs indicated that the expression of ER17 and neurotransmitter signatures primarily came from the microenvironment (as opposed to tumor region) in the case of high colocalization (Supplementary Fig. 5c). In low colocalization cases, the neurotransmitter signature expression primarily came from the tumor region whereas the ER17 signature expression was again predominantly in the microenvironment (Supplementary Fig. 5d). Focusing on colocalization features involving immune system signatures, the T effector cell vs Cysteine colocalization was identified as the most consistent spatial biomarker of risk in TCGA. This colocalization feature was associated with poor survival and also showed significant tumor enrichment in breast, squamous lung, and ovarian cancers (Fig. 6a, 6b). T effector and cysteine signatures were not individually found to be associated with risk in any of these indications (Supplementary Fig. 5a, 6b). Visual inspection of WSIs indicated that the expression of T effector and cysteine signatures primarily came from the microenvironment (as opposed to tumor region) in the case of high colocalization (Supplementary Fig. 6c). As high colocalization is a risk factor, this expression pattern may be suggestive of a cysteine-associated immunosuppressive TME. In low colocalization cases, the T effector signature expression primarily came from the microenvironment whereas the cysteine signature expression was predominantly in the tumor region (Supplementary Fig. 6d). The strongest survival effect for T effector vs cysteine colocalization was observed in breast cancer, where we investigated other immune cell types and found that immune vs cysteine colocalization was a general negative prognosis biomarker in this indication. Significant survival associations were observed for both lymphocytes/NK cells (Fig. 6c), and myeloid populations (Fig. 6d). The same immune vs cysteine colocalization features were also found to be negative prognostics in the Atezolizumab arm of Impower110 non-squamous cohort (Fig. 6e, 6f). Of note, most of these features were not prognostic in the chemotherapy arm (Supplementary Fig. 6e, 6f), however did not qualify as predictive biomarkers since the survival association differences between Atezolizumab and chemotherapy arms were not significant. DISCUSSION The MOSBY workflow achieves prediction of bulk omic profiles from H&E WSI features. Our results showed that, compared to ImageNet-based pretraining, self-supervised pretraining in large histological datasets allows creation of inference engines (e.g. RetCCL) that enable a more accurate mapping from image features to gene, signature, protein, and DNA-based measurements. We demonstrated that the most accurately predicted features by MOSBY involved processes such as proliferation, immune/stromal infiltration, differentiation, and epithelial-to-mesenchymal transition. This finding suggests that self-supervised features learned from pan-cancer histological datasets run the risk of accentuating biological processes and pathways that show the highest variation across different cancer indications. Training feature extractors on images from a single indication may be required to capture biological processes that play an important role in one or only a few cancer indications. RetCCL-based features showed promise by capturing angiogenesis in hepatocellular carcinoma, and fatty acid biology in low grade glioma. Yet, the accumulation of even larger histological datasets in the future have potential to allow refinement of image features relevant for indication-specific pathways, thus making possible the discovery of a greater number of clinically relevant spatial biomarkers. The utility of H&E images to predict high-throughput expression data has been largely explored in the context of individual genes 1,2,12 . The RNA-Seq noise level associated with individual gene transcript counts can be alleviated by utilizing focused gene sets (i.e. signatures) aiming to measure particular pathway or cell type levels. Also, proteins and posttranslational modification levels (e.g. phosphorylation, acetylation) capture a more faithful representation of tissue phenotype compared to gene transcript levels. We demonstrated with pathologist tumor annotation validations that training the model with signature or protein expression levels led to potentially more interpretable results compared to single gene-trained models. Moreover, the MOSBY workflow has the benefit of employing both RNA and protein level data to interrogate the spatial distribution of biological processes, each modality functioning as an independent validation platform for the other. Therefore, training MOSBY with signature and protein levels have the highest potential of uncovering spatial patterns relevant for cancer biology and clinical outcome. In the case of tumor cellularity, our results indicated that DNA-based measures outperformed RNA and protein-based epithelial markers, and importantly, highlighted that published epithelial signatures may be misleading due to being expressed also by normal cells. MOSBY, as with many other deep learning models, adopts a weakly-supervised approach by first making tile-level predictions and then aggregating tiles to obtain a prediction at the WSI level 4,45 . Although an intermediate output of the model, tile-level predictions enable in silico spatialization of whole slide-level annotations, opening the way to inferring intratumor heterogeneity for the slide-level information used as ground truth 12,44 . Spatial intratumor heterogeneity patterns learned from a cohort of patients subsequently allow investigation of clinically relevant biomarkers. In MOSBY, spatial patterns are captured by pairwise colocalization features. For a given signature pair, the colocalization value on a slide is defined as the Pearson correlation across all tiles. Thus, MOSBY colocalization features capture slide-level but not local spatial patterns. Local processes such as tertiary lymphoid structures are known to affect patient survival and response to cancer immunotherapy 46,47 . We demonstrated in this study that slide-level spatial patterns also carried survival signals, and increased predictive power of gene signature-based survival models in most TCGA indications. Moreover, we noted that the concordance indices of joint models (both colocalization features and gene signature levels) either surpassed or were comparable to those of multimodal deep learning models that employed the complete set of WSI and omic (RNA-seq, mutation status, copy number variation) data in TCGA 44 . This finding indicated that the 175-signature panel we defined was sufficient to capture most biological processes important for clinical outcome. End-to-end neural networks trained to predict survival have been lacking in terms of direct biological interpretation of image regions important for the model 44,48 . These models may incorporate mechanisms such as attention heatmaps and spatial credit assignment to increase interpretability 44,49 , yet still require pathologist efforts to examine important image regions whereby interpretation remains limited to phenotypes visible by human eye. In contrast, MOSBY spatial features are biologically interpretable by design, which is an advantage of this approach over end-to-end neural networks trained to predict survival 44,48 . In this study, we showed that the colocalization of an endoplasmic reticulum (ER) stress-related signature and a neurotransmitter signature is both elevated in tumors and associated with poor overall survival in four TCGA indications. The poor survival association in lung adenocarcinoma was also validated in the chemotherapy arm of an independent NSCLC cohort (nonsquamous samples in Impower110), suggesting this colocalization may be a chemotherapy-specific risk factor. Moreover, we identified the T effector cell vs cysteine signature colocalization as a TME-related risk factor in multiple TCGA indications, as well as in the immunotherapy arm of Impower110 nonsquamous cohort. These results showcase the high utility of colocalization maps for discovering biologically interpretable clinically relevant spatial biomarkers. Large oncology data sets have enabled the routine development of machine learning models trained simultaneously in multiple cancer indications from the same organ ( i.e. pan-tissue training), or from different organs across the body ( i.e. pan-cancer training). Such models have shown success in mapping molecular tumor features to patient outcomes as well as in mapping WSI features to transcriptomic data 12,50,51 . However, in the clinical setting where the machine learning model will be used for a single patient, it becomes critical to understand whether multi-indication training enables more robust models with higher generalization capacity or introduces bias into the model due to distribution differences across indications. Hence, we trained single-indication, pan-tissue, and pan-cancer versions of the MOSBY model to investigate the relative merits and potential clinical utility of multi- vs. single-indication models. We discovered that better generalization (test set performance) is achieved with single-indication training in indications where the single-indication performance is moderate to high. In contrast, multi-indication training has potential to increase generalization capacity in indications where single-indication performance is relatively poor. While these conclusions can be made for the median performance across tested signatures and protein/phosphoproteins, significant heterogeneity was observed among these features in terms of showing better performance in the single- vs. multi-indication setting. To prevent potential bias from multi-indication training, we recommend the use of single-indication models where the feature set in the model reflects an unbiased discovery approach. However, an explicit single- vs. multi-indication comparison may be required for settings where the studied feature set reflects a specific biology of interest. A limitation of our study is that MOSBY colocalization maps are not able to capture local spatial patterns. Future work involves the investigation of a graph neural network-based model on tile-level MOSBY predictions where we can capture local as well as slide-level patterns. Moreover, transformer-based architectures may increase the expressiveness of our model to allow a more accurate mapping from image features to bulk omic profiles. METHODS Data TCGA : Batch-normalized RNA sequencing, RPPA datasets as well as clinical and DNA-based data were obtained from the PanCanAtlas publications page of the Genomic Data Commons website ( https://gdc.cancer.gov/about-data/publications/pancanatlas ). H&E-stained slide images (i.e. Tissue Slides) were downloaded from GDC Data Portal ( https://portal.gdc.cancer.gov/ ). Spatially-resolved transcriptomic data : Publicly available breast cancer tissue slides and spatial transcriptomic assays processed by the Spatial Transcriptomics method 39,52 were downloaded from https://data.mendeley.com/datasets/29ntw7sh4r/5 . The dataset contained 68 WSIs from a total of 23 patients along with spot coordinates and respective RNA-seq expression values from the spatial transcriptomics assay. WSIs were tiled into 224x224 patches as input for the MOSBY model. Each tile was centered around the pixel coordinate of an assay spot to represent the 100µm region of the spot. Log-normalized gene expression values were used as ground truth. Signature scores were calculated by first creating an AnnData object using the anndata package then applying the score_genes function from scanpy 53 . All analyses were conducted using Python 3.10. MOSBY preprocessing Image tiling On whole slide images, foreground elements (tissue) and background (glass) were isolated through luminosity-based segmentation. The Python library OpenSlide was used as a backend for generating 224x224px foreground tiles at 0.5mpp resolution. The same tiling protocol was leveraged for both training and inference, allowing tiles to be mapped back to the original WSI positions for visualization (Supplementary Fig. 7). Feature Extraction Contrastive self-supervised learning-based RetCCL 8 was used to extract image features. RetCCL employs a ResNet-50 architecture to extract 2048 features for each image tile. Model Training TCGA A maximum of 8000 tiles were selected randomly for each slide to yield an unbiased representation of the slide, and all tiles were concurrently used for training. A 2-layer perceptron (512 and 256 nodes per layer) was used to map image features to omic variables. Number of epochs was set to a maximum of 300, with early stopping allowed with a patience of 30 epochs. Ground truth RNA-seq data were log-transformed. Ground truth signature and protein levels had negative values, and thus were shifted to make all values nonnegative. The model was trained with MSE loss between prediction and ground truth levels, while Spearman correlation between these quantities was used as early termination criterion in the validation set. A batch size of 64, and AdamW optimizer with 1e-3 weight decay were used. Learning rate scheduler was implemented with step size 5 and gamma 0.9. 5-fold cross-validation was performed (64% training, 16% validation, and 20% test set in each fold) to assess model performance. Information leak was prevented by assigning all WSIs from the same patient to the same partition. For inference, full models were trained with 80% of the data, with the remaining 20% used as validation set to check early stopping criterion. IMvigor210 and IMvigor211 A full model was trained using both IMvigor210 and IMvigor211 images, and largely following the parameters used in TCGA training. Patients were split into 80% training and 20% validation set, stratified by trial and preventing information leak. The signature and gene model consisted of 175 signatures (Supplementary Table 3) and 73 genes (Supplementary Table 2) respectively. Different from TCGA runs, a maximum of 4000 randomly selected tiles per WSI were used during training, and WSIs having less than 200 tiles were filtered out. Ground truth RNA-seq data were log-transformed and standardized to have 0 mean and unit variance. Both gene and signature expression levels were shifted to make all values nonnegative. Computational hardware and software MOSBY was built with the PyTorch library (v1.11.0) in Python (v3.10.8) as a novel implementation of the HE2RNA 12 model. Python libraries used for data processing included NumPy (v1.23.5), Pandas (v1.2.4), Scikit-learn (v0.24.1), OpenSlide (v1.1.1), Zarr (v2.12.0), Tifffile (v2020.10.1), and OpenCV-cv2 (v4.5.5). Whole slide image tiling, RetCCL feature extraction and MOSBY model training were implemented in NVIDIA Tesla V100 Tensor Core GPUs (graphics processing units). Deep learning models were trained with NVIDIA Cuda compiler (v12.1.105). Data visualization in Python was implemented with Matplotlib (v3.3.4) and Seaborn (v0.11.1) libraries. Python statistical analyses such as Spearman and Pearson correlation were implemented with the SciPy library (v1.10.1). Data processing in R (v4.1.1) was implemented with dplyr (v1.0.8), magrittr (v2.0.2), and reshape2 (v1.4.4) libraries. Data visualization in R was performed using ggplot2 (v3.3.5), ggpubr (v0.4.0), and ggsci (v2.9) libraries. Spatial transcriptomics validation analysis Derivation of an empirical null distribution for the distance between ground truth and predicted correlation structures For each slide, the gene-gene correlation matrix was computed for both ground truth and model predictions. The Frobenius distance between these two matrices was used as a test statistic. Subsequently, the MOSBY correlation matrix rows and columns were randomly shuffled, and the Frobenius distance to the ground truth matrix was computed. Random shuffling was repeated 1000 times to obtain an empirical null distribution for MOSBY vs ground truth distance. The test statistic was then compared with the empirical null distribution to obtain a p-value for model predictions. The same sequence of steps was followed for signature-signature correlations. CD8 IHC and H&E whole slide image alignment IMvigor210 and IMvigor211 slides were scanned by an external party (CellCarta, Montreal, QC) on a Pannoramic 250 (3DHistech, Budapest, Hungary) with a 20x or 40x objective. Digital files were transferred to Genentech and converted to the Aperio SVS image format for all viewing and analysis. Hematoxylin and Eosin (HE) slides were aligned with a CD8 stained section from the same block (1017 total sample pairs). These sections were almost never serial, and displayed variable degrees of adjacency characteristic of the block being refaced between sections. Image alignment was performed in Matlab (R2022a, Mathworks, Natick, MA.) by downsampling to ~ 20 um per pixel and converting them to normalized grayscale images, before calculating an affine transformation using mutual information as the underlying matching metric. The transformation was then upscaled before being applied to the original CD8 high magnification image data bringing it into alignment with the HE image. Various alignment metrics were produced (intersection over union, normalized cross correlation) to select for correctly aligned images (454 sample pairs) before a final manual QC check at low magnification (at least 95% of tissue present on both slides, with the majority of visible structures in close proximity) that was not exhaustive and resulted in 42 sample pairs subjected to further analysis. A binary mask was then produced from the aligned CD8 image using both HSV thresholding and a blue-normalized “brownness” algorithm 54 . CD8 IHC tile-level quantification Tile size used for H&E images was set to 224x224 pixels at 0.5mpp. In cases where the native resolution was different from 0.5mpp, tiling was performed with an adjusted tile size at the native resolution, and tiles were subsequently downsampled or upsampled to arrive at 224 pixels at 0.5mpp. In most IMvigor slides, the native resolution was 0.243mpp resulting in an adjusted tile size of approximately 460x460 pixels. CD8 IHC images were tiled using the ‘adjusted’ tile size to match with the corresponding H&E tiles. The CD8 IHC count for a tile was computed using a convolution approach: Each tile was split into 30x30p subtiles, for which the number of 1s (brown pixels) was counted. The counts across subtiles were averaged to obtain the value for the tile. This average value was log-transformed (log(x + 1)) in comparisons with MOSBY model predictions. Computation of signature scores from bulk RNA-seq data In TCGA and IMvigor RNA-seq datasets, normalized gene expression values were log-transformed (log(x + 1)), z-transformed across samples to have 0 mean, unit variance, and subsequently averaged across genes to arrive at a single signature score for each sample. In TCGA, signature scoring was performed on the batch-normalized pan-cancer RNA-seq dataset, and hence signature scores were comparable across cancer indications. Survival analysis TCGA concordance index analysis : MOSBY signature models (full models trained with 80%-20% training/validation split) were used to make inference on all slides in a cancer indication. For a given slide, the colocalization for a signature pair was computed with a Pearson correlation across tile-level predictions of the two signatures. The collection of all pairwise correlation values for a slide formed the ‘slide colocalization matrix’ (N = 175x175). For patients with multiple existing H&E slides, the patient-level colocalization matrix was computed as the average of all pertinent slide-level colocalization matrices. Three sets of cross-validation runs were implemented for L1-regularized Cox regression models. The inputs to these survival models consisted of: Model 1) Flattened patient-level colocalization maps (N = 15225 features). Model 2, patient-level signature scores (N = 175 signatures). Model 3) all tested signatures (N = 175) and lasso-selected colocalization features from Model 1. The input features were regressed against overall survival in all models and for all tested indications. In all three models, the sequence of possible values for lambda (shrinkage parameter) was internally determined in the cv.glmnet function from the glmnet 55 R library (v4.1.3) prior to cross-validation runs. The lambda maximizing the mean Harrel’s concordance measure across 10 cross-validation folds was chosen as optimal, and used to determine the concordance index estimate for the model. The standard error estimate for the model concordance index was calculated across the 10 cross-validation folds. Kaplan-Meier plots The survminer R library (v0.4.9) was used with a median cutoff to generate Kaplan-Meier plots. Log-rank test p-values were obtained internally in survminer using the survdiff function in the survival 56 library (v3.3.1). Data availability The accession URLs for publicly available data analyzed in this study (TCGA, spatial transcriptomics) are listed in the Data section of Methods. Datasets from clinical trials IMpower150 (H&E image data), IMpower110 (H&E image and clinical data), IMvigor210 (RNA-seq, H&E image and CD8 IHC image data), IMvigor211 (RNA-seq, H&E image and CD8 IHC image data) were also analyzed in the current study. IMvigor210 RNA-seq data is available at the European Genome-phenome archive (EGA) under the accession number EGAS00001002556, and was also published as an R package ( http://research-pub.gene.com/IMvigor210CoreBiologies/ ). IMpower110 and IMpower150 datasets are not publicly available as data release is designated for the pending primary biomarker manuscripts. IMvigor210 and IMvigor211 H&E and CD8 IHC image data as well as IMvigor211 RNA-seq data that support the findings of this study are available from Roche, but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of Roche. For up to date details on Roche's Global Policy on the Sharing of Clinical Information and how to request access to related clinical study documents, see here: https://go.roche.com/data_sharing . Declarations Funding: This work was supported by Genentech, Inc. Acknowledgements We would like to thank all of the study participants and their families, and all of the site investigators, study coordinators, and staff. We also would like to thank Brandon Kayser, Robert Johnston, Aïcha Bentaieb, Dan Ruderman, Hector Corrada Bravo, Jason Hackney, and colleagues from the Oncology Reverse Translation team for providing critical feedback on the manuscript. This work was supported by Genentech, Inc. Author contributions K.L., A.K, Y.S. devised the study. Y.S., V.P, A.K, J.E, E.L, E.W., B.N., M.S., M.B. acquired, analyzed, and interpreted the data. Y.S. wrote the manuscript with input from the remaining authors. Competing interests Y.S. , J.E. , E.L , B.N. , M.S. , M.B. , K.L. are employees of Genentech, Inc. and shareholders in F. Hoffmann La Roche, Ltd. V.P. is an external partner at Genentech, Inc. A.K. was an intern at Genentech, Inc. E.W. was an intern at Genentech, In References Comiter, C. et al. Inference of single cell profiles from histology stains with the Single-Cell omics from Histology Analysis Framework (SCHAF) . http://biorxiv.org/lookup/doi/10.1101/2023.03.21.533680 (2023) doi:10.1101/2023.03.21.533680. Alsaafin, A., Safarpoor, A., Sikaroudi, M., Hipp, J. D. & Tizhoosh, H. R. Learning to predict RNA sequence expressions from whole slide images with applications for search and classification. Commun. Biol. 6 , 304 (2023). Yang, K. D. et al. Multi-domain translation between single-cell imaging and sequencing data using autoencoders. Nat. Commun. 12 , 31 (2021). Tsai, P.-C. et al. Histopathology images predict multi-omics aberrations and prognoses in colorectal cancer patients. Nat. Commun. 14 , 2102 (2023). Haviv, D., Gatie, M., Hadjantonakis, A.-K., Nawy, T. & Pe’er, D. The covariance environment defines cellular niches for spatial inference . http://biorxiv.org/lookup/doi/10.1101/2023.04.18.537375 (2023) doi:10.1101/2023.04.18.537375. Chen, X., Fan, H., Girshick, R. & He, K. Improved Baselines with Momentum Contrastive Learning. Preprint at http://arxiv.org/abs/2003.04297 (2020). Oquab, M. et al. DINOv2: Learning Robust Visual Features without Supervision. Preprint at http://arxiv.org/abs/2304.07193 (2023). Wang, X. et al. RetCCL: Clustering-guided contrastive learning for whole-slide image retrieval. Med. Image Anal. 83 , 102645 (2023). Fremond, S. et al. Interpretable deep learning model to predict the molecular classification of endometrial cancer from haematoxylin and eosin-stained whole-slide images: a combined analysis of the PORTEC randomised trials and clinical cohorts. Lancet Digit. Health 5 , e71–e82 (2023). Schirris, Y., Gavves, E., Nederlof, I., Horlings, H. M. & Teuwen, J. DeepSMILE: Contrastive self-supervised pre-training benefits MSI and HRD classification directly from H&E whole-slide images in colorectal and breast cancer. Med. Image Anal. 79 , 102464 (2022). Ciga, O., Xu, T. & Martel, A. L. Self supervised contrastive learning for digital histopathology. Mach. Learn. Appl. 7 , 100198 (2022). Schmauch, B. et al. A deep learning model to predict RNA-Seq expression of tumours from whole slide images. Nat. Commun. 11 , 3877 (2020). Nahhas, O. S. M. E. et al. Regression-based Deep-Learning predicts molecular biomarkers from pathology slides. Preprint at http://arxiv.org/abs/2304.05153 (2023). Hoadley, K. A. et al. Cell-of-Origin Patterns Dominate the Molecular Classification of 10,000 Tumors from 33 Types of Cancer. Cell 173 , 291-304.e6 (2018). Bagaev, A. et al. Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. Cancer Cell 39 , 845-865.e7 (2021). Patil, N. S. et al. Intratumoral plasma cells predict outcomes to PD-L1 blockade in non-small cell lung cancer. Cancer Cell 40 , 289-300.e4 (2022). Liberzon, A. et al. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1 , 417–425 (2015). Miranda, A. et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc. Natl. Acad. Sci. U. S. A. 116 , 9020–9029 (2019). Şenbabaoğlu, Y. et al. Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. Genome Biol. 17 , 231 (2016). Taube, J. H. et al. Core epithelial-to-mesenchymal transition interactome gene-expression signature is associated with claudin-low and metaplastic breast cancer subtypes. Proc. Natl. Acad. Sci. U. S. A. 107 , 15449–15454 (2010). Masiero, M. et al. A core human primary tumor angiogenesis signature identifies the endothelial orphan receptor ELTD1 as a key regulator of angiogenesis. Cancer Cell 24 , 229–241 (2013). Harnoss, J. M. et al. IRE1α Disruption in Triple-Negative Breast Cancer Cooperates with Antiangiogenic Therapy by Reversing ER Stress Adaptation and Remodeling the Tumor Microenvironment. Cancer Res. 80 , 2368–2379 (2020). Gene Ontology Consortium et al. The Gene Ontology knowledgebase in 2023. Genetics 224 , iyad031 (2023). Hobert, O., Carrera, I. & Stefanakis, N. The molecular and gene regulatory signature of a neuron. Trends Neurosci. 33 , 435–445 (2010). Robertson, A. G. et al. Comprehensive Molecular Characterization of Muscle-Invasive Bladder Cancer. Cell 171 , 540-556.e25 (2017). Tsai, H. K. et al. Gene expression signatures of neuroendocrine prostate cancer and primary small cell prostatic carcinoma. BMC Cancer 17 , 759 (2017). Gillespie, M. et al. The reactome pathway knowledgebase 2022. Nucleic Acids Res. 50 , D687–D692 (2022). Xu, Q. et al. Pan-cancer transcriptome analysis reveals a gene expression signature for the identification of tumor tissue origin. Mod. Pathol. Off. J. U. S. Can. Acad. Pathol. Inc 29 , 546–556 (2016). Ben-Porath, I. et al. An embryonic stem cell-like gene expression signature in poorly differentiated aggressive human tumors. Nat. Genet. 40 , 499–507 (2008). Bhattacharya, B., Puri, S. & Puri, R. K. A review of gene expression profiling of human embryonic stem cell lines and their differentiated progeny. Curr. Stem Cell Res. Ther. 4 , 98–106 (2009). Shats, I. et al. Using a stem cell-based signature to guide therapeutic selection in cancer. Cancer Res. 71 , 1772–1780 (2011). Kim, J. et al. A Myc network accounts for similarities between embryonic stem and cancer cell transcription programs. Cell 143 , 313–324 (2010). Possemato, R. et al. Functional genomics reveal that the serine synthesis pathway is essential in breast cancer. Nature 476 , 346–350 (2011). Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12 , 453–457 (2015). Mariathasan, S. et al. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554 , 544–548 (2018). Böttcher, J. P. & Reis e Sousa, C. The Role of Type 1 Conventional Dendritic Cells in Cancer Immunity. Trends Cancer 4 , 784–792 (2018). Yoshihara, K. et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4 , 2612 (2013). Akbani, R. et al. A pan-cancer proteomic perspective on The Cancer Genome Atlas. Nat. Commun. 5 , 3887 (2014). Stenbeck, L., Bergenstråhle, L., Lundeberg, J. & Borg, Å. Human breast cancer in situ capturing transcriptomics. (2021) doi:10.17632/29ntw7sh4r.5. Balar, A. V. et al. Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. Lancet Lond. Engl. 389 , 67–76 (2017). Powles, T. et al. Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma (IMvigor211): a multicentre, open-label, phase 3 randomised controlled trial. Lancet Lond. Engl. 391 , 748–757 (2018). Diao, J. A. et al. Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes. Nat. Commun. 12 , 1613 (2021). Socinski, M. A. et al. Atezolizumab for First-Line Treatment of Metastatic Nonsquamous NSCLC. N. Engl. J. Med. 378 , 2288–2301 (2018). Chen, R. J. et al. Pan-cancer integrative histology-genomic analysis via multimodal deep learning. Cancer Cell 40 , 865-878.e6 (2022). Lu, M. Y. et al. Data-efficient and weakly supervised computational pathology on whole-slide images. Nat. Biomed. Eng. 5 , 555–570 (2021). Helmink, B. A. et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature 577 , 549–555 (2020). Trüb, M. & Zippelius, A. Tertiary Lymphoid Structures as a Predictive Biomarker of Response to Cancer Immunotherapies. Front. Immunol. 12 , 674565 (2021). Katzman, J. L. et al. DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network. BMC Med. Res. Methodol. 18 , 24 (2018). Javed, S. A. et al. Additive MIL: Intrinsically Interpretable Multiple Instance Learning for Pathology. (2022) doi:10.48550/ARXIV.2206.01794. Chowell, D. et al. Improved prediction of immune checkpoint blockade efficacy across multiple cancer types. Nat. Biotechnol. 40 , 499–506 (2022). Chen, R. J. et al. Whole Slide Images are 2D Point Clouds: Context-Aware Survival Prediction using Patch-based Graph Convolutional Networks. (2021) doi:10.48550/ARXIV.2107.13048. Ståhl, P. L. et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science 353 , 78–82 (2016). Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19 , 15 (2018). Brey, E. M. et al. Automated selection of DAB-labeled tissue for immunohistochemical quantification. J. Histochem. Cytochem. Off. J. Histochem. Soc. 51 , 575–584 (2003). Simon, N., Friedman, J., Tibshirani, R. & Hastie, T. Regularization Paths for Cox’s Proportional Hazards Model via Coordinate Descent. Journal of Statistical Software 39 , 1--13 (2011). Therneau, T. M. & Grambsch, P. M. Modeling Survival Data: Extending the Cox Model . (Springer, 2000). Additional Declarations Competing interest reported. Y.S. , J.E. , E.L , B.N. , M.S. , M.B. , K.L. are employees of Genentech, Inc. and shareholders in F. Hoffmann La Roche, Ltd. V.P. is an external partner at Genentech, Inc. A.K. was an intern at Genentech, Inc. E.W. was an intern at Genentech, Inc. Supplementary Files Supplementarylegends.docx SUPPFIGURE1.pdf SUPPFIGURE2.jpg SUPPFIGURE3.pdf SUPPFIGURE4.pdf SUPPFIGURE5.pdf SUPPFIGURE6.pdf SUPPFIGURE7.pdf SupplementaryTable1.xlsx SupplementaryTable2.xlsx SupplementaryTable3.xlsx SupplementaryTable4.xlsx Cite Share Download PDF Status: Published Journal Publication published 06 Aug, 2024 Read the published version in Scientific Reports → Version 1 posted Editorial decision: Revision requested 17 May, 2024 Reviews received at journal 15 May, 2024 Reviewers agreed at journal 05 May, 2024 Reviews received at journal 17 Mar, 2024 Reviewers agreed at journal 08 Mar, 2024 Reviewers agreed at journal 07 Mar, 2024 Reviewers invited by journal 05 Mar, 2024 Editor assigned by journal 05 Mar, 2024 Editor invited by journal 28 Feb, 2024 Submission checks completed at journal 28 Feb, 2024 First submitted to journal 07 Feb, 2024 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-3938444","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":275455079,"identity":"62705584-dbf6-444d-bf90-aec741fe3529","order_by":0,"name":"Yasin Senbabaoglu","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA1ElEQVRIiWNgGAWjYDACCQYGZgaGAwz8IE5CASlaJBtAWgxI0WJwAMQjRov87OZjjwtq7sgbn1+d+OGBAYM8v9gB/FoM7hxLN55x7JnhthtvN0sAHWY4c3YCAS0SOWbSPGyHGbfdOLsBpCXB4DYBLfIz8r9J8/w7bL95xtnNP4jSwnAjh02at+1w4gb+3m3E2WJwI81Membfs+QZN3i3WSQYSBD2i/yM5GfSBd/u2Pb3n91880eFjTy/NCGHwYEEWKUEscpBgP8AKapHwSgYBaNgJAEAKJ5ILMZXkCQAAAAASUVORK5CYII=","orcid":"","institution":"Genentech","correspondingAuthor":true,"prefix":"","firstName":"Yasin","middleName":"","lastName":"Senbabaoglu","suffix":""},{"id":275455080,"identity":"93077c03-bc68-4dc2-99c5-863f764372bb","order_by":1,"name":"Vignesh Prabhakar","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Vignesh","middleName":"","lastName":"Prabhakar","suffix":""},{"id":275455081,"identity":"042b6438-99ed-42c2-9e88-948521a9c279","order_by":2,"name":"Aminollah Khormali","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Aminollah","middleName":"","lastName":"Khormali","suffix":""},{"id":275455082,"identity":"b88fafe3-46ab-45b8-88d3-4f1d84d8528a","order_by":3,"name":"Jeff Eastham","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Jeff","middleName":"","lastName":"Eastham","suffix":""},{"id":275455083,"identity":"3e79cbaf-25e3-4c69-bf35-bebc00a601e9","order_by":4,"name":"Elisa Warner","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Elisa","middleName":"","lastName":"Warner","suffix":""},{"id":275455084,"identity":"40250b5c-8c6c-4493-8676-bdb58cbadf33","order_by":5,"name":"Evan Liu","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Evan","middleName":"","lastName":"Liu","suffix":""},{"id":275455085,"identity":"cd5c959c-9f28-4bf9-a78a-173c25cfbbb3","order_by":6,"name":"Barzin Nabet","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Barzin","middleName":"","lastName":"Nabet","suffix":""},{"id":275455086,"identity":"222ca1f5-c82f-4dc7-886b-5f3b9b65afa4","order_by":7,"name":"Minu Srivastava","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Minu","middleName":"","lastName":"Srivastava","suffix":""},{"id":275455087,"identity":"f8192157-7380-4e8a-8c4c-201c58c3b06c","order_by":8,"name":"Marcus Ballinger","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Marcus","middleName":"","lastName":"Ballinger","suffix":""},{"id":275455088,"identity":"4719403b-461c-4b64-b31b-791c28babde3","order_by":9,"name":"Kai Liu","email":"","orcid":"","institution":"Genentech","correspondingAuthor":false,"prefix":"","firstName":"Kai","middleName":"","lastName":"Liu","suffix":""}],"badges":[],"createdAt":"2024-02-08 01:44:06","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3938444/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3938444/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1038/s41598-024-69198-6","type":"published","date":"2024-08-06T15:56:54+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":51973021,"identity":"f0ea093d-9740-4e5b-9cfd-3860f5d2b0f1","added_by":"auto","created_at":"2024-03-04 18:59:03","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":2846614,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ea) \u003c/strong\u003eMOSBY flowchart: H\u0026amp;E-stained whole slide images were downloaded from TCGA and partitioned into tiles. Image features were extracted using RetCCL-pretrained weights and the ResNet-50 architecture. A 2-layer perceptron was trained to learn the mapping from image features to bulk omic profiles. A by-product of the process, tile level predictions were utilized to achieve virtual spatialization of omic profiles. Pairwise correlation of omic features resulted in a colocalization map for each patient. Colocalization maps were then flattened, and used as covariates in survival analysis with the goal of spatial biomarker discovery. \u003cstrong\u003eb-c-d-e) \u003c/strong\u003eTest set Spearman correlation between omic feature predictions and ground truth, averaged over 5 cross-validation folds. Each data point is an omic feature, and the x-axis shows different single- and multi-indication runs. Red and turquoise boxes indicate results from RetCCL vs. ImageNet-pretrained ResNet-50\u003cstrong\u003e \u003c/strong\u003earchitecture respectively. Results are shown for \u003cstrong\u003eb)\u003c/strong\u003e 55 TME-related genes, \u003cstrong\u003ec) \u003c/strong\u003e175 tumor cell and microenvironment-related signatures, \u003cstrong\u003ed) \u003c/strong\u003e191 reverse phase protein array antibodies, and \u003cstrong\u003ee) \u003c/strong\u003e4 DNA-based features. Mann-Whitney hypothesis tests for (b), (c), (d) were implemented with the ggpubr compare_means function to compare RetCCL and ImageNet based results (Hypothesis tests were not performed for (e) due to the small number of features, N=4). Significance levels are provided in the figure, original P-values are provided in Supplementary Table 4. ****p\u0026lt;0.0001, ***p\u0026lt;0.001, **p\u0026lt;0.01, *p\u0026lt;0.05, ns: p\u0026gt;0.05. Lower and upper hinges in box plots correspond to first and third quartiles, while whiskers extend to 1.5 times interquartile range. \u003cstrong\u003ef-g) \u003c/strong\u003eConcordance between ST ground truth and MOSBY-predicted levels of select \u003cstrong\u003ec)\u003c/strong\u003e signatures, \u003cstrong\u003ed)\u003c/strong\u003e genes.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/1c8564f8d0be1f935bdb1e4c.png"},{"id":51974661,"identity":"e3442804-a742-4d9e-bf3e-ec9c1e7300d1","added_by":"auto","created_at":"2024-03-04 19:15:04","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":3629138,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eValidation of MOSBY \u003c/strong\u003e\u003cem\u003e\u003cstrong\u003ein silico \u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003espatialization with CD8 IHC whole slide\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eImages in IMvigor211. a) \u003c/strong\u003eWorkflow for computational alignment of serially sectioned CD8 IHC and H\u0026amp;E whole slide images, and correlation of CD8 IHC-based ground truth values with MOSBY tile-level predictions. The derivation of tile level ground truth values by applying DAB mask on CD8 IHC images is described in \u003cstrong\u003eMethods\u003c/strong\u003e. \u003cstrong\u003eb)\u003c/strong\u003e Pearson correlation between tile-level ground truth (CD8 IHC) and MOSBY-predicted values for the 42 slides that satisfied quality control criteria in \u003cstrong\u003e(a)\u003c/strong\u003e. The x-axis shows gene (\u003cem\u003eCD8A, CD8B\u003c/em\u003e) and signature (Cibersort CD8 T cell, T effector cell) features compared with CD8 IHC. For a given slide, the feature with the highest CD8 IHC correlation is denoted with a red dot, while the other three features are denoted with blue dots. The number of slides where each feature had the highest CD8 IHC correlation is shown in parenthesis under the x-axis labels. Lower and upper hinges in box plots correspond to first and third quartiles, while whiskers extend to 1.5 times interquartile range. \u003cstrong\u003ec-j) \u003c/strong\u003eDescription of results for the slide where MOSBY predictions are the most concordant with CD8 IHC-based ground truth levels: \u003cstrong\u003ec) \u003c/strong\u003eDensity plot showing the correlation between tile-level ground truth data (DAB mask applied on CD8 IHC WSI) and MOSBY \u003cem\u003eCD8B\u003c/em\u003epredictions (Pearson r = 0.79, N=16859 tiles, two-sided p ≈ 0 using the exact beta distribution of r). \u003cstrong\u003ed)\u003c/strong\u003e Computationally aligned CD8 IHC and H\u0026amp;E whole slide images. \u003cstrong\u003ee)\u003c/strong\u003e Heat map showing CD8 IHC tile level quantification. DAB mask is applied on CD8 IHC images to call ‘brown’ pixels. Brown pixels are counted within each 100x100px window. Count values are log-transformed, and then clipped at 50th and 98th percentiles for contrast. Background regions are denoted with gray. \u003cstrong\u003ef)\u003c/strong\u003e Heat map showing \u003cem\u003eCD8B\u003c/em\u003e tile-level model predictions. Predicted values are clipped at 50th and 98th percentiles for contrast. Background regions are denoted with gray. \u003cstrong\u003eg)\u003c/strong\u003eHigh magnification example of a CD8 T cell-rich region on the CD8 IHC image. \u003cstrong\u003eh)\u003c/strong\u003e DAB mask classification of ‘brown’ pixels (denoted with green) overlaid on CD8 IHC image. \u003cstrong\u003ei)\u003c/strong\u003e Corresponding region on the H\u0026amp;E image as identified by slide-level computational alignment. \u003cstrong\u003ej) \u003c/strong\u003e\u003cem\u003eCD8B \u003c/em\u003emodel predictions inferred from the H\u0026amp;E image in \u003cstrong\u003e(i)\u003c/strong\u003e overlaid on the CD8 IHC image in \u003cstrong\u003e(g). \u003c/strong\u003eOrange and blue show high and low values respectively.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/c27401695b58ddb8316844b8.png"},{"id":51973026,"identity":"a95fd5ac-7b58-45f6-a168-03e99ede08ec","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":2885981,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMOSBY prediction of cancer cell regions using single gene, signature, and protein models. a)\u003c/strong\u003eH\u0026amp;E-stained WSI for three representative slides with pathologist penmark annotations indicating cancer cells. \u003cstrong\u003eb-c-d-e)\u003c/strong\u003eSpatialization of MOSBY-predicted expression levels of \u003cstrong\u003eb)\u003c/strong\u003e tumor purity, \u003cstrong\u003ec)\u003c/strong\u003eepithelial markers \u003cem\u003eEPCAM\u003c/em\u003e, \u003cstrong\u003ed)\u003c/strong\u003e E-cadherin, and \u003cstrong\u003ee)\u003c/strong\u003e Taube \u003cem\u003eet al.\u003c/em\u003eepithelial signature.\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/212d0032c73a31318b94a834.png"},{"id":51973029,"identity":"9a93c41e-eb57-4293-a8a0-a1e63a2f3a39","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":1368680,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eBest performing features in MOSBY signature, protein, gene, and DNA panels. a) \u003c/strong\u003eTop 25 signatures with best prediction accuracy according to cross-validated test set performance in single-indication models. \u003cstrong\u003eb)\u003c/strong\u003e Test set correlations of select signatures in a single split run (70% training, 15% validation, and 15% test sets). \u003cstrong\u003ec) \u003c/strong\u003ePercentage of signatures predicted with significant ground truth correlation (adjusted p-value \u0026lt; 0.05).\u003cstrong\u003e d) \u003c/strong\u003eTop 25 proteins/phosphoproteins with best prediction accuracy according to cross-validated test set performance in single-indication models. \u003cstrong\u003ee)\u003c/strong\u003e Test set correlations of select proteins/phosphoproteins in a single split run (70% training, 15% validation, and 15% test sets). \u003cstrong\u003ef) \u003c/strong\u003ePercentage of proteins/phosphoproteins predicted with significant ground truth correlation (adjusted p-value \u0026lt; 0.1). \u003cstrong\u003eg) \u003c/strong\u003eTop 25 genes with best prediction accuracy according to cross-validated test set performance in single-indication models. \u003cstrong\u003eh)\u003c/strong\u003eCross-validated test set performance of all tested DNA-based features. Lower and upper hinges in box plots correspond to first and third quartiles, while whiskers extend to 1.5 times interquartile range.\u003c/p\u003e","description":"","filename":"floatimage5.png","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/eea71045bc5c9a55982d025a.png"},{"id":51973030,"identity":"63ad26c2-b47d-47a3-a0b9-08a40b64b457","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1343420,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSurvival predictive power comparison between gene signatures and colocalization features.\u003c/strong\u003e \u003cstrong\u003ea) \u003c/strong\u003eConcordance indices of survival models involving only colocalization features (denoted with C), only signatures (denoted with S), or signatures and lasso-selected colocalization features (denoted with C+S) in 21 TCGA indications. Indications are ordered according to increasing concordance index in signature-only models. Error bars denote +/-1 standard error. \u003cstrong\u003eb) \u003c/strong\u003eColocalization features most consistently associated with poor overall survival and elevated tumor levels in multiple TCGA indications. \u003cstrong\u003ec-d)\u003c/strong\u003eSurvival and differential tumor colocalization evidence for the ER17 vs. Neurotransmitter signature pair in TCGA colon, liver, lung, and ovarian cancer cohorts. \u003cstrong\u003ec)\u003c/strong\u003e Median-split Kaplan-Meier plots and logrank test p-values, \u003cstrong\u003ed)\u003c/strong\u003e tumor vs normal box plots and Wilcoxon test p-values. \u003cstrong\u003ee-f)\u003c/strong\u003e Median-split Kaplan-Meier plots and logrank test p-values for the colocalization between ER17 and Neurotransmitter signatures in IMpower110 non-squamous cohort: \u003cstrong\u003ee)\u003c/strong\u003e Chemotherapy arm,\u003cstrong\u003e f) \u003c/strong\u003eAtezolizumab arm. Lower and upper hinges in box plots correspond to first and third quartiles, while whiskers extend to 1.5 times interquartile range. Red and blue color denote high and low levels respectively.\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/0ea60b194a14d5d8a54380e0.png"},{"id":51973024,"identity":"152713cd-8acd-4c58-8511-5bf28c116f26","added_by":"auto","created_at":"2024-03-04 18:59:03","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":1938504,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eImmune system-related spatial biomarkers consistently associated with risk. a-b) \u003c/strong\u003eSurvival and differential tumor colocalization evidence for T effector cell and Sabatini Cysteine signature pair in TCGA breast, lung squamous, and ovarian cancer cohorts: \u003cstrong\u003ea)\u003c/strong\u003e Median-split Kaplan-Meier plots and logrank test p-values, \u003cstrong\u003eb)\u003c/strong\u003e tumor vs normal box plots and Wilcoxon test p-value. Lower and upper hinges in box plots correspond to first and third quartiles, while whiskers extend to 1.5 times interquartile range. \u003cstrong\u003ec-d-e-f) \u003c/strong\u003eMedian-split Kaplan-Meier plots and logrank test p-values for Sabatini Cysteine colocalization with various immune cell populations: \u003cstrong\u003ec-d) \u003c/strong\u003eTCGA breast cancer cohort: \u003cstrong\u003ec)\u003c/strong\u003e Lymphocytes and NK cells, \u003cstrong\u003ed)\u003c/strong\u003emyeloid cell populations. \u003cstrong\u003ee-f) \u003c/strong\u003eAtezolizumab arm of IMpower110 non-squamous cohort (N=137): \u003cstrong\u003ee)\u003c/strong\u003e Lymphocytes and NK cells, \u003cstrong\u003ef)\u003c/strong\u003emyeloid cell populations. Red and blue color denote high and low levels respectively.\u003c/p\u003e","description":"","filename":"floatimage7.png","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/d24a364b225b260b75b68204.png"},{"id":62298110,"identity":"72c922be-f0f0-4a77-ad62-2ddaf8318742","added_by":"auto","created_at":"2024-08-12 16:08:36","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":18704779,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/b827b1d3-4f06-4d76-b4b0-38196ddc7d7d.pdf"},{"id":51973023,"identity":"a4b764d0-95c5-4c1d-8952-e0e109965941","added_by":"auto","created_at":"2024-03-04 18:59:03","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":29458,"visible":true,"origin":"","legend":"","description":"","filename":"Supplementarylegends.docx","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/f042c5d6baa1f5bd82524dfd.docx"},{"id":51973022,"identity":"c82ed044-2015-491f-95f1-54a38a1b7689","added_by":"auto","created_at":"2024-03-04 18:59:03","extension":"pdf","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":4864992,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/293b8cc59ad7af7b9144f9db.pdf"},{"id":51973028,"identity":"c041a717-52e7-47df-ae43-18596b634dab","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"jpg","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":10187607,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/50e606b808882c5c5e2b71f6.jpg"},{"id":51973040,"identity":"92be4f21-fee8-437b-8a50-f23e56f0dd52","added_by":"auto","created_at":"2024-03-04 18:59:05","extension":"pdf","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":14510787,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE3.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/c5d08d3ad9ca51bc3ad4fff1.pdf"},{"id":51974344,"identity":"48e17735-dbb4-43cd-be95-6cc12b4ff33f","added_by":"auto","created_at":"2024-03-04 19:07:04","extension":"pdf","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":4628910,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE4.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/141c6145ea02adc6fe610e8f.pdf"},{"id":51973025,"identity":"1273d755-b14c-485f-927d-0c1c228079e3","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"pdf","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":1589875,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE5.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/83fa3a0b2e0e43615dcd5ee1.pdf"},{"id":51973037,"identity":"4a297432-9e6c-4081-87e6-139719c2a852","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"pdf","order_by":7,"title":"","display":"","copyAsset":false,"role":"supplement","size":1799861,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE6.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/b0030ff189d238a1194c13c8.pdf"},{"id":51973033,"identity":"43826a8a-8428-4b8f-a2db-c74a6a6aa4b3","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"pdf","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":794165,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPFIGURE7.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/535f06f774a454782336084d.pdf"},{"id":51973031,"identity":"d85bde05-c4cc-4ed5-950e-c13b05ccfe41","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"xlsx","order_by":9,"title":"","display":"","copyAsset":false,"role":"supplement","size":9598,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/ad04c5a74a5cea5faaacc306.xlsx"},{"id":51973035,"identity":"2f4e4c7f-9e30-4a7b-a96e-b87ff77eafce","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"xlsx","order_by":10,"title":"","display":"","copyAsset":false,"role":"supplement","size":11722,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/4af47a394e779b0264d3fbaa.xlsx"},{"id":51973038,"identity":"73d42bcb-55ff-4bf7-9622-1f19bc0f2852","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"xlsx","order_by":11,"title":"","display":"","copyAsset":false,"role":"supplement","size":564121,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable3.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/7c790a038b928d9b47110c18.xlsx"},{"id":51973036,"identity":"17334471-2559-4549-9598-67ca0ba64ab8","added_by":"auto","created_at":"2024-03-04 18:59:04","extension":"xlsx","order_by":12,"title":"","display":"","copyAsset":false,"role":"supplement","size":17553,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTable4.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-3938444/v1/c3ac3ce62c428017f1bbb8ca.xlsx"}],"financialInterests":"Competing interest reported. Y.S. , J.E. , E.L , B.N. , M.S. , M.B. , K.L. are employees of Genentech, Inc. and shareholders in F. Hoffmann La Roche, Ltd.\nV.P. is an external partner at Genentech, Inc. \nA.K. was an intern at Genentech, Inc. \nE.W. was an intern at Genentech, Inc.","formattedTitle":"MOSBY enables multi-omic inference and spatial biomarker discovery from whole slide images","fulltext":[{"header":"Introduction","content":"\u003cp\u003eThe developments in high-throughput RNA, DNA, and protein assays as well as in digital pathology have enabled multiple different high-resolution perspectives on a given tissue. Sampled from the same tissue, these data modalities are inherently linked, due to arising from highly similar cell populations and potentially the same biological state. Successfully inferring this common biology is key for relating different data modalities, and predicting one modality from another. Deep neural nets have emerged as an effective and flexible framework to represent the common biology underlying different readouts, opening the way to \u0026lsquo;translate\u0026rsquo; any data modality into another\u003csup\u003e1\u0026ndash;5\u003c/sup\u003e. A high utility use case for this translation task involves digital pathology where relatively easy-to-acquire hematoxylin and eosin (H\u0026amp;E) stained whole slide images (WSI) can be used to infer high-throughput molecular data that is time consuming or expensive to acquire.\u003c/p\u003e \u003cp\u003eThe success of inferring the joint biology (i.e. underlying latent distribution) between WSIs and molecular data relies heavily on extracting informative features from WSIs. Contrastive self-supervised learning has made a breakthrough in computer vision by learning high quality image representations in no-annotation settings as shown in curated benchmark datasets such as ImageNet \u003csup\u003e6,7\u003c/sup\u003e. Contrastive learning aims to discriminate between positive and negative images where positive pairs are obtained from augmentations from the same image and negative pairs are augmentations from different images \u003csup\u003e8\u003c/sup\u003e. Contrastive self-supervised learning has also shown promise in identifying cancer-specific morphological features in large-scale and heterogeneous histology datasets and training feature extractor models that outperform ImageNet-pretrained networks \u003csup\u003e9\u0026ndash;11\u003c/sup\u003e. However, the superior performance of self-supervised features has largely been demonstrated for classification tasks, with limited attention on regression tasks. For instance, Wang \u003cem\u003eet al.\u003c/em\u003e\u003csup\u003e8\u003c/sup\u003e developed an oncology-focused feature extractor RetCCL that was trained in more than 34,000 histology images and achieved state-of-the-art performance in cancer subtype classification. Regression tasks, such as inferring gene expression levels from WSI features, have been addressed with multi-output regression networks \u003csup\u003e1,12\u003c/sup\u003e; yet, self-supervised learning-based image features remain largely underexplored in this context\u003csup\u003e13\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eAmong molecular data modalities, the inference of bulk and single-cell transcriptomics from H\u0026amp;E WSI features has received considerable attention\u003csup\u003e1,2,12\u003c/sup\u003e; however, an underexplored extension of this question is whether mapping the same image features to proteomic data holds greater promise due to proteins being more proximally associated with phenotype. The full power and limitations of image to molecular data \u0026lsquo;translations\u0026rsquo; will become more evident as these deep learning models are expanded to learn multiple -omic data modalities such as transcriptomics, proteomics, metabolomics, and whole exome or SNP array-based DNA measurements.\u003c/p\u003e \u003cp\u003eHere, we present the MOSBY (\u003cb\u003eM\u003c/b\u003eulti-\u003cb\u003eOmic\u003c/b\u003e translation of whole slide images for \u003cb\u003eS\u003c/b\u003epatial \u003cb\u003eB\u003c/b\u003eiomarker discover\u003cb\u003eY\u003c/b\u003e) model that employs contrastive self-supervised learning-based (RetCCL) features from H\u0026amp;E images to infer high-throughput molecular data such as transcriptomics (single genes or gene signatures), proteomics, and selected DNA-based features (e.g. cancer DNA fraction, average DNA methylation). MOSBY adopts a multiple instance learning approach where it partitions gigapixel WSIs to small tiles, makes tile-level predictions as an intermediate step, and then aggregates those to obtain slide-level predictions. At the inference stage, MOSBY enables \u003cem\u003ein silico\u003c/em\u003e spatialization of bulk omic profiles by reconstructing tiles from a given WSI. Correlation between two omic features across all tiles of a slide reveals patient-specific colocalization or spatial exclusion patterns. Expanded to the cohort, these spatial features allow identification of clinically relevant spatial biomarkers. We demonstrate the performance of MOSBY in pan-cancer TCGA data, and spatially validate results with 1) spatially resolved transcriptomic data in breast cancer, 2) serially-sectioned CD8 IHC images in bladder cancer, and 3) tumor microenvironment (TME) cell type predictions in nonsmall cell lung cancer from an orthogonal deep learning model. In contrast to a whole transcriptome approach, we focus on TME cell type marker genes for MOSBY gene models, whereby model predictions could be contrasted with validation data in (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e). Signature models, on the other hand, aim to encompass a wider biology including TME cell type, metabolic, and oncogenic pathway signatures. We finally utilize signature models to showcase the derivation of biologically interpretable colocalization features consistently associated with risk and disease state in human cancers.\u003c/p\u003e"},{"header":"RESULTS","content":"\u003cp\u003eWe developed the MOSBY model that learned a mapping from RetCCL (contrastive self-supervised learning)-based WSI features to bulk transcriptome, proteome, whole exome, SNP-array, and DNA methylome based profiles (Fig.\u0026nbsp;1a). We limited our study to 21 TCGA indications\u003csup\u003e14\u003c/sup\u003e with at least 200 paired RNA-seq and H\u0026amp;E whole slide images, resulting in a total of 12592, 10192, and 12090 images matching with transcriptome, proteome and DNA-based data respectively (breakdown by indication in Supplementary Table\u0026nbsp;1). Analyzed transcriptomic features consisted of 55 TME-related genes\u003csup\u003e15,16\u003c/sup\u003e (Supplementary Table\u0026nbsp;2) and 175 gene sets that covered tumor-related processes\u003csup\u003e17\u0026ndash;32\u003c/sup\u003e, metabolic pathways\u003csup\u003e33\u003c/sup\u003e, and TME cell type or process signatures\u003csup\u003e34\u0026ndash;37\u003c/sup\u003e (Supplementary Table\u0026nbsp;3). The proteome model involved 191 total and phosphoprotein antibodies from the TCGA reverse-phase protein array (RPPA) panel that focused on tumor-intrinsic and oncogenic processes \u003csup\u003e38\u003c/sup\u003e. Tested DNA-based features were limited to tumor purity, cancer DNA fraction, subclonal genome fraction, and average DNA methylation (all continuous measures bound to the interval [0,1]).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn addition to single indication models, we also trained MOSBY in pan-tissue, pan-cancer, pan-squamous and pan-adenocarcinoma settings. Pan-tissue models consisted of LUNG (lung adenocarcinoma and squamous cell carcinoma), KIDNEY (clear cell and papillary renal cell carcinoma), BRAIN (glioblastoma and low grade glioma), and COADREAD (colon and rectal adenocarcinoma). The pan-adenocarcinoma (ADENO) model consisted of pancreatic, lung, stomach, colon, rectal, prostate, and ovarian adenocarcinomas; while the pan-SQUAMOUS model included squamous cell carcinomas of the lung, cervix, and head and neck. The ADENO and SQUAMOUS models allowed us to investigate whether histological similarities beyond tissue architecture enabled MOSBY to learn a better mapping from WSI features to multi-omic data.\u003c/p\u003e \u003cp\u003eSimilar to the HE2RNA\u003csup\u003e12\u003c/sup\u003e model, a 2-layer perceptron was adopted as the multi-output regression network that mapped image features to omic variables, and a maximum of 8000 image tiles per slide were used. Separate models were trained for gene, signature, protein, and DNA variables. In contrast to HE2RNA, image tiles were randomly selected from each WSI to capture an unbiased representation of the entire slide, and were all used in training without clustering. In addition, batch-normalized transcriptomic and proteomic data were used as ground truth to enable across-indication comparisons with resulting MOSBY predictions. To obtain slide-level predictions, tile-level predictions of omic features were aggregated by averaging.\u003c/p\u003e \u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eContrastive self-supervised pretraining benefits prediction of omic data from H\u0026amp;E whole slide images\u003c/h2\u003e \u003cp\u003eThe RetCCL feature extractor utilized TCGA H\u0026amp;E images during contrastive training, however it was not supervised by gene expression information. Hence, there are no \u003cem\u003ea priori\u003c/em\u003e guarantees for RetCCL-based image features to predict gene expression with higher accuracy than ImageNet-based features. The success in this task hinges upon the diversity of semantic (in this case, biological) features in the representations learned by RetCCL. Thus we investigated the performance differences between MOSBY models trained with RetCCL- or ImageNet-based image features. Models were trained with 5-fold cross-validation (80/20 percent training/test set split). A random one fifth of the training set was also allocated as validation set to determine an early stopping criterion using the Spearman correlation between slide-level model prediction and ground truth omic data. Spearman correlation was preferred to mitigate the effects of outlier samples in relatively smaller cohorts. Correlation coefficients from test sets were subsequently averaged across five folds to obtain the final performance score for a given feature.\u003c/p\u003e \u003cp\u003eRetCCL-based image features consistently led to higher cross-validated test set averages in all four omic data types compared to features extracted with ImageNet-pretrained ResNet-50 architecture (Fig.\u0026nbsp;1b-1e). The pan-cancer model (PANCAN) with RetCCL-based features achieved the highest performance for all tested data types with median cross-validated Spearman correlation of 0.722 in single genes (0.673 with ImageNet-features) (Fig.\u0026nbsp;1b), 0.693 in signatures (0.647 with ImageNet-features) (Fig.\u0026nbsp;1c), 0.549 in proteomic data (0.512 with ImageNet-features) (Fig.\u0026nbsp;1d), and 0.6 in DNA features (0.533 with ImageNet-features) (Fig.\u0026nbsp;1e). In terms of single indication models, the thyroid cancer model (THCA) achieved best performance for single gene and protein expression data sets with RetCCL-based features. Median cross-validated Spearman correlations in these models were 0.59 vs 0.524 for single gene, and 0.443 vs 0.385 for protein expression data with RetCCL and ImageNet-based features respectively (Fig.\u0026nbsp;1b, 1d). The single indication models achieving best performance for signature and DNA data were the liver and bladder cancer models respectively (LIHC and BLCA) (Fig.\u0026nbsp;1c, 1e). Across tested signatures, the LIHC model showed a median cross-validated Spearman correlation of 0.532 with RetCCL-based vs 0.45 with ImageNet-based features. For tested DNA features, the BLCA model achieved a median cross-validated Spearman correlation of 0.553 with RetCCL-based vs 0.417 with ImageNet-based features. Taken together, these results indicated that contrastive self-supervised pretraining in large-scale histology datasets benefits prediction of multiomic data from H\u0026amp;E-stained whole slide images.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003eMOSBY tile level predictions are validated with spatially resolved transcriptomic data\u003c/h2\u003e \u003cp\u003eFor inference tasks, MOSBY \u0026lsquo;full models\u0026rsquo; were trained with data from 80% of patients, with the remaining 20% used as the validation set to determine an early stopping criterion. An MOSBY full model was trained in the TCGA breast cancer (BRCA) cohort (N\u0026thinsp;=\u0026thinsp;1576 WSIs), and subsequently deployed on H\u0026amp;E image tiles from a publicly available spatially resolved transcriptomic breast cancer dataset\u003csup\u003e39\u003c/sup\u003e (referred to as ST from here on, N\u0026thinsp;=\u0026thinsp;68 WSIs) that served as an independent validation cohort. Image tiles were centered around ST spot coordinates to enable one-to-one comparison between spot level ground truth data and tile level model predictions (\u003cb\u003eMethods\u003c/b\u003e). For individual slides, the MOSBY signature model predictions showed the highest concordance with ground truth for \u0026lsquo;poor differentiation\u0026rsquo; (Stemness_Kim_Myc, Spearman r\u0026thinsp;=\u0026thinsp;0.71) and stromal features (Stroma_Estimate, Spearman r\u0026thinsp;=\u0026thinsp;0.63) (Fig.\u0026nbsp;1f, Supplementary Fig.\u0026nbsp;1a). Across 68 slides, concordance was highest for a monocyte signature with a median Spearman correlation of 0.238 (Supplementary Fig.\u0026nbsp;1a), and a maximum of 0.611 (Fig.\u0026nbsp;1f). This large variation across slides was observed for all tested signatures suggesting that the quality of spatially resolved data is critical in validating MOSBY predictions. Of note, CD8 T cell infiltration and proliferation-related model predictions also showed strong concordance with ground truth, demonstrating the variety in phenotypes captured successfully by the model (Supplementary Fig.\u0026nbsp;1b).\u003c/p\u003e \u003cp\u003eSpot level concordance for single gene expression values was overall lower than that for signatures (Supplementary Fig.\u0026nbsp;1d), potentially driven by the large degree of zero reads (i.e. dropouts) in ground truth data (Supplementary Fig.\u0026nbsp;1e). However, model predictions for genes associated with stroma, plasma cells, and epithelial features showed good performance for individual slides (\u003cem\u003eCOL1A2, MZB1, EPCAM\u003c/em\u003e respectively) (Fig.\u0026nbsp;1g). \u003cem\u003eCD68\u003c/em\u003e, a myeloid marker, showed the highest median concordance in gene models (Spearman r\u0026thinsp;=\u0026thinsp;0.127), which was low in magnitude but consistent with the high concordance of myeloid cells in signature models (Supplementary Fig.\u0026nbsp;1d).\u003c/p\u003e \u003cp\u003eAs MOSBY aims to discover the joint latent space between H\u0026amp;E WSIs and bulk omic profiles, the correlation structure in ground truth data was also utilized to validate MOSBY model predictions (Methods). This analysis revealed that the MOSBY signature-signature correlation matrix was significantly different from random for 67 out of 68 slides (empirical p-value\u0026thinsp;\u0026lt;\u0026thinsp;0.05) (Supplementary Fig.\u0026nbsp;1c). For gene-gene correlation matrices, this value dropped to only 28 out of 68 slides (Supplementary Fig.\u0026nbsp;1f), again highlighting that computing gene signature scores is an effective strategy to deal with the dropout issue in spatially resolved transcriptomic datasets. Despite reduced overall performance in gene models, predicted gene expression levels still successfully captured specific highly correlated gene blocks derived from ground truth data, such as those among cancer-associated fibroblast (CAF) and immune cell marker genes (Supplementary Fig.\u0026nbsp;1g). An analysis across 68 slides indicated that MOSBY model predictions were particularly good at distinguishing CAF-related features from other TME population marker genes. For 63 out of 68 slides, the MOSBY gene-gene correlation matrix had a cluster significantly enriched for fibroblast marker genes, and this enrichment was higher than 89% for the large majority of the slides (Supplementary Fig.\u0026nbsp;1h).\u003c/p\u003e \u003cp\u003e \u003cb\u003eSerially sectioned CD8 IHC whole slide images validate MOSBY\u003c/b\u003e \u003cb\u003ein silico\u003c/b\u003e \u003cb\u003espatialization\u003c/b\u003e\u003c/p\u003e \u003cp\u003eMOSBY tile-level predictions were further validated with CD8 antibody-stained immunohistochemistry (IHC) WSIs. To this end, gene and signature models were trained with data from urothelial carcinoma immune checkpoint inhibitor (ICI) trials IMvigor210\u003csup\u003e40\u003c/sup\u003e and IMvigor211\u003csup\u003e41\u003c/sup\u003e (Methods, N\u0026thinsp;=\u0026thinsp;1460 paired H\u0026amp;E and RNA-seq samples). H\u0026amp;E-CD8 double-stained slides were not available, thus serially-sectioned single-stained H\u0026amp;E and CD8 IHC whole slide images were computationally aligned (Methods, N\u0026thinsp;=\u0026thinsp;1017 paired H\u0026amp;E and CD8 IHC WSIs). 454 slides that were correctly aligned based on multiple alignment metrics were subjected to a manual QC check and a stringent 95% tissue overlap threshold (Methods, Fig.\u0026nbsp;2a). The remaining 42 IHC WSIs constituted the analysis set, and were used to categorize \u0026ldquo;brownness\u0026rdquo; of pixels based on diaminobenzidine staining (\u003cem\u003ei.e.\u003c/em\u003e DAB mask)(Methods). Each DAB mask was split into tiles, and \u0026ldquo;brown\u0026rdquo; pixels in each tile were counted to determine tile-level CD8 IHC quantification (Methods). Pearson correlation coefficient was then computed between CD8 IHC and MOSBY-predicted tile-level values (Fig.\u0026nbsp;2a).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo compare with CD8 IHC quantification, gene features \u003cem\u003eCD8A, CD8B\u003c/em\u003e, and signature features Cibersort_CD8_T_cells and T_effector_cells were used from MOSBY model predictions. Across 42 slides, \u003cem\u003eCD8B\u003c/em\u003e and Cibersort_CD8_T_cells showed the strongest concordance with CD8 IHC, having the highest correlation in 22 and 15 slides respectively (Fig.\u0026nbsp;2b). Overall, \u003cem\u003eCD8B\u003c/em\u003e model predictions had the highest median correlation with CD8 IHC tile-level values (Pearson r\u0026thinsp;=\u0026thinsp;0.352), with a range of [-0.006, 0.789]. The density plot for r\u0026thinsp;=\u0026thinsp;0.789 correlation (N\u0026thinsp;=\u0026thinsp;16859 tiles, p\u0026thinsp;\u0026asymp;\u0026thinsp;0, two-tailed exact test) indicates stronger concordance for values closer to the higher and lower ends of the distribution (Fig.\u0026nbsp;2c). Visual inspection of CD8 IHC and H\u0026amp;E WSIs confirmed successful alignment, but also revealed minor differences between tissue sections potentially due to these being non-adjacent sections (Fig.\u0026nbsp;2d). At low magnification, visualization of DAB mask tile level quantification and \u003cem\u003eCD8B\u003c/em\u003e model predictions indicated that the model successfully captures CD8 T cell infiltration patterns in the tumor (Fig.\u0026nbsp;2e, 2f). High magnification inspection of model predictions and DAB mask (as well as the input H\u0026amp;E and CD8 IHC images) confirmed that MOSBY, despite being weakly supervised on a patient level bulk RNA-seq dataset, effectively learned a CD8 T cell-specific representation in H\u0026amp;E images (Fig.\u0026nbsp;2g\u0026ndash;2j).\u003c/p\u003e \u003cp\u003eIn order to gain insight into both high and low concordance cases, four representative slides were selected based on the distribution of predicted \u003cem\u003eCD8B\u003c/em\u003e vs DAB mask correlations (near the minimum, and 25th, 50th, and 75th percentiles). Selected slides had Pearson correlation values 0.48 (N\u0026thinsp;=\u0026thinsp;30032 tiles, p\u0026thinsp;\u0026asymp;\u0026thinsp;0, two-tailed exact test), 0.35 (N\u0026thinsp;=\u0026thinsp;22942 tiles, p\u0026thinsp;\u0026asymp;\u0026thinsp;0, two-tailed exact test), 0.25 N\u0026thinsp;=\u0026thinsp;11642 tiles, p\u0026thinsp;\u0026asymp;\u0026thinsp;0, two-tailed exact test, and 0.011 (N\u0026thinsp;=\u0026thinsp;16595 tiles, p\u0026thinsp;=\u0026thinsp;0.16, two-tailed exact test) respectively (Supplementary Fig.\u0026nbsp;2a\u0026ndash;2d). Low magnification inspection of DAB mask and MOSBY \u003cem\u003eCD8B\u003c/em\u003e predictions indicated that model predictions captured overall CD8 T cell infiltration patterns in the first three slides despite the correlation being as low as 0.25 (Supplementary Fig.\u0026nbsp;2a\u0026ndash;2c). For the last slide (r\u0026thinsp;=\u0026thinsp;0.011), spatial patterns in DAB mask and model predictions were conspicuously discordant (Supplementary Fig.\u0026nbsp;2d). High magnification inspection of CD8 IHC and DAB mask images from this last slide revealed that DAB mask captured brown staining artifacts in regions that did not show any expected CD8-specific staining patterns (Supplementary Fig.\u0026nbsp;2e, 2f). In addition, CD8 IHC and H\u0026amp;E tissue sections showed misalignment in small substructures, despite an overall high score in tissue alignment (Supplementary Fig.\u0026nbsp;2e-2h). Taken together, these results showed that MOSBY successfully learned a CD8 T cell-specific representation on H\u0026amp;E images, and correlations as low as 0.25 may be sufficient to capture overall CD8 T cell infiltration patterns in the tumor microenvironment, suggesting the potential of this tool for patient stratification in ICI trials.\u003c/p\u003e \u003cp\u003e \u003cb\u003eValidation with pathologist annotations of cancer regions indicates epithelial gene and signature markers may be misleading in accounting for tumor-specific expression\u003c/b\u003e \u003c/p\u003e \u003cp\u003eVarious RNA- and protein-based epithelial markers such as \u003cem\u003eEPCAM\u003c/em\u003e (GENE) and E-cadherin (PROTEIN) are frequently used in the literature to mark cancer cells in tissue specimens. However, epithelial markers can be expressed in both transformed and non-transformed epithelial cells. We utilized MOSBY tile level predictions in order to compare virtual spatialization of epithelial markers with pathologist tumor region annotations and assess the tumor-specificity of commonly used tumor markers.\u003c/p\u003e \u003cp\u003eMOSBY full models were trained in TCGA lung adenocarcinoma (LUAD), and deployed on WSIs in IMpower150\u003csup\u003e43\u003c/sup\u003e, an independent validation cohort in nonsquamous nonsmall cell lung cancer with pathologist annotations (i.e. pen marks) demarcating cancer epithelia.\u003c/p\u003e \u003cp\u003ePathologist annotations were compared with MOSBY features of tumor purity (DNA-based) and epithelial markers \u003cem\u003eEPCAM\u003c/em\u003e (GENE), E-cadherin (PROTEIN), and an epithelial gene signature by Taube \u003cem\u003eet al\u003c/em\u003e\u003csup\u003e20\u003c/sup\u003e. WSI level comparison indicated that MOSBY models trained with DNA-based quantification of cancer cells (tumor purity) achieves a stronger concordance with pathologist tumor region annotations (Fig.\u0026nbsp;3a, 3b). E-cadherin closely followed DNA-based features in terms of concordance with pathologist annotations, albeit noticeably lower specificity for the cancer region in Slide 2 (Fig.\u0026nbsp;3d vs 3b). Interestingly, predicted spatialization of \u003cem\u003eEPCAM\u003c/em\u003e transcript levels showed relatively poor concordance with cancer epithelium annotation (Fig.\u0026nbsp;3c); and the performance of the Taube \u003cem\u003eet al\u003c/em\u003e. epithelial signature was further reduced with highest predicted tiles corresponding to non-tumor regions (Fig.\u0026nbsp;3e). Collectively, these results suggest that employing DNA-based tumor cellularity estimates is the most reliable approach for accounting for tumor-specific expression in bulk datasets; and published epithelial signatures may not be fit-for-purpose in this task.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cb\u003eMOSBY enables investigation of intratumor heterogeneity with prediction of single gene, signature, and protein expression levels\u003c/b\u003e \u003c/p\u003e \u003cp\u003eThe MOSBY framework employs RNA-seq and RPPA data to predict gene (or signature) and protein expression levels respectively (Methods). The utility of the model in predicting single gene, signature, and protein expression levels and inferring intratumor heterogeneity is demonstrated in Impower150 in the context of cell proliferation features (Supplementary Fig.\u0026nbsp;3). Three representative H\u0026amp;E slides are shown in this figure together with pathologist annotations of cancer epithelia (Supplementary Fig.\u0026nbsp;3a), and MOSBY-predicted spatialization of cell proliferation markers Cyclin B1 (Supplementary Fig.\u0026nbsp;3b), Hallmark G2M checkpoint signature (Supplementary Fig.\u0026nbsp;3c) and \u003cem\u003eMKI67\u003c/em\u003e (Supplementary Fig.\u0026nbsp;3d). Even though tile-level predictions of all three cell cycle markers are highest in cancer epithelia-annotated regions, the spatialization of protein and signature markers (Cyclin B1 and Hallmark G2M checkpoint signature) demonstrate a stronger localization to cancer epithelia (Supplementary Fig.\u0026nbsp;3b, 3c, 3d). In addition, the spatialization of \u003cem\u003eMKI67\u003c/em\u003e expression in all three slides appears to suffer from checkerboard-like artifacts (Supplementary Fig.\u0026nbsp;3d). These results suggest that molecular data prediction from WSIs may be more successful with proteins or gene signatures, as opposed to single genes. Potential reasons for this observation are 1) proteins being closer to phenotype than genes and 2) gene signatures increasing the signal-to-noise ratio from RNA-seq data. Biologically, the predicted spatial pattern of Cyclin B1 and Hallmark G2M checkpoint signature levels show differential proliferation levels across cancer epithelia tiles (Slide 3 in Supplementary Fig.\u0026nbsp;3b, 3c), indicating MOSBY also has potential to infer intratumor heterogeneity.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eMOSBY predicts stroma, immune, and proliferation features with highest accuracy in TCGA\u003c/h2\u003e \u003cp\u003eCross-validated test set performance in single-indication signature and protein models was used to determine biological features with the highest prediction accuracy, and those with 25 highest performance are highlighted in Fig.\u0026nbsp;4a, 4d. The best-performing signatures ESTIMATE Stroma and ESTIMATE Immune had median correlations 0.627 and 0.616 across all tested indications (Fig.\u0026nbsp;4a). In randomly split test sets, the Spearman correlation of ESTIMATE Stroma and ESTIMATE Immune signatures reached as high as 0.787 and 0.781 in skin cutaneous melanoma and thyroid cancer respectively (Fig.\u0026nbsp;4b). Other best-performing signatures were again largely enriched in stroma and immune features as well as general mesenchymal characteristics (such as Hallmark EMT and Taube \u003cem\u003eet al\u003c/em\u003e mesenchymal signatures). (Fig.\u0026nbsp;4a). Of note, particular signatures known to play an important role in specific cancer indications also showed strong prediction accuracy in those pertinent settings. For instance, the Hallmark Angiogenesis signature had Spearman r\u0026thinsp;=\u0026thinsp;0.746 in the liver cancer model (LIHC), and a fatty acid elongation signature had Spearman r\u0026thinsp;=\u0026thinsp;0.741 in the low grade glioma setting (LGG) (Fig.\u0026nbsp;4b). After adjusting for multiple hypothesis testing with false discovery rate (adj. p\u0026thinsp;\u0026lt;\u0026thinsp;0.05), 14 out of 21 tested indications showed significant performance (as measured by Spearman r) for more than 90% of tested signatures (total N\u0026thinsp;=\u0026thinsp;175 signatures) (Fig.\u0026nbsp;4c).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eIn contrast with our signature set that represented both tumor-intrinsic pathways and cell populations in the TME, the TCGA RPPA antibody set was heavily focused on tumor-intrinsic pathways. Therefore, the best-performing proteins had a diversity of representation from tumor-intrinsic characteristics, such as proliferation (Cyclin-B1, FOXM1), DNA repair (MSH6, PARP1), and apoptosis (cleaved Caspase7) (Fig.\u0026nbsp;4d). MSH6 and Cyclin B1 were the 2 best-performing proteins that were tested in at least 20 indications (Fig.\u0026nbsp;4d). These proteins had median Spearman r\u0026thinsp;=\u0026thinsp;0.423 and 0.417 respectively across tested indications, but correlations in individual indications were as high 0.676 in PRAD for Cyclin B1, and 0.593 in LGG for MSH6. Overall, lower prediction accuracy in protein models was expected due to the lower signal-to-noise ratio in the RPPA technology compared to RNA-seq.\u0026nbsp;However, in specific settings such as the PRAD model, both total and phosphoproteins showed strong prediction accuracy in test sets. Progesterone receptor (PR) and cMET_pY1235 reached Spearman r\u0026thinsp;=\u0026thinsp;0.614 and 0.631 respectively in this indication (Fig.\u0026nbsp;4e). Also, the AKT/mTOR pathway showed evidence of strong prediction accuracy in the sarcoma setting where S6 and RICTOR_pT1135 antibodies showed Spearman r\u0026thinsp;=\u0026thinsp;0.731 and 0.695 respectively (Fig.\u0026nbsp;4e). Despite paucity of representation in the RPPA panel, stromal and immune features were also found among the best-performing proteins such as ECM-associated Collagen-VI, Fibronectin and T cell-associated Lck (Fig.\u0026nbsp;4d). Overall, 11 out of 20 tested indications showed significant performance (adj. p\u0026thinsp;\u0026lt;\u0026thinsp;0.1) for more than 90% of tested antibodies (total N\u0026thinsp;=\u0026thinsp;191 antibodies) (Fig.\u0026nbsp;4f).\u003c/p\u003e \u003cp\u003eIn terms of single gene models, best performing features confirmed the strong prediction accuracy associated with stromal features observed above; the highest-ranking genes were marker genes of fibroblasts (e.g. LUM, COL5A1) (Fig.\u0026nbsp;4g). Top-ranking genes also included markers of macrophages (e.g. CSF1R), and T cells (e.g. CD3E). In DNA models, tumor cellularity measures (tumor purity, cancer DNA fraction) achieved better performance than subclonal genome fraction, and average DNA methylation features (Fig.\u0026nbsp;4h).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eMulti-indication runs may show exaggerated performance\u003c/h2\u003e \u003cp\u003eThe results above demonstrated that MOSBY models trained with gene signature and protein expression data generally perform better than those trained with single genes. Thus, using signature and protein expression data, we next addressed whether single- vs multi-indication models (introduced in Fig.\u0026nbsp;1b) had more utility from a clinical perspective. We focused on performance in the unseen test set to evaluate the generalization capacity of each approach in a scenario emulating the clinical setting; \u003cem\u003ei.e.\u003c/em\u003e predicting gene and protein expression levels for a single patient. Performance was measured, as above, with the Spearman correlation between MOSBY slide-level predictions and ground truth levels of signature or (phospho)protein expression (N\u0026thinsp;=\u0026thinsp;175 for signatures, 191 for total and phosphoproteins).\u003c/p\u003e \u003cp\u003eMulti-indication models (PANCAN, ADENO, BRAIN, KIDNEY, SQUAMOUS, LUNG) showed superior performance for both signature and protein prediction tasks, compared to single-indication counterparts (Supplementary Fig.\u0026nbsp;4a). The PANCAN model had a median cross-validated Spearman correlation of 0.693 for tested signatures (Supplementary Fig.\u0026nbsp;4a), with epithelial, proliferation, and stemness signatures showing the best prediction accuracy (cross-validated Spearman r\u0026thinsp;=\u0026thinsp;0.826, 0.824, and 0.82 respectively). We focused on the Benporath_ES2 stemness signature to demonstrate the bias associated with multi-indication models. In a randomly split test set, this signature had Spearman r\u0026thinsp;=\u0026thinsp;0.84 between predicted and ground truth levels in the PANCAN dataset (Supplementary Fig.\u0026nbsp;4b, left). However, the correlations in individual indications varied between 0.244\u0026ndash;0.645 for this signature, suggesting that the PANCAN model performance may be spuriously inflated by combining multiple indications with different distributions; reminiscent of the Simpson\u0026rsquo;s paradox (e.g. OV and READ slides concentrated on the higher end, and THCA and KIRC concentrated on the lower end of the distribution; Spearman r\u0026thinsp;=\u0026thinsp;0.248, 0.326, 0.41, and 0.482 for OV, READ, THCA and KIRC respectively) (Supplementary Fig.\u0026nbsp;4b, right).\u003c/p\u003e \u003cp\u003eIn cross-validated protein models, Cyclin B1, Claudin-7 and progesterone receptor (PR) showed the best prediction accuracy in the PANCAN model (Spearman r\u0026thinsp;=\u0026thinsp;0.775, 0.772, and 0.749 respectively). To ask whether this performance was spuriously inflated, we again compared PANCAN and single indication models. In a randomly split test set, Cyclin B1 had Spearman correlation of 0.76 in the PANCAN model (Supplementary Fig.\u0026nbsp;4c, left), while correlations in individual indications were lower ranging from \u0026minus;\u0026thinsp;0.10 to 0.597. Simpson\u0026rsquo;s paradox was again evident in protein model predictions with Cyclin B1 levels in individual indications localizing to smaller parts of the pan-cancer distribution and showing lower prediction accuracy (e.g. Spearman r = \u0026minus;\u0026thinsp;0.103, 0.198, 0.214, and 0.263 for LIHC, THCA, READ and CESC respectively) (Supplementary Fig.\u0026nbsp;4c, right).\u003c/p\u003e \u003cp\u003eWe then investigated the potential bias in multi-indication models more systematically by comparing median test set performance (i.e. median correlation for signatures or proteins) in single \u003cem\u003evs\u003c/em\u003e multi-indication models (Supplementary Fig.\u0026nbsp;4e, 4g). The train/validate/test set splits in multi-indication models were implemented with indication stratification, enabling comparison with single-indication training. Results showed that neither single-indication nor multi-indication training was universally better (Supplementary Fig.\u0026nbsp;4e, 4g). However, evidence also showed that multi-indication training generally did not improve performance if a single-indication model reached a median performance of 0.45 for signatures (Supplementary Fig.\u0026nbsp;4e), and 0.35 for proteins (Supplementary Fig.\u0026nbsp;4g), with LGG as a potential exception benefiting from joint training with GBM.\u003c/p\u003e \u003cp\u003ePan-tissue training appeared to improve median test set accuracy in certain settings such as protein models in kidney, brain, lung, colon/rectum (Supplementary Fig.\u0026nbsp;4g), and signature models in brain (Supplementary Fig.\u0026nbsp;4e). However, at the individual signature/protein level, further heterogeneity in model performance was observed (Supplementary Fig.\u0026nbsp;4f, 4h). A subset of signatures/proteins benefited from multi-indication training whereas another subset did not. This observation prevented us from making general recommendations as to the deployment of single \u003cem\u003evs\u003c/em\u003e multi-indication models in a clinical setting. In circumstances where there is an interest in specific signatures/proteins, it may be suboptimal to make the model decision based on median performance. However, results in Supplementary Fig.\u0026nbsp;4e, 4g can be used as a guideline in other settings that lack a focus on specific biological features. We utilize single indication models in subsequent parts of this study in order to prevent potential bias from multi-indication training.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eSpatial patterns inferred from tile-level MOSBY predictions increase survival predictive power of gene signature-based models\u003c/h2\u003e \u003cp\u003eMOSBY tile-level predictions enable \u003cem\u003ein silico\u003c/em\u003e spatialization for a tested omic feature, as well as assess spatial correlation between two tested features (positive and negative correlations indicating colocalization and spatial exclusion respectively). We define a \u0026lsquo;colocalization feature\u0026rsquo; as the Pearson correlation between tile-level predictions of two omic features. As correlation coefficients were computed across all tiles on a slide, colocalization features represent slide-level as opposed to \u0026lsquo;local\u0026rsquo; spatial patterns. We focused on our signature panel (N\u0026thinsp;=\u0026thinsp;175) as the omic features to derive colocalization features and investigate survival associations, as the signature panel covered both tumor and non-tumor TME components. Moreover, using signatures as opposed to single genes enable the discovery of \u0026lsquo;biologically interpretable\u0026rsquo; spatial biomarkers as well-designed signatures capture pathway and cell type-related gene expression with higher fidelity. For a slide, correlations from all pairwise combinations of signatures (N\u0026thinsp;=\u0026thinsp;15225) (i.e. the collection of all colocalization features) are referred to as the \u0026lsquo;colocalization map\u0026rsquo; from here on (Fig.\u0026nbsp;1a).\u003c/p\u003e \u003cp\u003eA survival analysis was performed to ask whether slide-level patterns in colocalization maps harbored survival signals that could not be captured by the mere magnitude of signature expression. Three L1-regularized Cox proportional hazards regression models were fit to address this question, and the ability of the models to predict survival was computed with concordance indices (c-index) (Fig.\u0026nbsp;5a) (\u003cb\u003eMethods\u003c/b\u003e). The first model only used flattened colocalization features (MOSBY predictions, N\u0026thinsp;=\u0026thinsp;15225). The second Cox model only used signature expression levels (N\u0026thinsp;=\u0026thinsp;175) with the goal of assessing the survival predictive power of gene expression \u0026lsquo;magnitude\u0026rsquo;. The third model was a joint model combining all signatures but only lasso-selected colocalization features from the first model in order to prevent colocalization features from dominating the model. Each model was run across 10 cross-validation folds to optimize the shrinkage parameter and also obtain mean and standard error estimates for the relevant c-index (\u003cb\u003eMethods\u003c/b\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThis process was implemented separately for all tested TCGA indications, and c-index mean and standard error estimates were plotted (Fig.\u0026nbsp;5a). We observed that the joint (i.e. third) model had a higher c-index than the signature-only model in most indications. This finding revealed that slide-level spatial patterns, discovered with an inference engine such as MOSBY, have survival predictive power that could not be captured by gene expression alone. As MOSBY colocalization features are biologically interpretable, this opens the door to discovering potentially \u0026lsquo;actionable\u0026rsquo; colocalization or spatial exclusion patterns that predict clinical outcomes. Moreover, the c-indices achieved in the joint model were highly competitive with or higher than those reported in the literature from end-to-end survival-trained multimodal models utilizing genomic, transcriptomic, and image datasets\u003csup\u003e44\u003c/sup\u003e. Of note, comparing colocalization-only and signature-only models, we also found that the total survival predictive power of slide-level spatial patterns was not as high as that of signature levels in most TCGA indications. Ovarian and rectal cancer results were an exception to this general pattern (Fig.\u0026nbsp;5a), suggesting spatial biomarkers discovered in these indications may have the greatest potential to lead to novel insights.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eColocalization maps enable discovery of biologically interpretable spatial biomarkers\u003c/h2\u003e \u003cp\u003eWe next investigated consistent spatial predictors of risk that were supported across multiple indications and also showed evidence of tumor specificity. In each indication, survival effects and tumor specificity of colocalization features were explored with two tests: 1) A univariate Cox regression model for survival, and 2) a Mann-Whitney test to compare tumor-normal levels of colocalization value. A potential spatial biomarker of risk was defined as a colocalization feature that was significantly associated with poor overall survival (p\u0026thinsp;\u0026lt;\u0026thinsp;0.05) and also had elevated levels of colocalization in the tumor (adj. p\u0026thinsp;\u0026lt;\u0026thinsp;0.05). Given these criteria, four colocalization features had evidence in four different TCGA indications to be a spatial biomarker of risk (Fig.\u0026nbsp;5b). Of these four, the colocalization between an ER stress signature (XBP1s targets ER17\u003csup\u003e22\u003c/sup\u003e) and a neurotransmitter signature was associated with poor survival and malignant state in colon adenocarcinoma, lung adenocarcinoma, liver hepatocellular and ovarian cancers (Fig.\u0026nbsp;5c, 5d). In an independent non-squamous lung cancer study involving both immune checkpoint blockade (atezolizumab) and chemotherapy arms (IMpower110), this colocalization feature was also found to be associated with poor survival in the chemotherapy, but not immunotherapy arm, suggesting higher relevance as a resistance factor in chemotherapy (Fig.\u0026nbsp;5e, 5f). Of note, the ER17 and neurotransmitter signatures were not individually found to be associated with risk in any of the four mentioned indications (Supplementary Fig.\u0026nbsp;5a, 5b). Visual inspection of WSIs indicated that the expression of ER17 and neurotransmitter signatures primarily came from the microenvironment (as opposed to tumor region) in the case of high colocalization (Supplementary Fig.\u0026nbsp;5c). In low colocalization cases, the neurotransmitter signature expression primarily came from the tumor region whereas the ER17 signature expression was again predominantly in the microenvironment (Supplementary Fig.\u0026nbsp;5d).\u003c/p\u003e \u003cp\u003eFocusing on colocalization features involving immune system signatures, the T effector cell vs Cysteine colocalization was identified as the most consistent spatial biomarker of risk in TCGA. This colocalization feature was associated with poor survival and also showed significant tumor enrichment in breast, squamous lung, and ovarian cancers (Fig.\u0026nbsp;6a, 6b). T effector and cysteine signatures were not individually found to be associated with risk in any of these indications (Supplementary Fig.\u0026nbsp;5a, 6b). Visual inspection of WSIs indicated that the expression of T effector and cysteine signatures primarily came from the microenvironment (as opposed to tumor region) in the case of high colocalization (Supplementary Fig.\u0026nbsp;6c). As high colocalization is a risk factor, this expression pattern may be suggestive of a cysteine-associated immunosuppressive TME. In low colocalization cases, the T effector signature expression primarily came from the microenvironment whereas the cysteine signature expression was predominantly in the tumor region (Supplementary Fig.\u0026nbsp;6d).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eThe strongest survival effect for T effector vs cysteine colocalization was observed in breast cancer, where we investigated other immune cell types and found that immune vs cysteine colocalization was a general negative prognosis biomarker in this indication. Significant survival associations were observed for both lymphocytes/NK cells (Fig.\u0026nbsp;6c), and myeloid populations (Fig.\u0026nbsp;6d). The same immune vs cysteine colocalization features were also found to be negative prognostics in the Atezolizumab arm of Impower110 non-squamous cohort (Fig.\u0026nbsp;6e, 6f). Of note, most of these features were not prognostic in the chemotherapy arm (Supplementary Fig.\u0026nbsp;6e, 6f), however did not qualify as predictive biomarkers since the survival association differences between Atezolizumab and chemotherapy arms were not significant.\u003c/p\u003e \u003c/div\u003e"},{"header":"DISCUSSION","content":"\u003cp\u003eThe MOSBY workflow achieves prediction of bulk omic profiles from H\u0026amp;E WSI features. Our results showed that, compared to ImageNet-based pretraining, self-supervised pretraining in large histological datasets allows creation of inference engines (e.g. RetCCL) that enable a more accurate mapping from image features to gene, signature, protein, and DNA-based measurements. We demonstrated that the most accurately predicted features by MOSBY involved processes such as proliferation, immune/stromal infiltration, differentiation, and epithelial-to-mesenchymal transition. This finding suggests that self-supervised features learned from pan-cancer histological datasets run the risk of accentuating biological processes and pathways that show the highest variation across different cancer indications. Training feature extractors on images from a single indication may be required to capture biological processes that play an important role in one or only a few cancer indications. RetCCL-based features showed promise by capturing angiogenesis in hepatocellular carcinoma, and fatty acid biology in low grade glioma. Yet, the accumulation of even larger histological datasets in the future have potential to allow refinement of image features relevant for indication-specific pathways, thus making possible the discovery of a greater number of clinically relevant spatial biomarkers.\u003c/p\u003e \u003cp\u003eThe utility of H\u0026amp;E images to predict high-throughput expression data has been largely explored in the context of individual genes\u003csup\u003e1,2,12\u003c/sup\u003e. The RNA-Seq noise level associated with individual gene transcript counts can be alleviated by utilizing focused gene sets (i.e. signatures) aiming to measure particular pathway or cell type levels. Also, proteins and posttranslational modification levels (e.g. phosphorylation, acetylation) capture a more faithful representation of tissue phenotype compared to gene transcript levels. We demonstrated with pathologist tumor annotation validations that training the model with signature or protein expression levels led to potentially more interpretable results compared to single gene-trained models. Moreover, the MOSBY workflow has the benefit of employing both RNA and protein level data to interrogate the spatial distribution of biological processes, each modality functioning as an independent validation platform for the other. Therefore, training MOSBY with signature and protein levels have the highest potential of uncovering spatial patterns relevant for cancer biology and clinical outcome. In the case of tumor cellularity, our results indicated that DNA-based measures outperformed RNA and protein-based epithelial markers, and importantly, highlighted that published epithelial signatures may be misleading due to being expressed also by normal cells.\u003c/p\u003e \u003cp\u003eMOSBY, as with many other deep learning models, adopts a weakly-supervised approach by first making tile-level predictions and then aggregating tiles to obtain a prediction at the WSI level\u003csup\u003e4,45\u003c/sup\u003e. Although an intermediate output of the model, tile-level predictions enable \u003cem\u003ein silico\u003c/em\u003e spatialization of whole slide-level annotations, opening the way to inferring intratumor heterogeneity for the slide-level information used as ground truth\u003csup\u003e12,44\u003c/sup\u003e. Spatial intratumor heterogeneity patterns learned from a cohort of patients subsequently allow investigation of clinically relevant biomarkers. In MOSBY, spatial patterns are captured by pairwise colocalization features. For a given signature pair, the colocalization value on a slide is defined as the Pearson correlation across all tiles. Thus, MOSBY colocalization features capture \u003cem\u003eslide-level\u003c/em\u003e but not \u003cem\u003elocal\u003c/em\u003e spatial patterns. Local processes such as tertiary lymphoid structures are known to affect patient survival and response to cancer immunotherapy\u003csup\u003e46,47\u003c/sup\u003e. We demonstrated in this study that slide-level spatial patterns also carried survival signals, and increased predictive power of gene signature-based survival models in most TCGA indications. Moreover, we noted that the concordance indices of joint models (both colocalization features and gene signature levels) either surpassed or were comparable to those of multimodal deep learning models that employed the complete set of WSI and omic (RNA-seq, mutation status, copy number variation) data in TCGA\u003csup\u003e44\u003c/sup\u003e. This finding indicated that the 175-signature panel we defined was sufficient to capture most biological processes important for clinical outcome.\u003c/p\u003e \u003cp\u003eEnd-to-end neural networks trained to predict survival have been lacking in terms of direct biological interpretation of image regions important for the model\u003csup\u003e44,48\u003c/sup\u003e. These models may incorporate mechanisms such as attention heatmaps and spatial credit assignment to increase interpretability\u003csup\u003e44,49\u003c/sup\u003e, yet still require pathologist efforts to examine important image regions whereby interpretation remains limited to phenotypes visible by human eye. In contrast, MOSBY spatial features are biologically interpretable by design, which is an advantage of this approach over end-to-end neural networks trained to predict survival\u003csup\u003e44,48\u003c/sup\u003e. In this study, we showed that the colocalization of an endoplasmic reticulum (ER) stress-related signature and a neurotransmitter signature is both elevated in tumors and associated with poor overall survival in four TCGA indications. The poor survival association in lung adenocarcinoma was also validated in the chemotherapy arm of an independent NSCLC cohort (nonsquamous samples in Impower110), suggesting this colocalization may be a chemotherapy-specific risk factor. Moreover, we identified the T effector cell vs cysteine signature colocalization as a TME-related risk factor in multiple TCGA indications, as well as in the immunotherapy arm of Impower110 nonsquamous cohort. These results showcase the high utility of colocalization maps for discovering biologically interpretable clinically relevant spatial biomarkers.\u003c/p\u003e \u003cp\u003eLarge oncology data sets have enabled the routine development of machine learning models trained simultaneously in multiple cancer indications from the same organ (\u003cem\u003ei.e.\u003c/em\u003e pan-tissue training), or from different organs across the body (\u003cem\u003ei.e.\u003c/em\u003e pan-cancer training). Such models have shown success in mapping molecular tumor features to patient outcomes as well as in mapping WSI features to transcriptomic data \u003csup\u003e12,50,51\u003c/sup\u003e. However, in the clinical setting where the machine learning model will be used for a single patient, it becomes critical to understand whether multi-indication training enables more robust models with higher generalization capacity or introduces bias into the model due to distribution differences across indications. Hence, we trained single-indication, pan-tissue, and pan-cancer versions of the MOSBY model to investigate the relative merits and potential clinical utility of multi- \u003cem\u003evs.\u003c/em\u003e single-indication models. We discovered that better generalization (test set performance) is achieved with single-indication training in indications where the single-indication performance is moderate to high. In contrast, multi-indication training has potential to increase generalization capacity in indications where single-indication performance is relatively poor. While these conclusions can be made for the median performance across tested signatures and protein/phosphoproteins, significant heterogeneity was observed among these features in terms of showing better performance in the single- vs. multi-indication setting. To prevent potential bias from multi-indication training, we recommend the use of single-indication models where the feature set in the model reflects an unbiased discovery approach. However, an explicit single- vs. multi-indication comparison may be required for settings where the studied feature set reflects a specific biology of interest.\u003c/p\u003e \u003cp\u003eA limitation of our study is that MOSBY colocalization maps are not able to capture local spatial patterns. Future work involves the investigation of a graph neural network-based model on tile-level MOSBY predictions where we can capture local as well as slide-level patterns. Moreover, transformer-based architectures may increase the expressiveness of our model to allow a more accurate mapping from image features to bulk omic profiles.\u003c/p\u003e "},{"header":"METHODS","content":"\u003ch2\u003eData\u003c/h2\u003e\u003cp\u003e \u003cb\u003eTCGA\u003c/b\u003e: Batch-normalized RNA sequencing, RPPA datasets as well as clinical and DNA-based data were obtained from the PanCanAtlas publications page of the Genomic Data Commons website (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://gdc.cancer.gov/about-data/publications/pancanatlas\u003c/span\u003e\u003cspan address=\"https://gdc.cancer.gov/about-data/publications/pancanatlas\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). H\u0026amp;E-stained slide images (i.e. Tissue Slides) were downloaded from GDC Data Portal (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://portal.gdc.cancer.gov/\u003c/span\u003e\u003cspan address=\"https://portal.gdc.cancer.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003cspan type=\"Underline\" class=\"Underline\" name=\"Emphasis\"\u003e).\u003c/span\u003e\u003c/p\u003e\u003cp\u003e \u003cb\u003eSpatially-resolved transcriptomic data\u003c/b\u003e: Publicly available breast cancer tissue slides and spatial transcriptomic assays processed by the Spatial Transcriptomics method\u003csup\u003e39,52\u003c/sup\u003e were downloaded from \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://data.mendeley.com/datasets/29ntw7sh4r/5\u003c/span\u003e\u003cspan address=\"https://data.mendeley.com/datasets/29ntw7sh4r/5\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. The dataset contained 68 WSIs from a total of 23 patients along with spot coordinates and respective RNA-seq expression values from the spatial transcriptomics assay. WSIs were tiled into 224x224 patches as input for the MOSBY model. Each tile was centered around the pixel coordinate of an assay spot to represent the 100µm region of the spot. Log-normalized gene expression values were used as ground truth. Signature scores were calculated by first creating an AnnData object using the anndata package then applying the score_genes function from scanpy\u003csup\u003e53\u003c/sup\u003e. All analyses were conducted using Python 3.10.\u003c/p\u003e\u003ch2\u003eMOSBY preprocessing\u003c/h2\u003e\u003cp\u003e \u003cstrong\u003eImage tiling\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eOn whole slide images, foreground elements (tissue) and background (glass) were isolated through luminosity-based segmentation. The Python library OpenSlide was used as a backend for generating 224x224px foreground tiles at 0.5mpp resolution. The same tiling protocol was leveraged for both training and inference, allowing tiles to be mapped back to the original WSI positions for visualization (Supplementary Fig.\u0026nbsp;7).\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eFeature Extraction\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eContrastive self-supervised learning-based RetCCL\u003csup\u003e8\u003c/sup\u003e was used to extract image features. RetCCL employs a ResNet-50 architecture to extract 2048 features for each image tile.\u003c/p\u003e\u003ch2\u003eModel Training\u003c/h2\u003e\u003cp\u003e \u003cstrong\u003eTCGA\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eA maximum of 8000 tiles were selected randomly for each slide to yield an unbiased representation of the slide, and all tiles were concurrently used for training. A 2-layer perceptron (512 and 256 nodes per layer) was used to map image features to omic variables. Number of epochs was set to a maximum of 300, with early stopping allowed with a \u003cem\u003epatience\u003c/em\u003e of 30 epochs. Ground truth RNA-seq data were log-transformed. Ground truth signature and protein levels had negative values, and thus were shifted to make all values nonnegative. The model was trained with MSE loss between prediction and ground truth levels, while Spearman correlation between these quantities was used as early termination criterion in the validation set. A batch size of 64, and AdamW optimizer with 1e-3 weight decay were used. Learning rate scheduler was implemented with step size 5 and gamma 0.9.\u003c/p\u003e\u003cp\u003e5-fold cross-validation was performed (64% training, 16% validation, and 20% test set in each fold) to assess model performance. Information leak was prevented by assigning all WSIs from the same patient to the same partition. For inference, full models were trained with 80% of the data, with the remaining 20% used as validation set to check early stopping criterion.\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eIMvigor210 and IMvigor211\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eA full model was trained using both IMvigor210 and IMvigor211 images, and largely following the parameters used in TCGA training. Patients were split into 80% training and 20% validation set, stratified by trial and preventing information leak. The signature and gene model consisted of 175 signatures (Supplementary Table\u0026nbsp;3) and 73 genes (Supplementary Table\u0026nbsp;2) respectively. Different from TCGA runs, a maximum of 4000 randomly selected tiles per WSI were used during training, and WSIs having less than 200 tiles were filtered out. Ground truth RNA-seq data were log-transformed and standardized to have 0 mean and unit variance. Both gene and signature expression levels were shifted to make all values nonnegative.\u003c/p\u003e\u003ch2\u003eComputational hardware and software\u003c/h2\u003e\u003cp\u003eMOSBY was built with the PyTorch library (v1.11.0) in Python (v3.10.8) as a novel implementation of the HE2RNA\u003csup\u003e12\u003c/sup\u003e model. Python libraries used for data processing included NumPy (v1.23.5), Pandas (v1.2.4), Scikit-learn (v0.24.1), OpenSlide (v1.1.1), Zarr (v2.12.0), Tifffile (v2020.10.1), and OpenCV-cv2 (v4.5.5). Whole slide image tiling, RetCCL feature extraction and MOSBY model training were implemented in NVIDIA Tesla V100 Tensor Core GPUs (graphics processing units). Deep learning models were trained with NVIDIA Cuda compiler (v12.1.105). Data visualization in Python was implemented with Matplotlib (v3.3.4) and Seaborn (v0.11.1) libraries. Python statistical analyses such as Spearman and Pearson correlation were implemented with the SciPy library (v1.10.1).\u003c/p\u003e\u003cp\u003eData processing in R (v4.1.1) was implemented with dplyr (v1.0.8), magrittr (v2.0.2), and reshape2 (v1.4.4) libraries. Data visualization in R was performed using ggplot2 (v3.3.5), ggpubr (v0.4.0), and ggsci (v2.9) libraries.\u003c/p\u003e\u003ch2\u003eSpatial transcriptomics validation analysis\u003c/h2\u003e\u003cp\u003e \u003cstrong\u003eDerivation of an empirical null distribution for the distance between ground truth and predicted correlation structures\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eFor each slide, the gene-gene correlation matrix was computed for both ground truth and model predictions. The Frobenius distance between these two matrices was used as a test statistic. Subsequently, the MOSBY correlation matrix rows and columns were randomly shuffled, and the Frobenius distance to the ground truth matrix was computed. Random shuffling was repeated 1000 times to obtain an empirical null distribution for MOSBY vs ground truth distance. The test statistic was then compared with the empirical null distribution to obtain a p-value for model predictions. The same sequence of steps was followed for signature-signature correlations.\u003c/p\u003e\u003ch2\u003eCD8 IHC and H\u0026amp;E whole slide image alignment\u003c/h2\u003e\u003cp\u003eIMvigor210 and IMvigor211 slides were scanned by an external party (CellCarta, Montreal, QC) on a Pannoramic 250 (3DHistech, Budapest, Hungary) with a 20x or 40x objective. Digital files were transferred to Genentech and converted to the Aperio SVS image format for all viewing and analysis. Hematoxylin and Eosin (HE) slides were aligned with a CD8 stained section from the same block (1017 total sample pairs). These sections were almost never serial, and displayed variable degrees of adjacency characteristic of the block being refaced between sections. Image alignment was performed in Matlab (R2022a, Mathworks, Natick, MA.) by downsampling to ~ 20 um per pixel and converting them to normalized grayscale images, before calculating an affine transformation using mutual information as the underlying matching metric. The transformation was then upscaled before being applied to the original CD8 high magnification image data bringing it into alignment with the HE image. Various alignment metrics were produced (intersection over union, normalized cross correlation) to select for correctly aligned images (454 sample pairs) before a final manual QC check at low magnification (at least 95% of tissue present on both slides, with the majority of visible structures in close proximity) that was not exhaustive and resulted in 42 sample pairs subjected to further analysis.\u003c/p\u003e\u003cp\u003eA binary mask was then produced from the aligned CD8 image using both HSV thresholding and a blue-normalized “brownness” algorithm\u003csup\u003e54\u003c/sup\u003e.\u003c/p\u003e\u003ch2\u003eCD8 IHC tile-level quantification\u003c/h2\u003e\u003cp\u003eTile size used for H\u0026amp;E images was set to 224x224 pixels at 0.5mpp. In cases where the native resolution was different from 0.5mpp, tiling was performed with an adjusted tile size at the native resolution, and tiles were subsequently downsampled or upsampled to arrive at 224 pixels at 0.5mpp. In most IMvigor slides, the native resolution was 0.243mpp resulting in an adjusted tile size of approximately 460x460 pixels. CD8 IHC images were tiled using the ‘adjusted’ tile size to match with the corresponding H\u0026amp;E tiles. The CD8 IHC count for a tile was computed using a convolution approach: Each tile was split into 30x30p subtiles, for which the number of 1s (brown pixels) was counted. The counts across subtiles were averaged to obtain the value for the tile. This average value was log-transformed (log(x + 1)) in comparisons with MOSBY model predictions.\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eComputation of signature scores from bulk RNA-seq data\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eIn TCGA and IMvigor RNA-seq datasets, normalized gene expression values were log-transformed (log(x + 1)), z-transformed across samples to have 0 mean, unit variance, and subsequently averaged across genes to arrive at a single signature score for each sample. In TCGA, signature scoring was performed on the batch-normalized pan-cancer RNA-seq dataset, and hence signature scores were comparable across cancer indications.\u003c/p\u003e\u003ch2\u003eSurvival analysis\u003c/h2\u003e\u003cp\u003e \u003cb\u003eTCGA concordance index analysis\u003c/b\u003e: MOSBY signature models (full models trained with 80%-20% training/validation split) were used to make inference on all slides in a cancer indication. For a given slide, the colocalization for a signature pair was computed with a Pearson correlation across tile-level predictions of the two signatures. The collection of all pairwise correlation values for a slide formed the ‘slide colocalization matrix’ (N = 175x175). For patients with multiple existing H\u0026amp;E slides, the patient-level colocalization matrix was computed as the average of all pertinent slide-level colocalization matrices. Three sets of cross-validation runs were implemented for L1-regularized Cox regression models. The inputs to these survival models consisted of: Model 1) Flattened patient-level colocalization maps (N = 15225 features). Model 2, patient-level signature scores (N = 175 signatures). Model 3) all tested signatures (N = 175) and lasso-selected colocalization features from Model 1. The input features were regressed against overall survival in all models and for all tested indications.\u003c/p\u003e\u003cp\u003eIn all three models, the sequence of possible values for lambda (shrinkage parameter) was internally determined in the cv.glmnet function from the glmnet\u003csup\u003e55\u003c/sup\u003e R library (v4.1.3) prior to cross-validation runs. The lambda maximizing the mean Harrel’s concordance measure across 10 cross-validation folds was chosen as optimal, and used to determine the concordance index estimate for the model. The standard error estimate for the model concordance index was calculated across the 10 cross-validation folds.\u003c/p\u003e\u003cp\u003e \u003cstrong\u003eKaplan-Meier plots\u003c/strong\u003e \u003c/p\u003e\u003cp\u003eThe survminer R library (v0.4.9) was used with a median cutoff to generate Kaplan-Meier plots. Log-rank test p-values were obtained internally in survminer using the survdiff function in the survival\u003csup\u003e56\u003c/sup\u003e library (v3.3.1).\u003c/p\u003e\u003ch2\u003eData availability\u003c/h2\u003e\u003cp\u003eThe accession URLs for publicly available data analyzed in this study (TCGA, spatial transcriptomics) are listed in the \u003cspan refid=\"Sec11\" class=\"InternalRef\"\u003eData\u003c/span\u003e section of Methods. Datasets from clinical trials IMpower150 (H\u0026amp;E image data), IMpower110 (H\u0026amp;E image and clinical data), IMvigor210 (RNA-seq, H\u0026amp;E image and CD8 IHC image data), IMvigor211 (RNA-seq, H\u0026amp;E image and CD8 IHC image data) were also analyzed in the current study. IMvigor210 RNA-seq data is available at the European Genome-phenome archive (EGA) under the accession number EGAS00001002556, and was also published as an R package (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://research-pub.gene.com/IMvigor210CoreBiologies/\u003c/span\u003e\u003cspan address=\"http://research-pub.gene.com/IMvigor210CoreBiologies/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). IMpower110 and IMpower150 datasets are not publicly available as data release is designated for the pending primary biomarker manuscripts.\u003c/p\u003e\u003cp\u003eIMvigor210 and IMvigor211 H\u0026amp;E and CD8 IHC image data as well as IMvigor211 RNA-seq data that support the findings of this study are available from Roche, but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with permission of Roche. For up to date details on Roche's Global Policy on the Sharing of Clinical Information and how to request access to related clinical study documents, see here: \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://go.roche.com/data_sharing\u003c/span\u003e\u003cspan address=\"https://go.roche.com/data_sharing\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003eFunding: This work was supported by Genentech, Inc.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgements\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe would like to thank all of the study participants and their families, and all of the site investigators, study coordinators, and staff. We also would like to thank Brandon Kayser, Robert Johnston, A\u0026iuml;cha Bentaieb, Dan Ruderman, Hector Corrada Bravo, Jason Hackney, and colleagues from the Oncology Reverse Translation team for providing critical feedback on the manuscript. This work was supported by Genentech, Inc.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;\u003cstrong\u003eAuthor contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eK.L., A.K, Y.S. devised the study.\u003c/p\u003e\n\u003cp\u003eY.S., V.P, A.K, J.E, E.L, E.W., B.N., M.S., M.B. acquired, analyzed, and interpreted the data.\u003c/p\u003e\n\u003cp\u003eY.S. wrote the manuscript with input from the remaining authors.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;\u003cstrong\u003eCompeting interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eY.S. , J.E. , E.L , B.N. , M.S. , M.B. , K.L. are employees of Genentech, Inc. and shareholders in F. Hoffmann La Roche, Ltd.\u003c/p\u003e\n\u003cp\u003eV.P. is an external partner at Genentech, Inc.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eA.K. was an intern at Genentech, Inc.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003eE.W. was an intern at Genentech, In\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eComiter, C. \u003cem\u003eet al. Inference of single cell profiles from histology stains with the Single-Cell omics from Histology Analysis Framework (SCHAF)\u003c/em\u003e. http://biorxiv.org/lookup/doi/10.1101/2023.03.21.533680 (2023) doi:10.1101/2023.03.21.533680.\u003c/li\u003e\n\u003cli\u003eAlsaafin, A., Safarpoor, A., Sikaroudi, M., Hipp, J. D. \u0026amp; Tizhoosh, H. R. Learning to predict RNA sequence expressions from whole slide images with applications for search and classification. \u003cem\u003eCommun. Biol. \u003c/em\u003e\u003cstrong\u003e6\u003c/strong\u003e, 304 (2023).\u003c/li\u003e\n\u003cli\u003eYang, K. D. \u003cem\u003eet al.\u003c/em\u003e Multi-domain translation between single-cell imaging and sequencing data using autoencoders. \u003cem\u003eNat. Commun. \u003c/em\u003e\u003cstrong\u003e12\u003c/strong\u003e, 31 (2021).\u003c/li\u003e\n\u003cli\u003eTsai, P.-C. \u003cem\u003eet al.\u003c/em\u003e Histopathology images predict multi-omics aberrations and prognoses in colorectal cancer patients. \u003cem\u003eNat. Commun. \u003c/em\u003e\u003cstrong\u003e14\u003c/strong\u003e, 2102 (2023).\u003c/li\u003e\n\u003cli\u003eHaviv, D., Gatie, M., Hadjantonakis, A.-K., Nawy, T. \u0026amp; Pe\u0026rsquo;er, D. \u003cem\u003eThe covariance environment defines cellular niches for spatial inference\u003c/em\u003e. http://biorxiv.org/lookup/doi/10.1101/2023.04.18.537375 (2023) doi:10.1101/2023.04.18.537375.\u003c/li\u003e\n\u003cli\u003eChen, X., Fan, H., Girshick, R. \u0026amp; He, K. Improved Baselines with Momentum Contrastive Learning. Preprint at http://arxiv.org/abs/2003.04297 (2020).\u003c/li\u003e\n\u003cli\u003eOquab, M. \u003cem\u003eet al.\u003c/em\u003e DINOv2: Learning Robust Visual Features without Supervision. Preprint at http://arxiv.org/abs/2304.07193 (2023).\u003c/li\u003e\n\u003cli\u003eWang, X. \u003cem\u003eet al.\u003c/em\u003e RetCCL: Clustering-guided contrastive learning for whole-slide image retrieval. \u003cem\u003eMed. Image Anal. \u003c/em\u003e\u003cstrong\u003e83\u003c/strong\u003e, 102645 (2023).\u003c/li\u003e\n\u003cli\u003eFremond, S. \u003cem\u003eet al.\u003c/em\u003e Interpretable deep learning model to predict the molecular classification of endometrial cancer from haematoxylin and eosin-stained whole-slide images: a combined analysis of the PORTEC randomised trials and clinical cohorts. \u003cem\u003eLancet Digit. Health \u003c/em\u003e\u003cstrong\u003e5\u003c/strong\u003e, e71\u0026ndash;e82 (2023).\u003c/li\u003e\n\u003cli\u003eSchirris, Y., Gavves, E., Nederlof, I., Horlings, H. M. \u0026amp; Teuwen, J. DeepSMILE: Contrastive self-supervised pre-training benefits MSI and HRD classification directly from H\u0026amp;E whole-slide images in colorectal and breast cancer. \u003cem\u003eMed. Image Anal. \u003c/em\u003e\u003cstrong\u003e79\u003c/strong\u003e, 102464 (2022).\u003c/li\u003e\n\u003cli\u003eCiga, O., Xu, T. \u0026amp; Martel, A. L. Self supervised contrastive learning for digital histopathology. \u003cem\u003eMach. Learn. Appl. \u003c/em\u003e\u003cstrong\u003e7\u003c/strong\u003e, 100198 (2022).\u003c/li\u003e\n\u003cli\u003eSchmauch, B. \u003cem\u003eet al.\u003c/em\u003e A deep learning model to predict RNA-Seq expression of tumours from whole slide images. \u003cem\u003eNat. Commun. \u003c/em\u003e\u003cstrong\u003e11\u003c/strong\u003e, 3877 (2020).\u003c/li\u003e\n\u003cli\u003eNahhas, O. S. M. E. \u003cem\u003eet al.\u003c/em\u003e Regression-based Deep-Learning predicts molecular biomarkers from pathology slides. Preprint at http://arxiv.org/abs/2304.05153 (2023).\u003c/li\u003e\n\u003cli\u003eHoadley, K. A. \u003cem\u003eet al.\u003c/em\u003e Cell-of-Origin Patterns Dominate the Molecular Classification of 10,000 Tumors from 33 Types of Cancer. \u003cem\u003eCell \u003c/em\u003e\u003cstrong\u003e173\u003c/strong\u003e, 291-304.e6 (2018).\u003c/li\u003e\n\u003cli\u003eBagaev, A. \u003cem\u003eet al.\u003c/em\u003e Conserved pan-cancer microenvironment subtypes predict response to immunotherapy. \u003cem\u003eCancer Cell \u003c/em\u003e\u003cstrong\u003e39\u003c/strong\u003e, 845-865.e7 (2021).\u003c/li\u003e\n\u003cli\u003ePatil, N. S. \u003cem\u003eet al.\u003c/em\u003e Intratumoral plasma cells predict outcomes to PD-L1 blockade in non-small cell lung cancer. \u003cem\u003eCancer Cell \u003c/em\u003e\u003cstrong\u003e40\u003c/strong\u003e, 289-300.e4 (2022).\u003c/li\u003e\n\u003cli\u003eLiberzon, A. \u003cem\u003eet al.\u003c/em\u003e The Molecular Signatures Database (MSigDB) hallmark gene set collection. \u003cem\u003eCell Syst. \u003c/em\u003e\u003cstrong\u003e1\u003c/strong\u003e, 417\u0026ndash;425 (2015).\u003c/li\u003e\n\u003cli\u003eMiranda, A. \u003cem\u003eet al.\u003c/em\u003e Cancer stemness, intratumoral heterogeneity, and immune response across cancers. \u003cem\u003eProc. Natl. Acad. Sci. U. S. A. \u003c/em\u003e\u003cstrong\u003e116\u003c/strong\u003e, 9020\u0026ndash;9029 (2019).\u003c/li\u003e\n\u003cli\u003eŞenbabaoğlu, Y. \u003cem\u003eet al.\u003c/em\u003e Tumor immune microenvironment characterization in clear cell renal cell carcinoma identifies prognostic and immunotherapeutically relevant messenger RNA signatures. \u003cem\u003eGenome Biol. \u003c/em\u003e\u003cstrong\u003e17\u003c/strong\u003e, 231 (2016).\u003c/li\u003e\n\u003cli\u003eTaube, J. H. \u003cem\u003eet al.\u003c/em\u003e Core epithelial-to-mesenchymal transition interactome gene-expression signature is associated with claudin-low and metaplastic breast cancer subtypes. \u003cem\u003eProc. Natl. Acad. Sci. U. S. A. \u003c/em\u003e\u003cstrong\u003e107\u003c/strong\u003e, 15449\u0026ndash;15454 (2010).\u003c/li\u003e\n\u003cli\u003eMasiero, M. \u003cem\u003eet al.\u003c/em\u003e A core human primary tumor angiogenesis signature identifies the endothelial orphan receptor ELTD1 as a key regulator of angiogenesis. \u003cem\u003eCancer Cell \u003c/em\u003e\u003cstrong\u003e24\u003c/strong\u003e, 229\u0026ndash;241 (2013).\u003c/li\u003e\n\u003cli\u003eHarnoss, J. M. \u003cem\u003eet al.\u003c/em\u003e IRE1\u0026alpha; Disruption in Triple-Negative Breast Cancer Cooperates with Antiangiogenic Therapy by Reversing ER Stress Adaptation and Remodeling the Tumor Microenvironment. \u003cem\u003eCancer Res. \u003c/em\u003e\u003cstrong\u003e80\u003c/strong\u003e, 2368\u0026ndash;2379 (2020).\u003c/li\u003e\n\u003cli\u003eGene Ontology Consortium \u003cem\u003eet al.\u003c/em\u003e The Gene Ontology knowledgebase in 2023. \u003cem\u003eGenetics \u003c/em\u003e\u003cstrong\u003e224\u003c/strong\u003e, iyad031 (2023).\u003c/li\u003e\n\u003cli\u003eHobert, O., Carrera, I. \u0026amp; Stefanakis, N. The molecular and gene regulatory signature of a neuron. \u003cem\u003eTrends Neurosci. \u003c/em\u003e\u003cstrong\u003e33\u003c/strong\u003e, 435\u0026ndash;445 (2010).\u003c/li\u003e\n\u003cli\u003eRobertson, A. G. \u003cem\u003eet al.\u003c/em\u003e Comprehensive Molecular Characterization of Muscle-Invasive Bladder Cancer. \u003cem\u003eCell \u003c/em\u003e\u003cstrong\u003e171\u003c/strong\u003e, 540-556.e25 (2017).\u003c/li\u003e\n\u003cli\u003eTsai, H. K. \u003cem\u003eet al.\u003c/em\u003e Gene expression signatures of neuroendocrine prostate cancer and primary small cell prostatic carcinoma. \u003cem\u003eBMC Cancer \u003c/em\u003e\u003cstrong\u003e17\u003c/strong\u003e, 759 (2017).\u003c/li\u003e\n\u003cli\u003eGillespie, M. \u003cem\u003eet al.\u003c/em\u003e The reactome pathway knowledgebase 2022. \u003cem\u003eNucleic Acids Res. \u003c/em\u003e\u003cstrong\u003e50\u003c/strong\u003e, D687\u0026ndash;D692 (2022).\u003c/li\u003e\n\u003cli\u003eXu, Q. \u003cem\u003eet al.\u003c/em\u003e Pan-cancer transcriptome analysis reveals a gene expression signature for the identification of tumor tissue origin. \u003cem\u003eMod. Pathol. Off. J. U. S. Can. Acad. Pathol. Inc \u003c/em\u003e\u003cstrong\u003e29\u003c/strong\u003e, 546\u0026ndash;556 (2016).\u003c/li\u003e\n\u003cli\u003eBen-Porath, I. \u003cem\u003eet al.\u003c/em\u003e An embryonic stem cell-like gene expression signature in poorly differentiated aggressive human tumors. \u003cem\u003eNat. Genet. \u003c/em\u003e\u003cstrong\u003e40\u003c/strong\u003e, 499\u0026ndash;507 (2008).\u003c/li\u003e\n\u003cli\u003eBhattacharya, B., Puri, S. \u0026amp; Puri, R. K. A review of gene expression profiling of human embryonic stem cell lines and their differentiated progeny. \u003cem\u003eCurr. Stem Cell Res. Ther. \u003c/em\u003e\u003cstrong\u003e4\u003c/strong\u003e, 98\u0026ndash;106 (2009).\u003c/li\u003e\n\u003cli\u003eShats, I. \u003cem\u003eet al.\u003c/em\u003e Using a stem cell-based signature to guide therapeutic selection in cancer. \u003cem\u003eCancer Res. \u003c/em\u003e\u003cstrong\u003e71\u003c/strong\u003e, 1772\u0026ndash;1780 (2011).\u003c/li\u003e\n\u003cli\u003eKim, J. \u003cem\u003eet al.\u003c/em\u003e A Myc network accounts for similarities between embryonic stem and cancer cell transcription programs. \u003cem\u003eCell \u003c/em\u003e\u003cstrong\u003e143\u003c/strong\u003e, 313\u0026ndash;324 (2010).\u003c/li\u003e\n\u003cli\u003ePossemato, R. \u003cem\u003eet al.\u003c/em\u003e Functional genomics reveal that the serine synthesis pathway is essential in breast cancer. \u003cem\u003eNature \u003c/em\u003e\u003cstrong\u003e476\u003c/strong\u003e, 346\u0026ndash;350 (2011).\u003c/li\u003e\n\u003cli\u003eNewman, A. M. \u003cem\u003eet al.\u003c/em\u003e Robust enumeration of cell subsets from tissue expression profiles. \u003cem\u003eNat. Methods \u003c/em\u003e\u003cstrong\u003e12\u003c/strong\u003e, 453\u0026ndash;457 (2015).\u003c/li\u003e\n\u003cli\u003eMariathasan, S. \u003cem\u003eet al.\u003c/em\u003e TGF\u0026beta; attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. \u003cem\u003eNature \u003c/em\u003e\u003cstrong\u003e554\u003c/strong\u003e, 544\u0026ndash;548 (2018).\u003c/li\u003e\n\u003cli\u003eB\u0026ouml;ttcher, J. P. \u0026amp; Reis e Sousa, C. The Role of Type 1 Conventional Dendritic Cells in Cancer Immunity. \u003cem\u003eTrends Cancer \u003c/em\u003e\u003cstrong\u003e4\u003c/strong\u003e, 784\u0026ndash;792 (2018).\u003c/li\u003e\n\u003cli\u003eYoshihara, K. \u003cem\u003eet al.\u003c/em\u003e Inferring tumour purity and stromal and immune cell admixture from expression data. \u003cem\u003eNat. Commun. \u003c/em\u003e\u003cstrong\u003e4\u003c/strong\u003e, 2612 (2013).\u003c/li\u003e\n\u003cli\u003eAkbani, R. \u003cem\u003eet al.\u003c/em\u003e A pan-cancer proteomic perspective on The Cancer Genome Atlas. \u003cem\u003eNat. Commun. \u003c/em\u003e\u003cstrong\u003e5\u003c/strong\u003e, 3887 (2014).\u003c/li\u003e\n\u003cli\u003eStenbeck, L., Bergenstr\u0026aring;hle, L., Lundeberg, J. \u0026amp; Borg, \u0026Aring;. Human breast cancer in situ capturing transcriptomics. (2021) doi:10.17632/29ntw7sh4r.5.\u003c/li\u003e\n\u003cli\u003eBalar, A. V. \u003cem\u003eet al.\u003c/em\u003e Atezolizumab as first-line treatment in cisplatin-ineligible patients with locally advanced and metastatic urothelial carcinoma: a single-arm, multicentre, phase 2 trial. \u003cem\u003eLancet Lond. Engl. \u003c/em\u003e\u003cstrong\u003e389\u003c/strong\u003e, 67\u0026ndash;76 (2017).\u003c/li\u003e\n\u003cli\u003ePowles, T. \u003cem\u003eet al.\u003c/em\u003e Atezolizumab versus chemotherapy in patients with platinum-treated locally advanced or metastatic urothelial carcinoma (IMvigor211): a multicentre, open-label, phase 3 randomised controlled trial. \u003cem\u003eLancet Lond. Engl. \u003c/em\u003e\u003cstrong\u003e391\u003c/strong\u003e, 748\u0026ndash;757 (2018).\u003c/li\u003e\n\u003cli\u003eDiao, J. A. \u003cem\u003eet al.\u003c/em\u003e Human-interpretable image features derived from densely mapped cancer pathology slides predict diverse molecular phenotypes. \u003cem\u003eNat. Commun. \u003c/em\u003e\u003cstrong\u003e12\u003c/strong\u003e, 1613 (2021).\u003c/li\u003e\n\u003cli\u003eSocinski, M. A. \u003cem\u003eet al.\u003c/em\u003e Atezolizumab for First-Line Treatment of Metastatic Nonsquamous NSCLC. \u003cem\u003eN. Engl. J. Med. \u003c/em\u003e\u003cstrong\u003e378\u003c/strong\u003e, 2288\u0026ndash;2301 (2018).\u003c/li\u003e\n\u003cli\u003eChen, R. J. \u003cem\u003eet al.\u003c/em\u003e Pan-cancer integrative histology-genomic analysis via multimodal deep learning. \u003cem\u003eCancer Cell \u003c/em\u003e\u003cstrong\u003e40\u003c/strong\u003e, 865-878.e6 (2022).\u003c/li\u003e\n\u003cli\u003eLu, M. Y. \u003cem\u003eet al.\u003c/em\u003e Data-efficient and weakly supervised computational pathology on whole-slide images. \u003cem\u003eNat. Biomed. Eng. \u003c/em\u003e\u003cstrong\u003e5\u003c/strong\u003e, 555\u0026ndash;570 (2021).\u003c/li\u003e\n\u003cli\u003eHelmink, B. A. \u003cem\u003eet al.\u003c/em\u003e B cells and tertiary lymphoid structures promote immunotherapy response. \u003cem\u003eNature \u003c/em\u003e\u003cstrong\u003e577\u003c/strong\u003e, 549\u0026ndash;555 (2020).\u003c/li\u003e\n\u003cli\u003eTr\u0026uuml;b, M. \u0026amp; Zippelius, A. Tertiary Lymphoid Structures as a Predictive Biomarker of Response to Cancer Immunotherapies. \u003cem\u003eFront. Immunol. \u003c/em\u003e\u003cstrong\u003e12\u003c/strong\u003e, 674565 (2021).\u003c/li\u003e\n\u003cli\u003eKatzman, J. L. \u003cem\u003eet al.\u003c/em\u003e DeepSurv: personalized treatment recommender system using a Cox proportional hazards deep neural network. \u003cem\u003eBMC Med. Res. Methodol. \u003c/em\u003e\u003cstrong\u003e18\u003c/strong\u003e, 24 (2018).\u003c/li\u003e\n\u003cli\u003eJaved, S. A. \u003cem\u003eet al.\u003c/em\u003e Additive MIL: Intrinsically Interpretable Multiple Instance Learning for Pathology. (2022) doi:10.48550/ARXIV.2206.01794.\u003c/li\u003e\n\u003cli\u003eChowell, D. \u003cem\u003eet al.\u003c/em\u003e Improved prediction of immune checkpoint blockade efficacy across multiple cancer types. \u003cem\u003eNat. Biotechnol. \u003c/em\u003e\u003cstrong\u003e40\u003c/strong\u003e, 499\u0026ndash;506 (2022).\u003c/li\u003e\n\u003cli\u003eChen, R. J. \u003cem\u003eet al.\u003c/em\u003e Whole Slide Images are 2D Point Clouds: Context-Aware Survival Prediction using Patch-based Graph Convolutional Networks. (2021) doi:10.48550/ARXIV.2107.13048.\u003c/li\u003e\n\u003cli\u003eSt\u0026aring;hl, P. L. \u003cem\u003eet al.\u003c/em\u003e Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. \u003cem\u003eScience \u003c/em\u003e\u003cstrong\u003e353\u003c/strong\u003e, 78\u0026ndash;82 (2016).\u003c/li\u003e\n\u003cli\u003eWolf, F. A., Angerer, P. \u0026amp; Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. \u003cem\u003eGenome Biol. \u003c/em\u003e\u003cstrong\u003e19\u003c/strong\u003e, 15 (2018).\u003c/li\u003e\n\u003cli\u003eBrey, E. M. \u003cem\u003eet al.\u003c/em\u003e Automated selection of DAB-labeled tissue for immunohistochemical quantification. \u003cem\u003eJ. Histochem. Cytochem. Off. J. Histochem. Soc. \u003c/em\u003e\u003cstrong\u003e51\u003c/strong\u003e, 575\u0026ndash;584 (2003).\u003c/li\u003e\n\u003cli\u003eSimon, N., Friedman, J., Tibshirani, R. \u0026amp; Hastie, T. Regularization Paths for Cox\u0026rsquo;s Proportional Hazards Model via Coordinate Descent. \u003cem\u003eJournal of Statistical Software \u003c/em\u003e\u003cstrong\u003e39\u003c/strong\u003e, 1--13 (2011).\u003c/li\u003e\n\u003cli\u003eTherneau, T. M. \u0026amp; Grambsch, P. M. \u003cem\u003eModeling Survival Data: Extending the Cox Model\u003c/em\u003e. (Springer, 2000).\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-3938444/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3938444/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe utility of deep neural nets has been demonstrated for mapping hematoxylin-and-eosin (H\u0026amp;E) stained image features to expression of individual genes. However, these models have not been employed to discover clinically relevant spatial biomarkers. Here we develop MOSBY (\u003cb\u003eM\u003c/b\u003eulti-\u003cb\u003eOmic\u003c/b\u003e translation of whole slide images for \u003cb\u003eS\u003c/b\u003epatial \u003cb\u003eB\u003c/b\u003eiomarker discover\u003cb\u003eY\u003c/b\u003e) that leverages contrastive self-supervised pretraining to extract improved H\u0026amp;E whole slide images features, learns a mapping between image and bulk omic profiles (RNA, DNA, and protein), and utilizes tile-level information to discover spatial biomarkers. We validate MOSBY gene and gene set predictions with spatial transcriptomic and serially-sectioned CD8 IHC image data. We demonstrate that MOSBY-inferred colocalization features have survival-predictive power orthogonal to gene expression, and enable concordance indices highly competitive with survival-trained multimodal networks. We identify and validate 1) an ER stress-associated colocalization feature as a chemotherapy-specific risk factor in lung adenocarcinoma, and 2) the colocalization of T effector cell vs cysteine signatures as a negative prognostic factor in multiple cancer indications. The discovery of clinically relevant biologically interpretable spatial biomarkers showcases the utility of the model in unraveling novel insights in cancer biology as well as informing clinical decision-making.\u003c/p\u003e","manuscriptTitle":"MOSBY enables multi-omic inference and spatial biomarker discovery from whole slide images","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-03-04 18:58:58","doi":"10.21203/rs.3.rs-3938444/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2024-05-17T09:25:10+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2024-05-15T17:18:15+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"9ca3427e-c3e0-4f0b-a6cc-7c5ffb915972","date":"2024-05-05T09:14:57+00:00","index":"hide","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2024-03-17T18:29:55+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"746d119e-23c9-4d9b-ab76-55bf7fb99f05","date":"2024-03-08T07:56:12+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"49910fae-921a-4bc2-b95f-aa8981894a7e","date":"2024-03-07T14:30:53+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2024-03-05T11:20:56+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2024-03-05T11:07:36+00:00","index":"","fulltext":""},{"type":"editorInvited","content":"","date":"2024-02-28T18:49:48+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2024-02-28T18:43:45+00:00","index":"","fulltext":""},{"type":"submitted","content":"Scientific Reports","date":"2024-02-08T01:32:13+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"scientific-reports","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"scirep","sideBox":"Learn more about [Scientific Reports](http://www.nature.com/srep/)","snPcode":"","submissionUrl":"","title":"Scientific Reports","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"stoa","reportingPortfolio":"Scientific Reports","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"02a07227-2746-4893-a3ba-50c746660419","owner":[],"postedDate":"March 4th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2024-08-12T15:58:52+00:00","versionOfRecord":{"articleIdentity":"rs-3938444","link":"https://doi.org/10.1038/s41598-024-69198-6","journal":{"identity":"scientific-reports","isVorOnly":false,"title":"Scientific Reports"},"publishedOn":"2024-08-06 15:56:54","publishedOnDateReadable":"August 6th, 2024"},"versionCreatedAt":"2024-03-04 18:58:58","video":"","vorDoi":"10.1038/s41598-024-69198-6","vorDoiUrl":"https://doi.org/10.1038/s41598-024-69198-6","workflowStages":[]},"version":"v1","identity":"rs-3938444","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3938444","identity":"rs-3938444","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}