Deep Learning and Explainable AI: New Pathways to Genetic Insights

Input interpretability aims to reveal key patterns within the input features or to assess the importance of each feature for the prediction of the model. Unlike traditional k-mer-based feature engineering methods that directly reveal feature importance, deep neural networks employ multi-layer nonlinear transformations for automated feature extraction, resulting in opaque decision processes. To address this, there exist many input interpretability methods, which we further categorize into three classes: convolutional kernel visualization, gradient-based methods, and perturbation methods. Figure 2 illustrates the underlying principles of each method.

Convolutional kernel visualization. With the advancement of deep learning technology, convolutional neural networks (CNNs) have been widely applied in genomic research, such as enhancer prediction [24], transcription factor binding site prediction [41], and chromatin accessibility prediction [14]. CNNs can automatically extract features from DNA sequences. Specifically, CNNs slide multiple filters across the DNA sequence in the convolutional layers, obtaining activation values for each position relative to each filter. These activation values are then used as features and passed to fully connected layers for prediction. As shown in Figure 2(a), the convolutional kernel weights in the first layer are first converted to a position frequency matrix (PFM), followed by log-scaling to produce a standard position weight matrix (PWM) [1, 28], which captures a short sequence motif. Therefore, visualizing these filters can help us understand which sequence patterns the model has extracted that are crucial for prediction. A potential issue here is the constraints in weight learning: unconstrained or inappropriate constraints during weight learning may lead to scaling problems, causing the importance of certain sequences to be disproportionately emphasized or underestimated. In most studies, a commonly used visualization strategy involves computing the activation values of multiple sequences for each trained filter, indicating the degree of match between each sequence and the filter. The sequences with the highest activation values are then selected to generate a corresponding PWM, representing the motif learned by the filter. Then we can obtain motifs learned by multiple filters and match them with existing motif databases to verify whether the model has captured biologically relevant and important features. For example, DeepEnhancer [24] uses a CNN to predict enhancers from DNA sequences. Its first convolutional layer applies 128 filters of length 8 to extract features. By visualizing these filters, the model identifies motifs and matches them with the known transcription factor binding site database JASPAR [23], successfully recognizing biologically relevant motifs. Furthermore, Basset [14] recovered a large number of known DNA-binding protein motifs by analyzing the parameters of the 300 convolutional filters in the convolutional layer.

Gradient-based methods. Gradient-based importance analysis methods achieve global quantitative evaluation of feature importance by computing the partial derivatives of model outputs relative to input features. As shown in Figure 2(c), this approach utilizes backpropagation to compute the gradients of the loss function with respect to sequence features, where gradient magnitudes directly measure the contribution of individual nucleotide in the DNA sequence. For instance, in genomic analyses, positive/negative gradient peaks correspond to enhancer and silencer regions, respectively [15]. However, gradient-based methods are susceptible to vanishing gradients or gradient saturation. In deep neural networks or when using saturating activation functions, excessively small gradients may introduce bias in feature importance estimation. This limitation can be mitigated by integrated gradients, which correct importance scores through path integration [33]. Moreover, an alternative approach, DeepLIFT [31], provides local explanations by comparing the prediction differences between test instances and reference sequences (e.g., background nucleotide frequencies), thereby avoiding the saturation effects of gradient methods. DeepLIFT has been employed to identify critical nucleotides in splice sites [44]. The application of DeepLIFT algorithm successfully identified expression-predictive motifs (EPMs) in both 5’UTR and 3’UTR regions, which exhibit distinct sequence-specific positional preferences [29]. Notably, the attribution rules of DeePLIFT match Shapley values, providing a robust feature importance framework [21]. Both integrated gradients and DeepLIFT require a predefined reference baseline, and the choice of this baseline can significantly influence the accuracy of attribution results. However, there is currently no consensus or established guideline on how to select an appropriate reference.

