A Multimodal Fusion-Based Method for Predicting Chromatin Topological Association Domain Boundaries

Abstract

This invention discloses a method for predicting the boundaries of chromatin topological association domains (TADs) based on multimodal fusion. Combining word embedding concepts, it employs k-mer embedding encoding on DNA sequence data and fuses epigenetic modification data and Hi-C interaction frequency data to construct a multi-omics dataset. This dataset is preprocessed, and a deep neural network model is constructed for multimodal feature extraction. The extracted results are flattened, fused, and used for boundary prediction. The model is trained on a training set using the backpropagation algorithm to obtain a model for TAD boundary prediction. This method improves the accuracy of TAD boundary prediction, provides biologically meaningful predictions, and reduces time and financial costs.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, specifically to a method for predicting the boundaries of chromatin topological association domains based on multimodal fusion. Background Technology

[0002] Chromatin is the form of the genome in eukaryotes, composed of coiled and folded DNA strands. Genomic regulatory elements that are linearly distant may be close together in three-dimensional space. The three-dimensional structure of chromatin influences DNA replication, transcription, and regulatory processes. Existing research has shown that changes in chromosome structure in specific genomic regions and under specific conditions are highly correlated with various human diseases, including cancer, limb deformities, and brain developmental abnormalities. Therefore, studying the three-dimensional structure of chromatin is of great value in explaining cellular processes such as gene expression and regulation, thereby guiding research on genetic diseases and medical issues such as gene therapy.

[0003] Topologically Associated Domains (TADs) are structural units on chromatin, discovered through chromatin interaction data. A TAD is a contiguous region of chromatin where chromatin interactions occur much more frequently than between different regions. Three-dimensional chromatin structural modeling studies have shown that TADs are independent structures on chromatin and participate as a whole in the formation of higher-order chromatin structures.

[0004] The TAD boundary is a chromatin region between different TADs, acting as an insulating "fence" that restricts the interaction between cis-regulatory sequences and target genes. TAD boundaries are rich in complex trait heritability, evolutionary constraints, CTCF protein binding sites, and housekeeping genes. Studies have shown that TAD boundaries are closely related to the pathogenesis of genetic diseases, and the development of many rare diseases is caused by disruption of TAD boundaries.

[0005] Due to the uncertainties surrounding TAD length and formation mechanisms, predicting TAD boundaries using traditional biological experiments would be extremely costly in terms of manpower, financial resources, and time. Therefore, to better determine whether a chromatin fragment represents a TAD boundary, researchers need reliable computational models to guide biological experiments in predicting TAD boundaries.

[0006] Currently, based on the different data used in the models, the main computational methods in this field can be divided into two categories: methods based on epigenetic modification data and prediction methods based on DNA sequences. The basic principle of methods based on epigenetic modification data is to design corresponding computational models for prediction by utilizing the enrichment or absence of CTCF protein binding sites, histone modifications, and other gene regulatory elements present at the TAD boundary; however, these methods suffer from difficulties in feature selection, leading to limitations in prediction accuracy. In recent years, researchers have discovered that DNA sequences can be used to accurately predict chromatin interaction frequencies and epigenetic modification signal intensity, thus DNA sequences are increasingly used for TAD boundary prediction. DNA sequence-based prediction methods input DNA sequence features near the TAD boundary into the model, but the main problem is that human DNA sequences are too long, making feature extraction difficult, and DNA data lacks cell line specificity. In addition, Hi-C data can also be used for TAD-related research; it is a high-resolution chromatin regulatory element interaction map obtained by studying the spatial relationships of DNA across the entire chromatin genome using high-throughput sequencing technology.

[0007] Existing methods for predicting the boundaries of chromosome topological association domains (TADs) have some shortcomings. Regarding data acquisition, TAD prediction methods only include epigenetic data or DNA sequence data, and lack Hi-C interaction frequency data that is highly correlated with TADs, thus the model accuracy could be further improved. In terms of model construction, DNA coding data is missing and is not integrated with Hi-C interaction frequency features, which to some extent affects the performance of existing methods. Summary of the Invention

[0008] To address the aforementioned problems, the present invention aims to provide a method for predicting the boundaries of chromatin topological association domains (TADs) based on multimodal fusion. This method constructs a multi-omics dataset containing DNA sequence data, epigenetic modification data, and Hi-C interaction frequency data, extracts data features, and inputs them into a multimodal fusion model to improve the accuracy of TAD boundary prediction. The technical solution is as follows:

[0009] A method for predicting the boundary of chromatin topological association structural domains based on multimodal fusion includes the following steps:

[0010] S1: Data Acquisition: Acquire multi-omics datasets, including Hi-C dataset, epigenetic modification dataset, and DNA sequence dataset, and obtain TAD boundary label data through models;

[0011] S2: Data preprocessing: Based on the TAD boundary label data, feature pre-extraction is performed on the multi-omics dataset, dataset cleaning is completed, and the DNA dataset is encoded and embedded.

[0012] S3: Model Building: Construct a deep neural network model and extract multimodal features, including DNA sequence features, epigenetic modification features, and Hi-C data features;

[0013] S4: Model prediction: Flatten and splice the DNA sequence features, epigenetic modification features and Hi-C data features to complete multimodal feature fusion, and predict the boundary results based on the fusion results;

[0014] S5: Model Training: Train the model on the training set using the backpropagation algorithm to obtain a model for TAD boundary prediction.

[0015] Furthermore, S1 specifically includes:

[0016] S11: Download Hi-C data, epigenetic modification data, and DNA sequence data for specific types of human cell lines;

[0017] S12: Apply the TopDom model to the above Hi-C data to obtain positive sample label data of TAD boundary; obtain negative sample label data of TAD boundary by randomly selecting the same number of non-boundary sites as the positive sample label data.

[0018] Furthermore, step S2 specifically includes:

[0019] S21: Construct a TAD boundary label data file based on the TAD boundary label data, extract DNA sequence fragments of a specified region and resolution from the file, and construct a DNA sequence dataset;

