Method for enhancer-promoter interaction prediction based on multi-omics data

Abstract

This invention proposes a method for predicting enhancer promoter interactions based on multi-omics data, belonging to the field of bioinformatics technology. It addresses the problems of commonly used techniques failing to adequately consider nonlinearity when integrating features and neglecting global analysis of EPIs from a multi-cell line perspective. The method includes: Step 1: Constructing a deep neural network model; Step 2: Acquiring multi-omics data and constructing a multi-omics dataset; Step 3: Preprocessing the multi-omics dataset; Step 4: Training the deep neural network model based on the preprocessed multi-omics dataset; Step 5: Predicting enhancer promoter interactions based on the trained deep neural network model.

Description

Technical Field

[0001] This invention relates to a method for predicting enhancer promoter interactions based on multi-omics data, and belongs to the field of bioinformatics technology. Background Technology

[0002] Enhancers and promoters are two important DNA sequence elements in cells. Enhancers play a crucial role in the transcriptional regulation of genes, enhancing the expression of target genes. Promoters, acting as "switches" for gene expression, are typically defined as a region 1500 bp upstream and 500 bp downstream of the transcription start site (TSS). As cis-regulatory elements in the genome, enhancers regulate gene expression through interactions with the promoters of target genes. Recent studies of the three-dimensional genome have revealed that distal enhancers can interact with proximal promoters to regulate the expression of target genes. Accurate identification of EPIs is crucial for understanding gene regulation, cell differentiation, and disease mechanisms. Several variants associated with cancer susceptibility exist in enhancer regions. For example, the 8q24 region contains multiple functional enhancers associated with various tumor types. However, accurate identification of EPIs remains challenging due to the often large distances between enhancers and their target promoters.

[0003] Previous studies inferred EPIs by expressing quantitative trait loci (eQTLs), but eQTL mapping typically requires large sample sizes and primarily identifies short-range EPIs. With the development of high-throughput sequencing technologies, researchers can now utilize techniques such as paired-end-tag sequencing (ChIA-PET) and high-throughput chromosome conformation capture (Hi-C) to detect long-range chromatin interactions, which can be used to identify EPIs. For example, Hi-C technology has revealed multiple topologically related domains (TADs), showing that mammalian chromatin interactions occur significantly more frequently within these TADs than outside them. However, while these high-throughput technologies demonstrate powerful capabilities in detecting EPIs, they are time-consuming, and sample preparation, library construction, and sequencing equipment are expensive.

[0004] To address these issues, numerous computational methods have been proposed for large-scale, rapid, and accurate identification of EPIs. These methods apply data that can be broadly categorized into two types: DNA sequences and epigenetic data. Most of these computational methods use single-cell-line prediction models, only identifying EPIs specific to that cell line. Identifying EPIs from other tissues requires model retraining. Therefore, constructing a global EPI prediction model—that is, a global analysis of EPIs from a multi-cell-line perspective—is essential. Summary of the Invention

[0005] To address the issues of insufficient consideration of nonlinearity and lack of global analysis of EPIs from a multi-cell line perspective when integrating features, this invention proposes a method for predicting enhancer promoter interactions based on multi-omics data.

[0006] The technical solution adopted by the present invention to solve the above problems is as follows: The present invention specifically includes:

[0007] Step 1: Construct a deep neural network model;

[0008] Step 2: Acquire multi-omics data and construct a multi-omics dataset;

[0009] Step 3: Preprocess the multi-omics dataset;

[0010] Step 4: Train the deep neural network model based on the preprocessed multi-omics dataset;

[0011] Step 5: Predict enhancer promoter interactions based on the trained deep neural network model.

[0012] Preferably, the deep neural network model in step 1 consists of a CNN network, a BiLSTM network, and a Kansformer encoder;

[0013] CNN networks are convolutional neural networks used to identify local patterns and key signal features in input multi-omics datasets;

[0014] BiLSTM networks are bidirectional long short-term memory networks used for feature processing from both the forward and backward directions of multi-omics datasets, enhancing the ability to capture temporal dependencies.

[0015] The Kansformer encoder is used to capture long-distance dependencies in multi-omics datasets.

[0016] Preferably, the Kansformer encoder is composed of a fusion of KAN layers and Transformer, including alternating stacks of multi-head self-attention layers and KAN layers. The multi-head self-attention layers are used to capture the relationships between features, and the KAN layers are used to enhance the prediction of sub-promoter interaction relationships.

