Knowledge-guided artificial intelligence technologies for decoding complex multiomics interactions in cells

1. Learning the sequence preference of proteins on DNA and RNA binding

The efficiency of the interaction between proteins and DNA is often determined by the compatibility between protein structure and DNA sequence [11]. In particular, many DNA-binding proteins, including transcription factors (TFs), have an intrinsic preference for short stretches of DNA with unique ordering of nucleotides, or sequence motifs [12]. Sequence motifs bound by TFs, namely TF-binding sites (TFBSs), are conventionally determined by high-throughput experiments such as chromatin immunoprecipitation combined with DNA sequencing (ChIP-seq) [13] or SELEX-seq [14] followed by motif-finding algorithms including MEME [15] and HOMER [16]. The resulting TFBSs are often represented as position weight matrices (PWMs) that encode the relative base preference for each position within the motif.

DeepBind [17] pioneered the application of deep learning for binding prediction of TF or RNA-binding proteins (RBPs) using high-throughput experimental data. For each TF or RBP, individual DeepBind models were trained to classify whether the corresponding protein would bind to the given nucleotide sequence in a supervised manner. More specifically, a number of short nucleotide sequences that were experimentally validated to be bound by the protein were prepared along with random negative sequences. The model weights were optimized to make the model produce correct predictions for either of the 2 classes.

Although no explicit prior knowledge was incorporated in DeepBind training, it is worth noting that the choice of model architecture provides minimal guidance or inductive bias for the model to successfully detect motifs. DeepBind adopts a convolutional neural network (CNN) architecture [18], which was originally developed for computer vision applications such as image recognition. Key operators at the CNN core are convolutional filters, which iteratively slide throughout the image and produce a scalar for each of its patches. The scalar value is computed as a dot product between the weight of the filter and the pixel intensity of the image patch, thus producing high values at the patch that are preferred by the filter. A simple but powerful analogy can be derived from this formulation when we consider a convolutional filter as a learnable PWM that slides through the sequence and detects the short substring that is compatible with the PWM. The PWM denoting TF-binding motifs can be learned by optimizing the model to respond positively to TF-binding sequences.

DeepBind models showed improved accuracy of TF-binding specificity prediction compared to conventional models. Furthermore, it could predict mutations that would disrupt TF binding, and many have already been reported as disease-associated mutations. Since the successful introduction of deep learning models for the task of revealing protein-DNA/RNA interactions, many variants of DeepBind models have been developed. DeeperBind [19] appended recurrent neural network (RNN) architecture [20] after convolutional layers to capture dependencies between motifs as well as positional bias in probes used for high-throughput screens such as protein-binding microarrays. DanQ [21] is another model that uses the hybrid structure of CNN and RNN, but it differs from the aforementioned models in that it is trained by multitask learning [22]. In other words, it is a single unified model that predicts the binding specificity of about 1,000 TFs. The utility of multitask learning is especially highlighted in biological applications because the most fundamental dogmas are shared across every cell. More recently, a CNN-based model without pooling, called BPNet [23], along with a motif discovery method TF-MoDISco [24], further captured the detailed syntaxes of TF-binding motifs, such as helical periodicity of motifs and multiple motifs repeatedly occurring at a fixed distance because of the cooperative DNA binding of interacting proteins. The fine resolution of the analysis could be achieved using high-resolution TF-binding signals from ChIP-nexus [25] experiments, underscoring the importance of training data quality.

2. Learning DNA sequence determinant of gene expression and RNA splicing from scratch

As discussed above, TFBSs mark hotspots of protein-DNA interactions and potentiate the activation of downstream genes upon binding to the corresponding TFs. In other words, TFBSs can serve as major sequence determinants of gene expression. If a deep learning model is powerful enough, by training the model to predict whether a gene is highly or lowly expressed solely based on the nearby DNA sequence (especially promoters), the model may learn the sequence motifs of TFBSs and localize them. At first glance, this seems infeasible because the model is trained without any information on TFs, but surprisingly, deep learning models can successfully fill the missing links of TF binding between DNA sequences and gene expression.