[0020] S22: Based on the TAD boundary label data file, extract the signal intensity in specific regions upstream and downstream of each boundary label location as the epigenetic modification features of the boundary.

[0021] S23: Based on the Hi-C data, extract the interaction frequency data in specific areas upstream and downstream of each boundary label position to obtain the Hi-C interaction frequency feature vector for each sample;

[0022] S24: Remove invalid samples from the DNA sequence dataset;

[0023] S25: Perform k-mer encoding, one-hot encoding, and embedding operations on the DNA sequence data to obtain a DNA embedding coding matrix, which reflects the relationship between the current sequence and the structural characteristics of its location.

[0024] Furthermore, the specific process of obtaining the DNA embedding coding matrix in S25 is as follows:

[0025] S251: For a DNA sequence of length n, S = (s1, s2, ..., s...). nK-mer encoding is performed, which involves moving a sliding window from the beginning to the end of the sequence and extracting a subsequence of length k each time.

[0026] S252: Perform one-hot encoding on the k-mer subsequence, concatenate them in window order, and pad with zeros at the end to align the lengths to obtain the one-hot encoded matrix.

[0027] S253: An embedding layer is introduced after the one-hot encoding matrix, using the embedding matrix. With one-hot encoding matrix Multiplying these components reduces the feature dimension, resulting in the DNA embedding coding matrix D∈R. n×d_embed , where d_embed represents the embedding dimension.

[0028] Furthermore, step S3, multimodal feature extraction, includes:

[0029] S31: DNA sequence feature extraction: Input the DNA embedding coding matrix into a deep neural network for DNA feature extraction to extract DNA sequence features;

[0030] S32: Epigenetic modification feature extraction: Input the epigenetic modification feature data into the deep neural network for epigenetic modification feature extraction to extract the epigenetic modification features, which are used to learn the enrichment or absence of eight histone modification signals and CTCF binding factors at several sites upstream and downstream of the current chromatin fragment.

[0031] S33: Hi-C Data Feature Extraction: Input Hi-C data into a Hi-C feature extraction deep neural network to extract Hi-C signal features, which are used to learn the interaction frequency features between the current chromatin fragment and several upstream and downstream sites.

[0032] Furthermore, the DNA feature extraction deep neural network in S31 consists of 3 layers;

[0033] A1: Layers 1 and 2 are convolutional neural network layers, each containing two sub-layers: a one-dimensional convolutional sub-layer and a pooling sub-layer. Features are input into the one-dimensional convolutional sub-layer, and the calculation formula is as follows:

[0034]

[0035] In this context, the superscript "1" indicates a feature related to the DNA sequence data; This is the output of the one-dimensional convolutional sublayer in this layer; x 1 The input features for the convolutional neural network layers represent epigenetic modification features in layer 1 and the output of layer 1 in layer 2; W 1 ∈R out_c_1×k_size_1The convolution kernel parameters are automatically updated through network training, where out_c_1 is the channel dimension of the output features, and k_size_1 is the convolution kernel size parameter; "*" indicates a convolution operation; b 1 This is a one-dimensional convolution bias.

[0036] Each one-dimensional convolutional sublayer is followed by a batch normalization process, calculated using an activation function, as follows:

[0037]

[0038] in, This represents the output after passing through the activation function; ReLU indicates activation via the ReLU activation function; BatchNorm is the implementation of batch normalization.

[0039] The pooling sublayer implements one-dimensional max pooling to reduce feature length and preserve effective features within the window. The calculation formula is as follows:

[0040]

[0041] Among them, maxpool is the implementation of one-dimensional max pooling functionality. This is the output of the pooling sublayer.

[0042] A2: The output of the first pooling sub-layer is input into the second layer, and then output after passing through the same one-dimensional convolutional sub-layer and pooling sub-layer in sequence;

[0043] A3: The third layer includes a BiLSTM sublayer and a Dropout sublayer; the output of the second layer is input into the BiLSTM sublayer to learn long-range relationships in the DNA sequence, and the calculation formula is as follows:

[0044]

[0045] Among them, BiLSTM is the implementation of the BiLSTM sublayer function, which adopts a bidirectional long short-term memory network; This is the output of the second layer; The output of the BiLSTM sublayer;

[0046] Obtain the output of the BiLSTM sublayer Then through the Dropout sublayer:

[0047]

[0048] Here, dropout is the dropout function, which randomly deactivates nodes in the network to reduce overfitting; For dropout probability, This is the result of DNA feature extraction.

[0049] Furthermore, the deep neural network for extracting epigenetic modification features in S32 consists of three layers:

[0050] B1: Layers 1 and 2 are convolutional neural network layers; features are input into the one-dimensional convolutional neural network, and the calculation formula is as follows:

[0051]

[0052] In this context, the superscript "2" indicates a characteristic related to epigenetic modification data; This is the output of the one-dimensional convolutional neural network in this layer; x 2 The input features for the convolutional neural network layers represent epigenetic modification features in layer 1 and the output of layer 1 in layer 2; W 2 ∈R out_c_2×k_size_2 The parameters for the convolution kernels are automatically updated during network training. `out_c_2` represents the channel dimension of the output features, i.e., the number of convolution kernels, and `k_size_2` is the kernel size parameter. "*" indicates a convolution operation. 2 This is a one-dimensional convolution bias.

[0053] Each one-dimensional convolutional network includes a batch normalization process, calculated using an activation function, as follows:

[0054]

[0055] in, The output is the result of batch normalization after passing through the activation function.

[0056] B2: The output of layer 1 is input into layer 2, and then passed through the same convolutional neural network and batch normalized before being output.

[0057] B3: Layer 3 includes a BiLSTM sublayer and a Dropout sublayer. The output of layer 2 is input into the BiLSTM sublayer to learn the long-range relationships of epigenetic modification data. The calculation formula is as follows:

[0058]

[0059] in, This is the output of the second layer; The output of the BiLSTM sublayer;

