A Fair Experimental Comparison of Neural Network Architectures for Latent Representations of Multi-Omics for Drug Response Prediction

For our experiments, we used a publicly available drug response data set containing responses for six drugs: Docetaxel, Cisplatin, Gemcitabine, Paclitaxel, Erlotinib, and Cetuximab \autociteHossein:2020:dataset. The data set was chosen, because it had patient-derived xenograft (PDX) or patient data necessary to test the method’s translatability \autociteHossein:2019. Two external test sets are available for Gemcitabine.

The data set contains the drug response as target and data for three omics: somatic point mutations, somatic copy number profiles and gene expression profiles. Gene expressions are standardized and pairwise homogenized. Gene-level copy number estimates are binarized, representing copy-neutral genes as zeroes and genes with overlapping deletions or amplifications as ones. Somatic point mutations are also in a binary format: one for mutated genes and zero for genes without a mutation \autociteHossein:2019. For a fair comparison, the same data processing and training procedures were used for all methods. We used a variance threshold to filter genes with a small variance, as they hardly provide additional information. We set the variance thresholds in the same way as Park et al. \autocitesuperfelt.

Sharifi-Noghabi et al. acquired the data from PDX mice models \autocitepdx, The Cancer Genome Atlas (TCGA) \autociteWeinstein:2013 and Genomics of Drug Sensitivity in Cancer (GDSC) \autocitegdsc. GDSC consists of cell line data and was used for training, validation and testing because of its high number of samples. The trained neural networks were additionally tested on either PDX or TCGA to validate the algorithms’ translatability.

The characteristics of the data set per drug are summarized in Table 1.

To compare the algorithms fairly, different precautions were taken: First, the same preprocessing was performed in all experiments to provide the same input data, and second, the hyperparameters of the algorithms were optimized with an equal number of iterations by random search from a fixed grid. All algorithms draw parameters from the same grids (Table 2).

A $5\times 5$ stratified cross-validation was performed to reduce the dependence on the data splitting. 200 hyperparameter sets were created for each iteration of the outer cross-validation and each set was used for training in the inner cross-validation. The mean area under the receiver operating characteristic (AUROC) was used as performance measure. The algorithm was retrained with the best hyperparameter set on the combined train and validation sets. The trained network was used to compute the final results on the test sets from cross-validation and the external test set.

The data sets are imbalanced as they contain only few responders. In the drug discovery process, researchers are mainly interested in the positive samples to find an effective drug. To account for this requirement, the area under precision recall curve (AUPRC) \autociteauprc was additionally computed.

One requirement to use these methods is for them to work on patient data sets which are normally too small for transfer learning and still detect positive samples. This was covered by using the final model of an outer cross-validation iteration on an external patient or PDX data set.

The comparison of multiple algorithms on multiple data sets can lead to ambiguous results, if none of them performs significantly better. We compute the mean ranking as single value for a better comparison. Additionally, we compute the critical difference (CD) \autocitecritical_difference with the Nemenyi significance test with $\alpha =0.05$.

Most deep learning algorithms for multi-omics integration are based upon the concept of encoding the feature space into a lower-dimensional latent space. The encoded representation of features in the latent space is commonly called latent features. The subnetwork that computes the latent features is called encoder and computes a non-linear dimensionality reduction. The smaller dimension of the latent representation assures that the encoder does not simply learn the identity function. An autoencoder transforms an input first into the latent representation and then decodes it back into the input dimension. The reconstructed sample should resemble the input as far as possible. The difference between input and reconstruction is measured and used as loss function \autociteautoencoder.

Next, we explain briefly the six multi-omics integration architectures that were used. The first one is Early integration (EI), which serves as baseline in our experiments. EI concatenates the omics before they serve as single input of a neural network, which consists of an encoding subnetwork and a classifier layer. The network is trained by minimizing the binary cross-entropy loss function, given that our target is to classify subjects into responders and non-responders.

The schematic architecture for three omics in visualized in Figure 1.

The next method is Multi-Omics Late Integration (MOLI) developed by Sharifi-Noghabi et al. \autociteHossein:2019, which is, notwithstanding its name, an intermediate integration method. MOLI uses an individual encoder for each omics to compute the latent representations, which are concatenated and used as input for the classifier (Figure 2).

Due to the small sample size, the neural network is prone to overfitting, so MOLI uses the triplet loss \autocitetriplet_loss for regularization. The rationale behind triplet loss is that instances of the same class should have shorter distance between them than instances of different classes. To calculate the loss value, the triplet loss uses the embedding $f\left(x\right)\in R^d$ of three samples. The first one is the anchor sample $x_i^a$, the second $x_i^p$ belongs to the same class, and at last a sample $x_i^n$ is of the opposite class.

With them, the following loss function is minimized:

where $\alpha$ is a margin that defines the minimum distance between pairs of different classes. In MOLI the concatenated latent representations are used as embeddings. The final loss function is the sum of classification loss and triplet loss:

where the influence of the triplet loss is weighted by the hyperparameter $\gamma$ and $L_{Classification}$ is the binary cross-entropy. Supervised Feature Extraction Learning using Triplet Loss (Super.FELT) is a variation of MOLI, in which the encoding and classification is not performed jointly, but in two different phases \autocitesuperfelt. The first phase trains a supervised encoder with triplet loss for each omics to extract latent features of the high-dimensional omics data to avoid overfitting. The latent features are concatenated and used to train a classifier. Figure 3 (a) shows the encoding phase and (b) the classification phase of Super.FELT.

Super.FELT added a variance threshold in the feature selection to remove features with a low variance \autocitesuperfelt. It is based on the assumption that genes with a low variance might contain less important information and can be safely removed to reduce the dimensionality, hence increasing the focus on the more variable genes.

Additionally, we developed a novel extension of MOLI, which we call Omics Stacking. It is inspired by stacking and a combination of intermediate and late integration. The method stacks the results of different neural network classifier layers that use different latent features as input, with a meta-learner. The meta-learner can be any classifier, but we opted for a fully connected layer to train the neural network end-to-end.

The omics are transformed to the latent space with individual encoders, but instead of only classifying the concatenated features, Omics Stacking trains a separate classifier for each omics and the concatenated omics (Figure 4). The triplet loss is still used on the concatenated embeddings to regularize. The outputs of the classifiers are combined by a meta-learner that enables the weighting of different results to emphasize the most accurate one.

Omics Stacking combines the advantages of intermediate and late integration. It models the interaction between omics and retains the weak signals of individual omics. This method relies not only on the combined omics, but also on the individual omics and so fosters generalization. We performed an ablation study to validate different versions of Omics Stacking in Appendix B.

The former methods were developed and tested on drug response prediction, but in fact can be used for other tasks as well. Next we present two methods that were validated for the prediction of a tumor’s type or stage.

The first of these methods, Multi-task Attention Learning Algorithm for Multi-omics Data (MOMA) \autocitemoma, uses a geometrical representation and an attention mechanism \autociteattention. It is composed of three components: First, it builds a module for each omics using a module encoder. Each omics has its own module encoder consisting of two fully connected layers that convert omics features to modules. A module is represented by a normalized two-dimensional vector.

Second, it focuses on important modules between omics using a module attention mechanism. This mechanism is designed to act as a mediator to identify modules with high similarity among multiple omics. The relevance between modules is measured by the cosine similarity and is converted to a probability distribution with the softmax function. The distributions are then used to create an attention matrix that stores the relationship information between modules of different omics. To highlight important modules, the module vectors are multiplied by the attention matrix (Figure 5).

Subsequently, fully connected layers and the logistic function are applied to flatten the multidimensional vectors and to compute the final probabilities for each omics \autocitemoma.

MOMA is trained with cross-entropy loss between the true label and the omics specific outputs. After the training of the neural network, a logistic regression is fit on the omics specific outputs to generate the combined prediction.

The last architecture we tested OmiEmbed \autociteomiEmbed, which is based on a supervised variational autoencoder (VAE) \autocitevae. OmiEmbed was developed as a unified end-to-end multi-view multitask deep learning framework for high-dimensional multi-omics data. The method learns the latent features with the auxiliary unsupervised task of reconstructing the input omics. The trained latent features can be used for one or more supervised tasks. The overall architecture of OmiEmbed comprises a deep embedding module and one or multiple downstream task modules (Figure 6).

The loss function of OmiEmbed is the sum of two parts: a reconstruction loss and the loss of the downstream task. $L_{embed}$ is the unsupervised loss function of the VAE:

$BCE$ is the binary cross-entropy to measure the differences between the input $x$ and the reconstruction $x^{\prime }$ and is computed individually for each of the $M$ omics. $D_{KL}$ is the Kullback–Leibler divergence between the learned distribution and a standard normal distribution.

The embedding loss function is used together with the loss of the downstream task, which is in our case classification:

where $L_{CE}$ is the cross-entropy loss and $\lambda$ a balancing weight.

Instead of training all layers at the same time, the network learns in three phases: First, only the VAE is trained. In the second phase, the VAE weights are fixed and only the downstream network is trained and in the last phase the complete network is fine-tuned.

Table 3 summarizes the components and architectures of the described methods.

This work was funded by the German Federal Ministry for Education and Research as part of the DIASyM project under grant number [031L0217A].

TH implemented the methods and experiments and wrote the manuscript.
SK supervised the study, developed the concept and wrote and previewed the manuscript.