Attention-based Interpretable Regression of Gene Expression in Histology

Interpretability of deep learning is widely used to evaluate the reliability of medical imaging models and reduce the risks of inaccurate patient recommendations. For models exceeding human performance, e.g. predicting RNA structure from microscopy images, interpretable modelling can be further used to uncover highly non-trivial patterns which are otherwise imperceptible to the human eye. We show that interpretability can reveal connections between the microscopic appearance of cancer tissue and its gene expression profiling. While exhaustive profiling of all genes from the histology images is still challenging, we estimate the expression values of a well-known subset of genes that is indicative of cancer molecular subtype, survival, and treatment response in colorectal cancer. Our approach successfully identifies meaningful information from the image slides, highlighting hotspots of high gene expression. Our method can help characterise how gene expression shapes tissue morphology and this may be beneficial for patient stratification in the pathology unit. The code is available on GitHub.


Introduction
The wide variability of existing interpretability methods, in the form of posthoc explanations [1,2] or models with ad-hoc transparency constraints [3], has been mainly dedicated to ensuring the safety and reliability of opaque deep neural networks. In medical applications such as digital pathology, saliency maps highlighted the relevance of anomalous nuclei in the detection of tumorous tissue [4], and concept-based analyses confirmed that clinically relevant measures on nuclei area and appearance are learned as intermediate features [5]. Undesired hidden biases and behaviours were detected and corrected to improve model performance [6,7,8]. Only recently, interpretability techniques were proposed to uncover unknown insights about models with super-human performance, e.g. the algorithm AlphaZero [9] defeating the world chess master. Similarly, deep learning models can predict biomarkers invisible to the human eye [10,11]. For instance, DNA mutations [10] and gene expression profiles were inferred from hematoxylin and eosin (H&E) stained whole slide images (WSIs) [11], demonstrating that complex transcriptomic patterns can be captured from images without the support of sophisticated sequencing machines.
While ensuring patient safety remains the top priority in high-risk clinical environments [12,13], scientific discovery may benefit from unravelling the complex association between RNA expression and tissue microscopy. Identifying the histological patterns that are predictive of gene expression could have an important impact on diagnostic routines, facilitating early patient stratification and aiding the identification of molecular subtypes that are informative of prognosis, survival, and treatment response [14]. For example, clinical trials show that some colorectal (CRC) cancer patients with hypermutated microsatellite instability (MSI) may respond better to immunotherapy than chemotherapy [15].
Our central question is whether we can clarify which regions in CRC H&E slides are the most informative about RNA transcriptomics by using a trainable attention mechanism [16]. It is yet unclear whether the existing algorithms [11,17] infer the bulk average expression of RNA from hotspots of high expression or as a uniform distribution on the entire slide. To obtain interpretable insights, we learn the expression of individual genes rather than the entire transcriptomic profile at once, training several attention-based multiple instance regression models independently. This enables a fine grained analysis of each gene individually and reduces the need for extremely large dataset sizes [11,17] In our results, the attention mechanism brings transparency to the morphological patterns in the tissue that are learned by the model. Hotspots of high expression are highlighted in most of our visualisations and hypotheses can be formulated to relate the histological appearance of tissue to varying gene expression. We show that meaningful information is successfully filtered from the WSIs by the attention, leading to more accurate gene expression estimates and patient stratification. A reduction of 10% in the regression error is seen across genes and patients.

Datasets
We use images of colon adenocarcinoma (COAD) and rectal adenocarcinoma (READ) that are publicly available together with matched transcriptomic profiles at The Cancer Genome Atlas (TCGA) 5 . Each biopsy sample is split into three adjacent portions, two of which are used to generate H&E-stained frozen tissue slides, called top-section (TS) and bottom-section (BS), used to verify the presence of sufficient tumour content before sequencing. The central section is used for RNA sequencing. Differently from [11,17], we focus on frozen tissue sections rather than on diagnostic slides, since they are directly adjacent to the sequenced tissue, and thus constitute the best available representation of the RNA profiles [18]. The WSIs are preprocessed with HistoQC [19] to mask out the background and blurred locations.
Gene expression profiles are obtained from the UCSC Xena Browser [20], which links to the Genomic Data Commons [21] version of the TCGA COAD and READ projects. The High Throughput-Sequencing (HT-Seq) raw counts are log2-based Fragments per Kilobase of transcript per Million mapped reads (FPKM) normalised. Because of a distributional shift among institutions, we retain all the patient measurements from a single institution, ignoring duplicates for 23 patients. We focus on the 45 genes in the ColoType signature [22], a gene set that is predictive of CRC prognosis [23] and clinico-pathological variables [24]. The list of genes is in the Appendix Table A.1. Single-cell RNA sequencing data are also used for validation 6 . Namely, we use the profiling of 969 single-cells from the CRC resected primary tumours of 11 patients in [25].
As in [26], we exclude patients in preoperative therapy and the rare subtypes of neuroendocrine and signet cell tumours. In total, we use 774 WSIs at 20X magnification from 364 patients. The test set is built by selecting randomly 82 patients and the remaining 282 patients are used for five-fold cross validation. Our model is an adaptation of the attention-based multiple instance learning approach in [16], which follows the implementation in [27]. We consider a patient p with associated H&E-stained scans of the top and bottom frozen tissue sections x t and x b , both images in R w×h , where w is the image width and h is the height. The associated bulk-RNA expression for a given gene g is y g ∈ R. To deal with the large input sizes of WSIs, x t and x b are represented as a single collection of N non-overlapping adjacent patches obtained by a sliding window {x i } N i=0 [28]. The task is to learn a permutation invariant mapping to the real-valued gene label y g from this joint collection. We assume that each collection contains at least one instance associated with the label [29]. The gene expression y g can be obtained as a linear-weighted combination of the y i,g predictions for each of the x i patches, where the weights are given by the trainable attention. Differently from [16,27], the attention weights {a i,g } N i=0 are optimised to predict the continuous gene expression label rather than a binary outcome.