[0060] Obtain the output of the BiLSTM sublayer Then through the Dropout sublayer:

[0061]

[0062] in, The results of epigenetic modification feature extraction; This represents the dropout probability.

[0063] Furthermore, the Hi-C feature extraction deep neural network in S33 consists of four layers:

[0064] C1: Batch Normalization Layer 1: Batch normalization is performed on the Hi-C interaction frequency feature vector. The calculation formula is as follows:

[0065]

[0066] The superscript "3" indicates a feature related to Hi-C data; x 3 Data representing epigenetic modification characteristics; This represents the output result of batch normalization layer one.

[0067] C2: BiLSTM layer: Learns the long-distance relationships between interaction frequency data, calculated using the following formula:

[0068]

[0069] in, This is the output of the BiLSTM layer.

[0070] C3: Batch Normalization Layer 2: Results of the BiLSTM Layer Output Batch normalization is performed using the following formula:

[0071]

[0072] in, This represents the output result of batch normalization layer 2.

[0073] C4: Dropout layer: Uses random deactivation to reduce network overfitting. The calculation formula is as follows:

[0074]

[0075] in, The results of feature extraction from Hi-C data; This represents the dropout probability.

[0076] Furthermore, step S4 specifically includes:

[0077] S41: Flatten the DNA sequence features, epigenetic modification features, and Hi-C data features obtained in S3 respectively, and calculate them using the following formulas;

[0078]

[0079] in, The result is the flattened result of DNA sequence features (t=1), epigenetic modification features (t=2), and Hi-C interaction frequency data features (t=3). "flatten" is the flattening process. The extracted DNA sequence features, epigenetic modification features, and Hi-C interaction frequency data features;

[0080] S42: Features after stitching and flattening are used to complete multimodal fusion. The calculation formula is as follows:

[0081]

[0082] Where, x fusion This is the result of multimodal fusion; concat indicates vector concatenation.

[0083] S43: Combine the fused result x fusion The input consists of m consecutive fully connected layers; each fully connected layer includes a feedforward neural network with a linear transformation, an activation function layer, and a Dropout layer, calculated as follows:

[0084]

[0085] Wherein, the superscript “(i)” represents the relevant features of the i-th fully connected layer, and “(0)” represents the relevant features of the initial input fully connected layer; Let the dropout probability be such that the value of a node in the i-th fully connected layer is increased by a probability. It becomes 0; GELU indicates activation processing via the GELU activation function; The output of the linear transformation in the i-th fully connected layer is calculated as follows:

[0086]

[0087] in, This represents the output of the (i-1)th fully connected layer. This represents the initial input of the fully connected layer, i.e., x. fusion W (i) b represents the weight parameters of the i-th fully connected layer; (i) W represents the bias parameters of the i-th fully connected layer. (i) b (i) All are learned automatically by the model;

[0088] S44: Pass the output of the fully connected layer through an activation function to output the prediction result. The calculation formula is as follows:

[0089]

[0090] Where, x predict For each sample, the predicted result is either 0 or 1; sigmoid indicates activation via the sigmoid activation function, as shown in the formula below:

[0091] sigmoid(x) = 1 / (1+e) -x )

[0092] Among them, e -x It represents the reciprocal of the power of the natural base e.

[0093] Furthermore, step S5 specifically includes:

[0094] S51: Use gradient descent to reduce model error; employ the Adam optimizer to accelerate model convergence; use the binary cross-entropy loss function for model training, as shown in the following formula:

[0095]

[0096] Among them, y j This represents the predicted value of the j-th sample; This represents the label value of the j-th sample; J represents the total number of samples.

[0097] S52: The learning rate is dynamically reduced as the model iterates using a constant step decay method, as shown in the following formula:

[0098]

[0099] Where, α decayed α represents the initial learning rate after decay; α represents the decay rate; decay_rate represents the decay coefficient; global_steps represents the global step size; decay_steps represents the decay step size.

[0100] S53: Employ an early stop strategy to end training prematurely.

[0101] The beneficial effects of this invention are:

[0102] 1) This invention proposes a prediction method for the TAD boundary prediction problem, which combines multi-omics data, including DNA sequence data, epigenetic modification data, and Hi-C interaction frequency data, to improve prediction accuracy:

[0103] 2) The method described in this invention fully considers the multimodal contextual information of TAD sites during the feature extraction process, providing more effective interactive features for subsequent models;

[0104] 3) This invention designs a deep neural network for feature extraction and combines it with multimodal fusion technology to fully learn data features;

[0105] 4) Compared with existing models that do not fuse multimodal data on the same dataset, it achieved the best predictive performance;

[0106] 5) Validate the method proposed in this invention on existing biological data to prove that the method has biological significance and feasibility, and achieves the goal of reducing time and financial resources. Attached Figure Description

[0107] Figure 1 This is a flowchart of the chromatin topological association domain boundary prediction method and system based on multimodal fusion of the present invention.

[0108] Figure 2 This is a performance comparison chart between the method proposed in this invention and methods using non-fused multimodal data (DNA sequence data, epigenetic modification data, and Hi-C interaction frequency data).

[0109] Figure 3 This is a graph showing the results of the analysis of the enrichment of histone modification signals at the TAD boundary using the method proposed in this invention.

[0110] Figure 4 This is a graph showing the results of the analysis of the motif enrichment level of the TAD boundary predicted by the method proposed in this invention. Detailed Implementation

[0111] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0112] This invention proposes a method and system for predicting chromatin topological association domain (TAD) boundaries based on multimodal fusion. A multi-omics dataset containing DNA sequence data, epigenetic modification data, and Hi-C interaction frequency data is constructed. Data features are extracted and input into a multimodal fusion deep neural network, outputting prediction results. Experiments demonstrate that compared to models without multimodal fusion, this invention improves the accuracy of TAD boundary prediction. This embodiment provides a method for predicting chromatin topological association domain (TAD) boundaries based on multimodal fusion. The process is based on Python 3.9 and the PyTorch framework, implementing some functions by calling interfaces provided by the PyTorch framework. The method includes:

[0113] S1: Data Acquisition: Acquire multi-omics datasets, including Hi-C dataset, epigenetic modification dataset, and DNA sequence dataset, and obtain TAD boundary label data through models;

[0114] S2: Data preprocessing: Based on the TAD boundary label data, perform feature pre-extraction on the multi-omics dataset and complete the dataset cleaning; perform encoding and embedding processing on the DNA dataset;

[0115] S3: Model Building: Build a deep neural network model and perform multimodal feature extraction;

[0116] S4: Model Prediction: Flatten and fuse the extracted results for boundary result prediction;

[0117] S5: Model Training: Train the model on the training set using the backpropagation algorithm to obtain a model for TAD boundary prediction.

[0118] Specifically, the detailed process of the above five steps is as follows:

[0119] I. Data Acquisition

[0120] The data in this invention are based on five human cell lines, from which Hi-C data, DNA sequence data, and TAD boundary data (tags) were extracted. The five cell lines include: erythroleukemia cells (K562), human B lymphocytes (GM12878), human microvascular endothelial cells (HMEC), human epidermal keratinocytes (NHEK), and human embryonic lung fibroblasts (IMR90).

[0121] 1. Hi-C Data: Downloaded from the GEO database of the National Center for Biotechnology Information (NCBI), access number GSE63525. This dataset includes data at resolutions of 1Mb, 100kb, 50kb, 25kb, 10kb, and 5kb. This invention selects to use the 10kb resolution.

[0122] 2. Epigenetic modification data: including eight histone modification signals and CTCF-binding proteins, all downloaded from the ENCODE project. The eight histone modification data are: H3k4me1, H3k4me2, H3k4me3, H3k9ac, H3k9me3, H3k27ac, H3k27me3, and H3k36me3.

[0123] 3. DNA Sequence Data: Downloaded from the UCSC Genome Browser, with hg19 as the reference genome. DNA sequence data is typically stored in FASTA and FASTQ file formats, containing chromatin labels and base sequences. The base sequence is expressed as a, g, c, and t, representing the four bases. Due to sequencing accuracy limitations, some bases in certain regions were not accurately detected; these are represented by N to indicate unknown bases. Therefore, the DNA sequence data contains five possible scenarios. Each FASTA file contains two lines: the first line stores the sequence location information in the format ">chrNum:Start-End", indicating which chromosome the fragment belongs to and its position on that chromosome; the second line stores the base information for that chromosome fragment.

[0124] 4. TAD Boundary Data (Labels): Positive sample label data is obtained by applying the TopDom model to the Hi-C dataset mentioned above, while negative sample label data is obtained by randomly selecting the same number of non-boundary sites.

[0125] II. Data Preprocessing

[0126] 1. TAD Boundary Data Preprocessing: This invention selects Hi-C data with a resolution of 10kb and uses KRnorm for preprocessing to eliminate data noise caused by distance. Considering the novelty of the TAD prediction model and its acceptance by researchers, this invention chooses to use the TopDom model to calculate the TAD boundary, obtaining 44,824 positive examples. Non-boundary regions of the same data and length on the same chromosome are randomly selected as negative examples.

[0127] 2. DNA Sequence Data Preprocessing: Based on the TAD boundary tag data obtained in the previous step, construct a BED format TAD boundary tag data file. Use the BEDtools tool to extract DNA sequence fragments from specified regions of the BED file. This invention uses a 10kb resolution, therefore the DNA sequence length corresponding to each TAD boundary site is 10,000.

[0128] 3. Based on the TAD boundary label data file, extract the signal intensity in specific regions upstream and downstream of each boundary label location as the epigenetic modification signal features of the boundary.

[0129] 4. Hi-C interaction data preprocessing: Extract the interaction frequency data within 500kb upstream and downstream of the TAD boundary site obtained in the previous step. This invention uses a 10kb resolution, and the final Hi-C interaction frequency feature dimension for each sample is [1, 100].

[0130] 5. Dataset Cleaning: Samples with an N value greater than 300 were removed from the constructed dataset to reduce the negative impact of excessive N values ​​on training performance. The final dataset size was 89,244 samples, including 44,449 positive samples and 44,795 negative samples, with a positive-to-negative sample ratio of approximately 1:1. The dataset was then randomly divided into training and test sets at an 8:2 ratio and saved in HDF5 format.

[0131] 6. DNA sequence data encoding: K-mer one-hot encoding and embedding operations are performed on DNA sequence data to obtain DNA embedding codes, which reflect the relationship between the current sequence and the structural characteristics of its location.

[0132] In a DNA sequence of length n, S = (s1, s2, ..., snc) nK-mer encoding is performed, where a sliding window moves from the beginning to the end of the sequence, extracting a sequence of length k each time, resulting in (n-k+1) k-mer subsequences. In this invention, k is set to 3 and n to 10000, therefore the number of subsequences obtained is (10000-3+1).

[0133] One-hot encoding is a fundamental method for converting DNA sequences into feature vectors. For each feature, if it has m possible values, then after one-hot encoding, it becomes m binary features, where each bit has only two values: 0 or 1, as shown in the formula below:

[0134]

[0135] Because there are 4 3 With +1 possible arrangements of sequences within a window, the k-mer encoding result is one-hot encoded. After concatenating the sequences according to the window order and padding with zeros at the end to align the lengths, the one-hot encoded matrix can be obtained.

[0136] To incorporate semantic relationships between subsequences, an embedding layer is introduced after one-hot encoding. An embedding matrix is ​​used to convert the one-hot encoded vector into a low-dimensional vector. Embedding refers to learning a mapping f that transforms the k-mer encoding into another vector representation.

[0137] embed = f(k_mer)

[0138] Input the mapped embedded vector into the model. Specifically, use the embedding matrix. Multiplying the feature matrix by the one-hot encoding matrix reduces the feature dimension, yielding the DNA embedding encoding matrix D∈R. 10000×d_embed Where d_embed represents the embedding dimension, which is set to 10 in this invention. The embedding matrix is ​​updated during model training.