Basenji [26] is one of the first CNN-based attempts to predict gene expression levels from DNA sequences. It adopts a multitask learning scheme that predicts read coverage representing various modalities in 128-bp genomic bins using genomic sequence information. Gene expression measurements were obtained by cap analysis of gene expression followed by sequencing (CAGE-seq [27]), which allows the precise quantification of gene expression specifically in the vicinity of the transcription start sites. Moreover, the models were trained to predict chromatin features measured by DNase-seq and ChIP-seq. Interestingly, Basenji was able to predict expression quantitative trait loci and disease-associated variants without any prior knowledge of them. Similarly, a multitask CNN-based model named Expecto [28] allows tissue-specific stratification of variants in terms of their transcriptional effect; using the genomic information of other species, especially mice, it boosts the accuracy of gene expression prediction based on deep learning models [29]. Another CNN-based model, Xpresso [30], trained solely based on promoter sequences, can explain about 60% of the variation in human gene expression and clearly reveals general sequencebased features dictating gene expression, such as promoter CpG contents. More recently, modern deep learning architectures, apart from CNNs, have started being applied to genomic tasks. Enformer [31] is one of the forefront applications of transformer architecture [32] for gene expression prediction. Transformer layers allow the model to capture the long-range interactions of distal elements at most 100 kb away from the gene, whereas Basenji only captures a 20-kb window.

Meanwhile, mRNA splicing is a complex biological process that confers the functional diversity and plasticity of cells. Since the process involves sequential binding of RBPs to pre-mRNAs, it is conceivable that the signals governing the efficiency of splicing lie within DNA sequences. Several well-known grammars of mRNA splicing have been firmly established, including the consensus sequence of donor and acceptor sites [33] and the existence of branch points [34]. However, there are numerous combinations of donor and acceptor sites within a gene, only a subset of which is actually spliced out. In other words, the current knowledge on the grammars of mRNA splicing is insufficient to specifically determine the splice site. SpliceAI [35] filled this gap by training a deep learning model with reference genome sequence and exon annotations, and it almost perfectly predicted whether a base is a splice donor or acceptor. One interesting observation drawn from SpliceAI training is that the performance steadily increased as the model was allowed to see larger genomic windows (up to 10 kbp), implying the long-range sequence determinant of mRNA splicing.

3. Learning the correlation between genome and epigenome

The intimate relationship between genomic sequence and epigenomic features has been of great research interest because it enabled the high-throughput measurement of epigenomes. Revealing the dependency of the epigenome on genomic sequences is especially important in clinical applications since it may reveal uncharted pathologic roles of noncoding variants. However, characterizing their quantitative relationship is challenging because there are many complex types of nonlinear interactions across different epigenomic features involving diverse protein machineries. Accordingly, there are great opportunities for the application of deep learning approaches.

DeepSEA [36] is a CNN-based model that predicts allele-specific chromatin profiles from a 1,000-bp sequence window. Models were trained to simultaneously predict 125 DNase I-hypersensitive site profiles and 104 histone mark profiles along with 690 TF-binding profiles. The training data were compiled from ENCODE (Encyclopedia of DNA Elements) [37] and Roadmap Epigenomics projects [38]. The model not only accurately predicted chromatin features based on DNA sequences, it also showed the capability to functionally prioritize single nucleotide variants. The functional impact of each variant was determined by mutating a single base from the original input and propagating the mutated input through the model to obtain the perturbed prediction. The discrepancy between the perturbed and original predictions represents the functional impact of the variant. The power of DeepSEA for functional prioritization of variants implies that deep learning models can capture the general rules of sequence-based regulation of the epigenome while not being overfitted to the training data.