[0017] Preferably, the multi-omics data in step 2 includes CTCF binding sites, DNase-I signals, five histone modification features, and DNA methylation information. The multi-omics datasets include the Hi-C dataset, the CHIA-PET dataset, and feature signals.

[0018] Preferably, step 3 specifically includes:

[0019] Step 3.1: For the Hi-C and CHIA-PET datasets, enhancer-gene pairs are converted into enhancer-promoter pairs. Following the negative sample generation strategy, the measured enhancers are paired with genes that have not been experimentally verified to interact with them. The gene selection criterion is that the distance between the gene and the corresponding enhancer is within 95% of the enhancer-gene distance in the positive samples. This completes the preprocessing of the Hi-C and CHIA-PET datasets.

[0020] Step 3.2: For the feature signals, store the feature signals in the bigwig file and the narrowPeak file, convert the bigwig file and the narrowPeak file into pt files that store tensors, and normalize the features to complete the preprocessing of the feature signals.

[0021] Preferably, step 4 specifically includes:

[0022] Step 4.1: The preprocessed multi-omics dataset is used as the input sequence, and features are extracted from the input sequence through a CNN network;

[0023] Step 4.2: The BiLSTM network receives the extracted features and enhances the ability of the extracted features to capture temporal dependencies through forward and backward LSTM networks, outputting a matrix Z, Z∈R. l×h R is a matrix containing enhancer and promoter site features, where R is the real number field, l is the number of rows in the matrix, corresponding to the number of sequences, and h is the number of columns in the matrix, corresponding to the dimension of the features.

[0024] Step 4.3: Input the output matrix Z of the BiLSTM network into the Kansformer encoder, and capture the relationships between the features in the output matrix Z through a multi-head self-attention layer to obtain the enhanced features.

[0025] Step 4.4: Enhance the features The input to the KAN layer is used for feature integration to obtain the output H0∈R of the Kansformer encoder. l×h ;

[0026] Step 4.5: Use the self-attention embedding method to learn the low-dimensional representation of H0, and obtain the aggregated features H∈R. 4h ;

[0027] Step 4.6: Input the aggregated features H into the KAN classifier to calculate the probability p used to predict enhancer-promoter interaction and the distance d between enhancer and promoter. p ;

[0028] Step 4.7: Based on the calculated enhancer-promoter interaction probability p and the distance d between the enhancer and the promoter. p Constructing a binary cross-entropy loss function for EPI prediction and the mean squared error loss function for predicting enhancer promoter distance. Based on binary cross-entropy loss function Mean squared error loss function Construct the total loss function

[0029] Step 4.8: Based on the total loss function The deep neural network model is trained to obtain the trained deep neural network model;

[0030] The formula for calculating the enhancer-promoter interaction probability p is:

[0031] p=σ(KAN1(H))* (1);

[0032] In formula (1), σ is the sigmoid function, p is the probability of enhancer-promoter interaction, and the value of p ranges from 0 to 1;

[0033] The distance d between enhancer and promoter p The calculation formula is:

[0034] d p =KAN2(H)(2);

[0035] Binary cross-entropy loss function The expression is:

[0036]

[0037] In formula (3), y i p represents the true label (0 or 1) of the i-th sample. i The probability of prediction;

[0038] Mean squared error loss function The expression is:

[0039]

[0040] In formula (4), d p For the predicted enhancer promoter distance, d t For the actual distance;

[0041] Total loss function The expression is:

[0042]

[0043] In formula (5), This is the total loss function.

[0044] Preferably, step 4.3 specifically includes:

[0045] Step 4.3.1: Calculate the query Q, key Q, and value vector V based on matrix Z and the weight matrix used for mapping transformation;

[0046] Step 4.3.2: Calculate the global relation matrix M based on query Q and key Q;

[0047] Step 4.3.3: Perform feature enhancement on the value vector V using the global relation matrix M to obtain the enhanced features.

[0048] The formulas for calculating query Q, key Q, and value vector V are as follows:

[0049] Q = ZW q K = ZW k V=ZW v (6);

[0050] In formula (6), W q W k and W v These are three weight matrices used for mapping transformation, and K, Q, and V are three feature codes of matrix Z after linear transformation;

[0051] The formula for calculating the global relation matrix M is:

[0052]

[0053] In formula (7), d is the feature dimension of matrix Z;

[0054] Enhanced features The calculation formula is:

[0055]

[0056] In formula (8), b is the deviation.

[0057] Preferably, in step 4.4, the KAN layer is a K-layer KAN network, which includes multiple nested KAN units. The feature integration steps include:

[0058] Step 4.4.1: For each KAN layer, through the activation function For input The activation value is obtained after transformation.