[0139] III. Model Construction

[0140] 1. DNA Sequence Feature Extraction: The DNA embedding encoding matrix obtained from the model preprocessing is passed through a feature extraction network to obtain DNA sequence features. The input feature vector dimension is [4]. 3 +1, 10000]. Specifically, the DNA feature extraction network consists of three layers. The first and second layers are convolutional neural network layers, each containing two sub-layers: a one-dimensional convolutional sub-layer and a pooling sub-layer. The features are input into the one-dimensional convolutional sub-layer, and the calculation formula is as follows:

[0141]

[0142] In this context, the superscript "1" indicates a feature related to the DNA sequence data; This is the output of the one-dimensional convolutional sublayer in this layer; x 1 The input features for the convolutional neural network layers are represented in layer 1 as the epigenetic modification feature data obtained from S22, and in layer 2 as the output result of layer 1; W 1 ∈R out_c_1×k_size_1 The convolution kernel parameters are automatically updated through network training, where out_c_1 is the channel dimension of the output features, and k_size_1 is the convolution kernel size parameter; "*" indicates a convolution operation; b 1 This is a one-dimensional convolutional bias. In this invention, the size of the DNA feature extraction network convolutional kernel is in. channel ×k_size_1, where in channel The dimension of the input features is k_size_1 = 3; padding = same; the number of convolutional kernels in the first layer is 64, and the number of convolutional kernels in the second layer is 32.

[0143] Each one-dimensional convolutional sublayer is followed by a batch normalization process, calculated using an activation function, as follows:

[0144]

[0145] in, This represents the output after passing through the activation function; ReLU indicates activation via the ReLU activation function; BatchNorm is the implementation of batch normalization.

[0146] The pooling sublayer implements one-dimensional max pooling, which reduces the feature length and preserves effective features within the window. The calculation formula is as follows:

[0147]

[0148] Among them, maxpool is the implementation of one-dimensional max pooling functionality. This is the output of the pooling sublayer. The first dimension is always the same as the dimension of the input data to the pooling layer, the other dimension is 4, and the stride is 4. After the first layer, the output tensor length is 2500, and after the second layer, the output tensor length is 625.

[0149] The output of the first pooling sub-layer is input into the second layer, and then output after passing through the same one-dimensional convolutional sub-layer and pooling sub-layer in sequence.

[0150] The third layer contains two sublayers: a BiLSTM sublayer and a Dropout sublayer. The output of the above process is input into the BiLSTM sublayer to learn long-range relationships in the DNA sequence, using the following calculation formula:

[0151]

[0152] Among them, BiLSTM is the implementation of the BiLSTM sublayer function, which adopts a bidirectional long short-term memory network; This is the output of the second layer. The hidden layer parameter dimension is 4, resulting in an output feature dimension of [8, 625].

[0153] After obtaining x lstm Then through the Dropout sublayer:

[0154]

[0155] Where dropout is the dropout function, which randomly deactivates nodes in the network to reduce overfitting; The dropout probability is used to make the value of a node drop out with probability. It becomes 0. This invention takes... It is 0.25.

[0156] 2. Epigenetic Modification Feature Extraction: Epigenetic modification feature data is input into a deep neural network for epigenetic modification feature extraction to extract epigenetic modification features. The purpose is to learn the enrichment or absence of eight histone modification signals and CTCF binding factors at 50 sites upstream and downstream of the current chromatin fragment. The input feature vector dimension is [9, 100]. Specifically, the epigenetic modification feature extraction network consists of three layers, with the first and second layers being convolutional neural network layers. The features are input into a one-dimensional convolutional neural network, and the calculation formula is as follows:

[0157]

[0158] In this context, the superscript "2" indicates a characteristic related to epigenetic modification data; This is the output of the one-dimensional convolutional neural network in this layer; x 1 The input features for the convolutional neural network layers represent epigenetic modification features in layer 1 and the output of layer 1 in layer 2; W 2 ∈R out_c_2×k_size_2 The parameters for the convolution kernels are automatically updated during network training. `out_c_2` represents the channel dimension of the output features, i.e., the number of convolution kernels, and `k_size_2` is the kernel size parameter. "*" indicates a convolution operation. 2 This is a one-dimensional convolutional bias. In this invention, the size of the convolutional kernel of the epigenetic modification feature extraction network is in. channel ×k_size_2, where in channel The dimension of the input features is k_size_2 = 1; padding = same; the number of convolutional kernels in the first layer is 32, and the number of convolutional kernels in the second layer is 9.

[0159] Each one-dimensional convolutional network includes a batch normalization process, calculated using an activation function, as follows:

[0160]

[0161] in The output is the result of batch normalization after passing through the activation function.

[0162] The output of the first layer is input into the second layer, and then passed through the same convolutional neural network and batch normalized before being output.

[0163] The third layer contains two sublayers: a BiLSTM sublayer and a Dropout sublayer. The output of the above process is input into the BiLSTM sublayer to learn the long-range relationships of epigenetic modification data, using the following calculation formula:

[0164]

[0165] in, This is the output of the second layer. The hidden layer parameter dimension is set to 9, resulting in an output feature dimension of [18, 100].

[0166] In obtaining Then through the Dropout sublayer:

[0167]

[0168] in, The results of epigenetic modification feature extraction; The dropout probability is used to make the value of a node drop out with probability. It becomes 0. This invention takes... It is 0.25.

[0169] 3. Hi-C Data Feature Extraction: Hi-C data is input into a Hi-C feature extraction deep neural network to extract Hi-C signal features. The purpose is to learn the interaction frequency features between the current chromatin fragment and 50 upstream and downstream sites. The input feature dimension is [1, 100]. Specifically, the Hi-C feature extraction network consists of 4 layers. The first layer is batch normalization layer 1, which performs batch normalization on the Hi-C feature vector. The calculation formula is as follows:

[0170]

[0171] The superscript "3" indicates a feature related to Hi-C data; x 3 Data representing epigenetic modification characteristics; This represents the output result of batch normalization layer one.

[0172] The second layer is a BiLSTM layer, which learns the long-distance relationships between interaction frequency data. The calculation formula is as follows:

[0173]

[0174] in, This is the output of the BiLSTM layer.

[0175] Layer 3 is batch normalization layer 2, for Batch normalization is performed using the following formula:

[0176]

[0177] in, This represents the output of the second batch normalization layer. The output feature dimension is [2, 100].

[0178] The fourth layer is the Dropout layer: it uses random deactivation to reduce overfitting in the network, and the calculation formula is as follows:

[0179]

[0180] in, The results of feature extraction from Hi-C data; The dropout probability is used to make the value of a node drop out with probability. It becomes 0. This invention takes... It is 0.25.

[0181] IV. Model Prediction

[0182] The obtained DNA sequence features, epigenetic modification features, and Hi-C data features were flattened, and the calculation formulas are as follows;

[0183]

[0184] in, The result is the flattened version of DNA sequence features, epigenetic modification features, and Hi-C data features. "flatten" is the flattening process. These are DNA sequence characteristics, epigenetic modification characteristics, and Hi-C data characteristics.

[0185] The features after stitching and flattening are used to complete multimodal fusion. The calculation formula is as follows:

[0186]

[0187] Where, x fusion This is the result of multimodal fusion; concat indicates vector concatenation.

[0188] The result after fusion x fusion The input consists of m consecutive fully connected layers, where m is set to 4 in this invention. Specifically, each fully connected layer includes a feedforward neural network with a linear transformation, an activation function layer, and a Dropout layer, calculated as follows:

[0189]

[0190] Wherein, the superscript “(i)” represents the relevant features of the i-th fully connected layer, and “(0)” represents the relevant features of the initial input fully connected layer; Let the dropout probability be such that the value of a node in the i-th fully connected layer is increased by a probability. The invention takes 0. Both are 0.25; GELU indicates activation processing via the GELU activation function; The output of the linear transformation in the i-th fully connected layer is calculated as follows:

[0191]

[0192] in, This represents the output of the (i-1)th fully connected layer. This represents the initial input of the fully connected layer, i.e., x. fusion .

[0193] The output of the fully connected layer is passed through an activation function to output the prediction result, calculated as follows:

[0194]

[0195] Where, x predict For each sample, the value is either 0 or 1 to predict the result; sigmoid indicates activation via the sigmoid activation function, as shown in the following formula:

[0196] sigmoid(x) = 1 / (1+e) -x )

[0197] Among them, e -x It represents the reciprocal of the power of the natural base e.

[0198] V. Model Training

[0199] Gradient descent is used to reduce model error; the Adam optimizer is used to accelerate model convergence; and binary cross-entropy loss (BCE) is used for model training, as shown in the following formula:

[0200]

[0201] Among them, y j This represents the predicted value of the j-th sample; Let J represent the label value of the j-th sample; J represents the total number of samples.

[0202] A constant-step decay method is used to dynamically decrease the learning rate as the model iterates, thus mitigating the problem of the model oscillating around the global optimum due to a large learning rate in the later stages of training. The formula is as follows:

[0203]

[0204] Where, α decayed The decay rate is represented by α; the initial learning rate is represented by α; the decay rate is represented by decay_rate; the global step size is represented by global_steps; and the decay step size is represented by decay_steps.

[0205] By employing an early stopping strategy to terminate training in advance, the model can achieve better generalization ability.

[0206] The following experiments will verify the beneficial effects of the present invention.

[0207] The experimental data of this invention are based on five human cell lines and are used to verify the effectiveness of the method of this invention.

[0208] This invention can be abstracted as a binary classification supervised learning problem, and is evaluated using the following metrics:

[0209] 1) Accuracy

[0210] Accuracy is used as a statistical measure to correctly identify samples in a classification problem. It is equal to the proportion of accurate predictions (TP and TN) out of the total number of predicted samples, and is calculated as follows:

[0211]

[0212] The drawback of accuracy evaluation is that it cannot objectively evaluate the performance of a model when the sample classes are imbalanced.

[0213] 2) Precision

[0214] Precision is equal to the proportion of accurately predicted positive samples (TP) out of all samples predicted as positive (TP and FP), and is calculated as follows:

[0215]

[0216] 3) Recall

[0217] Recall is equal to the proportion of accurately predicted positive samples (TP) out of all true positive samples (TP and FN), and is calculated as follows:

[0218]

[0219] 4) F1-Score

[0220] The weighted harmonic mean of precision and recall is calculated as follows:

[0221]

[0222] 5) AUC

[0223] The area under the ROC curve, with the false positive rate (FPR) on the x-axis and the true positive rate (TPR) on the y-axis.

[0224] The accuracy, F1-Score, and AUC value of the method proposed in this invention are compared with those of methods using unfused multimodal data. For all three metrics, higher values ​​indicate better method performance. Figure 2 It is evident that the method of this invention produces the best prediction of TAD boundaries, exhibiting high prediction accuracy. Figure 2 The study confirmed that DNA sequence data, epigenetic modification data, and Hi-C interaction frequency data all contribute to the final prediction results, with Hi-C data contributing the most significantly. Although all three types of data can effectively predict TAD boundaries, the multimodal fusion technology employed in this invention significantly improves the numerical values ​​of the evaluation indicators, demonstrating the effectiveness of the proposed method.

[0225] In addition to comparing the performance of different methods, this invention also verifies whether the TAD boundaries predicted by the method of this invention have biological significance. Figure 3 This shows the predicted enrichment of histone modification signals at the TAD region boundary using the method of this invention. The TAD boundary exhibits characteristic phenomena of specific histone modification signal enrichment or absence. Figure 3 The results showed that histones H3k4me2, H3k4me3, H3k9a, and H3k36me3 were significantly enriched at the TAD boundary, while showing no trend of enrichment or loss at non-TAD boundaries. Histone H3k9me3 showed significant loss at the TAD boundary, but enrichment at non-TAD boundaries. These results are similar to those summarized in existing literature, demonstrating from the perspective of regulatory factors that the boundary sites predicted by the model are closely related to the physical TAD boundary.