Perturbation-based methods. Perturbation-based methods infer the importance of features by manipulating the input features and observing the changes in the output of the model. Perturbation-based interpretability methods were first introduced in the field of computer vision, where specific regions of an image are masked to observe changes in predictions, thereby assessing the importance of those regions [38]. This intuitive approach aligns well with human reasoning and has since been widely adopted in other domains. Similarly to image perturbation, one can consider altering certain specific segments in biological gene sequences to determine the importance of the features corresponding to these segments. For example, we can perturb parts of the sequence and observe changes in the output, as seen in Figure 2(b). The discrepancy between the predictions for these new sequences and the original sequence is quantified as the attribution score, serving as a critical metric for evaluating the functional significance of individual nucleotides in the context of the trained model. It is also feasible to mutate a specific nucleotide within the sequence into each of the other three nucleotides and observe the resultant changes in the output. By repeating this operation for all nucleotides, a matrix of dimensions $4\times L$ can be obtained, which is commonly referred to as the attribution map [32]. iMAP [20] is used to knock out genes for the discovery of therapeutic targets for diseases. The perturbation model based on $\beta$-VAES [4] conducts perturbation analysis on single-cell RNA datasets to predict gene expression changes during gene knockout, toxin response, and embryonic development. A perturbation map [6] was developed to provide a scalable approach for evaluating how specific genetic alterations impact the local, proximal, and distal tumor microenvironment (TME) states. However, many neural networks are trained in a way that resists Dropout [12]. As a result, multiple neurons may be associated with the importance of the same feature, which can lead to the failure of perturbation methods.

For input interpretability, while weight matrices corresponding to convolutional kernels can reflect the importance of sequences, the lack of proper constraints during their learning process may distort the perceived importance of these sequences. Gradient-based methods can evaluate feature importance by computing the partial derivatives of model output with respect to input features, while they may also face the potential issue of gradient vanishing. Perturbation-based methods assess the importance of altered features by introducing minor perturbations to the input and observing the resulting changes in model outputs, yet they may lead to an underestimation of feature importance due to neuronal redundancy.

Model interpretability aims to facilitate human understanding of internal mechanisms and decision-making processes through architectural design. For example, traditional linear models are inherently interpretable, as we can directly understand the importance of each feature through its assigned weight. In contrast, deep neural networks typically consist of multiple hidden layers and numerous parameters, making it difficult to comprehend their decision-making processes solely based on their architecture or parameters. Therefore, many researchers are currently exploring model-based interpretability techniques to enhance the credibility of the model and uncover new biological insights. These methods can be mainly categorized into two types: attention mechanisms and transparent models based on biological prior knowledge. The principles of the various methods are illustrated in Figure 3.

Attention mechanisms. The attention mechanism is a widely adopted technique in neural networks. A representative example is the Transformer model [37], which leverages self-attention to capture dependencies between tokens within a sequence. The attention weights inherently provide a form of interpretability [27], offering insight into which input tokens the model focuses on when making predictions and revealing interactions between them. The visualization results are shown in Figure 3(a). EpiBERT [13] introduces a masking technique for genomic loci and enhancer regions, extracts query and key matrices from each attention layer, averages the resulting weights across all layers and attention heads, and visualizes these weights to indicate the degree of correlation between different positions. Enformer [2] extracts the query row at the transcription start site (TSS), where the keys denote different spatial positions and the attention values reflect the focus of the model on these positions during the prediction of TSS. In addition, C.Origami [34] exploits attention scores from Transformer models to identify cis-regulatory elements that play critical roles in 3D chromatin architecture. However, in DNA sequences, sequences at multiple positions may be associated with the same biological feature, which is mathematically termed multicollinearity, and this characteristic may lead to instability in the weights of attention mechanism matrices.