[0059] Step 4.4.2: Sum the post-activation values ​​of all inputs to obtain the activation values ​​of the neuron and the output H0∈R of the Kansformer encoder. l×h ;

[0060] The expression for the KAN network is:

[0061] KAN(Z)=(Φ K-1 ,Φ K-2 ,…,Φ1,Φ0)(Z)(9);

[0062] In formula (9), Φ i Let be the i-th layer of the KAN network. For each KAN layer, the input dimension is n. in The output dimension is n out Each layer of the KAN network is defined as follows:

[0063]

[0064] Post-activation value The calculation formula is:

[0065]

[0066] Preferably, step 4.5 specifically includes:

[0067] Step 4.5.1: Input H0 into a two-layer fully connected network to obtain matrix A. s ;

[0068] Step 4.5.2: Connect H0 with matrix A s Multiplying them together yields the weighted embedding H1∈R r×l ;

[0069] Step 4.5.3: Calculate the mean and maximum value of the second dimension of H1, and concatenate the mean and maximum value of the second dimension of H1 as a low-dimensional embedding into H1;

[0070] Step 4.5.4: Concatenate the low-dimensional embedding with the hidden states at the locations of the enhancer and promoter to obtain the aggregated feature H∈R. 4h ;

[0071] Matrix A s The calculation formula is:

[0072]

[0073] In formula (12), A s ∈R r×l Matrix A s The sum of each row of elements is 1, and each row of elements is a set of attention coefficients at each position of H0;

[0074] The aggregated features H∈R 4h The calculation formula is:

[0075] H=AvgPool(H1)||MaxPool(H1)||h e ||h p (13).

[0076] The beneficial effects of this invention are:

[0077] 1. This invention proposes the KansformerEPI model, which uses a fusion encoder of KAN and Transformer to predict enhancer-promoter interactions using sequence features and epigenetic features. By using KAN to integrate features, the nonlinear relationships between features are fully considered, strengthening the important connections between features and significantly improving the accuracy of the model in predicting enhancer-promoter interactions.

[0078] 2. The KansformerEPI model achieves excellent predictive performance on various mantle cell lines. Compared with other state-of-the-art predictive models, the KansformerEPI model demonstrates superior performance across multiple evaluation metrics. Attached Figure Description

[0079] Figure 1 A flowchart illustrating the method for predicting enhancer promoter interactions based on multi-omics data provided by the present invention;

[0080] Figure 2 A schematic diagram of the structure of the deep neural network model provided by this invention;

[0081] Figure 3 The performance comparison chart of KansformerEPI with other methods for predicting EPIs on the HMEC, IMR90, K562, and NHEK datasets provided in this invention; Figure 3 In the diagram, a is a schematic diagram showing the difference in AUC index between the KansformerEPI model and the TransEPI model constructed in this invention on four cell lines, and b is a schematic diagram showing the difference in AUPR index between the KansformerEPI model and the TransEPI model constructed in this invention on four cell lines.

[0082] Figure 4 A comparison chart showing the impact of adding a methylation feature set on EPIs prediction results provided by this invention;

[0083] Figure 5 A comparison chart of the HMEC, IMR90, K562, and NHEK datasets provided by this invention after data balancing and after imbalanced data balancing;

[0084] Figure 6 Performance comparison chart of KansformerEPI model and TransEPI model after balancing data provided in this invention;

[0085] Figure 7 This invention provides feature importance analysis and window importance analysis of the KansformerEPI model on the GM12878 dataset. Figure 7 In the figure, a is a schematic diagram of the relative importance of the seven features of the KansformerEPI model on the GM12878 cell line, and b is a schematic diagram of the differences in importance scores of the seven features in different windows.

[0086] Figure 8 Box plots showing the predicted probability distribution of the connection modes between enhancers and promoters provided by this invention in four cell lines. Detailed Implementation

[0087] Specific implementation method one: Combining Figure 1 and Figure 2 This implementation method is described as follows: Figure 1 As shown, the steps of the method for predicting enhancer promoter interactions based on multi-omics data described in this embodiment include:

[0088] S1: Acquire multi-omics data and construct a multi-omics dataset;

[0089] The multi-omics data includes CTCF binding sites, DNase-I signals, five histone modification features, and DNA methylation information. The multi-omics datasets include the Hi-C dataset, the CHIA-PET dataset, and feature signals.

[0090] S2: Preprocessing multi-omics datasets;

