Extractable DNA sequence motif extraction and visualization method based on deep learning

Abstract

The invention belongs to the technical field of deep learning, and particularly relates to an interpretable DNA sequence motif extraction and visualization method based on deep learning, and a motif extraction method based on combination of baseline input constructed based on a data set AGCT background and DeepLIFT attribution. A target motif length is parameterized into a normalized cut-off frequency, and a length controllable positioning frequency domain noise suppression scheme is realized through adjustable low-pass filtering; a screened motif set is obtained according to a dual screening standard formed by Fisher single-side inspection, BH-FDR correction and an effect quantity threshold value in a combined mode; an integrated link from site-by-site attribution to frequency domain noise suppression and length adjustability to statistics double-threshold robust screening to database comparison verification is established while the capability of capturing long-short motifs with different scales is ensured; the motif extraction not only can reflect the internal discrimination basis of the model, but also can maintain repeatable, comparable and traceable biological reliability in different data fields and application scenes.

Description

Technical Field

[0001] This invention belongs to the field of deep learning technology, and in particular relates to a method for extracting and visualizing interpretable DNA sequence motifs based on deep learning. Background Technology

[0002] Sequence motifs are core clues in epigenetic and transcriptional regulatory mechanisms, directly related to key processes such as promoter / enhancer activity, transcription factor binding, and chromatin accessibility. With the rapid accumulation of high-throughput omics (such as ChIP-seq, ATAC-seq, and single-cell sequencing) and reference databases, the demand from research and industry for "stable and reproducible extraction of biologically reliable motifs in long-sequence scenarios" continues to rise. Especially in tasks involving essential genes, longer or complex regulatory fragments are often significantly associated with functional necessity, requiring algorithms that can suppress noise and provide tunable control over length scales to meet the requirements of generalization and reproducibility across datasets and platforms. Existing work has gradually formed a general paradigm of "attribution—signal processing—candidate localization—statistical screening—database comparison," using the site-by-site contribution of deep models as the signal source for motif discovery, and then improving the interpretability and comparability of results through statistical and knowledge base validation. However, bottlenecks still exist in unified modeling across long and short scales, robustness of boundary localization, and transferability across scenarios.

[0003] In the area of ​​"motif discovery based on importance scores from deep models," a mature technical framework has been established in academia. This framework involves first training a sequence model capable of characterizing regulatory signals, then performing base-by-base attribution on the input sequence, using recurring, high-contribution local fragments from the attribution curves to abstract comparable motif representations (PWM / PFM / CWM), and finally performing similarity matching with authoritative databases to complete biological annotation. Key steps in this approach include: calculating the contribution of input features to the model output using integrated gradients or similar methods; identifying and clustering high-contribution fragments (seqlets) using representative algorithms for "motif discovery from importance scores" (such as TF-MoDISco), summarizing them into visualized motifs; and comparing and statistically evaluating candidate motifs with standard libraries such as HOCOMOCO using TomTom on the results side, thereby grounding the "discrimination criteria learned by the model" in "interpretable regulatory patterns." This process has been validated in numerous genome modeling and transcription factor studies and has become the mainstream implementation path closest to the goals of this invention.

[0004] The aforementioned approximation schemes are highly effective in identifying peak regions of short to medium length, but their performance and efficiency are still constrained by several factors: attribution operators and hyperparameters such as thresholds / windows have a significant impact on seqlet detection and clustering stability, and cross-dataset transfer often requires readjustment; in scenarios with long sequences, weak signals, or complex / nested structures, local segment-driven clustering and alignment are more susceptible to noise and distribution differences, thus increasing computational and manual verification costs; furthermore, while database comparisons provide standardized similarity and significance measures, different database versions and scoring criteria can lead to differences in comparability. These limitations are common pain points faced by existing approximation schemes in practical applications, and they constitute the direct motivation and baseline for this invention to further improve the robustness, scale adjustability, and cross-scenario transferability of long sequences.

[0005] Therefore, this invention aims to provide a method for interpretable DNA sequence motif extraction and visualization based on deep learning; a motif extraction method that combines baseline input constructed based on the AGCT dataset background with DeepLIFT attribution; a frequency domain noise reduction scheme that parameterizes the target motif length as a normalized cutoff frequency and achieves length-controllable localization through adjustable low-pass filtering; and a dual screening criterion composed of Fisher's one-sided test, BH-FDR correction, and effect size thresholds (OR and 95% CI). Summary of the Invention

[0006] The purpose of this invention is to provide a deep learning-based method for interpretable DNA sequence motif extraction and visualization, in order to solve the technical problems existing in the prior art.

[0007] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows: A deep learning-based method for interpretable DNA sequence motif extraction and visualization includes the following steps: S1: Obtain a publicly available human gene-dependent dataset containing DNA sequences, gene necessity tags, and continuous fitness scores, and preprocess the dataset to obtain a training subset and an evaluation subset; S2: Gene sequence encoding is performed on each base in the DNA sequences of the training subset and the evaluation subset to obtain an initial matrix. The initial matrix is ​​then converted into an embedding tensor E to align the encoding dimension with the number of channels in subsequent network layers. S3: Based on the embedding tensor E, a sequentially connected feature extraction network is constructed to extract features, resulting in the CNN feature matrix Hcnn, the attention feature matrix H, and the bidirectional hidden state Hbi at each base position; S4: Globally pool the bidirectional hidden state Hbi output by the BiLSTM network to obtain the sequence-level representation z. Construct a classification or regression output layer through a fully connected layer. Train the deep learning model to complete the gene necessity prediction with the gene necessity label or fitness mean as the target. S5: Motif extraction and screening: Perform site-by-site attribution calculation, frequency domain noise reduction and candidate localization, and dual-threshold joint screening to obtain the screened motif set; S6: Compare and verify the selected motif set with a standard database to generate a complete motif list containing motif sequences, statistical test results, and database comparison evidence; draw a color fluctuation curve, mark the peak and valley nodes, and annotate the motif regions after comparison and verification with orange to achieve visual output of the motifs.