Transparent models based on biological prior knowledge. Unlike traditional neural networks, where neurons in the hidden layer typically lack clear biological meaning, transparent neural network models are explicitly designed so that their hidden nodes correspond directly to specific biological units at a fine-grained level, thereby significantly enhancing interpretability [22, 9, 36]. To construct models with intrinsically interpretable units, it is necessary to incorporate prior knowledge into the network architecture design. For instance, by mapping bottom-layer nodes to molecular entities (e.g., genes or proteins) and deeper-layer nodes to functional modules (e.g., metabolic pathways or organelles), a hierarchical biological representation system can be established, as illustrated in Figure 3(b). Deep neural networks build higher-order biological representations through progressive integration of lower-level features, for example, the second layer may capture cooperative interactions between transcription factor binding motifs, while deeper layers can map to complex systems such as complete biological pathways [28]. As a breakthrough in biological applications, DCell [22] pioneered a modeling approach that combines both neural network computational power and model transparency. The input layer of DCell directly corresponds to gene nodes, while its second layer constructs functional group nodes based on the hierarchical structure of Gene Ontology, establishing an association model between genotype and yeast growth rate. The GenNet [36] model utilizes genetic variants (SNPs) as input layer nodes and maps them to gene nodes in the second layer through NCBI RefSeq gene annotations, thereby constructing predictive association models between genotypes and complex phenotypes (e.g., schizophrenia, hair color). Through transparent variable construction based on multi-omics data, P-NET [9] has identified MDM4 as a biomarker in metastatic prostate cancer. Although transparent models offer significant interpretability advantages, their construction relies on domain-specific prior knowledge, which to some extent limits their applicability. Furthermore, the technical approaches used to ensure model transparency may adversely affect predictive performance.

For model interpretability, the attention mechanism reveals the importance of input tokens and reflects the correlations between different tokens. Multicollinearity among input tokens can lead to numerical instability in attention weights. Transparent neural networks can fully incorporate biologically relevant prior knowledge to reveal specific biological functional relationships, yet they may lack universality and cause models to overlook not clearly defined relationships.

Analysis reveals three fundamental limitations in input interpretability approaches for deep neural networks: First, unconstrained weight optimization in convolutional layers provably causes scaling instability, mathematically demonstrating disproportionate feature amplification. Second, the inherent neuronal redundancy from dropout mechanisms obscures node-level importance attribution. Third, the null gradient of ReLU for negative inputs, combined with the multiplicative accumulation of small weights through successive layers, results in progressive gradient decay during backpropagation.

The unconstrained learning of weights. CNNs for DNA sequence prediction extract local features by sliding multiple convolutional kernels across the input sequence. The learned kernels can ultimately be interpreted as PWM, representing the key motifs identified by the model from the input sequences. However, it should be noted that the learning of weights should be regulated by constraints. Otherwise, it is highly likely to lead to scaling issues.

Suppose we study a classification problem where a DNA sequence is input to determine whether each position is related to a certain structure. The input sequences are typically represented as a $L\times 4$ matrix, where $s_{ij}$ represents the entry at position (i, j) and the convolutional kernel is a matrix of size $k\times 4$. To facilitate the analysis of the problem, the network structure is simplified. Specifically, the fully connected layers are abstracted into a single layer, and the binary cross-entropy loss function is adopted as the loss function. The sigmoid function $L=\sigma(z)$ is used as the activation function for the last layer of the neural network. The weights of the convolutional kernel are updated through gradient descent. Then in each learning iteration, the weights will be updated in the following manner:

To further analyze the scaling problem of the weights, the partial derivatives will be expanded in accordance with the chain rule. Taking $w_{11}$, which is located at the first row and the first column of the convolutional kernel as an example:

The sequence feature values $y_{1},y_{2},...,y_{L}$ are obtained by convolutional operation (using $y_{1}$ as an example):

The output $z$ is obtained by multiplying the fully-connected weight matrix with the feature sequence.

Furthermore, the update formula for the weights of the convolutional kernel can be specifically written.