[0091] S201: For the Hi-C and CHIA-PET datasets, enhancer-gene pairs are converted into enhancer-promoter pairs. Following the negative sample generation strategy, the measured enhancers are paired with genes that have not been experimentally verified to interact with them. The gene selection criterion is that the distance between the gene and the corresponding enhancer is within 95% of the enhancer-gene distance in the positive samples. This completes the preprocessing of the Hi-C and CHIA-PET datasets.

[0092] S202: For the feature signal, store the feature signal in the bigwig file and the narrowPeak file, convert the bigwig file and the narrowPeak file into a pt file to store tensors, and normalize the features to complete the preprocessing of the feature signal.

[0093] S3: Construct a deep neural network model;

[0094] like Figure 2 As shown, the deep neural network model used in this embodiment is the KansformerEPI model, which consists of a CNN network, a BiLSTM network, and a Kansformer encoder.

[0095] CNN networks are convolutional neural networks used to identify local patterns and key signal features in input multi-omics datasets;

[0096] BiLSTM networks are bidirectional long short-term memory networks used for feature processing from both the forward and backward directions of multi-omics datasets, enhancing the ability to capture temporal dependencies.

[0097] The Kansformer encoder is used to capture long-distance dependencies in multi-omics datasets. It consists of a fusion of KAN layers and Transformer layers, including alternating stacks of multi-head self-attention layers and KAN layers. The multi-head self-attention layers are used to capture the relationships between features, while the KAN layers are used to enhance the prediction of sub-promoter interaction relationships.

[0098] S4: Train the deep neural network model using the preprocessed multi-omics dataset;

[0099] S401: Features are extracted from the input signal using a Convolutional Neural Network (CNN). The convolutional operations of the CNN layers can effectively identify local patterns and key signal features in the input sequence, providing high-quality feature representations for subsequent layers; BiLSTM enhances the ability to capture temporal dependencies through its forward and backward LSTM layers, outputting matrix Z.

[0100] S402: The output Z of the BiLSTM is fed into the Kansformer encoder, which uses a Transformer-like structure to capture long-range dependencies. Multi-Head Self-Attention is used in the encoder.

[31] The mechanism can capture the relationships between features and calculate the query, key, and value vectors. The calculation formula is as follows:

[0101] Q = ZW q K = ZW k V=ZW v (1);

[0102] In formula (1), matrix Z∈R l×h W represents a matrix containing features of enhancer and promoter sites. q W k and W vThese represent three weight matrices used for mapping transformation, while K, Q, and V represent the three feature codes of Z after linear transformation.

[0103] Then, the global relation matrix M is calculated using K and Q:

[0104]

[0105] In formula (2), d represents the feature dimension of Z. Next, feature enhancement is performed on V using M to obtain the enhanced features.

[0106]

[0107] S403: By using KAN layers, the non-linear relationship between features is taken into account when integrating different features. A K-layer KAN network can be represented as a nested network of multiple KAN units:

[0108] KAN(Z)=(Φ K-1 ,Φ K-2 ,…,Φ1,Φ0)(Z)(4);

[0109] In formula (4), Φ i Let be the i-th layer of the KAN network. For each KAN layer, the input dimension is n. in The output dimension is n out Each layer of the KAN network is defined as follows:

[0110]

[0111] For each KAN layer, through the activation function For input The activation value is obtained after transformation. Summing the post-activation values ​​of all inputs yields the activation values ​​of the neuron and the output H0∈R of the Kansformer encoder. l×h ;

[0112] Post-activation value The calculation formula is:

[0113]

[0114] S404: Convert the output H0∈R from the Kansformer encoder l×h Low-dimensional representation learning of H0 is performed using a self-attention embedding method:

[0115] Input H0 into a two-layer fully connected network:

[0116]

[0117] In formula (7), A s ∈R r×l Matrix A s The sum of each row of elements is 1, and each row of elements is a set of attention coefficients at each position of H0;

[0118] Connect H0 and A s Multiplying them together yields the weighted embedding H1∈R r×l Then, the mean and maximum values ​​of the second dimension of H1 are calculated and concatenated as a low-dimensional embedding to H1. Finally, this low-dimensional embedding is concatenated with the hidden states at the enhancer and promoter positions to obtain the final aggregated feature H∈R. 4h ;

[0119] The aggregated features H∈R 4h The calculation formula is:

[0120] H=AvgPool(H1)||MaxPool(H1)||h e ||h p (8);

[0121] S405: The aggregated H is input into the KAN classifier to predict the enhancer-promoter interaction probability p. Additionally, to make the model sensitive to the positions of enhancers and promoters, a KAN module is used to predict the distance d between enhancers and promoters. p :