[0008] Preferably, the specific process of preprocessing the dataset in step S1 to obtain the training subset and the evaluation subset is as follows: S11: Data format standardization: unify the base encoding format of DNA sequences, retain only the four valid bases A, T, C, and G, mark unknown bases as N, and associate each sequence with the corresponding gene ID, cell line origin, tissue type, fitness score, and necessity label to form a structured data table; S12: The IsolationForest algorithm is used for anomaly detection. The algorithm contamination rate is set to 0.01. The algorithm models the length characteristics of all DNA sequences, automatically identifies and removes abnormally long sequences that exceed the normal length distribution range, so that the DNA sequence length is concentrated within 10000bp. S13: Select core essential genes using majority voting rules: If a gene is determined to be an essential gene in more than half of the cell lines, then the gene is included in the core essential gene set. S14: Based on the processed sequences, the set of core essential genes, and the mean fitness value, the training subset and the evaluation subset are randomly divided according to a specified ratio.

[0009] Preferably, the specific process of step S2 is as follows: S21: Take any DNA sequence s of length L from the training subset or evaluation subset, in the form s=( s 1 , s 2 ,...,s L ); S22: For each base in the sequence s t Perform one-hot encoding to generate a vector of dimension 4. xt This corresponds to the binary identifier of the four effective bases; S23: Generate the one-hot vectors corresponding to all bases in the sequence. x 1, x 2,..., x L The sequences are vertically concatenated in sequence to form an initial matrix X with dimensions of "sequence length × base coding dimension", expressed as: S24: Transform the initial matrix X, with a fixed dimension of 4, into an embedding tensor that matches the number of channels in subsequent network layers. E .

[0010] Preferably, the specific process of step S24 is as follows: Introducing the linear projection formula: ; in, The projection weight matrix is... This represents the number of feature channels for subsequent model layers; X is the bias vector used to compensate for the offset of the linear projection; X is the initial matrix. E This is the embedded tensor of the output.

[0011] S241: For the initial matrix X The Middle t row vector x t , with weight moments Perform matrix multiplication to obtain the intermediate vector. x t ; S242: Intermediate vector x t Superimposed bias This yields the embedding vector of the t-th base. e t = x t ; S243: Embedding vector of all bases e 1, e 2,..., e L By vertically concatenating the sequences, we obtain the embedded tensor E=[ e 1; e 2; ...; e L ] .

[0012] Preferably, the specific process of step S3 is as follows: S31: One-dimensional convolution extracts local n-qram patterns in the time dimension using multi-scale kernels while maintaining the same length. The convolution output is the CNN feature matrix. ; S32: By stacking Transformer encoders consisting of self-attention and feedforward networks, long-range dependencies are aggregated, and an attention feature matrix is ​​output. ; S33: Use a bidirectional LSTM to model the time dependency in both forward and backward directions to obtain the bidirectional hidden state at each base position. .

[0013] Preferably, the specific process of performing global pooling on the bidirectional hidden state Hbi output by the BiLSTM network in step S4 to obtain the sequence-level representation z is as follows: S41: The cells of BiLSTM at each time step Receive vectors from the attention block And the previous moment Calculate the forward and reverse chains separately. And spliced ​​together at each position as follows: ; Let the input of the attention layer be... , No. The query, key, and value of each head are: in ; Scaling Dot Product Attention: ; S42: By concatenating all the heads and linearly mapping them, we obtain: ; S43: Perform residual summation and layer normalization, then pass through a feedforward network to form the output of an encoder layer; stacking The results are then fed into a BiLSTM layer. S44: When the input is a DNA sequence, the tensor shape of the data stream changes gradually along the pipeline; when , 256. Number of heads LSTM hidden state dimension : The original onehot was Embedded / First Layer After convolution, we get One-dimensional convolution maintains the same length and outputs... ; When entering a single-layer multi-head attention, each head's The shapes are respectively The attention matrix is After splicing Back Stacking The layer retains its shape; BiLSTM maps sequences to S45: The sequence-level vector is obtained after global pooling. .

[0014] Preferably, the specific process of motif extraction in step S5 is as follows: S51: Standardize the length of DNA sequences in the training subset and evaluation subset, uniformly adjust them to a fixed length of 1800bp, fill the part that is less than 1800bp with the unknown base "N", and cut off the part that exceeds 1800bp. S52: Based on the background base probability distribution of AGCT in the dataset, construct the corresponding one-hot form of baseline input; S53: Input the constructed baseline along with the actual input into the deep learning model, and use DeepLIFT to calculate the attribution contribution of each location to the prediction result.

[0015] Preferably, the specific process of step S53 is as follows: S531: The length-normalized real DNA sequence and the baseline sequence are used as input pairs and fed into the deep learning model. S532: By comparing the model's output with the baseline sequence input, backpropagation is used to calculate the contribution of each base position to the output difference. S533: Obtain a "site-by-site attribution contribution score sequence" of length 1800, denoted as: S attr =(s1, s2, ..., s 1800 ); in, s t Indicates the first t The attribution contribution score for each base position.