Similarly, Basset [39] is another CNN-based model trained to predict sequence-level grammars by controlling genome accessibility. As mentioned above, its modified version, Basenji [26], accepts a much larger genomic window as an input and predicts more diverse profiles, including ChIP-seq and CAGE profiles. Strikingly, Basenji successfully captured the importance of distal regulatory elements. By computing the gradient of the model prediction with respect to each of the input base positions, one can obtain the effect of small perturbation or mutation in each position in the sequence on the model prediction. Accordingly, a base is considered important when the magnitude of the corresponding effect size, or saliency, is large. Peaks in the saliency map derived from Basenji revealed that the model focuses on distal regulatory elements, especially enhancers, in addition to promoters. This observation underscores the importance of using sufficiently large genomic windows as inputs to allow the model to capture as many unbiased sequence features as possible.

DNA methylation is another major epigenetic feature that regulates gene expression. It is defined as the covalent attachment of a methyl group at the fifth carbon of cytosine bases, and it commonly refers to the methylation of cytosine bases within CG dinucleotides. Since DNA methylation is tightly associated with the regulation of nearby chromatin states, the link between aberrant patterning of DNA methylation and diseases has long been studied. However, systematic characterization of the role of genomic sequences in the regulation of DNA methylation patterns remains challenging.

DeepCpG [40] predicts a single-cell-level DNA methylation state based on the DNA sequence context and nearby methylation states of multiple cells observed by single-cell bisulfite sequencing. It consists of 2 modules, in which the CpG module summarizes the neighboring methylation states of cells with bidirectional RNN and the DNA module summarizes the genomic sequence with the CNN. The outputs of the 2 modules were combined in a joint module. Investigating the filters of the first convolutional layer revealed de novo sequence motifs associated with cell-to-cell methylation variability as well as average methylation levels. While DeepCpG conducts binary classification of DNA methylation states at the single-cell level, MRCNN [41] is a CNN-based model that aims to regress cell population–level methylation.

Finally, deeply learned correlative relationships between genomic and epigenomic features allow us to predict or impute the missing observations for certain epigenomic features from the other present features. The correct imputation of missing epigenomic profiles is important because it may significantly reduce experimental costs and thus facilitate the large-scale measurement of epigenomic features. Avocado [42] is a clever deep factorization-based approach that decomposes the signal value into multiscale genomic position factors, assay factors, and cell-type factors. By jointly learning the 3-factor embeddings and weights for their nonlinear combinations to produce accurate signals, the model could correctly predict epigenomic signals that were not observed during training. As the learned latent embedding enhanced the performance of downstream predictive tasks, such as the prediction of gene expression, promoter-enhancer interaction, and replication timing, Avocado successfully encoded the underlying correlative relationship between epigenetic features.

4. Integrating multiple multiomics modalities at once

The methods discussed so far mainly focus on binary interactions between the 2 omics layers. However, the actual landscape of intracellular multiomics interactions is far more complex because it involves multiple interactions among several omics modalities.

Although we still have a long way to go to achieve clear deep learning model understanding and explaining of the multiomics interaction landscape, there have been several initial attempts to show the potential of deep learning-based multiomics integration. One approach utilized a bottlenecked autoencoder to integrate mRNA expression, DNA methylation, and miRNA expression levels and showed that deeply integrated features perform well as a biomarker for predicting the prognosis of hepatocellular carcinoma [43]. Here, a bottlenecked autoencoder is defined as a model that is trained to emit an output that is identical to the input, where the dimension of the bottleneck layer is far smaller than that of the input. The expression and methylation levels were concatenated into a single vector and fed into the model. The rationale behind this study is that the model should learn a compact representation that captures the interactions involving 2 or more omics layers. On the other hand, MOLI [44] adopts the late integration of mutation, copy number, and gene expression profiles with deep learning. Each omics feature was separately encoded by an individual neural network called encoding subnetwork, and the representations were concatenated and fed into the fully connected layers to predict the drug response of the corresponding sample.

1. Using biological interactions as a guidance