[0122] The formula for calculating the enhancer-promoter interaction probability p is:

[0123] p=σ(KAN1(H))* (9);

[0124] In formula (1), σ is the sigmoid function, p is the probability of enhancer-promoter interaction, and the value of p ranges from 0 to 1;

[0125] The distance d between enhancer and promoter p The calculation formula is:

[0126] d p =KAN2(H)(10);

[0127] S406: Based on the calculated enhancer-promoter interaction probability p and the distance d from enhancer to promoter. p Constructing a binary cross-entropy loss function for EPI prediction and the mean squared error loss function for predicting enhancer promoter distance. Based on binary cross-entropy loss function Mean squared error loss function Construct the total loss function

[0128] Binary cross-entropy loss function The expression is:

[0129]

[0130] In formula (11), y i p represents the true label (0 or 1) of the i-th sample. i The probability of prediction;

[0131] Mean squared error loss function The expression is:

[0132]

[0133] In formula (12), d p For the predicted enhancer promoter distance, d t For the actual distance;

[0134] Total loss function The expression is:

[0135]

[0136] In formula (13), This is the total loss function.

[0137] S407: Based on the total loss function The deep neural network model is trained to obtain the trained deep neural network model;

[0138] S5: Predict enhancer promoter interactions based on a trained deep neural network model.

[0139] In summary, this implementation uses an encoder that integrates KAN and Transformer to predict enhancer-promoter interactions using sequence features and epigenetic features. By using KAN to integrate features, the nonlinear relationships between features are fully considered, strengthening the important connections between features and significantly improving the accuracy of the model in predicting enhancer-promoter interactions.

[0140] Specific Implementation Method Two: Combining Figures 3-8This embodiment will be described in detail. In order to fully evaluate the effectiveness of the prediction method in this embodiment, the constructed KansformerEPI model was compared with five state-of-the-art deep learning models (TransEPI, TargetFinder, SPEID, 3DPredictor, and DeepTACT). This embodiment uses the same chromosome segmentation scheme to train all models, thereby eliminating errors that may be introduced by differences in data segmentation. The GM12878 and HeLa-S3 datasets were used for training, and the HMEC, IMR90, K562, and NHEK datasets were used for testing.

[0141] In this implementation, 30 experiments were run on each cell line, and the model was trained and tested identically according to the model parameters provided by TransEPI. To further evaluate the predictive performance of the KansformerEPI model, this implementation compared it with the existing TransEPI model, and statistically analyzed the mean differences between the two methods in the four cell lines using a two-pronged t-test. This implementation set a significance level of 0.05. Under this condition, if the p-value was less than 0.05, the predictive performance of the two methods on the corresponding cell lines was considered significantly different; conversely, if the p-value was greater than 0.05, the difference was considered insignificant. Figure 3 As shown in Figure a, the differences in AUC between the two methods across four cell lines are illustrated. The results demonstrate that the KansformerEPI method of this embodiment outperforms TransEPI in all four cell lines, and the differences are statistically significant. Specifically, in the IMR90 and NHEK cell lines, KansformerEPI improved the AUC by approximately 7% and 8%, respectively, while in the HMEC and K562 cell lines, the improvement was approximately 4%. Figure 3 As shown in b, the differences in AUPR index between the two methods on four cell lines are illustrated. The effects were comparable on the K562 cell line, but there were significant differences on the HMEC, IMR90, and NHEK cell lines, with all showing significant improvement.

[0142] Table 1 shows the performance comparison of the model of this embodiment with other existing models in terms of AUROC and AUPR. In terms of AUROC, the model of this embodiment improved the performance of the best-performing models in the HMEC, IMR90, and K562 cell lines by 0.54%, 5.7%, and 3.3%, respectively. However, in the NHEK cell line, the model of this embodiment was 2.7% lower than the best-performing model. In terms of AUPR, this embodiment found that in the K562 cell line, the model of this embodiment performed comparably to the best-performing models; while in the HMEC, IMR90, and NHEK cell lines, the model of this embodiment outperformed the runner-up model, improving performance by 11.7%, 11.1%, and 4.2%, respectively. Overall, despite a slight difference in the NHEK cell line, the model of this embodiment outperformed other existing models in most cell lines.

[0143] Table 1

[0144]