[0016] Preferably, the specific process of frequency domain noise reduction, candidate localization, and dual-threshold joint screening in step S5 is as follows: S54: Assigning contribution scores to each location sequentially S attr Input a low-pass filter, output a smoothed contribution score sequence S smooth ; S55: Contribution score sequence after smoothing S smooth In this process, the start and end positions of the signal peak are located, and the continuous regions of the signal peak correspond to potential motif segments. S56: Based on the preset target motif length, remove segments whose length exceeds the range; For candidate segments with an overlap of more than 80%, non-maximum suppression is used to remove duplicates, retaining the segment with the highest peak contribution score and deleting the remaining overlapping segments; Output a set of candidate motifs, each containing information such as start position, end position, length, corresponding DNA subsequence, and mean contribution score after smoothing; S57: By using a dual threshold screening method of statistical significance and effect size, a reliable set of selected motifs is finally obtained.

[0017] The beneficial effects of this invention include: 1. To address the problem that existing methods often rely on short-window peaks or fixed windows for local detection of base-by-base importance curves, resulting in insufficient capture capability for long / complex motifs spanning tens of bp and easily broken boundaries, which in turn affects downstream functional interpretation and experimental site selection, this invention aims to achieve continuous and consistent extraction of motifs of varying sizes (10–50 bp) under the same process through length-adjustable noise suppression and candidate localization that adaptively maps to the target length. This significantly improves the detection stability and reproducibility of long sequence scenarios.

[0018] 2. Since the importance curve contains a lot of high-frequency noise and is highly sensitive to thresholds / windows, the same model may exhibit unstable candidates, duplicate detections, and an increase in spurious peaks on different datasets. This invention aims to use frequency domain noise reduction and data-driven thresholding as the core, map the target length to the cutoff frequency and adaptively set the peak threshold, smoothly improve the signal-to-noise ratio and reduce hyperparameter sensitivity, thereby maintaining the consistency of the candidate set when migrating across datasets.

[0019] 3. To address the problem that relying solely on the significance p-value for screening leads to the inclusion of fragments that are "significant but have weak effects," and that while the FDR is controllable, it is difficult to guarantee the substance and stability of the biological effect, this invention proposes to adopt a dual-threshold joint criterion of "Fisher + FDR + effect size (OR and 95% confidence interval)" to retain statistically significant and clearly defined candidates, reduce false positives, and improve the reproducibility of cross-batch verification (maintaining a higher proportion of true positives under a unified q threshold). Attached Figure Description

[0020] Figure 1 This is a flowchart illustrating the deep learning-based interpretable DNA sequence motif extraction and visualization method of the present invention.

[0021] Figure 2 This is a schematic diagram of the sequence extraction process of the present invention. Detailed Implementation

[0022] The following is in conjunction with the appendix Figures 1-2 The present invention will be further described in detail below: Example 1 See appendix Figure 1 As shown, the deep learning-based method for interpretable DNA sequence motif extraction and visualization includes the following steps: S1: Obtain a publicly available human gene-dependent dataset containing DNA sequences, gene essentiality tags, and continuous fitness scores. The dataset covers multiple tissue types and cancer cell lines. Use the IsolationForest algorithm with a contamination rate of 0.01 to remove abnormally long sequences from the dataset, so that 99% of the DNA sequences after processing are concentrated within 10,000 bp. Define a set of core essential genes in the processed dataset based on the majority voting rule, and set the regression target as the fitness mean across cancer cell lines to construct training and evaluation subsets.

[0023] S2: Sequence Encoding: The DNA sequences in the training and evaluation subsets are encoded. Each base is one-hot encoded to obtain an initial matrix. If an unknown base N exists in the sequence, it is encoded as an all-zero vector or a single UNK vector. Linear projection is then used to... E=XW emb +b emb Where X is the one-hot encoding initial matrix, Wemb is the projection weight matrix, and bemb is the bias vector. The initial matrix is ​​converted into an embedding tensor E to align the encoding dimension with the number of channels in subsequent network layers. S3: Construct a sequentially connected feature extraction network based on the embedded tensor E. The feature extraction network includes a one-dimensional CNN, a multi-head self-attention Transformer encoder with N stacked layers, and a BiLSTM network. The one-dimensional CNN extracts local n-qram patterns in the sequence dimension through multi-scale kernels while keeping the sequence length unchanged, and outputs a CNN feature matrix Hcnn. The Transformer encoder performs long-range dependency aggregation on Hcnn and outputs an attention feature matrix H. The BiLSTM network models forward and backward time dependencies on H and outputs the bidirectional hidden state Hbi at each base position.

[0024] S4: Globally pool the bidirectional hidden state Hbi output by the BiLSTM network to obtain the sequence-level representation z. Construct a classification or regression output layer through a fully connected layer. Train the deep learning model to complete the gene necessity prediction with the gene necessity label or fitness mean as the objective.

[0025] S5: Site-by-site attribution calculation: Based on the AGCT background base probability distribution of DNA sequences in the training subset and evaluation subset of step 1, construct a one-hot baseline input; feed the baseline input and the real DNA sequence inputs in the training subset and evaluation subset into the trained deep learning model, and use the DeepLIFT algorithm to calculate the attribution contribution score of each base position to the model output result.

[0026] Frequency domain noise suppression and candidate localization: The obtained attribution contribution scores are input into an adjustable low-pass filter, which is an eighth-order Butterworth filter with the following transfer function: in, f c The cutoff frequency, f s The sampling frequency is used; the filter parameters are dynamically configured according to the mapping relationship between the target motif length and the normalized cutoff frequency (30bp corresponds to 0.1, 10bp corresponds to 0.2, and 5bp corresponds to 0.4) to suppress high-frequency noise in the attribution contribution score; the peak node of the filtered signal is located by extreme value detection to achieve stable positioning of continuous segments of different scales from 5 to 50bp, and a candidate motif set is obtained.