Among them, $s_{ij}$ $\in$ $\{0,1\}$. If the true label is 1, then $\sigma(z)-y<0$. If the subsequent part is greater than 0, it will cause the weight to increase. At this time, if there is no constraint, there will be a risk that the weight will continue to increase. Eventually, a certain weight will become too large.

Then it is common practice to transform the PWM into a PFM[28] through the softmax operation, and then obtain the standard PWM matrix via a logarithmic transformation.

The softmax operation on the convolutional kernel is shown in the following equation (still taking $w_{11}$ as an example):

It can be seen that if $w_{11}$ is too large while the other weight values are relatively small, it leads to $w^{\prime}_{11}\approx\frac{e^{w_{11}}}{e^{w_{11}}+0+0+0}=1$.

And if $w_{11}$ is small and the other weight is relatively large, it will lead to $w^{\prime}_{11}\approx\frac{0}{0+e^{w_{12}}+e^{w_{13}}+e^{w_{14}}}=0$.

This will lead to the irrationality of the weights. It will cause the functions of certain bases in the sequence to be overemphasized, while the significance of some other bases is neglected.

For the purpose of ensuring the positive and negative nature of the values, when performing the logarithmic transformation, it is usually carried out by $w_{ij_{log}}=log(w^{\prime}_{ij}+\epsilon)$, which $\epsilon$ is always determined according to the condition to satisfy the stability of the values. Due to $\epsilon$ being always small, when $w^{\prime}_{ij}$ approaches 0, the value after logarithmic transformation will tend to negative infinity. When $w^{\prime}_{ij}$ approaches 1, the value after the logarithmic transformation will tend to 0. Both scenarios will lead to the distortion of $w_{ijlog}$.

This may lead to the effect of the bases in certain parts of the input DNA sequence being overly amplified, while the roles of the bases at other positions are ignored. This process theoretically analyzes why unconstrained learning can lead to scaling issues, thus rendering the interpretation of the DNA sequence lacking in reliability.

The redundancy of neurons. Perturbation-based methods can evaluate the contribution of individual neurons by selectively removing them and observing the resulting changes in model predictions. However, deep neural networks are usually trained in a way that resists Dropout. The mechanism of Dropout is to deactivate some neurons during each training session. This ensures that certain specific features in the input DNA sequence do not rely on a fixed neuron, but can be captured by multiple neurons. More specifically, it is $r_{j}^{l}\sim Bernouli(p)$, which means that the $j$-th neuron in the $l$-th layer has the probability of $p$ being retained during training process.

Suppose that there are two different neurons, $N_{i}$ and $N_{j}$, in the first convolutional layer, $y_{i}$ and $y_{j}$ are the final contributions of neuron $i$ and neuron $j$ in the entire neural network.

In standard CNNs, adjacent layer neurons are tightly coupled through weight updates, enabling collaborative feature learning. Dropout randomly disconnects the interneuron pathways during training, forcing independent feature acquisition. This causes redundant node formation: multiple subnetworks repeatedly learn identical DNA sequence features, meaning minimal output changes when removing specific neurons do not necessarily indicate low predictive importance of these nodes.

Since a neural network contains numerous neurons, as well as parameters such as weights and biases, it can be regarded as a high-dimensional manifold $\mathcal{M}$. The number of neurons in each neural network as a whole, the initial weights and other factors together constitute an abstract parameter $\theta$. From this perspective, the neural network training process is to select a minimum point of the loss function $L(\theta)$ on $\mathcal{M}$. Suppose that the overall parameter settings under the initial conditions are $\theta$. At this time, for any input $x$, the output of the model can be expressed as $f_{\theta}(x)$. $N_{i}$ and $N_{j}$ can be regarded as two distinct minor parameters within the parameter space $\theta$. Removing such a neuron is equivalent to perturbing in a certain parameter direction. If these two neurons are redundant with respect to each other, that is, there exist parameter directions $v_{i},v_{j}\in\mathbb{R^{d}}$ such that for any input $x$:

This equation indicates that the output results after small perturbations in two directions are approximately equal, where $O(\epsilon^{2})$ denotes an extremely small error related to $\epsilon$.

And these two parameter directions can be regarded as symmetric, and the corresponding covariant derivatives commute.

Then the components of the Riemann curvature tensor $\mathcal{R}$ of the parameter space in the directions of $v_{1}$ and $v_{2}$ are:

It represents that the manifold is locally flat in the two parameter directions represented by the removal of redundant nodes. Assume that $N_{i}$ is removed. Then it can be regarded as a perturbation in the direction of $v_{i}$ in the parameter space, with $\theta^{\prime}=\theta-\alpha v_{i}$. Subsequently, the variations in the output can be calculated:

Due to the large number and complexity of parameters in the entire neural network, the removal of a single neuron can be regarded as a minor change in a certain direction, and $\alpha$ is an infinitesimal quantity. At the same time, since the manifold $\mathbb{R}$ is locally flat in the subspace spanned by $v_{1}$ and $v_{2}$, $\nabla v_{1}f_{\theta}(x)$ and $\nabla v_{1}^{2}f_{\theta}(x)$ approach 0.

In summary, we obtain that $\Delta f$ approaches 0. That is, when redundancy occurs in neurons due to the Dropout mechanism, the output of the model does not change significantly after removing a single neuron, which makes it hard to access the importance of the node.

Vanishing gradient in DNA sequences. Vanishing gradient is a common issue in gradient-based model interpretation methods, and is closely tied to the choice of activation function. For instance, the ReLU activation function outputs zero for negative inputs while propagating positive inputs unchanged. However, its constant unity gradient for positive inputs may still lead to vanishing gradients during backpropagation in deep networks. Consider a single-layer neural network with the output given by:

Here, $W$ represents the weights matrix, $X$ is the input, and $b$ is the bias term. The normal input and output relationship can be expressed as $Y=\text{ReLU}(W\cdot X+b)=W\cdot X+b$. The gradient of the output $Y$ with respect to the input $X$ can be calculated as:

According to the derivative definition of ReLU, if $W\cdot X+b>0$, then $\text{ReLU}^{\prime}(W\cdot X+b)=1$, and the gradient is $W$. Conversely, if $W\cdot X+b\leq 0$, then $\text{ReLU}^{\prime}(W\cdot X+b)=0$, and the gradient is 0. Extending this to a multi-layer neural network, assume output of each layer passes through a ReLU activation function, with the output of $l$-th layer given by:

Here, $Z^{(0)}=X$ is the input, and $Z^{(L)}=Y$ is the output. To calculate the gradient of the output $Y$ with respect to the input $X$, the chain rule must be applied:

The gradient of each layer is:

If the input to a layer $W^{(l)}\cdot Z^{(l-1)}+b^{(l)}\leq 0$, then $\text{ReLU}^{\prime}(W^{(l)}\cdot Z^{(l-1)}+b^{(l)})=0$, and thus, the gradient for that layer is 0. This means that the gradient at this layer vanishes and cannot continue to propagate backward. Additionally, the gradient is also influenced by the weight matrices, as shown by:

When $W^{(1)},W^{(L-1)},\ldots,W^{(L)}$ are very small, the value of the gradient will also be very small, leading to vanishing gradient.

This section systematically demonstrates the inherent limitations of model interpretability approaches in DNA sequence analysis. First, attention mechanisms intrinsically exhibit numerical instability, which can be mathematically proven to induce matrix ill-conditioning under multicollinearity. Second, biologically-constrained transparent models may yield higher loss values compared to their unconstrained counterparts.

Instability of in the estimation of the weights. The computation of attention weights typically involves two vectors, namely the query vector and the key vector. The attention score is calculated by taking the dot product of the query matrix ${Q}=[q_{1},q_{2},...,q_{n}]$ and the key matrix ${K}=[k_{1},k_{2},...,k_{n}]$:

Since the input features often exhibit multicollinearity, some of the vectors in the key matrix are linearly correlated. That is, there exists an index $i$ such that $k_{i}=\sum_{j\neq i}\alpha_{j}k_{j}+\beta$.

Consider the covariance matrix $C$ of $K$, and $C=K^{T}K$, due to the multicollinearity, the rank of the matrix $K$ will satisfy $r(K)=r<n$, then we have:

The matrix $C$ also satisfies the condition of partial linear correlation. Therefore, $|C|=0$. Since the matrix $C$ is a positive semi-definite matrix, all the eigenvalues of $C$ are greater than or equal to 0. At this time, there will be:

So there must exist $\lambda_{i}$ which satisfies that $\lambda_{i}=0$, i.e., $\lambda_{min}=0$.

To demonstrate that multicollinearity leads to numerical instability, we first present the following lemma and theorem.

This means that multicollinearity will lead to an excessively large condition number of the matrix, which in turn poses a potential risk of numerical instability [40].

Constructing transparent models with prior knowledge. Transparent models constructed based on prior knowledge have attracted widespread attention due to their computational efficiency and interpretability. However, transparent models also have their limitations. Biological constraints within transparent models may diminish their predictive performance. For example, Hard-coding TF-TF interactions may underperform unconstrained models. (although the use of milder regularization techniques can mitigate this decrease to some extent). The following serves as an illustration of this. Suppose that the DNA samples input into two models are identical, with all other confounding factors excluded. The loss function for the hard-coded transparent model is defined as follows:

There, $N$ denotes the number of input samples, $m_{i}$ is the true label of the $i$-th sample, and $p_{i}^{\text{hard}}$ represents the predicted probability of the $i$-th sample by the hard-coded model. This prediction is subject to the constraint $p_{i}^{\text{hard}}=f(W_{\text{prior}},W_{\text{adapt}},m_{i})$. The regularization coefficient $\lambda$ is associated with the prior constraint term, while $\gamma$ corresponds to the regularization coefficient for the adjustable parameters in the hard-coded model, which usually has a smaller value. $W_{\text{prior}}$ signifies the fixed weights in the hard-coded model, $W_{\text{ref\_prior}}$ is the preset baseline weight reference value, and $W_{\text{adapt}}$ is the weight in the hard-coded model that can be freely adjusted. The loss function of the hard-coded model is primarily determined by the first term, which is the cross-entropy loss. The loss function for the unconstrained model is defined as follows:

In comparison, the predicted probability $p_{i}^{\text{free}}=f(W_{\text{free}},m_{i})$ of the unconstrained model is determined by the free parameters $W_{\text{free}}$. Here, $\beta$ is the regularization coefficient for the unconstrained model (with a smaller value). Similarly, the loss function of the unconstrained model is primarily determined by the first term, which is the cross-entropy loss. Due to the limitations imposed by hard coding, it holds that $p_{i}^{\text{hard}}\leq p_{i}^{\text{free}}$. To demonstrate that the performance of the hard-coded model may potentially decrease compared to the unconstrained model, it is only necessary to prove that the difference between the loss functions of the hard-coded model and the unconstrained model is not less than zero.

Where $\sum_{i=1}^{N}m_{i}\left(\log\frac{p_{i}^{\text{free}}}{p_{i}^{\text{hard}}}\right)\geq 0$, $\lambda\|W_{\text{prior}}-W_{\text{ref\_prior}}\|_{2}\geq 0$, and $\gamma\|W_{\text{adapt}}\|_{2}-\beta\|W_{\text{free}}\|_{2}\approx 0$. From this, it can be concluded that $L_{\text{hard}}-L_{\text{free}}\geq 0$, which proves that the hard-coded model may perform worse than the unconstrained model.