[0145] In this embodiment, based on information from CTCF binding sites, DNase-I signals, and five histone modifications, DNA methylation data was further introduced. After extracting the methylation signals, all features were normalized. Referring to the research of Whalen et al., this embodiment obtained methylation feature sets for four cell lines: GM12878, HeLa-S3, IMR90, and K562. These four cell lines were selected as the dataset for studying the impact of methylation on EPI prediction results. For training and testing, this embodiment used the previous method to fuse the data from the GM12878 and HeLa-S3 cell lines to construct a training set. Independent testing was then conducted on the IMR90 and K562 cell lines. During the testing process, this embodiment conducted 30 rounds of experiments, and the results are as follows: Figure 4 As shown, in the IMR90 dataset, the AUROC index is comparable, but the AUPR index is approximately 1.67% higher. In another dataset, K562, both AUROC and AUPR are approximately 1% lower. While methylation shows a certain promoting effect on EPI prediction on some datasets, the opposite result was obtained on other datasets. This discrepancy suggests that the effect of methylation on EPI prediction may not be universally applicable. Currently, this embodiment only tested data from these two methylated datasets, therefore, there is insufficient evidence to prove that methylation has a significant promoting effect on EPI prediction. Although some studies have shown that methylation has a promoting effect on EPI prediction, this phenomenon was not observed in this embodiment.

[0146] In the task of predicting enhancer-promoter interactions, the original dataset is extremely imbalanced, with some ratios reaching 1:20. Positive samples are determined through wet-lab experiments, while negative samples are generated by pairing measured enhancers with non-interacting genes. To observe the impact of data balance on predicted EPIs, this paper uses an undersampling strategy in the experimental design to adjust the positive-to-negative sample ratio to 1:1. This implementation suggests that using a 1:1 balanced dataset can effectively reduce class bias in the model. Since the model may be biased towards the majority class (negative samples) in an imbalanced dataset, a balanced dataset allows the model to focus more on positive samples, thereby improving its ability to identify positive samples. This data balancing strategy helps improve the stability of the model's performance on positive samples, such as... Figure 5 As shown, this embodiment performed a two-pronged t-test analysis on the balanced and imbalanced datasets. It can be seen that the AUPR values ​​of all four cell lines on the balanced dataset were significantly higher than those on the imbalanced dataset. Through this dataset adjustment, this embodiment demonstrates a significant improvement in the AUPR index on the balanced dataset. It also helps evaluate the model's generalization ability in detecting minority class samples, allowing for better observation of the model's performance differences on balanced and imbalanced data. This verifies whether the model can accurately identify enhancer-promoter interactions when the positive and negative sample distributions change. Furthermore, this embodiment compares its performance with other state-of-the-art models on the balanced dataset; the comparison results are as follows. Figure 6 As shown, the models in this embodiment are far superior to TransEPI, which demonstrates that the models in this embodiment are also better than TransEPI models on balanced datasets, and are also more stable.

[0147] This implementation utilizes a random forest model to conduct an in-depth analysis of key features of enhancer-promoter interactions, focusing on the importance of features and their contributions within different windows. To this end, this implementation first calculates the importance score of each feature across multiple windows to assess its impact on predicted EPIs. To ensure the systematic nature of the analysis, this implementation selects some windows from 5000 windows where enhancers and promoters are relatively important, and calculates the contribution of seven features in each window. Through this windowed analysis method, this implementation can observe the changes in the importance of each feature within different gene regulatory regions, thereby identifying which features have a significant impact on predicted EPIs within certain specific windows, such as... Figure 7As shown in Figure a, this embodiment calculated the importance scores of these features within the window and took their average value to obtain the relative importance of the seven features in the GM12878 cell line. The analysis results showed that the DNase-I signal had the highest score, indicating its most prominent importance in prediction. In contrast, the CTCF binding site had the lowest score, while the H3K4me1 and H3K4me3 features were relatively important. This may be because these two features typically appear in enhancer regions, and enhancers are relatively active in EPI prediction.

[0148] Subsequently, this implementation uses heatmaps to visualize the importance score distribution of features within each window, such as... Figure 7 As shown in b, the heatmap illustrates the differences in importance scores for seven features (including CTCF, DNase-I, and five histone modification signals) within different windows. The color intensity in the heatmap reflects the feature importance score; darker colors indicate higher importance for predicting EPIs, while lighter colors indicate lower importance. Through heatmap analysis, this implementation method can more intuitively understand the contribution of each feature in different genomic regions. e and p represent enhancer and promoter windows, while e-1 and p-1 windows represent nearby windows. Within enhancer windows, this implementation method observes that H3K4me1 has the darkest color; this feature frequently appears in enhancer regions. In the window near the enhancer, the features H3K9me3, H3K27me3, and CTCF are found to be of higher importance. Similarly, in the promoter window, the feature H3K9me3 is the darkest and more important. In the window near the promoter, the features H3K4me1 and H3K4me3 are more important. This observation provides new biological insights into the regulatory mechanism of EPI and further verifies the key role of these features in gene regulation. In addition, these results provide important reference for subsequent in-depth biological research and help this implementation method better understand EPI and its biological functions.