Large-scale networks for functional interactions between proteins and genes include BioGRID [48], STRING [49], HumanNet [50], and Reactome [51]. Because edges in such networks represent paths through which biological information flows, its modeling will be straightforward if we locally propagate information only through the neighboring nodes. To this end, graph neural networks (GNNs) [52] have been widely applied for networkguided modeling of multiomics interactions (Fig. 2A). An early study [53] utilizing graph convolutional filtering on GNNs to capture localized patterns of gene expressions showed promising increase in performance of predicting subtypes of breast cancer samples, and another work using the Kyoto Encyclopedia of Genes and Genomes pathways [54] along with attention-based interpretation revealed subtype-specific aberrations in the biological pathways [55].

Fig. 2.

Strong biological guidance of deep learning models. (A) Graph neural networks (GNNs) are suitable for the modeling of interaction networks. Gene expression values for each sample are assigned to the corresponding nodes in the network to instantiate the network as an input for GNNs. Information of each gene is propagated to its neighbors by graph convolution. After a few iterations of graph convolution and pooling, information of the whole node is aggregated through readout function. Aggregated information is used to predict output values. (B) Knowledge-primed neural network. Nodes in a knowledge-primed neural network directly correspond to genes or proteins, and edges represent the interaction and transcriptional regulation between them. After training the network to predict the observed biological outcome upon certain stimuli, the model is clearly interpretable by edge weights and, thus, the core regulators of the process can be identified. (C) DCell and DrugCell incorporate hierarchical representations of biological knowledge to their network structure called a visible neural network (VNN). While input nodes denote the mutational states of genes, the nodes in hidden layers correspond to the biological concepts. Note that the nodes close to the output layer represent the broader concept. The VNN output, an embedding of the genotype, is subsequently used for phenotype prediction in DCell and drug response prediction in DrugCell.

On the other hand, networks can also be indirectly used as a regularizer for the latent features learned by the model. The multiview factorization autoencoder [56] adopts a general scheme called graph Laplacian regularizer [57] to incorporate network information in the training of autoencoders. The graph Laplacian regularizer measures the overall discrepancy of the learned features between neighboring nodes (e.g., genes) in terms of Euclidean distance, and its value is added to the loss term. Thus, a pair of genes connected by an edge in the network is forced to have similar latent representations.

Knowledge-primed neural networks (KPNNs) [58] are examples of biologically transparent and interpretable model architectures that use biological networks (Fig. 2B). Each node in the KPNN represents a gene or a protein, and the edges between nodes denote the known interactions between them. Given that the model achieves the desired performance for modeling a biological phenomenon, the model itself can be interpreted as a quantitative hypothesis for the phenomenon. The weight of a connection between 2neurons in the network represents the importance of the regulatory connection between them, thus allowing the prioritization of genes or proteins in the biological process. Notably, a KPNN trained for single-cell experimental results of T-ell receptor (TCR) stimulation showed prediction accuracy comparable to that of a generic deep neural network. Understandably, the topology of KPNN resembles that of a biological network in terms of 4 properties: (1) shortcuts between layers, (2) scale-freeness, (3) modularity, and (4) restricted reachability of hidden neurons to input neurons. As a result, the key mediators of TCR signaling could be identified by analyzing the activation of the edge weights connected to each gene upon stimulation compared to control inputs. Since the model could be clearly explained in terms of genes and proteins while accurately predicting the system-level outcome, it would be a good starting point for system-level characterization and simulation of a cell. Recent work on estimating transcriptomic age using KPNN supports this as known associations between the core pathways and aging were captured, but virtual knockdown of genes accurately recapitulated the effects on aging that were experimentally validated.

2. Using the biological knowledge hierarchy as guidance

Many biomedical concepts can be represented hierarchically. Therefore, these concepts are naturally organized as tree-like data structures. For example, the concept “cell cycle” encompasses more specific concepts such as “mitotic cell cycle” and “meiotic cell cycle,” and “mitotic cell cycle” includes concepts like “G1/S transition” and “DNA replication.” The most prominent example of the hierarchical compilation of biomedical concepts is the gene ontology (GO) terms [59], where the concepts are organized into 3 trees rooted in 3 categories: biological processes, molecular functions, and cellular components. MeSH (medical subject headings) [60] is another exemplary database that uses a conceptual hierarchy to organize biomedical terms.