[0027] Dual-threshold joint screening: The candidate motif set was statistically screened. A 2×2 contingency table was constructed based on the occurrence frequency of each candidate motif in the positive / negative sample groups of the training subset and evaluation subset and the total number of sequences in each group. Fisher's one-sided exact test was used to evaluate the significance of motif enrichment, and FDR correction was performed using the Benjamini–Hochberg procedure. The odds ratio (OR) and its 95% confidence interval for each candidate motif were calculated. Candidate motifs that met the criteria of "corrected p<0.05, OR>2 and lower bound of confidence interval>1" were retained to obtain the screened motif set.

[0028] S6: The selected motif set is compared and verified with the HOCOMOCO standard database using the Tomtom tool to generate a complete motif list containing motif sequences, statistical test results, and database comparison evidence; a color fluctuation curve is plotted based on the attribution contribution score, peak and trough nodes are marked, and the motif regions after comparison and verification are annotated in orange to achieve visual output of the motifs.

[0029] In this embodiment, the specific process of preprocessing the dataset to obtain the training subset and the evaluation subset in step S1 is as follows: S11: Standardize the data format: unify the base encoding format of DNA sequences, retain only the four valid bases A, T, C, and G, mark unknown bases as N, and associate each sequence with the corresponding gene ID, cell line origin, tissue type, fitness score, and necessity label to form a structured data table.

[0030] Anomaly identification: The original data contains DNA sequences with extremely long tails, with some sequences longer than 40kb. These excessively long sequences can lead to significant differences in the input dimensions of subsequent deep learning models, increasing training difficulty and reducing model convergence stability. Therefore, it is necessary to remove abnormally long sequences first.

[0031] S12: The Isolation Forest algorithm is used for anomaly detection. This algorithm is suitable for handling outliers in high-dimensional data and has high accuracy in identifying anomalies in long-tailed data. The algorithm pollution rate is set to 0.01, that is, the proportion of abnormal sequences in the preset dataset does not exceed 1%. The algorithm models the length characteristics of all DNA sequences, automatically identifies and removes abnormally long sequences that exceed the normal length distribution range, and ensures that 99% of the DNA sequence length is concentrated within 10,000 bp, which significantly optimizes the consistency of the input space of subsequent models and improves training efficiency and stability.

[0032] S13: To address the issue of differing essentiality labels for the same gene in different cell lines, a majority voting rule is used to screen core essential genes: if a gene is identified as essential in more than half of the cell lines, then that gene is included in the core essential gene set. This rule can reduce the impact of single cell line detection errors on gene essentiality determination, ensure that the essentiality of genes within the core set is consistent across cell lines, and improve the reliability of labels for subsequent model training.

[0033] Regression target setting: The target value of the regression task is defined as the fitness mean across cancer cell lines. For each gene, the arithmetic mean of its fitness score in all corresponding cell lines is calculated. This is used as a quantitative indicator of the "necessity" of the gene, avoiding the interference of fitness score fluctuations in a single cell line on model training. S14: Based on the processed sequences (length compliant), core essential gene set (reliable labels), and mean fitness (regression target), the ratio of 7:3 or 8:2 can be flexibly adjusted according to the amount of data. It is necessary to ensure that the tissue type and gene type distribution of the two sets of data are consistent, and randomly divide them into training subset and evaluation subset.

[0034] The partitioning process must adhere to the stratified sampling principle: ensuring that the ratio of core essential genes to non-core genes, the sequence proportion of each tissue type, and the distribution range of fitness scores are consistent between the training subset and the evaluation subset, so as to avoid distortion of evaluation results due to data distribution bias, and finally obtain training and evaluation subsets that can be directly used for deep learning model training and generalization ability evaluation.

[0035] Example 2 Based on Example 1, the specific process of step S2 is as follows: S21: Take any DNA sequence s of length L from the training subset or evaluation subset, in the form s=( s 1 , s 2 ,...,s L ), where the bases at each position s t , t =1,2,...,L, belonging to only one of {A,T,C,G,N}, where A,T,C,G are valid bases and N is an unknown base.

[0036] S22: For each base in the sequence s t Perform one-hot encoding to generate a vector of dimension 4. x t This corresponds to the binary identifier of the four effective bases; like s t =A, then x t = onehot (A)=[1,0,0,0] ; like s t =T, then x t = onehot (T)=[0,1,0,0] ; like s t =C, then x t = onehot (C)=[0,0,1,0] ; like s t =G, then x t = onehot (G)=[0,0,0,1] ; like s t =N (unknown base), then x t Set to a vector of all zeros [0,0,0,0] Alternatively, it can be defined separately as an UNK vector [0.25, 0.25, 0.25, 0.25]. Both processing methods must be set uniformly before model training to ensure consistency of input format.

[0037] S23: Generate the one-hot vectors corresponding to all bases in the sequence. x 1, x 2,..., x L The sequences are vertically concatenated in sequence to form an initial matrix X with dimensions of "sequence length × base coding dimension", expressed as: For example: If the DNA sequence length L=5 and the bases are [A,T,N,C,G], then the initial matrix X is: ; S24: The initial matrix X with a fixed dimension of 4 is converted into an embedding tensor E that matches the number of channels of subsequent network layers (multi-head self-attention, BiLSTM), thus solving the dimension adaptation problem between low-dimensional base encoding and high-dimensional feature extraction layers.