[0149] Table 2

[0150]

[0151] Enhancer-promoter interactions typically form subnetworks of various biological processes. To better analyze and understand these interaction patterns, this embodiment calculates the number of promoters regulated by each enhancer in the sample (out-degree) and the number of enhancers regulated by each promoter (in-degree). In the data of this embodiment, on average, each enhancer regulates approximately 1 to 2 promoters, but a promoter is regulated by 2 to 3 enhancers. This embodiment is influenced by Sushmita Roy.

[32] Inspired by others, the connection patterns between these enhancers and promoters were categorized into four types: (i) Single interaction: a one-to-one relationship where one enhancer interacts with only one promoter; (ii) Multiple-output component (MO): one enhancer regulates multiple promoters; (iii) Multiple-input component (MI): one promoter is simultaneously regulated by multiple enhancers; (iv) Multiple-input multiple-output component (MIMO): multiple enhancers regulate multiple promoters, forming a network. In analyzing four test cell lines, this embodiment statistically analyzed the number of these four regulatory relationships and their predicted probabilities in the model of this embodiment. The results showed that MI had the largest proportion among the four regulatory relationships, as shown in Table 2. Simultaneously, through... Figure 8 Box plots were used to observe the predicted value distributions of the four regulatory relationships in this implementation, revealing that the predicted probability of MI was generally optimal. This result indicates that in the BENGI dataset, the MI relationship is not only dominant in quantity but also the best in terms of prediction accuracy. These examples, combined, confirm the excellent predictive performance of KansformerEPI.

[0152] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent substitutions, and improvements made to the above embodiments without departing from the scope of the present invention, based on the technical essence of the present invention and within the spirit and principles of the present invention, shall still fall within the protection scope of the present invention.

Claims

1. A method for predicting enhancer promoter interactions based on multi-omics data, characterized in that, The steps of the method for predicting enhancer promoter interactions based on multi-omics data include: Step 1: Construct a deep neural network model; The deep neural network model in step 1 consists of a CNN network, a BiLSTM network, and a Kansformer encoder; The CNN network is a convolutional neural network used to identify local patterns and key signal features in the input multi-omics dataset; The BiLSTM network is a bidirectional long short-term memory network used for feature processing from both the forward and reverse directions of multi-omics datasets, enhancing the ability to capture temporal dependencies. The Kansformer encoder is used to capture long-distance dependencies in multi-omics datasets; The Kansformer encoder is composed of a fusion of KAN layers and Transformer layers, including alternating stacks of multi-head self-attention layers and KAN layers. The multi-head self-attention layers are used to capture the relationships between features, and the KAN layers are used to enhance the prediction of sub-promoter interaction relationships. Step 2: Acquire multi-omics data and construct a multi-omics dataset; Step 3: Preprocess the multi-omics dataset; Step 4: Train the deep neural network model based on the preprocessed multi-omics dataset; Step 4 specifically includes: Step 4.1: The preprocessed multi-omics dataset is used as the input sequence, and features are extracted from the input sequence through the CNN network; Step 4.2: The BiLSTM network receives the extracted features and enhances the ability of the extracted features to capture time dependencies through forward and backward LSTM networks, outputting a matrix. Z， Let R be a matrix containing features of enhancer and promoter sites, where R is a real number field. The number of rows in the matrix corresponds to the number of sequences. The number of columns in the matrix corresponds to the dimension of the feature. Step 4.3: Convert the output matrix of the BiLSTM network Z The input to the Kansformer encoder is captured by a multi-head self-attention layer to capture the output matrix. Z The relationship between the features is used to obtain the enhanced features. ; Step 4.4: Enhance the features The input to the KAN layer is used for feature integration to obtain the output of the Kansformer encoder. ; Step 4.5: Use the self-attention embedding method to... Perform low-dimensional representation learning to obtain aggregated features. H ; Step 4.6: Aggregate the features H The input is fed into the KAN classifier to calculate the probability used to predict enhancer-promoter interactions. And predict the distance between enhancers and promoters. ; Step 4.7: Based on the calculated enhancer-promoter interaction probability and the distance between enhancers and promoters Constructing a binary cross-entropy loss function for EPI prediction and the mean squared error loss function for predicting enhancer promoter distance. Based on the binary cross-entropy loss function Mean squared error loss function Construct the total loss function ; Step 4.8: Based on the total loss function The deep neural network model is trained to obtain the trained deep neural network model; enhancer-promoter interaction probability The calculation formula is: （1）； In formula (1), For the sigmoid function, To increase the probability of promoter-subpromoter interactions, The value ranges from 0 to 1; Distance between enhancer and promoter The calculation formula is: （2）； Binary cross-entropy loss function The expression is: （3）； In formula (3), For the first i The true label of each sample (0 or 1). The probability of prediction; Mean squared error loss function The expression is: （4）； In formula (4), For the predicted enhancer promoter distance, For the actual distance; Total loss function The expression is: （5）； In formula (5), This is the total loss function; Step 5: Predict enhancer promoter interactions based on the trained deep neural network model.