DCell [61] was the first approach to constructing a deep learning architecture that inherits the topology of the GO hierarchy (Fig. 2C). Its model, the visible neural network (VNN), consists of thousands of subsystems corresponding to each term in the GO hierarchy. Multiple neurons constitute a subsystem, and the connections between neurons are only allowed for the 2 neurons belonging to a pair of subsystems with hierarchical relationships. Since the leaf nodes of the GO hierarchy are genes, each node at the input layer of the VNN corresponds to a gene and its state represents the mutation status of the corresponding gene. Therefore, the input as a vector collectively denotes the cell’s genotype. The model was trained to predict the growth phenotype of yeast cells based on their genotypes. DCell could be used for mechanistic interpretation of perturbed biological processes due to mutations. Because the information flow originating from the mutations at the input layer is visible and interpretable throughout the network, a mechanistic explanation of genotype-phenotype association and their simulation is possible.

DrugCell [62] extends DCell to predict the drug responses of human cancer cells based on their genotype. The genotypes are embedded with VNN as in DCell, and the latent representation of drugs is computed using generic deep neural networks. The resulting 2 representations, genotypes and drugs, are concatenated and fed to a fully connected layer with nonlinearities to predict the cell’s response to the drug. The authors expressed this procedure as an in silico treatment of a cell with a drug. Similar to DCell, DrugCell was able to prioritize subsystems according to their predictive strength for the drug response of the cell. In other words, DrugCell revealed specific mechanisms that convey the effects of drug treatment. Since the model could identify the most crucial intracellular subsystem for each drug, strikingly, the principled design of synergistic drug combinations that maximizes the drug efficacy was demonstrated as possible.

3. Other approaches for biological knowledge-guided deep learning

Although much biological knowledge can be structured in the form of networks as discussed so far, much knowledge cannot be simply organized as networks. Therefore, it is becoming an active area for researchers to apply their own ideas to transform rich information into an organized form that can explicitly guide models.

As discussed in the previous section, the binding motif preference of TFs can be learned de novo through supervised learning. However, it would be desirable to utilize the precompiled TF-binding motif preference data such as JASPAR [63], TRANSFAC [64], or GTRD [65]. This idea was initially examined by Ploenzke and Irizarry [66], where the weights of the convolutional filters were initialized with PWMs of JASPAR motifs. Because not all TFs are important for a specific prediction task, the filters are discarded if their contribution to the model performance is negligible. New randomly initialized filters were added to the model to allow de novo motif learning. On the other hand, because the PWMs allow us to compute the putative binding sites and binding affinities of each subsequence, we can incorporate the precomputed TFBSs into the deep learning model. Kang et al. [67] realized this idea by adding the binding site information and the expression levels of corresponding TFs as extra information for the gene expression prediction task. Interestingly, guiding the model with putative TFBSs resulted in slightly improved model performance. Combined with attention weights, the model enabled the mechanistic interpretation of the interaction between TF binding and methylation level on gene expression.

Analysis of the 3-dimentional (3D) organization of chromatin has recently been of great interest. Accordingly, there are still many unexplored possibilities when the traditional 1-dimensional view of genomic sequences is modified to a 3D view. In particular, the prediction of gene expression levels may benefit from such new modeling of genomes because it allows the efficient incorporation of distal regulatory relationships. The 3D chromatin contacts are usually measured by Hi-C and are represented as a pairwise contact map representing the likelihood of contact between genomic fragments. Thus, utilizing the contact map in deep learning may increase its performance. ChromeGCN [68] was the first model to use 3D chromatin contacts to predict chromatin profiles. Using the 500,000 most likely Hi-C contacts between 1,000-bp fragments per chromosome, ChromeGCN views each chromosome as a graph in which the nodes and edges represent genomic fragments and contacts between them, respectively. As a result, the study showed that exchanging the information of genomic sequences through graph edges based on 3D contact using a graph convolutional network outperformed the state-of-the-art model that uses only a local sequence context.