==== Front Bioinformatics Bioinformatics bioinformatics Bioinformatics 1367-4803 1367-4811 Oxford University Press 37387182 10.1093/bioinformatics/btad232 btad232 Biomedical Informatics AcademicSubjects/SCI01060 AttOmics: attention-based architecture for diagnosis and prognosis from omics data https://orcid.org/0000-0001-7199-5242 Beaude Aurélien IBISC, Université Paris-Saclay, Univ Evry, 23 Boulevard de France, Evry-Courcouronnes 91020, France Artificial Intelligence & Deep Analytics, Omics Data Science, Sanofi R&D Data and Data Science, 1 Av. Pierre Brossolette, Chilly-Mazarin 91385, France https://orcid.org/0000-0002-1785-0234 Rafiee Vahid Milad Sanofi R&D Data and Data Science, Artificial Intelligence & Deep Analytics, Omics Data Science, 450 Water Street, Cambridge, MA 02142, United States https://orcid.org/0000-0002-5641-4152 Augé Franck Artificial Intelligence & Deep Analytics, Omics Data Science, Sanofi R&D Data and Data Science, 1 Av. Pierre Brossolette, Chilly-Mazarin 91385, France https://orcid.org/0000-0001-6278-1680 Zehraoui Farida IBISC, Université Paris-Saclay, Univ Evry, 23 Boulevard de France, Evry-Courcouronnes 91020, France https://orcid.org/0000-0002-5606-8296 Hanczar Blaise IBISC, Université Paris-Saclay, Univ Evry, 23 Boulevard de France, Evry-Courcouronnes 91020, France Corresponding author. E-mail: aurelien.beaude@universite-paris-saclay.fr (A.B.); blaise.hanczar@univ-evry.fr (B.H.) 6 2023 30 6 2023 30 6 2023 39 Suppl 1 ISMB/ECCB 2023 Proceedings i94i102 © The Author(s) 2023. Published by Oxford University Press. 2023 https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract Motivation The increasing availability of high-throughput omics data allows for considering a new medicine centered on individual patients. Precision medicine relies on exploiting these high-throughput data with machine-learning models, especially the ones based on deep-learning approaches, to improve diagnosis. Due to the high-dimensional small-sample nature of omics data, current deep-learning models end up with many parameters and have to be fitted with a limited training set. Furthermore, interactions between molecular entities inside an omics profile are not patient specific but are the same for all patients. Results In this article, we propose AttOmics, a new deep-learning architecture based on the self-attention mechanism. First, we decompose each omics profile into a set of groups, where each group contains related features. Then, by applying the self-attention mechanism to the set of groups, we can capture the different interactions specific to a patient. The results of different experiments carried out in this article show that our model can accurately predict the phenotype of a patient with fewer parameters than deep neural networks. Visualizing the attention maps can provide new insights into the essential groups for a particular phenotype. Availability and implementation The code and data are available at https://forge.ibisc.univ-evry.fr/abeaude/AttOmics. TCGA data can be downloaded from the Genomic Data Commons Data Portal. public–private partnership 2021-1047 ==== Body pmc1 Introduction The disruption of different biological processes (BPs) can negatively affect an organism and lead to a disease state. Early diagnosis plays an important role in precision medicine in order to improve clinical decision-making. The development of high-throughput methods influenced precision medicine by enabling easy access to a large amount of biological information for each patient, known as omics profile. Omics profiles are high-dimensional complex signatures resulting from interactions of many molecular entities. The first common step in most machine learning approaches used in precision medicine is a feature selection procedure that reduces the data’s size to construct a classifier from a single-omics (Kourou et al. 2015). As the selection procedure is decoupled from the prediction task, only selected features are used for downstream predictions (Liu et al. 2019). Using only selected features limit the model’s capacity to extract hidden information from the omitted features. Deep learning, on the contrary, can extract and exploit the complete information from all features and their interactions. This characteristic may be useful to achieve better-predicting performances. Following the recent successes of deep learning in computer vision or natural language processing (LeCun et al. 2015), different deep learning architectures were successfully applied to omics data. It allows for a high level of abstraction of features with nonlinear modeling and can handle complex dependencies in data to create informative representations. Assuming that omics data does not have any particular structure, unlike images or texts, multilayer perceptrons (MLP) were used to perform predictions (Yu et al. 2021) and autoencoders (AE) for dimension reduction (Gore and Azad 2022). Other approaches tried to integrate a structure in the model by embedding biological knowledge and applying convolutional neural networks (CNN) (Elbashir et al. 2019) or graph neural networks (GNN) (Ramirez et al. 2020). Cellular functions are governed by the combined action of multiple molecular entities which are specific to a patient. The expression of one gene may impact the expression of other genes differently in different patients. With classical deep learning approaches, these interactions which are learned during training, are assumed to be identical for all patients in the inference phase. It would be more beneficial to compute feature interactions that are specific to each patient. Self-attention can be used to improve the representation of the features vector by incorporating dynamically computed relationships between elements of the vector. It has been shown that the transformer architecture’s promising results extensively rely on attention mechanisms (Vaswani et al. 2017). Here, we propose a new method based on the self-attention mechanism to capture interactions between different molecular entities in order to predict the phenotype of patients, e.g. cancer types or the risk of death from omics data. Using a self-attention mechanism allows the model to capture feature interactions specific to each patient dynamically. Applying self-attention on high dimensional vectors such as omics profiles is challenging as the self-attention memory requirements scale quadratically with the number of elements. To overcome this problem, we propose to consider groups of features and apply self-attention to these groups. The architecture was tested on three different omics data: gene (mRNA) expression, methylation (DNAm), and micro-RNA (miRNA) expression and was compared with state-of-the-art deep learning-based methods. The results show that our proposed architecture better considers feature interactions in omics data and improves the model performance. 2 Related work Different deep learning approaches have already been tested on omics data: CNN, GNN, MLP, AE, and variational autoencoder (VAE). Hanczar et al. (2022) showed that MLP outperforms the classic machine learning methods on large gene expression datasets. Yu et al. explored different MLP architectures by varying the number of neurons in each layer and the number of layers. They showed that wider networks perform better than deeper ones. DeepCC (Gao et al. 2019) applies an MLP on biologically informed features by transforming gene expression data into a functional spectrum, i.e. a list of enrichment scores calculated by gene set enrichment analysis. There are also other approaches that included biological knowledge in the design of the neural networks by restricting connections between neurons to known biological relations, such as Gene Ontology (GO) (Bourgeais et al. 2021) or REACTOME (Hao et al. 2021). A VAE unsupervised training has been used as a pre-training for an MLP classifier (Levy et al. 2020), and the VAE latent space has also been directly passed to a classifier (Wang and Wang 2018). Another approach constrained the latent space to learn relevant features for the classification by end-to-end training a network on both the unsupervised and supervised tasks (Gore and Azad 2022). Promising results of CNN in computer vision inspired its application in precision medicine. Different strategies were developed to create a 2D image from an expression vector. Some approaches reshaped a 1D omics vector into a 2D image to exploit the capacity of CNN architectures to extract relevant visual patterns (Elbashir et al. 2019; Mostavi et al. 2020; Rukhsar et al. 2021). Ma and Zhang (2019) created an image by transforming an expression vector into a tree map based on the Kyoto encyclopedia of genes and genomes Brite structure. Instead of forcing a 2D representation, 1D convolution has been applied to an ordered expression vector (Mostavi et al. 2020; Zhao et al. 2020). In Zhao et al. (2020), the expression vector was reordered according to the chromosomal locations before applying a 1D inception architecture. Expression profiles have also been represented as graphs to represent the interactions of the different molecular entities. A graph convolutional network based on a co-expression network or a protein–protein interaction (PPI) network was used to predict cancer types from gene expression (Ramirez et al. 2020). Ramirez et al. (2021) explored the combination of a co-expression (CoExp) graph and a graph constructed from the GeneMania database. In Levy et al. (2021), an architecture inspired by the capsule network was used to predict the central nervous system tumors subtype. Methylation features have been grouped into capsules to create context-specific embeddings, and dynamic routing was then applied to make a prediction. Deep learning architectures using attention mechanisms have been little explored for the application of omics data. Some approaches helped the network to focus on relevant genes for the predicted phenotype by computing feature importance scores with a small neural network (Beykikhoshk et al. 2020; Lee 2022). Those methods were inspired by the attention mechanism but did not use the original dot product self-attention (Vaswani et al. 2017) layers directly on the omics data. Computing self-attention on a high dimensional vector is hardware limited as memory requirements scale quadratically with the number of elements. The Gene transformer (Khan and Lee 2021) was the first architecture to apply self-attention to mRNA data. The authors proposed to use 1D convolution layers combined with maximum pooling to reduce the dimension of the gene expression vector. Using a pooling layer is equivalent to a dimension reduction that does not consider all possible feature interactions. In our approach, AttOmics, we propose to embed groups of features in a lower dimension by considering all interactions inside this group and a new way of applying self-attention (Vaswani et al. 2017) to omics data that takes into account inter-group interactions. AttOmics can be applied to vectors of various sizes, and consequently, detect feature interactions in different omics data. 3 Model architecture 3.1 Architecture details The model includes a grouping module and an encoder followed by a predictor, illustrated in Fig. 1a. Instead of considering each feature individually, features are divided into different groups. The encoder is a stack of n blocks used to construct a new representation of the inputs. Each block is formed of a grouped fully connected network (gFCN) module where each group is projected into a lower dimension with a fully connected network (FCN). Segregating features in groups restrict the potential interactions between features to the ones inside the same group. Multihead self-attention (MHSA) is applied to the set of groups to recover all possible interactions between groups. Around the self-attention block, a residual connection is added before applying a normalization. The encoder output is transmitted to an FCN used as the predictor. Figure 1. The AttOmics architecture is composed of a grouping module, an encoder and a predictor (a). The grouping module transforms the input features into a set of different groups. In each of the n encoder blocks, each group is projected into a lower dimensional space (b). Interactions between the different groups is computed with the MHSA (c). In each of the h heads of the MHSA, a scaled dot product attention is computed between the different groups (d). A residual connection is added around the MHSA before applying a normalization. The new representation obtained with the encoder is transmitted to the predictor, an FCN followed by a normalization. 3.1.1 Grouped FCN Let X∈Rp be a training example where p is the number of features and Y is the associated label. The training example X is split into groups according to a grouping strategy (see Section 3.2), XG={Xgi}1≤i≤k, where k is the number of groups. For each group Xgi, a group embedding is independently computed as Xgi′ by projecting it into an s-dimensional space with an FCN, a succession of fully connected layers (FCL) (Fig. 1b). Each FCL is the composition of an affine transformation of its inputs with a rectified linear unit (ReLU) activation function: FCL(x)=ReLU(Wx+b)=max(0,Wx+b) After processing each group Xgi by the successive FCL, we obtain the set of group embeddings X′G: where Xgi′=FCN(Xgi). XG′={Xgi′∈Rs}1≤i≤k, Each group projection is only computed using elements from the same group. To create a representation of the expression vector based on all possible interactions, MHSA is then applied to X′G. 3.1.2 Multihead self-attention MHSA is applied to construct a new representation of the groups, U={Ugi}1≤i≤k, by allowing them to interact with each other (Fig. 1c). MHSA is performed with h different heads to learn different types of interactions. For each head j, self-attention is applied to each group gi (1≤i≤k), in order to obtain: U(j)={Ugi(j)∈Rl}1≤i≤k, where l=sh∈N. Ugi(j) is defined by where Agi(j) is the attention vector computed by the usual dot product attention (Fig. 1d) (Vaswani et al. 2017): Ugi(j)=Agi(j)·[Xg1′·WjV,…,Xgk′·WjV]T, Agi(j)=softmax([Agi,g1(j),…,Agi,gk(j)]),Agi,gk(j)=(Xgi′·WjQ)T·(Xgk′·WjK)s. Projection matrix WjQ (respectively WjK and WjV) maps the group Xgi′, from an s-dimensional space to an l-dimensional space. In the transformers formulation Xgi′·WjQ, Xgk′·WjK, and Xgi′·WjV are called, query, key, and value, respectively. Each element of Ugi is obtained by concatenating the representation of all groups in the different heads and projecting each group to an s-dimensional space using a projection matrix WO∈Rs×s as: Ugi=concat(Ugi(1),…,Ugi(h))·WO, 3.1.3 Residual connection and normalization The value of Xgi′ is added to Ugi, through a residual connection to prevent vanishing gradients. The last step in the encoder module consists of applying a normalization to obtain the final representation Zgi of group gi defined as Zgi=Norm(Xgi′+Ugi). The output of an encoder block is Z={Zgi∈Rs}1≤i≤k, which is a representation of the groups capturing their interactions. 3.1.4 Prediction module The vectors Zgi are concatenated into a new vector Z′∈Rks. The output of the encoder Z′ is then fed to a FCN followed by a normalization layer to predict the cancer type or the prognosis Y^ (Fig. 1a). For classification tasks, the output layer has one neuron per class, and a softmax activation function is applied to get the probability vector P=[pc]1≤c≤M, where M denotes the number of classes. For the survival analysis, the output is a single neuron with a linear activation function. 3.2 Grouping strategies The AttOmics architecture can be applied to any group specification. We explore different grouping strategies such as random groups, groups obtained with clustering, groups based on biological information like the GO (Gene Ontology Consortium 2021) or the hallmarks collection available in MSigDB (Liberzon et al. 2015). 3.2.1 Random With the random strategy, groups are formed by randomly sampling the input features in groups of similar sizes, p/k. 3.2.2 Gene ontology We used BP gene ontology as it groups different molecular activities in a shared process which are more likely linked to the same cancer phenotype. To avoid possible problems with selecting the GO terms (i.e. groups) of interest, we restrict ourselves to terms available in GO slims. Inside the BP slim ontology, a gene can belong to more than one group; on average, they belong to two groups. Before applying self-attention, each group must be projected into the same dimensional space. Each group is projected with a different number of layers to have the same reduction ratio across different groups. This grouping strategy can only be applied to mRNA data. 3.2.3 Hallmarks In the MSigDB hallmarks collection (Liberzon et al. 2015), there are 50 groups. Each one represents a well-defined BP. Each group is projected with a different number of layers to ensure identical reduction ratio across different groups. This grouping strategy can only be applied to mRNA data. 3.2.4 Clustering The clustering strategy groups features based on their expression levels. Traditional clustering methods, like K-Means or hierarchical clustering, can return sets of highly unbalanced clusters that may negatively affect the efficiency of our model. Large groups would require many parameters to be projected into a space with a dimension lower than the smallest group. Group unbalances would also imply a high compression of larger groups and almost no compression for the smallest group. To prevent this, we used constrained K-means clustering to ensure comparable group sizes (Bradley et al. 2000). 3.3 Model training For classification problems, our model is trained end-to-end with a weighted cross-entropy loss to account for class imbalance: where wc denotes the weight (inversely proportional to the size) of class c∈{1,…,M} and θ the model parameters. L(θ)=−∑c=1MwcYc log (pc), For survival analysis, our model is end-to-end trained with a partial log-likelihood loss, as proposed in DeepSurv (Katzman et al. 2018): where δi specifies if the event occurred for patient i, Ti represents the time associated to the event and Nδi=1 is the number of patients for which the event occurred (δi=1). ηi=eY^i is the predicted risk for patient i. R(Ti)={j:Tj>Ti} is the risk set, the set of patients who are still at risk of death at time Ti. L(θ)=1Nδi=1∑i:δi=1(Y^i−log ∑j∈R(Ti)ηj), 4 Experiments 4.1 Data TCGA data were used to evaluate our proposed approach AttOmics. We collected DNA methylation, gene expression, and miRNA expression data for 8416 patients of 19 different cancers and 361 normal samples from the GDC Data Portal (https://portal.gdc.cancer.gov/). FFPE samples and bad replicates were removed according to TCGA consortium recommendation. Methylation data was restricted to the probes common to both HumanMethylation27 and HumanMethylation450 platforms. No feature selection was applied, and data were standardized to a zero mean and unit variance. Patients with incorrect survival information were removed: 8349 patients were available for survival prediction. A total of 70% of the data are used as a training set, 15% forms the validation set, and the remaining 15% forms the test set while preserving the proportion of each cancer. The training set is used to perform two predicting tasks: phenotype prediction, 19 different cancers and normal, and survival risk prediction. 4.2 Comparative study For a comprehensive and comparative evaluation, we choose three deep learning architectures for comparison: CNN, GNN, and MLP. For the CNN (CNN1d), we ordered features based on their position on the genome, then used a 1D convolution, followed by a ReLU activation and a maximum pooling. For the GNN architecture, two graphs were used: PPI (GNN—PPI) and co-expression (GNN—CoExp) graphs. The PPI graph is based on data available in the STRING database (Szklarczyk et al. 2020) and was constructed by retaining only high-confidence links: edges with a score higher than 700. The CoExp graph was constructed similarly to Ramirez et al. (2020). The Spearman correlation matrix between gene expressions was computed. If the correlation was higher than a threshold and the associated P-value was lower than .05, then an edge between the two features was added to the graph. For mRNA and miRNA, the correlation threshold was set to 0.6. For DNAm, a 0.7 correlation threshold was used. Self-loops were not considered in the graph construction, and isolated nodes were removed. The PPI graph and the CoExp graph for mRNA have 9384 genes in common. Each graph is described in the Supplementary Table S2. MLP architecture has two hidden layers with ReLU activation and makes use of batch normalization. We also consider three state-of-the-art non-deep-learning models for comparison: support vector machine (SVM), random forest (RF), and extreme gradient boosting (XGBoost). For the non-deep-learning approaches, the 2000 most discriminative features are selected with a t-test-based selection. The hyperparameters of each approach are tuned on each omics data with a random search to achieve the best performances. The different values tested for each parameter are defined in the Supplementary Table S3. For each hyperparameter at each search iteration, a value is randomly drawn from the defined range. A model is constructed using these parameters, trained on the training set, and evaluated on the validation set. The selected hyper-parameters for each architecture are presented in the Supplementary Table S5. AttOmics is trained end-to-end using the Adam optimizer with a learning rate of 0.0001 and a batch size of 512. The maximum number of epochs was set to 100. An early stopping strategy is deployed to avoid over-fitting with a patience of 8 and a delta of 0.001 on the validation metric between two epochs. For the classification task, models were evaluated with the error rate. Prognosis prediction is evaluated with the concordance index (Harrell et al. 1996). It estimates that for a pair of individuals, the predicted risks, η, are concordant with their actual survival times. C-Index=∑i,j1Tjηiδj∑i,j1Tj