2. The method for predicting enhancer promoter interactions based on multi-omics data according to claim 1, characterized in that, Step 2 includes multi-omics data such as CTCF binding sites, DNase-I signals, five histone modification features, and DNA methylation information. The multi-omics datasets include the Hi-C dataset, the CHIA-PET dataset, and feature signals.

3. The method for predicting enhancer promoter interactions based on multi-omics data according to claim 1, characterized in that, Step 3 specifically includes: Step 3.1: For the Hi-C and CHIA-PET datasets, enhancer-gene pairs are converted into enhancer-promoter pairs. Following the negative sample generation strategy, the measured enhancers are paired with genes that have not been experimentally verified to interact with them. The gene selection criterion is that the distance between the gene and the corresponding enhancer is within 95% of the enhancer-gene distance in the positive samples. This completes the preprocessing of the Hi-C and CHIA-PET datasets. Step 3.2: For the feature signals, store the feature signals in the bigwig file and the narrowPeak file, convert the bigwig file and the narrowPeak file into pt files that store tensors, and normalize the features to complete the preprocessing of the feature signals.

4. The method for predicting enhancer promoter interactions based on multi-omics data according to claim 1, characterized in that, Step 4.3 specifically includes: Step 4.3.1: Based on matrix Z The query is obtained by calculating the weight matrix used for mapping changes. ,key Sum value vector ; Step 4.3.2: Based on query s and keys The global relation matrix is ​​calculated. M ; Step 4.3.3: Through the global relation matrix M value vector Perform feature enhancement to obtain enhanced features. ; Query ,key Sum value vector The calculation formula is: （6）； In formula (6), ， and These are three weight matrices used for mapping changes. K, Q, V For matrix Z The three feature codes after linear transformation; Global Relationship Matrix M The calculation formula is: （7）； In formula (7), d For matrix Z Feature dimensions; Enhanced features The calculation formula is: （8）； In formula (8), This is a deviation.

5. The method for predicting enhancer promoter interactions based on multi-omics data according to claim 1, characterized in that, In step 4.4, the KAN layer is a K-layer KAN network. A K-layer KAN network consists of multiple nested KAN units. The feature integration steps include: Step 4.4.1: For each KAN layer, through the activation function For input The activation value is obtained after transformation. ; Step 4.4.2: Sum the post-activation values ​​of all inputs to obtain the activation values ​​of the neurons and the output of the Kansformer encoder. ; The expression for the KAN network is: （9）； In formula (9), For the KAN network i For each KAN layer, the input dimension is... The output dimension is Each layer of the KAN network is defined as follows: （10）； Post-activation value The calculation formula is: （11）。 6. The method for predicting enhancer promoter interactions based on multi-omics data according to claim 1, characterized in that, Step 4.5 specifically includes: Step 4.5.1: [The text appears to be incomplete and contains several grammatical errors. A more accurate translation would require the full context.] The input is fed into a two-layer fully connected network to obtain a matrix. ; Step 4.5.2: [The text appears to be incomplete and contains several grammatical errors. A more accurate translation With matrix Multiply to obtain weighted embedding. ; Step 4.5.3: For The second dimension calculates the average and maximum values, and... The mean and maximum values ​​of the second dimension are concatenated as a low-dimensional embedding. ; Step 4.5.4: Concatenate the low-dimensional embedding with the hidden states at the locations of the enhancer and promoter to obtain the aggregated feature H. ; matrix The calculation formula is: （12）； In formula (12), ,matrix The sum of the elements in each row is 1, and the elements in each row are... A set of attention coefficients at each location; Aggregated features H The calculation formula is: （13）。

Images