[0038] The specific process of step S24 is as follows: Introducing the linear projection formula: ; in, The projection weight matrix has dimensions of . , The projection weight matrix is ​​iteratively optimized through backpropagation during model training to determine the number of feature channels for subsequent model layers (such as CNN and attention layers). Let be the bias vector, with dimension . This is used to compensate for the offset of the linear projection, and it is also optimized during model training; X is the initial matrix with dimension 1. , as input for linear projection; E The output embedding tensor must have dimensions that satisfy "sequence length × number of model channels" to obtain the embedding tensor. .

[0039] The essence of linear projection is to apply it to each row of the initial matrix X (a one-hot vector). x t A linear transformation is performed, and the specific process is as follows: S241: For the initial matrix X The Middlet row vector x t , with weight moments Perform matrix multiplication to obtain the intermediate vector. x t ; S242: Intermediate vector x t Superimposed bias This yields the embedding vector of the t-th base. e t = x t ; S243: Embedding vector of all bases e 1, e 2,..., e L By vertically concatenating the sequences, we obtain the embedded tensor E=[ e 1; e 2; ...; e L ] .

[0040] The specific process of step S3 is as follows: S31: One-dimensional convolution extracts local n-qram patterns using multi-scale kernels in the time dimension (sequence dimension) while maintaining the same length. The convolution output is the CNN feature matrix. ; S32: By stacking Transformer encoders consisting of self-attention and feedforward networks, long-range dependencies are aggregated, and an attention feature matrix is ​​output. ; S33: Use a bidirectional LSTM to model the time dependency in both forward and backward directions to obtain the bidirectional hidden state at each base position. .