[0226] Furthermore, HOMER software was used to perform motif enrichment analysis on the TAD boundaries identified by the method of this invention. A motif refers to a conserved local region in a DNA sequence; it is a small sequence pattern shared by a group of sequences and has a specific biological regulatory function. Figure 4As can be seen, CTCF, BORIS, Sp1, Sox17, HINFP, E2F4, and PGR motifs are significantly enriched in the TAD boundary, and Prop1, ZBTB12, and Znf263 are also enriched. The results of motif enrichment analysis are also highly correlated with gene regulation and chromatin topological constraints, indicating that the TAD boundary predicted by the multimodal fusion model has the corresponding biological function, thus proving the effectiveness of the method.

[0227] In summary, this invention presents a method and system for predicting chromatin topological association domain (TAD) boundaries based on multimodal fusion, achieving improved performance for TAD boundary prediction. The significance of this invention lies in providing researchers in the biomedical field with a novel method and system for identifying TAD boundaries, thereby improving prediction accuracy. This invention has significant value in elucidating cellular processes such as gene expression and regulation, and contributes to further research in disease mechanisms and precision medicine.

Claims

1. A method for predicting the boundary of chromatin topological association structural domains based on multimodal fusion, characterized in that, Includes the following steps: S1: Data Acquisition: Acquire multi-omics datasets, including Hi-C dataset, epigenetic modification dataset, and DNA sequence dataset, and obtain TAD boundary label data through models; S2: Data preprocessing: Based on the TAD boundary label data, feature pre-extraction is performed on the multi-omics dataset, dataset cleaning is completed, and the DNA dataset is encoded and embedded. S3: Model Building: Construct a deep neural network model and extract multimodal features, including DNA sequence features, epigenetic modification features, and Hi-C data features; S4: Model prediction: Flatten and splice the DNA sequence features, epigenetic modification features and Hi-C data features to complete multimodal feature fusion, and predict the boundary results based on the fusion results; S5: Model Training: Train the model on the training set using the backpropagation algorithm to obtain a model for TAD boundary prediction; Step S2 specifically involves: S21: Construct a TAD boundary label data file based on the TAD boundary label data, extract DNA sequence fragments of a specified region and resolution from the file, and construct a DNA sequence dataset; S22: Based on the TAD boundary label data file, extract the signal intensity in specific regions upstream and downstream of each boundary label location as the epigenetic modification features of the boundary. S23: Based on the Hi-C data, extract the interaction frequency data in specific areas upstream and downstream of each boundary label position to obtain the Hi-C interaction frequency feature vector for each sample; S24: Remove invalid samples from the DNA sequence dataset; S25: Perform k-mer encoding, one-hot encoding, and embedding operations on DNA sequence data to obtain a DNA embedding coding matrix, which reflects the relationship between the current sequence and the structural characteristics of its location; The specific process of obtaining the DNA embedding coding matrix in S25 is as follows: S251: For DNA sequences of length n Perform k-mer encoding, which involves moving a sliding window from the beginning to the end of the sequence and extracting a subsequence of length k each time; S252: Perform one-hot encoding on the k-mer subsequence, concatenate them in window order, and pad with zeros at the end to align the lengths to obtain the one-hot encoded matrix. ; S253: An embedding layer is introduced after the one-hot encoding matrix, using the embedding matrix. With one-hot encoding matrix Multiply to reduce feature dimensionality and obtain the DNA embedding coding matrix. ,in, Indicates the embedding dimension.

2. The chromatin topological association domain boundary prediction method based on multimodal fusion according to claim 1, characterized in that, S1 specifically includes: S11: Download Hi-C data, epigenetic modification data, and DNA sequence data for specific types of human cell lines; S12: Apply the TopDom model to the above Hi-C data to obtain positive sample label data of TAD boundary; obtain negative sample label data of TAD boundary by randomly selecting the same number of non-boundary sites as the positive sample label data.

3. The chromatin topological association domain boundary prediction method based on multimodal fusion according to claim 1, characterized in that, The multimodal feature extraction step S3 includes: S31: DNA sequence feature extraction: Input the DNA embedding coding matrix into a deep neural network for DNA feature extraction to extract DNA sequence features; S32: Epigenetic modification feature extraction: Input the epigenetic modification feature data into the deep neural network for epigenetic modification feature extraction to extract the epigenetic modification features, which are used to learn the enrichment or absence of eight histone modification signals and CTCF binding factors at several sites upstream and downstream of the current chromatin fragment. S33: Hi-C Data Feature Extraction: Input Hi-C data into a Hi-C feature extraction deep neural network to extract Hi-C signal features, which are used to learn the interaction frequency features between the current chromatin fragment and several upstream and downstream sites.