Multiple Instance Regression of Gene Expression
The model, illustrated in Figure 1, comprises: (i) a convolutional backbone to obtain low-dimensional representations of the input patches, i.e. a ResNet18 [30] pretrained on ImageNet [31]; (ii) a dense layer to predict the gene expression y i,g for each input x i ; (iii) an attention network to learn the attention weights.
The output of our model is the expression of a single gene, and we train a different model for each gene in our selection. This is opposed to the approach of the SoA [11], where regression values for over 30000 gene outputs are optimised in a single training. To allow for a fair comparison with the SoA, we adapt their model to output the expression of a single gene at a time. The adapted model predicts the gene expression for each patch, and the predictions are then aggregated by a weighted average with larger weights given to large-valued patch predictions as in [11]. Our trainable attention removes the need for this heuristic.

Attention-based Model Interpretability
The trainable attention provides transparency on how the model filters the information in the WSIs. Salient regions are identified without the need for an explicit localisation module [16,32]. However, recent debate argued that multiple plausible attention-based explanations may exist and that attention should be interpreted carefully [33,34]. The research in [35] addressed that debate, mentioning that attention weights learned by multiple models should be ensembled by either max or average pooling to reduce the risks of obtaining misleading interpretations. We thus average the attention weights of the models trained on the five folds, ensuring a trustworthy interpretation that focuses only on signalbearing instances.

Evaluation of performance and interpretability
Model performance is evaluated in multiple ways. First, the Median Average Percentage Error (MAPE) of our predictions is compared against the SoA. Both models are also compared to the lower bound of random guessing given by predicting the mean of the training labels. Finally, the significance of our results is verified by computing the non-parametric one-tailed Wilcoxon test between the patient-wise percentage error made by our model and the SoA. Successfully regressed genes are identified as in [11] by evaluating the Pearsons's correlation coefficient ρ between the true and the predicted labels.
The evaluation of model interpretability is non-trivial [17,35]. For morphological patterns that are known to be associated with RNA abnormal expression, e.g. mucinous tissue showing high levels of MUC2, we verify that our explanations confirm the existing knowledge. For the newly identified patterns, however, there is little to no ground truth available. In this case, we evaluate how well patch-wise attention weights correlate to the ground truth obtained from analysing single-cell data. We identify, for instance, genes that are co-expressed in single-cell data according to their Pearson's correlation, and we verify that the patch-wise attentions in overlapping WSI locations correlate for those models. Finally, the benefits of the attention are quantified by quantifying the performance improvement on an auxiliary task, i.e. inferring MSI status. The higher the quality of the attention-based localisation, the more accurate we expect to be the MSI-based stratification of the patients.

Network Training
The models are trained by optimising the mean squared error loss between the predicted gene expression valueŷ g and the label y g . Stochastic gradient descent (SGD) is used with standard hyperparameters (learning rate 0.0001, momentum 0.9, weight decay 0.01) and early stopping (12 epochs patience). The gradient updates affect only the dense layers and the attention mechanism. MSI status is learned in Sec 3.4 by two additional dense layers trained by SGD minimisation of the weighted binary cross-entropy. For a single model trained on a GPU Tesla V100 for 5 hours, we estimate a carbon footprint of 0.65 kgCO 2 [36].

Quantitative Model Evaluation
The gene expression estimates obtained by our method are more accurate than the SoA. Across all genes and patients, our model obtains a reduction of 10% in the MAPE, i.e. 0.65 ± 0.08 in our model against 0.72±0.12 in the SoA. Random guessing as described in Sec.   Figures 3 and 4 show the predictions and the distribution of the attention for our best gene models on the input WSI. The attention identifies hotspot regions, rather than being uniformly distributed over all patches. Figure 5 visualises more in detail the morphological patterns in the highlighted hotspots of high and low gene expression. We retrieve, for instance, the patches that received the highest normalised attention weights, namely the highest N p a i,g ∀i, p in the testing set, where i = 1, . . . , N p , and N p is the total number of patches for patient p 7 . Additional visualisations are in the Appendix and the code repository 8 .