[0041] Example 3 Based on Example 1 or Example 2, the specific process of global pooling the bidirectional hidden state Hbi output by the BiLSTM network in step S4 to obtain the sequence-level representation z is as follows: S41: The cells of BiLSTM at each time step Receive vectors from the attention block And the previous moment Calculate the forward and reverse chains separately. And spliced ​​together at each position as follows: ; Let the input of the attention layer be... , No. The query for the head, key, and value. in ; Scaling Dot Product Attention: ; S42: By concatenating all the heads and linearly mapping them, we obtain: ; S43: Perform residual summation and layer normalization, then pass through a feedforward network to form the output of an encoder layer; stacking The results are then fed into a BiLSTM layer. S44: When the input is a DNA sequence, the tensor shape of the data stream changes gradually along the pipeline; when , 256. Number of heads (thus LSTM hidden state dimension : The original onehot was Embedded / First Layer After convolution, we get One-dimensional convolution maintains the same length and outputs... When entering a single-layer multi-head attention, each head's The shapes are respectively The attention matrix is After splicing Back Stacking The layer retains its shape; then BiLSTM maps the sequence to... (Because bidirectional splicing is) S45: The sequence-level vector is obtained after global pooling. .

[0042] The specific process of motif extraction in step S5 is as follows: S51: For the DNA sequences in the training subset and the evaluation subset, first standardize the length and adjust them to a fixed length of 1800bp. Fill the part that is less than 1800bp with the unknown base "N" and cut off the part that exceeds the length to ensure that it is consistent with the input length during model training and avoid attribution bias. S52: Based on the background base probability distribution of AGCT in the dataset, construct the corresponding one-hot form of baseline input; For example, if the background probability of A is 25%, T is 25%, C is 25%, and G is 25%, then the one-hot vector at each position in the baseline sequence is [0.25, 0.25, 0.25, 0.25], or adjusted according to the actual statistical probability; the role of the baseline is to serve as a reference benchmark for non-specific signals, highlighting the specific contribution by comparing it with the real sequence. S53: Input the constructed baseline along with the actual input into the deep learning model, and use DeepLIFT to calculate the attribution contribution of each location to the prediction result.

[0043] In this embodiment, the specific process of step S53 is as follows: S531: The length-normalized real DNA sequence (one-hot encoded) and the baseline sequence (one-hot encoded) are used as input pairs and fed into the deep learning model. S532: By comparing the model output when the real sequence is input with the model output when the baseline sequence is input, backpropagation is used to calculate the contribution of each base position to the output difference. If the contribution of a certain base position is positive, it means that the presence of the base enhances the prediction tendency of "gene essentiality"; if it is negative, it inhibits the tendency; the larger the absolute value of the contribution, the more significant the influence of the base on the prediction result. S533: Obtain a "site-by-site attribution contribution score sequence" of length 1800, denoted as: S attr =(s1, s2, ..., s 1800 ); in, s t Indicates the first t The sequence is visualized using a "color fluctuation curve" (x-axis represents base position, y-axis represents contribution score, positive contribution is in red, negative contribution is in blue), which intuitively shows the distribution of high contribution regions.

[0044] The specific process of frequency domain noise reduction, candidate localization, and dual-threshold joint screening in step S5 is as follows: S54: Assigning contribution scores to each location sequentially S attr Input an 8th-order Butterworth low-pass filter, output a smoothed contribution score sequence. S smooth In this sequence, high-frequency noise is significantly suppressed, and the high-contribution region forms a continuous and smooth "signal peak", providing a clear signal for subsequent positioning; S55: Contribution score sequence after smoothing S smoothIn this context, the starting and ending positions of the "signal peak" are located, and the continuous regions of the signal peak correspond to potential motif segments (continuous sets of high-contribution bases). Local extremum detection algorithm is used: Traversal S smooth The contribution score of each position is determined as a "local peak" when the score of a position is greater than the scores of its two adjacent positions. When a value satisfies the condition of being "less than the fractions of its two adjacent positions", it is determined to be a "local valley value". Peak region definition: Using "two adjacent local valleys" as boundaries, the local peak region in the middle is enclosed to form a "continuous high contribution segment"; For example, if a valley is detected at position 100, the next valley is detected at position 200, and the peak is at position 150, then the continuous base fragment from 100 to 200 is identified as a "potential candidate motif fragment". S56: Based on the preset target motif length (5-50bp), remove fragments that exceed this range (too short may be noise residue, too long may contain irrelevant sequences); For candidate segments with an overlap of more than 80%, non-maximum suppression (NMS) is used to remove duplicates, retaining the segment with the highest peak contribution score and deleting the remaining overlapping segments to avoid duplicate candidates. Output: A set of candidate motifs is obtained. Each candidate motif contains information such as "start position, end position, length, corresponding DNA subsequence (e.g., CCAAGGAGG), and mean contribution score after smoothing". S57: By using a dual threshold screening method of "statistical significance + effect size", a reliable set of selected motifs is finally obtained; For each candidate motif, based on the positive and negative sample groups of the "training and evaluation subsets" (positive samples: core essential gene sequences; negative samples: non-core essential gene sequences), the following four key data points are collected, and a 2×2 contingency table is constructed, including: a The number of sequences in the positive sample group that contain this candidate motif; b The number of sequences in the positive sample group that do not contain the candidate motif; c The number of sequences in the negative sample group that contain this candidate motif; d The number of sequences in the negative sample group that do not contain the candidate motif; Statistical significance test: A one-sided Fisher's exact test (with the direction set to "positive sample group enrichment") was used to assess whether the frequency of candidate motifs in the positive sample group was significantly higher than that in the negative sample group. Validation logic: Based on the hypergeometric distribution, calculate the probability of "observing the current contingency table or a more extreme case", i.e., the original p-value; the smaller the p-value, the more significant the enrichment of the candidate motif in the positive sample group. FDR Correction: Since multiple candidate motifs are tested simultaneously (multiple tests), the Benjamini–Hochberg (BH) procedure is used to correct the false discovery rate (FDR) of the original p-values ​​to control the overall false positive rate. Threshold setting: retain candidate motifs with "corrected p-value < 0.05" and remove statistically insignificant fragments; Effect size and confidence interval constraints: Statistical significance alone cannot guarantee the substantiality of a biological effect; further screening using effect size (odds ratio, OR) and its confidence interval is necessary. Odds ratio (OR) calculation: OR quantifies the difference in frequency of candidate motifs between the positive and negative sample groups. The formula is: OR = a×d / b×c ; OR > 1: Indicates that the candidate motif is more common in the positive sample group; the larger the OR, the stronger the effect; the technical disclosure sets OR > 2 as the effect size threshold to ensure that the biological effect of the candidate motif is substantial; 95% confidence interval (CI) calculation: Calculate the 95% confidence interval of OR using the normal approximation method or the Wilson score interval method; the narrower the confidence interval, the more stable the OR estimate; Threshold setting: Retain candidate motifs with "OR > 2 and 95% CI lower bound > 1"—a CI lower bound > 1 ensures that at the 95% confidence level, the OR is significantly greater than 1, avoiding bias due to... 4. The output of the filtered motif set combines the above two screening steps and finally retains the candidate motifs that simultaneously meet the criteria of "corrected p < 0.05, OR > 2, and 95% CI lower bound > 1", forming the "filtered motif set". Each motif needs to be associated with complete information: sequence information: start position, end position, length, DNA subsequence; signal information: mean contribution score after smoothing, peak position; statistical information: original p value, corrected p value, OR value, 95% CI. This set provides reliable input for subsequent motif validation and visualization.

[0045] In another implementation, to achieve the goal of stably and reproducibly extracting and validating biologically significant DNA motifs from long sequences, a series of interchangeable or combinable alternative pathways can be adopted without changing the overall technical objectives and effects. At the attribution signal level, DeepLIFT can be replaced by methods such as integral gradient, DeepSHAP, Grad-CAM / Grad-CAM++ (one-dimensional adaptation), SmoothGrad, and Gradient⊙Input. Baseline construction, in addition to using one-hot encoding with AGCT background probability, can also utilize permutation sequences that maintain dinucleotide frequencies, Markov backgrounds with k-mer frequency matching, random masking and reverse complementary backgrounds, or generate multiple baselines through Monte Carlo sampling and take the desired attribution to enhance the robustness of baseline selection. Length normalization can be changed from fixed-length padding to variable-length batch processing or "position-relative" encoding (relative transcription start sites, exon boundaries, etc.), and attribution aggregation of ultra-long sequences can be achieved through sliding windows or block overlap.

[0046] In the frequency domain and noise suppression stages, the eighth-order Butterworth low-pass filter can be replaced with a Chebyshev type I / II or elliptic filter to obtain a steeper pass / stopband transition. Alternatively, Savitzky-Golay smoothing, Gaussian / bidirectional exponential moving average, LOESS / LOWESS local regression smoothing, and learnable window functions implemented by FFT convolution can be used. Discrete cosine transform (DCT), discrete wavelet transform (DWT, Daubechies, Symlet, Coiflet, etc.) and continuous wavelet transform (CWT) can be introduced for multi-scale decomposition and reconstruction, thereby extracting the main trend of the target length range in the manner of "band selection + inverse transform". The mapping relationship between the cutoff frequency and the target length can be changed from empirical formula to data-driven self-calibration (stacked grid search, Bayesian optimization) or learning-based parameterized mapping (SincNet / learnable filter kernel), and a robust threshold can be automatically set by generating an "empty distribution" in the domain through permutation test on different organizations or platforms.

[0047] Candidate detection and boundary determination can be extended from extremum detection to peak model fitting and segmentation methods, including zero-crossing localization based on second derivatives, region labeling of HMM / semi-Markov models, sequence labeling of conditional random fields (CRF), aggregation of high-attribution segments using density clustering (DBSCAN / HDBSCAN), and consensus peak merging after dynamic time warping (DTW) alignment. Redundancy can be controlled using non-maximum suppression (NMS), hierarchical merging based on segment IoU, topology-preserving interval pruning, and targeted minimum redundancy maximum correlation (mRMR) selection. Consensus motif generation can be achieved by directly summarizing into PWM / PFM, or by correcting through information weighting and position weighting, or by first aligning with Smith–Waterman / Needleman–Wunsch before constructing consensus. For composite or nested structures, a composite PWM group can be formed using "multi-template parallelism + cross-attention," or the co-occurrence structure of "motif topics" can be learned on the segment set using topic models (LDA / CTM).

[0048] In terms of statistical screening and multiple correction, Fisher's one-sided exact test and BH-FDR can be replaced by procedures such as Storeyq-value, Benjamini–Yekutieli (dependency robustness), and Holm–Bonferroni (family error rate control), or converted to a generalized linear model (logistic regression / Firth correction) to directly regress "whether it belongs to a positive sample", and correct for GC content, sequence complexity, and fragment length as covariates; in addition to the odds ratio (OR), the effect size can be combined with log2 fold-enrichment, Cohen'sh, or relative risk (RR), and combined with bootstrapping / permutation to generate confidence intervals; at the evidence integration layer, a Bayesian hierarchical model can be used to output posterior enrichment and Bayes factor to obtain more robust statistical conclusions in multi-data source or multi-batch scenarios.

[0049] In summary, the deep learning-based interpretable DNA sequence motif extraction and visualization method provided in this invention proposes a set of interpretable and transferable DNA motif extraction and verification systems for long sequences. It uses site-by-site attribution as the core signal source, and without altering the original model structure, constructs a one-hot baseline based on the AGC background distribution of the dataset, and obtains positional contribution scores along with the real input using DeepLIFT. To address the issues of weak signals and unstable boundaries in long sequences, an adjustable low-pass filter digital signal processing pipeline is introduced, establishing a mapping relationship between the "target motif length" and the "normalized cutoff frequency." A higher-order Butterworth algorithm is employed to enhance frequency domain selectivity, achieving stable localization of continuous fragments at different scales (5–50 bp) with the assistance of extremum detection. To ensure statistical robustness and biological relevance, a dual-threshold joint screening method combining Fisher one-sided filtering, Benjamini-Hochberg FDR, and effect size thresholding is proposed, and annotation loop closure is completed using standard library similarity comparison.

[0050] Meanwhile, by using length standardization and non-maximum suppression for deduplication, the burden of cross-dataset parameter tuning and manual verification costs are reduced, ultimately outputting a list of motifs with statistical and database evidence chains and visualization results. Compared with existing schemes that rely on fixed windows and local fragment clustering, the significant distinguishing feature of this invention is that it embeds frequency domain noise reduction and gating screening with a length-frequency reversible mapping into the attribution-driven extraction process, uniformly solving the problems of scale consistency and boundary stability of long / short and single / composite motifs. This significantly improves cross-scenario reproducibility and engineering usability without changing the downstream model and evaluation criteria.

Claims

1. A deep learning-based method for interpretable DNA sequence motif extraction and visualization, characterized in that, Includes the following steps: S1: Obtain a publicly available human gene-dependent dataset containing DNA sequences, gene necessity tags, and continuous fitness scores, and preprocess the dataset to obtain a training subset and an evaluation subset; S2: Gene sequence encoding is performed on each base in the DNA sequences of the training subset and the evaluation subset to obtain an initial matrix. The initial matrix is ​​then converted into an embedding tensor E to align the encoding dimension with the number of channels in subsequent network layers. S3: Based on the embedding tensor E, a sequentially connected feature extraction network is constructed to extract features, resulting in the CNN feature matrix Hcnn, the attention feature matrix H, and the bidirectional hidden state Hbi at each base position; S4: Globally pool the bidirectional hidden state Hbi output by the BiLSTM network to obtain the sequence-level representation z. Construct a classification or regression output layer through a fully connected layer. Train the deep learning model to complete the gene necessity prediction with the gene necessity label or fitness mean as the target. S5: Motif extraction and screening: Perform site-by-site attribution calculation, frequency domain noise reduction and candidate localization, and dual-threshold joint screening to obtain the screened motif set; S6: Compare and verify the selected motif set with the standard database to generate a complete motif list containing motif sequences, statistical test results and database comparison evidence; Draw a color fluctuation curve, mark the peak and valley nodes, and annotate the sequence regions after comparison and verification with orange to achieve sequence visualization output.

2. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 1, characterized in that, The specific process of preprocessing the dataset in step S1 to obtain the training subset and the evaluation subset is as follows: S11: Data format standardization: unify the base encoding format of DNA sequences, retain only the four valid bases A, T, C, and G, mark unknown bases as N, and associate each sequence with the corresponding gene ID, cell line origin, tissue type, fitness score, and necessity label to form a structured data table; S12: The Isolation Forest algorithm is used for anomaly detection. The algorithm contamination rate is set to 0.

01. The algorithm models the length characteristics of all DNA sequences, automatically identifies and removes abnormally long sequences that exceed the normal length distribution range, so that the DNA sequence length is concentrated within 10,000 bp. S13: Select core essential genes using majority voting rules: If a gene is determined to be an essential gene in more than half of the cell lines, then the gene is included in the core essential gene set. S14: Based on the processed sequences, the set of core essential genes, and the mean fitness value, the training subset and the evaluation subset are randomly divided according to a specified ratio.

3. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 1, characterized in that, The specific process of step S2 is as follows: S21: Take any DNA sequence s of length L from the training subset or evaluation subset, in the form s=( s 1 ,s 2 ,..., s L ); S22: For each base in the sequence s t Perform one-hot encoding to generate a vector of dimension 4. x t This corresponds to the binary identifier of the four effective bases; S23: Generate the one-hot vectors corresponding to all bases in the sequence. x 1, x 2,..., x L The sequences are vertically concatenated in sequence to form an initial matrix X with dimensions of "sequence length × base coding dimension", expressed as: S24: Transform the initial matrix X, with a fixed dimension of 4, into an embedding tensor that matches the number of channels in subsequent network layers. E .

4. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 3, characterized in that, The specific process of step S24 is as follows: Introducing the linear projection formula: ; in, The projection weight matrix is... This represents the number of feature channels for subsequent model layers; X is the bias vector used to compensate for the offset of the linear projection; X is the initial matrix. E For the output embedding tensor; S241: For the initial matrix X The Middle t row vector x t , with weight moments Perform matrix multiplication to obtain the intermediate vector. x t ; S242: Intermediate vector x t Superimposed bias This yields the embedding vector of the t-th base. e t = x t ; S243: Embedding vector of all bases e 1, e 2,..., e L By vertically concatenating the sequences, we obtain the embedded tensor E=[ e 1; e 2; ...; e L ] .

5. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 1, characterized in that, The specific process of step S3 is as follows: S31: One-dimensional convolution extracts local n-qram patterns in the time dimension using multi-scale kernels while maintaining the same length. The convolution output is the CNN feature matrix. ; S32: By stacking Transformer encoders consisting of self-attention and feedforward networks, long-range dependencies are aggregated, and an attention feature matrix is ​​output. ; S33: Use a bidirectional LSTM to model the time dependency in both forward and backward directions to obtain the bidirectional hidden state at each base position. .

6. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 1, characterized in that, The specific process of performing global pooling on the bidirectional hidden state Hbi output by the BiLSTM network in step S4 to obtain the sequence-level representation z is as follows: S41: The cells of BiLSTM at each time step Receive vectors from the attention block And the previous moment Calculate the forward and reverse chains separately. And spliced ​​together at each position as follows: ； Let the input of the attention layer be... , No. The query, key, and value of each head are: in ; Scaling Dot Product Attention: ； S42: By concatenating all the heads and linearly mapping them, we obtain: ； S43: Perform residual summation and layer normalization, then pass through a feedforward network to form the output of an encoder layer; stacking The results are then fed into a BiLSTM layer. S44: When the input is a DNA sequence, the tensor shape of the data stream changes gradually along the pipeline; when , 256. Number of heads LSTM hidden state dimension : The original onehot is Embedded / First Layer After convolution, we get One-dimensional convolution maintains the same length and outputs... ; When entering a single-layer multi-head attention, each head's The shapes are respectively The attention matrix is After splicing Back Stacking The layer retains its shape; BiLSTM maps sequences to S45: The sequence-level vector is obtained after global pooling. .

7. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 1, characterized in that, The specific process of motif extraction in step S5 is as follows: S51: Standardize the length of DNA sequences in the training subset and evaluation subset, uniformly adjust them to a fixed length of 1800bp, fill the part that is less than 1800bp with the unknown base "N", and cut off the part that exceeds 1800bp. S52: Based on the background base probability distribution of AGCT in the dataset, construct the corresponding one-hot form of baseline input; S53: Input the constructed baseline along with the actual input into the deep learning model, and use DeepLIFT to calculate the attribution contribution of each location to the prediction result.

8. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 7, characterized in that, The specific process of step S53 is as follows: S531: The length-normalized real DNA sequence and the baseline sequence are used as input pairs and fed into the deep learning model. S532: By comparing the model's output with the baseline sequence input, backpropagation is used to calculate the contribution of each base position to the output difference. S533: Obtain a "site-by-site attribution contribution score sequence" of length 1800, denoted as: S attr =(s1,s2,...,s 1800 ); in, s t Indicates the first t The attribution contribution score for each base position.

9. The method for extracting and visualizing interpretable DNA sequence motifs based on deep learning according to claim 7, characterized in that, The specific process of frequency domain noise reduction, candidate localization, and dual-threshold joint screening in step S5 is as follows: S54: Assigning contribution scores to each location sequentially S attr Input a low-pass filter, output a smoothed contribution score sequence S smooth ; S55: Contribution score sequence after smoothing S smooth In this process, the start and end positions of the signal peak are located, and the continuous regions of the signal peak correspond to potential motif segments. S56: Based on the preset target motif length, remove segments whose length exceeds the range; For candidate segments with an overlap of more than 80%, non-maximum suppression is used to remove duplicates, retaining the segment with the highest peak contribution score and deleting the remaining overlapping segments; Output a set of candidate motifs, each containing information such as start position, end position, length, corresponding DNA subsequence, and smoothed mean contribution score; S57: By using a dual threshold screening method of statistical significance and effect size, a reliable set of selected motifs is finally obtained.

Images