4. The chromatin topological association domain boundary prediction method based on multimodal fusion according to claim 3, characterized in that, The deep neural network for DNA feature extraction in S31 consists of three layers; A1: Layers 1 and 2 are convolutional neural network layers, each containing two sub-layers: a one-dimensional convolutional sub-layer and a pooling sub-layer. Features are input into the one-dimensional convolutional sub-layer, and the calculation formula is as follows: ； In this context, the superscript "1" indicates a feature related to the DNA sequence data; This is the output of the one-dimensional convolutional sublayer in this layer; The input features for the convolutional neural network layers represent epigenetic modification feature data in the first layer and the output results of the first layer in the second layer. The convolution kernel parameters are automatically updated and learned during network training, where... The channel dimension of the output feature. This is the kernel size parameter; "*" indicates the convolution operation; For one-dimensional convolution bias; Each one-dimensional convolutional sublayer is followed by a batch normalization process, calculated using an activation function, as follows: ； in, This represents the output after passing through the activation function; ReLU indicates activation via the ReLU activation function. This is for batch normalization. The pooling sublayer implements one-dimensional max pooling to reduce feature length and preserve effective features within the window. The calculation formula is as follows: ； in, This is for implementing one-dimensional max pooling functionality. This is the output of the pooling sublayer; A2: The output of the first pooling sub-layer is input into the second layer, and then output after passing through the same one-dimensional convolutional sub-layer and pooling sub-layer in sequence; A3: The third layer includes a BiLSTM sublayer and a Dropout sublayer; the output of the second layer is input into the BiLSTM sublayer to learn long-range relationships in the DNA sequence, and the calculation formula is as follows: ； in, for The sub-layer functionality is implemented using a bidirectional long short-term memory network. This is the output of the second layer; The output of the BiLSTM sublayer; Obtain the output of the BiLSTM sublayer Then through the Dropout sublayer: ； in, This is the dropout function, used to randomly deactivate nodes in the network to reduce overfitting; For dropout probability, This is the result of DNA feature extraction.

5. The chromatin topological association domain boundary prediction method based on multimodal fusion according to claim 3, characterized in that, The deep neural network for extracting epigenetic modification features in S32 consists of 3 layers. composition: B1: Layers 1 and 2 are convolutional neural network layers; The features are input into a one-dimensional convolutional neural network, and the calculation formula is as follows: ； The superscript "2" indicates a characteristic related to epigenetic modification data; This is the output of the one-dimensional convolutional neural network in this layer; The input features for the convolutional neural network layers represent epigenetic modification feature data in the first layer and the output results of the first layer in the second layer. The convolution kernel parameters are automatically updated and learned during network training, where... _2 represents the channel dimension of the output features, i.e., the number of convolutional kernels. This is the kernel size parameter; "*" indicates the convolution operation; For one-dimensional convolution bias; Each one-dimensional convolutional network includes a batch normalization process, calculated using an activation function, as follows: ； in, The output of the batch normalization result after passing through the activation function; B2: The output of layer 1 is input into layer 2, and after passing through the same convolutional neural network and batch normalization, it is output. B3: Layer 3 includes a BiLSTM sublayer and a Dropout sublayer. The output of layer 2 is input into the BiLSTM sublayer to learn the long-range relationships of epigenetic modification data. The calculation formula is as follows: ； in, This is the output of the second layer; The output of the BiLSTM sublayer; Obtain the output of the BiLSTM sublayer Then through the Dropout sublayer: ； in, The results of epigenetic modification feature extraction; This represents the dropout probability.

6. The chromatin topological association domain boundary prediction method based on multimodal fusion according to claim 3, characterized in that, The Hi-C feature extraction deep neural network in S33 consists of 4 layers: C1: Batch Normalization Layer 1: Batch normalization is performed on the Hi-C interaction frequency feature vector. The calculation formula is as follows: ； The superscript "3" indicates a feature related to Hi-C data; Data representing epigenetic modification characteristics; This represents the output result of batch normalization layer one; C2: BiLSTM layer: Learns the long-distance relationships between interaction frequency data, calculated using the following formula: ； in, The output of the BiLSTM layer; C3: Batch Normalization Layer 2: Results of the BiLSTM Layer Output Batch normalization is performed using the following formula: ； in, This represents the output of the second batch normalization layer. C4: Dropout layer: Uses random deactivation to reduce network overfitting. The calculation formula is as follows: ； in, The results of feature extraction from Hi-C data; This represents the dropout probability.

7. The chromatin topological association domain boundary prediction method based on multimodal fusion according to claim 1, characterized in that, Step S4 specifically involves: S41: Flatten the DNA sequence features, epigenetic modification features, and Hi-C interaction frequency data obtained in S3. The calculation formula is as follows: ； in, The results are after flattening DNA sequence features, epigenetic modification features, and Hi-C interaction frequency data. To achieve the flattening process; The extracted DNA sequence features, epigenetic modification features, and Hi-C interaction frequency data features; S42: Features after stitching and flattening are used to complete multimodal fusion. The calculation formula is as follows: ； in, This is the result after multimodal fusion; This indicates vector concatenation; S43: The fused result The input consists of m consecutive fully connected layers; each fully connected layer includes a feedforward neural network with a linear transformation, an activation function layer, and a Dropout layer, calculated as follows: ； Wherein, the superscript "(i)" represents the relevant features of the i-th fully connected layer, and "(0)" represents the relevant features of the initial input fully connected layer; Let the dropout probability be such that the value of a node in the i-th fully connected layer is increased by a probability. It becomes 0; GELU indicates activation processing via the GELU activation function; The output of the linear transformation in the i-th fully connected layer is calculated as follows: ； in, This represents the output of the (i-1)th fully connected layer. This represents the initial input of the fully connected layer, i.e. ; This represents the weight parameters of the i-th fully connected layer; This represents the bias parameters of the i-th fully connected layer; , All are learned automatically by the model; S44: Pass the output of the fully connected layer through an activation function to output the prediction result. The calculation formula is as follows: ； in, For each sample, the predicted result is either 0 or 1. This indicates activation via the sigmoid function, as shown in the following formula: ； in, It represents the reciprocal of the power of the natural base e.

8. The method for predicting chromatin topological association domain boundaries based on multimodal fusion according to claim 1, characterized in that, Step S5 specifically involves: S51: Use gradient descent to reduce model error; employ the Adam optimizer to accelerate model convergence; use the binary cross-entropy loss function for model training, as shown in the following formula: ； in, This represents the predicted value of the j-th sample; This represents the label value of the j-th sample; J represents the total number of samples. S52: The learning rate is dynamically reduced as the model iterates using a constant step decay method, as shown in the following formula: ； in, This represents the decayed learning rate; Indicates the initial learning rate; Indicates the attenuation coefficient; Indicates the global step size; Indicates the decay step size; S53: Employ an early stop strategy to end training prematurely.

Images