Attention-based Identification of Hotspots and Patterns
Genes that are co-expressed in single-cell data show high attention on the same WSI areas. For example, NRP2 and MUC2 show a correlation value of ρ = 0.40 at the single-cell gene expression level. Similarly, they display a correlation of 0.63 in the patch-wise attention values. Conversely, MUC2 is not co-expressed at the single-cell level with AOC3 and COL8A2, and the attention localise on 7 The normalisation enables the comparison across patients, since

Quantitative Evaluation of the Attention
The MSI status is inferred from the gene expression values regressed from the H&E slides. The input noise is reduced by selecting only the genes that were predicted with Pearson's p-value< 0.05 by each method. This leads to 37 genes with our method and 35 with the SoA from the initial pool of 45 genes. We perform 5x2 cross-validation to obtain more stable estimates through repetitions on non-overlapping training sets [37]. We observe an increase in the AUC from 0.60 ± 0.08 of the SoA to 0.71 ± 0.05 with our attention-based model. In comparison, using the ground-truth sequencing data yields AUC 0.86 ± 0.04.

Discussion
The experiments in Sec. 3.2 quantify the accuracy of our bulk-RNA predictions from WSIs in terms of the prediction error, i.e. MAPE, and correlation with the ground-truth labels, i.e. Pearson correlation ρ. Our trainable attention mechanism reduces the error for the majority of the analysed genes, achieving significant improvements over the SoA for the SPINK4, QPRT, MUC2, SDC2 and AOC3 genes. Besides accuracy improvements, learning the attention removes the use of a priori heuristic choices and hyperparameter tuning in [11]. The main contribution is the analysis of the signals that are picked by the models to obtain our performance improvements. Focusing on one gene at a time enables a fine-grained analysis of the model attention patterns since the entire model capacity is optimised to capture salient features of individual genes. The attention maps in Sec. 3.3 give interesting insights. For instance, Figure 3 shows that the attention localises over hotspot regions rather than being uniformly spread across the slide. This answers the questions raised in [17] regarding the distribution of the information related to gene expression. As expected, high attention is seen in overlapping WSI regions for genes that are co-expressed in the same cell and in differing regions for genes that are not co-expressed. Moreover, the highest attention weights point to visibly different histological patterns, as shown in Figure 5. For example, the high-expression patches of COL8A2 present higher concentrations of stroma than the low-expression ones.
An interesting example is that of MUC2, a particularly well-known gene that impacts tissue histology.MUC2 over-expression leads to high amounts of visible intra-and extra-cellular mucin, characterizing the tissue as mucinous adenocarcinoma. An example of mucinous adenocarcinoma is shown in Figure 4b. Our model predicts high values of MUC2 for the regions in the slides where the mucinous tissue is the most visible. When these are further weighted by the attention, the localisation is highly focused on a few hotspots of high expression. From a quantitative standpoint, the impact of the attention mechanism for this specific gene is remarkably significant. The SoA predictions for MUC2 are similar to random guessing (with MAPE 0.95 and ρ 0.13), whereas our method significantly reduces the MAPE to 0.73 and the correlation with ground truth values, i.e. ρ, to 0.46 with p-value< 0.0001. This result is also beneficial to the prediction of MSI, which is known to correlate with the presence of mucin [26].

Conclusion
We proposed an attention-based multiple instance regression model to infer bulk gene expression of tissue sections from H&E histology slides of colorectal adenocarcinoma. The automated regression of transcriptomics from H&E is not yet a feasible replacement for sequencing the tissue, but the analysis of interpretability patterns can highlight some histological patterns associated with specific gene expression levels. We show here how the attention mechanism successfully filters informative content from the vast amount of information in the WSIs, leading to significantly lower errors and interesting insights into how gene expression impacts tissue morphology. Further developments of this method should introduce pathologists in the loop to validate the identified patterns and further investigate their association to differentially expressed genes. In turn, this could help pathologists to identify cancer molecular subtypes from WSIs and consequently stratify patients for targeted treatment.

E. Alpaydmn. Combined 5 × 2 CV F Test for Comparing Supervised Classification
Learning Algorithms. Neural Computation, 11 (8):1885-1892, 1999.  Detailed results are reported in Table B.1, which was not included in the main paper because of the limited space. The model of the gene expression for NRP2 leads to the best results overall, with MAPE=0.47 ± 1.19 and ρ = 0.62 against the SoA MAPE=0.51 ± 1.17 and ρ = 0.61 of the baseline.

C Additional Visualizations
Additional visualizations of the spatialized prediction maps with and without attention-based weighting are given in Figure