Genome regulatory element analysis method based on deep learning

Abstract

The invention relates to a genome regulatory element analysis method based on deep learning. The method comprises the following steps: acquiring a genome DNA sequence and DNA shape information; on the basis of a shared hybrid neural network, transcription factor binding site features are extracted from a genome sequence, and binding specificity prediction is carried out aiming at the sequence level; aiming at the specificity of transcription factor binding, accurately positioning a transcription factor binding site on a nucleotide level by adopting a full convolutional neural network coupled with different models, and identifying a motif; for a chromosome three-dimensional genome, a chromatin ring is recognized in the three-dimensional whole genome based on a deep learning model, the relationship between the chromatin ring and regulatory elements in all cell lines is analyzed, and the chromatin ring represents cell type specificity and co-enrichment of transcription factor motifs among different cell lines and species. Compared with the prior art, the method has the advantages that the regulatory elements of the regulatory mechanism can be identified and analyzed by utilizing multiple information sources, and the accuracy of predicting the whole genome chromatin ring is effectively improved.

Description

Technical Field

[0001] The present invention relates to the field of gene transcription regulation, and in particular to a method for analyzing genomic regulatory elements based on deep learning. Background Art

[0002] The normal development and function of an organism's life depends on the expression of all its genes in a specific precise space and time. This precise expression pattern is mainly maintained through the transcriptional regulatory process. Dysregulation of gene expression often leads to disease. The transcriptional regulatory process is mediated by proteins binding to regulatory elements on DNA in a combinatorial manner, where a specific combination of transcription factor binding sites establishes a specific regulatory code. All functional elements in the genome, including genes and regulatory sequences, are identified and annotated as a fundamental challenge in genomics and computational biology. The rapid development of sequencing technology has provided researchers with a large amount of biological genome sequence data for the development of methods to detect transcriptional regulatory elements. Computational methods have made significant progress in identifying and detecting DNA regulatory elements, but since regulatory elements are usually short and variable, it is still difficult to identify and discover them using computational algorithms. Advances in the study of three-dimensional chromatin folding technology have highlighted the role of genome conformation in transcriptional regulation, and progress has also been made in identifying cis-regulatory modules and higher-order structures of regulatory sequences, which is crucial for understanding transcriptional regulation in metazoan genomes. Identifying and characterizing gene regulatory elements has far-reaching significance. It is not only crucial for understanding the molecular mechanisms that control specific gene expression patterns, but also an important step in understanding how gene expression changes under pathological conditions. However, there is currently a lack of an effective genomic regulatory element analysis method to model and analyze the identification, localization, and interaction of regulatory elements in the genome. Summary of the invention

[0003] The purpose of the present invention is to provide a genomic regulatory element analysis method based on deep learning.

[0004] The purpose of the present invention can be achieved by the following technical solutions:

[0005] A method for analyzing genomic regulatory elements based on deep learning, comprising the following steps:

[0006] S1, obtain genomic DNA sequence and DNA shape information;

[0007] S2, based on a shared hybrid neural network, extracts transcription factor binding site features from genomic sequences and predicts binding specificity at the sequence level;

[0008] S3, targeting the specificity of transcription factor binding, a fully convolutional neural network coupled with different models was used to accurately locate transcription factor binding sites and identify motifs at the nucleotide level;

[0009] S4, targeting the chromosome three-dimensional genome, identifies chromatin loops in the three-dimensional whole genome and analyzes their relationship with regulatory elements in all cell lines based on a deep learning model, wherein the chromatin loops characterize cell type specificity and co-enrichment of transcription factor motifs between different cell lines and species.

[0010] The shared hybrid neural network is constructed based on a CNN network and an RNN network, and includes an input layer, an encoding layer, a channel matching module, a feature extraction module and a feature integration module, wherein the input layer takes DNA sequence data and its shape features as input; the encoding layer converts the DNA sequence into a DNA sequence matrix through one-hot encoding, and pre-processes the DNA shape features through a sliding window; the channel matching module processes the shape features to match the number of channels of the hybrid shared neural network; the feature extraction module extracts features based on the DNA sequence matrix given by the encoding layer and the DNA shape information processed by the channel matching module; the feature integration module processes the features extracted by the feature extraction module through a batch normalization layer, and then merges the features through two fully connected layers to obtain the predicted binding specificity and outputs it.

[0011] The S3 comprises the following steps:

[0012] S31, according to the characteristics of the in vivo data set, using a fully convolutional neural network coupled with a global average pooling model to locate the transcription factor binding site and predict the motif, wherein the fully convolutional neural network coupled with a global average pooling model is trained with annotated transcription factor-DNA binding sequence data, and the annotation method is: collecting a position count matrix from a motif database, and using the position count matrix to annotate each nucleotide in the DNA sequence;

[0013] S32, according to the characteristics of the in vitro dataset, uses each nucleotide of the sequence as a putative binding site or non-binding site, and uses the FCN and GRU coupling model to achieve base pair resolution prediction, thereby directly locating the transcription factor binding site on the in vitro dataset. The FCN and GRU coupling model includes two branches, one branch is used to perform nucleotide-level classification tasks, predict the probability of each nucleotide and locate the transcription factor binding site; the other branch is used to perform sequence-level regression tasks to predict the strength of each sequence. This branch extracts the long-term dependency between sequence features and motifs through three convolutional blocks and a GRU layer, and inputs them into a fully connected layer for feature fusion. A fully connected layer with 1 neuron is used to predict the binding strength of each sequence.

[0014] The fully convolutional neural network coupled global average pooling model is a symmetrical architecture consisting of a top-down encoding process and a bottom-up decoding process, where:

[0015] The top-down encoding process includes three convolutional modules and a global average pooling layer. Each convolutional module consists of a convolutional layer, a ReLU layer, a maximum pooling layer and a dropout layer. The global average pooling layer is used to capture the global context information of the transcription factor-DNA binding sequence.

[0016] The bottom-up decoding process symmetrically has four deconvolution modules and four skip lines, where each deconvolution module consists of an upsampling layer, a batch normalization layer, a ReLU layer and a convolution layer, and each skip line represents a summation operation for combining high-level semantic information in the deep layer with detailed position information in the shallow layer.

[0017] The loss function of the fully convolutional neural network coupled with the global average pooling model is the difficult negative sample mining loss for locating transcription factor binding sites and predicting motifs. The calculation process is as follows: 1) Calculate the loss of all positive and negative samples; 2) Sort the losses of negative samples from high to low, and select the highest k losses, where k is determined by a specified ratio; 3) Take the average of the losses of all positive samples and the losses of the selected negative samples, and output the sum of the two to obtain the difficult negative sample mining loss. The calculation process is expressed as:

[0018]

[0019] Where S + It is the output of the last deconvolution block. The subscript pos represents positive samples, neg represents negative samples, ratio represents the ratio of negative samples, Crossentropy represents cross entropy loss, top-k represents the highest k loss operations, and mean represents the average value.

[0020] The output of the fully convolutional neural network coupled with the global average pooling model is a nucleotide-level prediction, which is post-processed to locate transcription factor binding sites and predict motifs:

[0021] Step 1) Locate transcription factor binding sites

[0022] A threshold is preset. If the predicted value of the nucleotide level output by the fully convolutional neural network coupled with the global average pooling model is greater than the threshold, the probability of all nucleotides is converted to 1, otherwise it is converted to 0; the number of label 1s is calculated using a sliding window of a preset short length, where the motif length < preset short length << sequence length; the region containing the largest number of label 1s is selected from the entire sequence as the location of the transcription factor binding site;

[0023] Step 2) Predict motifs

[0024] Repeat the above step 1) to locate the transcription factor binding site and obtain the localization region from the DNA sequence; score each sub-region of the localization region based on the training weight of the first convolutional layer, and select the sub-region with the highest score; align the selected sub-regions to calculate the position frequency matrix; match the position frequency matrix and the motifs in the database based on TOMTOM to achieve motif prediction.

[0025] The FCN and GRU coupling model includes three convolution blocks, one GRU layer and three deconvolution blocks. It takes the normalized DNA sequence and the corresponding position weight matrix as input and includes two stages: encoding stage and decoding stage.

[0026] The encoding stage is used to extract basic features of binding specificity from the input sequence and realize high-level encoding of features. It is based on three convolutional blocks and a GRU layer. Each convolutional block includes a convolutional layer, a ReLU layer, a maximum pooling layer, and a discard layer. The last convolutional block is connected to a GRU layer to capture the long-term dependencies between patterns and the spatial distances between DNA sequences.

[0027] The decoding stage is used to learn the position of TFBS in the input sequence and restore the size of the feature map. It is implemented based on three deconvolution blocks and three skipping layers, where each deconvolution block consists of an upsampling layer, a batch normalization layer, a ReLU layer and a convolution layer. After the skip operation and upsampling, the low-level position information is directly transmitted from the encoder to the high-level semantic information of the same height in the decoder, providing more refined features for locating transcription factor binding sites. A sigmoid layer is connected after the convolution layer to infer the probability of each nucleotide belonging to each category.

[0028] The output of the FCN and GRU coupled model is nucleotide-level predictions, which are post-processed to locate transcription factor binding sites and predict motifs:

[0029] A threshold is preset. If the predicted value of the nucleotide level output by the FCN and GRU coupling model is greater than the threshold, the probability of all nucleotides is converted to 1, otherwise it is converted to 0;

[0030] Use a sliding window of 20 bp to select the region containing the maximum number of labels 1 from the FCN and GRU coupling model output as the localization region; use the trained weights of the first convolutional layer to score each subregion of the localization region and select the subregion with the highest score; align and calculate the selected subregions to obtain the position frequency matrix; convert the position frequency matrix into the position weight matrix according to the following equation:

[0031]

[0032] Among them, fb,i represents the frequency of base b at position i, p b is the background frequency of base b in the genome;

[0033] The learned position weight matrix representing the predicted motif is compared with the known motifs in the standard database through TOMTOM to achieve motif prediction.

[0034] The input of the deep learning model includes Hi-C matrix, interactions involving target regions or proteins of interest, and corresponding accessible chromatin landscape data. A three-layer CNN model is applied to extract features from the Hi-C matrix and the accessible chromatin landscape matrix, and outputs the identified whole-genome chromatin loops. The deep learning model consists of three convolutional blocks, each of which consists of a convolutional layer, a ReLU layer, a dropout layer, and a global average pooling layer. The convolutional layer is used to directly capture local features, the global average pooling layer is used to capture global structural information from the Hi-C matrix and the chromatin landscape matrix, and the dropout layer is used to avoid overfitting and reduce complex adaptation relationships between neurons. Two fully connected layers of 64 neurons are used to fuse the features of the Hi-C matrix and the chromatin landscape matrix output by the last convolutional block, and the output of the fully connected layer is connected to a sigmoid layer to output the probability of the candidate chromatin loop.

[0035] The best deep learning model obtained through training is used to predict all potential chromatin loops in the corresponding chromosomes, which includes two stages. In the first stage, the best deep learning model is used to score all non-zero pixels, and the non-zero pixels indicate that there are potential chromatin loops on each chromosome. In the second stage, the candidate chromatin loops are merged, and only the candidate chromatin loops with a contact frequency greater than the average of all candidate chromatin loops are retained. A greedy algorithm is used to cluster all candidate chromatin loops, and the most representative pixels are selected as the identified chromatin loops.

[0036] Compared with the prior art, the present invention has the following beneficial effects:

[0037] (1) The present invention proposes a hybrid shared neural network model, which aims to improve the recognition accuracy of transcription factor binding sites by combining DNA sequence and DNA shape features. The innovation of this method lies in: 1) using a shared hybrid convolutional neural network and recurrent neural network, CNN is used to learn spatial information from a given DNA sequence and shape features, and RNN mainly captures the long-term dependencies between sequence motifs, which can effectively extract features from large-scale genome sequences obtained by high-throughput technology; 2) common patterns are found from DNA sequences and their corresponding DNA shape features, and then they are connected to calculate the predicted specificity; 3) the local structure of the DNA sequence can be captured without relying entirely on DNA shape features. 4) The identified transcription factor binding site motifs are highly matched with the verified motifs in the standard library.

[0038] (2) The present invention proposes a fully convolutional pooling network model and a fully convolutional gated recurrent unit network model, aiming to directly locate the binding sites of transcription factors. The innovation of this method lies in: 1) overcoming the limitations of most computational methods and deep learning methods aimed at predicting the specificity of transcription factor and DNA binding sequences, which is equivalent to converting the sequence-level regression task into a nucleotide-level classification problem; 2) in view of the characteristics of the human in vivo dataset, a model combining a fully convolutional module with a global average pooling is constructed, which can roughly locate TFBS and accurately identify motifs; 3) in view of the characteristics of the human in vitro dataset, a model combining a fully convolutional neural network with a gated recurrent unit is constructed. This method is a dual-task framework: one is a nucleotide-level classification task, predicting the probability of each nucleotide and locating the transcription factor binding site; the other is a sequence-level regression task, predicting the binding strength of each sequence. The dual-task framework has shown great advantages in both the nucleotide-level classification task and the sequence-level regression task.

[0039] (3) The present invention proposes a three-dimensional genome regulatory element analysis model based on convolutional neural networks, aiming to further analyze regulatory elements in high-dimensional spatial structures. The innovation of this method is: 1) by combining one-dimensional chromatin landscape data and three-dimensional chromatin conformation maps, it makes up for the deficiency that most current methods for analyzing regulatory elements ignore the real chromatin spatial structure information; 2) using some orthogonal data to generate positive samples with a wider contact matrix, it provides the possibility of finding more potential whole-genome chromatin loops; 3) using deep learning models to fully extract the characteristics of chromatin landscape and three-dimensional chromatin conformation maps, so as to fully mine chromatin loops in the whole genome and analyze the interaction of regulatory elements. This method effectively improves the accuracy of predicting whole-genome chromatin loops, and can identify some unique chromatin loops. In addition, the regulatory elements in the identified chromatin loops further reveal the phenomenon of cell type specificity and co-enrichment of transcription factor motifs between different cell lines and species. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 is a flow chart of the method of the present invention;

[0041] Figure 2 It is a schematic diagram of the structure of a shared hybrid neural network;

[0042] Figure 3 A performance comparison result of a shared hybrid neural network and the prior art in an embodiment;

[0043] Figure 4 This is a schematic diagram of the structure of a fully convolutional neural network coupled with a global average pooling model;

[0044] Figure 5 This is a schematic diagram of the workflow of the FCN and GRU coupling model;

[0045] Figure 6 This is the workflow diagram of the deep learning model;

[0046] Figure 7 Performance comparison chart with existing methods when using different parameters for deep learning models;

[0047] Figure 8 The figure shows the overall performance comparison of deep learning models and existing methods on different cell lines. DETAILED DESCRIPTION

[0048] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. This embodiment is implemented based on the technical solution of the present invention, and provides a detailed implementation method and specific operation process, but the protection scope of the present invention is not limited to the following embodiments.

[0049] This embodiment provides a method for analyzing genomic regulatory elements based on deep learning, such as Figure 1 As shown, the following steps are included:

[0050] S1, obtain genomic DNA sequence and DNA shape information.

[0051] S2, based on a shared hybrid neural network, extracts transcription factor binding site features from genomic sequences and predicts binding specificity at the sequence level.

[0052] Considering that the current shared CNN model that combines DNA shape features and DNA sequences has achieved certain results, there are still some defects that need to be improved, such as data bias, extracting local dependencies between motifs, and visualizing the binding motifs of transcription factors (TFs). Therefore, this step aims to combine DNA sequences and their corresponding local DNA shape features to improve the accuracy of predicting transcription factor binding sites, and proposes a new hybrid deep neural network architecture that shares a combination of CNN and RNN. In this model framework, CNN is used to capture low-level spatial information from a given DNA sequence and DNA shape features, and RNN is used to capture long-term dependencies between sequences. The framework mainly consists of data input, model training, and predicted binding strength output. The architecture diagram is shown in the figure. Figure 2 As shown, it includes an input layer, an encoding layer, a channel matching module, a feature extraction module and a feature integration module.

[0053] The input layer takes DNA sequence data and its shape features as input, where the DNA shape features are obtained from the nucleotide sequence by Monte Carlo simulation method, and the sliding window method is used to obtain various structural features of DNA (such as minor groove width, rolling, propeller twist and helical twist).

[0054] The encoding layer converts the DNA sequence into a DNA sequence matrix through one-hot encoding and preprocesses the DNA shape features through a sliding window.

[0055] The channel matching module processes the shape features to match the number of channels of the hybrid shared neural network. This module consists of a convolutional layer and an activation layer ReLU (Rectified Linear Unit). Since DNA sequences and DNA shape features are heterogeneous data, a reasonable integration strategy should be adopted to fuse two different types of features. In order to facilitate the integration of more shape features in future work, the same general channel matching module as the DLBSS model is used, that is, a convolutional layer with a convolution kernel size of 1 is used to process DNA shape features to match the number of channels of the hybrid network. More specifically, the model uses a convolutional layer to positionally correspond the preprocessed shape feature matrix 4*l with the DNA sequence to ensure that each position corresponds to 4 feature vectors. For each convolutional network layer, the output of the layer is calculated by the following formula:

[0056]

[0057] Where X is the input, i is the index of the output position, k is the index of the convolution kernel, W is the convolution weight tensor, which can be understood as a DNA shape feature motif detector; b is the bias term. And the convolution output is converted by the activation function, in which the element-by-element nonlinear function ReLU is used. ReLU is an activation function widely used in deep learning, which can alleviate the gradient vanishing problem in the back-propagation training process and has better convergence performance. ReLU is defined as follows:

[0058]

[0059] The feature extraction module performs feature extraction based on the DNA sequence matrix given by the encoding layer and the DNA shape information processed by the channel matching module. The module consists of a convolutional layer, a ReLU layer, a maximum pooling layer, a LSTM layer, and a dropout layer. The DNA shape information processed by the channel matching module and the original DNA sequence data are input into the shared hybrid neural network to extract features for predicting TFBS. The use of a shared hybrid model not only greatly reduces the parameters of the network but also allows for parallel training. First, a CNN layer with 16 convolution kernels + ReLU layer is designed to capture spatial information from the given DNA sequence and DNA shape features, which can not only easily integrate multiple types of input data and features, but also automatically discover effective information through representation learning. Then, a maximum pooling layer is used to select the maximum value based on the output of the convolutional layer, thereby reducing the dimension of the input and making the model more computationally efficient. The pooling operation is defined as:

[0060] pooling(X) ik =max({X iM+1,k ,X (iM+1),k ,…,X (iM+M-1),k})

[0061] Then, the output of the maximum pooling layer is input to the LSTM layer to capture the directional and spatial long-term dependencies between sequences. Secondly, a dropout layer with a probability of 0.2 is added to avoid overfitting by ignoring half of the feature detectors, making the model more generalizable. Next, the DNA sequence and DNA shape information result groups output by the dropout layer are concatenated into a feature vector and input to the output stage.

[0062] The feature integration module consists of two fully connected layers, a batch normalization layer, and a dropout layer. At the output stage, batch normalization is applied before it is input into the fully connected layer, which not only avoids the gradient problem during back-propagation, but also simplifies the initialization process of network parameters. The output of the batch normalization layer is then input into a fully connected layer with 32 hidden neurons to integrate features, passes through a dropout layer, and finally passes through a fully connected layer containing only one neuron for predicting TF-DNA binding specificity.

[0063] For each dataset, the model uses the mean squared error (MSE) loss function and minimizes it to train the proposed hybrid model. This example conducted a series of experiments on 66 in vitro datasets and evaluated the performance of the model using the R 2 The performance of the proposed model and competing methods were compared with PCC on 66 in vitro datasets. Figure 3 As shown, it is shown that the performance of the present invention is significantly better than that of the competing methods.

[0064] S3, targeting the specificity of transcription factor binding, a fully convolutional neural network coupled with different models was used to accurately locate transcription factor binding sites (TFBS) and identify motifs at the nucleotide level.

[0065] S3 includes the following steps:

[0066] S31, according to the characteristics of the in vivo dataset, a fully convolutional neural network coupled with a global average pooling model is used to locate transcription factor binding sites and predict motifs, wherein the fully convolutional neural network coupled with a global average pooling model is trained with annotated transcription factor-DNA binding sequence data, and the annotation method is as follows: a position count matrix is ​​collected from the motif database, and each nucleotide in the DNA sequence is annotated using the position count matrix.

[0067] like Figure 4 As shown in Figure 2, the fully convolutional neural network coupled global average pooling model is a symmetric architecture consisting of a top-down encoding process (such as Figure 4 Left part) and bottom-up decoding process (such as Figure 4 right part).

[0068] The top-down encoding process includes three convolution modules and a global average pooling layer. Each convolution module consists of a convolution layer, a ReLU layer, a maximum pooling layer, and a dropout layer. The calculation process is as follows:

[0069]

[0070] Where S is the DNA sequence matrix encoded by the one-hot method, X represents the convolution of the convolution kernel M and the matrix S, and ⊙ represents the convolution operation. The convolution layer is usually regarded as a motif scanner to score each fragment of the sequence. The maximum pooling layer is used to reduce the computational complexity and select the best response of the local neighboring area. The dropout layer is often used to alleviate the overfitting problem.

[0071] In the application of image segmentation, several studies have shown that global average pooling can capture the global context information of the image. Similarly, considering the global context of the DNA sequence is also important for motif discovery, for example, the regional information on both sides of the peak can affect the binding activity of TF-DNA. Therefore, a global average pooling layer is added on top of the last convolutional block to capture the global context information of the TF-DNA binding sequence.

[0072] The bottom-up decoding process symmetrically has four deconvolution modules and four skip lines, where each deconvolution module consists of an upsampling layer, a batch normalization (BN) layer, a ReLU layer and a convolution layer. Each skip line represents a summation operation. The calculation process is as follows:

[0073]

[0074] Where Z is the output of the same-level encoding process, and ⊙ represents the convolution operation. The upsampling layer is used to restore the size of the downsampled feature map. The skip line is used to combine the high-level semantic information of the deep layer with the detailed position information of the shallow layer, which helps to accurately predict the label of each position.

[0075] When the categories are unbalanced, useless easy negative samples will deviate from the overall learning direction of the model, resulting in invalid learning effects and inability to correctly distinguish specific target objects. In the field of target detection, some effective methods have been proposed to solve the problem of data imbalance. In particular, the hard negative mining loss method is usually used to process unbalanced positive and negative data. This embodiment designs a hard negative mining loss for locating TFBS and predicting motifs as a loss function of the fully convolutional neural network coupled global average pooling model, and its calculation process is: 1) Calculate the loss of all positive and negative samples; 2) Sort the losses of negative samples from high to low, and select the highest k losses, where k is determined by a specified ratio; 3) Take the average of the losses of all positive samples and the selected negative samples, and output the sum of the two to obtain the hard negative mining loss. The calculation process is expressed as:

[0076]

[0077] Where S + It is the output of the last deconvolution block. The subscript pos represents positive samples, neg represents negative samples, ratio represents the ratio of negative samples, which is set to 0.3 in this embodiment. Crossentropy represents cross entropy loss, top-k represents taking the highest k loss operations, and mean represents taking the average value.

[0078] Since the output of the fully convolutional neural network coupled global average pooling model is a nucleotide-level prediction, it is post-processed to locate transcription factor binding sites and predict motifs:

[0079] Step 1) Locate transcription factor binding sites

[0080] A threshold (such as 0.9) is preset. If the predicted value of the nucleotide level output by the fully convolutional neural network coupled with the global average pooling model is greater than the threshold, the probability of all nucleotides is converted to 1, otherwise it is converted to 0; a short sliding window of 50bp is used to calculate the number of label 1s, where the motif length <50bp<<sequence length; the region containing the largest number of label 1s is selected from the entire sequence as the location of the transcription factor binding site.

[0081] Step 2) Predict motifs

[0082] Repeat the above step 1) to locate the transcription factor binding site and obtain the localization region from the DNA sequence; score each sub-region of the localization region based on the training weight of the first convolutional layer, and select the sub-region with the highest score; align the selected sub-regions to calculate the position frequency matrix; match the position frequency matrix and the motifs in the database based on TOMTOM to achieve motif prediction.

[0083] S32, according to the characteristics of the in vitro dataset, uses each nucleotide of the sequence as a putative binding site or non-binding site, and utilizes the FCN and GRU coupling model based on the combination of FCN and GRU to achieve base pair resolution prediction, thereby directly locating the transcription factor binding site on the in vitro dataset. The FCN and GRU coupling model includes two branches, one branch is used to perform nucleotide level classification tasks, predict the probability of each nucleotide and locate the transcription factor binding site; the other branch is used to perform sequence level regression tasks to predict the strength of each sequence. This branch extracts the long-term dependency between sequence features and motifs through three convolution blocks and one GRU layer, and inputs them into a fully connected layer for feature fusion. A fully connected layer with 1 neuron is used to predict the binding strength of each sequence.

[0084] like Figure 5 The workflow of coupling FCN with GRU includes three steps: data processing, model training and motif query. (i) One-hot encoding is used after data normalization of DNA sequences; (ii) The classification task at the nucleotide level contains three convolutional layers, one GRU layer and three deconvolutional layers; the regression task at the sequence level contains one FC layer, followed by three convolutional layers and a GRU layer; (iii) TOMTOM is used to match the learned PWM with experimentally verified motifs in the standard database.

[0085] The FCN and GRU coupling model includes three convolution blocks, one GRU layer and three deconvolution blocks. It takes the normalized DNA sequence and the corresponding position weight matrix as input, and includes two stages: encoding stage and decoding stage.

[0086] The encoding stage is used to extract basic features of binding specificity from the input sequence and realize high-level encoding of features, which is implemented based on three convolutional blocks and a GRU layer. Each convolutional block includes a convolutional layer, a ReLU layer, a maximum pooling layer, and a dropout layer. The deeper the convolutional layer, the more abstract the TF-DNA features obtained from the DNA sequence. The low-level feature map is obtained after downsampling, which provides contextual information of TFBS in the entire DNA sequence and can be understood as a feature reflecting the relationship between the target and its background. The maximum pooling layer is a dimensionality reduction process that allows the model to pay more attention to global features rather than local positions, so that some important feature information can be retained. The dropout layer can effectively alleviate overfitting. Considering the shortcomings of FCN in image segmentation, such as classifying each pixel without fully considering the relationship between pixels. To make up for this shortcoming, a GRU layer is connected after the last convolutional block to capture the long-term dependencies between patterns and the spatial distances between DNA sequences.

[0087] The decoding stage is used to learn the position of TFBS in the input sequence and restore the size of the feature map. It is implemented based on three deconvolution blocks and three skipping layers, where each deconvolution block consists of an upsampling layer, a batch normalization layer, a ReLU layer and a convolution layer. After the skip operation and upsampling, the low-level position information is directly transmitted from the encoder to the high-level semantic information of the same height in the decoder, providing more refined features for locating transcription factor binding sites. A sigmoid layer is connected after the convolution layer to infer the probability of each nucleotide belonging to each category.

[0088] The output of the FCN coupled with GRU model is nucleotide-level predictions, which are post-processed to locate transcription factor binding sites and predict motifs:

[0089] A threshold is preset. If the predicted value of the nucleotide level output by the FCN and GRU coupling model is greater than the threshold, the probability of all nucleotides is converted to 1, otherwise it is converted to 0;

[0090] Use a sliding window of 20 bp to select the region containing the maximum number of labels 1 from the FCN and GRU coupling model output as the positioning region;

[0091] Use the trained weights of the first convolutional layer to score each sub-region of the localization region and select the sub-region with the highest score;

[0092] Align and calculate the selected sub-regions to obtain the position frequency matrix PFM;

[0093] The position frequency matrix is ​​converted into the position weight matrix PWM according to the following equation:

[0094]

[0095] Among them, f b,i represents the frequency of base b at position i, p b is the background frequency of base b in the genome, assuming that all four bases occur equally, i.e., p b =0.25;

[0096] Finally, the predicted motifs represented by the learned position weight matrix are compared with the known motifs in the standard database through TOMTOM to achieve motif prediction.

[0097] S4, targeting the chromosome three-dimensional genome, identifies chromatin loops in the three-dimensional whole genome based on a deep learning model and analyzes their relationship with regulatory elements in all cell lines. Among them, chromatin loops characterize cell type specificity and co-enrichment of transcription factor motifs between different cell lines and species.

[0098] like Figure 6 As shown in Figure 1, the input of the deep learning model includes Hi-C matrix, some interactions involving target regions or proteins of interest verified by biological enrichment experiments such as ChIA-PET / HiChIP and Capture Hi-C, and the corresponding accessible chromatin landscape data. Positive samples are generated based on the input data, and negative samples are generated based on the distance close to or large from the positive samples. Then, a three-layer CNN model is applied to extract features from the Hi-C matrix and the accessible chromatin landscape matrix, and output the identified genome-wide chromatin loops. Specifically, the deep learning model consists of three convolutional blocks, each of which consists of a convolutional layer, a ReLU layer, a dropout layer and a global average pooling layer. The convolutional layer is used to directly capture local features, the global average pooling layer is used to capture global structural information from the Hi-C matrix and the chromatin landscape matrix, and the dropout layer is used to avoid overfitting and reduce the complex adaptation relationship between neurons, which is set to 0.2; two 64-neuron fully connected layers are used to fuse the features of the Hi-C matrix and the chromatin landscape matrix output by the last convolutional block, and the output of the fully connected layer is connected to a sigmoid layer to output the probability of the candidate chromatin loop.

[0099] The prediction of genome-wide chromatin loops is regarded as a binary classification task, so the binary cross-entropy loss (BCELoss) is used to train the model.

[0100] The best deep learning model obtained through training is used to predict all potential chromatin loops in the corresponding chromosomes, which includes two stages. In the first stage, the best deep learning model is used to score all non-zero pixels, and the non-zero pixels indicate that there are potential chromatin loops on each chromosome. In the second stage, the candidate chromatin loops are merged, and only the candidate chromatin loops with a contact frequency greater than the average of all candidate chromatin loops are retained. A greedy algorithm is used to cluster all candidate chromatin loops, and the most representative pixels are selected as the identified chromatin loops.

[0101] In order to train and verify the performance of the proposed deep learning model, the deep learning model of the present invention was implemented on three human cell lines (GM12878, K562 (replicate sample 1), H1-ESC) and one mouse cell line (mESC) to verify its classification performance. F1-score, PRAUC, Precision and Recall indicators were used to verify the performance of the present invention and existing methods in identifying chromatin loops.

[0102] The model proposed in this embodiment uses a deep learning framework to automatically learn features, rather than the manually designed features used in Peakachu to identify chromatin loops, which is one of the innovations of the present invention. In order to fully extract the features of chromatin loops, a larger window (23*23) is used to generate the feature matrix of positive and negative samples. To better illustrate this issue, the training samples with a window of 11*11 were expanded and experiments were performed on the GM12878 cell line H3K27ac HiChIP and RAD21 ChIA-PET. Figure 7 As shown in (a), it can be seen that although the deep learning model of the present invention uses a window size of 11*11 to generate features for each interaction, its performance is generally better than peakachu. Figure 7 (a) and (b), the experimental results confirm our hypothesis that the deep learning model of the present invention with a larger window (23*23) performs better than Peakachu in identifying chromatin loops. The window size is just a parameter, and a larger window is more suitable for the deep learning model of the present invention to improve the recognition accuracy of chromatin loops. Although Peakachu and the deep learning model of the present invention use different window sizes, it is still relatively fair for experimental purposes. Another innovation of the deep learning model of the present invention is to effectively combine the one-dimensional open chromatin landscape with 3D genome data for chromatin loop prediction, which is not considered by Peakachu. Therefore, a larger window size (23*23) was used in the deep learning model of the present invention for subsequent experiments.

[0103] from Figure 8It can be seen that the average values ​​of F1-score and PRAUC of the deep learning model of the present invention are greater than Peakachu, which indicates that the classification performance of the deep learning model of the present invention is better than Peakachu on all cell lines. Figure 8 As shown, (a) and (b) are experimental results obtained by training five enrichment experimental data sets in the GM12878 cell line. It can be observed that the average F1-score and PRAUC values ​​in CTCF ChIA-PET, H3K27ac HiChIP, SMC1HiChIP, and RAD21 ChIA-PET are close to 0.95 or greater than 0.95, respectively, showing relatively excellent classification performance. It is worth noting that the F1-score of the deep learning model of the present invention is about 8% higher than Peakachu in the K562 (replicate 1) cell line, but compared with other cell lines, the F1-score and PRAUC values ​​of these two methods are relatively low.

[0104] The preferred specific embodiments of the present invention are described in detail above. It should be understood that a person skilled in the art can make many modifications and changes based on the concept of the present invention without creative work. Therefore, any technical solution that can be obtained by a person skilled in the art through logical analysis, reasoning, or limited experiments based on the concept of the present invention on the basis of the prior art should be within the scope of protection determined by the claims.

Claims

1. A method for analyzing genomic regulatory elements based on deep learning, characterized in that: The following steps are involved: S1, obtain genomic DNA sequence and DNA shape information; S2, based on a shared hybrid neural network, extracts transcription factor binding site features from genomic sequences and predicts binding specificity at the sequence level; S3, targeting the specificity of transcription factor binding, a fully convolutional neural network coupled with different models was used to accurately locate transcription factor binding sites and identify motifs at the nucleotide level; S4, targeting the chromosome three-dimensional genome, identifies chromatin loops in the three-dimensional whole genome and analyzes their relationship with regulatory elements in all cell lines based on a deep learning model, wherein the chromatin loops characterize cell type specificity and co-enrichment of transcription factor motifs between different cell lines and species.

2. A method for analyzing genomic regulatory elements based on deep learning according to claim 1, characterized in that: The shared hybrid neural network is constructed based on a CNN network and an RNN network, and includes an input layer, an encoding layer, a channel matching module, a feature extraction module and a feature integration module, wherein the input layer takes DNA sequence data and its shape features as input; the encoding layer converts the DNA sequence into a DNA sequence matrix through one-hot encoding, and pre-processes the DNA shape features through a sliding window; the channel matching module processes the shape features to match the number of channels of the hybrid shared neural network; the feature extraction module extracts features based on the DNA sequence matrix given by the encoding layer and the DNA shape information processed by the channel matching module; the feature integration module processes the features extracted by the feature extraction module through a batch normalization layer, and then merges the features through two fully connected layers to obtain the predicted binding specificity and outputs it.

3. The method for analyzing genomic regulatory elements based on deep learning according to claim 1, characterized in that: The S3 comprises the following steps: S31, according to the characteristics of the in vivo data set, using a fully convolutional neural network coupled with a global average pooling model to locate the transcription factor binding site and predict the motif, wherein the fully convolutional neural network coupled with a global average pooling model is trained with annotated transcription factor-DNA binding sequence data, and the annotation method is: collecting a position count matrix from a motif database, and using the position count matrix to annotate each nucleotide in the DNA sequence; S32, according to the characteristics of the in vitro dataset, uses each nucleotide of the sequence as a putative binding site or non-binding site, and uses the FCN and GRU coupling model to achieve base pair resolution prediction, thereby directly locating the transcription factor binding site on the in vitro dataset. The FCN and GRU coupling model includes two branches, one branch is used to perform nucleotide-level classification tasks, predict the probability of each nucleotide and locate the transcription factor binding site; the other branch is used to perform sequence-level regression tasks to predict the strength of each sequence. This branch extracts the long-term dependency between sequence features and motifs through three convolutional blocks and a GRU layer, and inputs them into a fully connected layer for feature fusion. A fully connected layer with 1 neuron is used to predict the binding strength of each sequence.

4. A method for analyzing genomic regulatory elements based on deep learning according to claim 3, characterized in that: The fully convolutional neural network coupled global average pooling model is a symmetrical architecture consisting of a top-down encoding process and a bottom-up decoding process, where: The top-down encoding process includes three convolutional modules and a global average pooling layer, each convolutional module consists of a convolutional layer, a ReLU layer, a maximum pooling layer and a dropout layer, and the global average pooling layer is used to capture the global context information of the transcription factor-DNA binding sequence; The bottom-up decoding process symmetrically has four deconvolution modules and four skip lines, where each deconvolution module consists of an upsampling layer, a batch normalization layer, a ReLU layer and a convolution layer, and each skip line represents a summation operation for combining high-level semantic information in the deep layer with detailed position information in the shallow layer.

5. The method for analyzing genomic regulatory elements based on deep learning according to claim 3, characterized in that: The loss function of the fully convolutional neural network coupled with the global average pooling model is the difficult negative sample mining loss for locating transcription factor binding sites and predicting motifs. The calculation process is as follows: 1) Calculate the loss of all positive and negative samples; 2) Sort the losses of negative samples from high to low, and select the highest k losses, where k is determined by a specified ratio; 3) Take the average of the losses of all positive samples and the losses of the selected negative samples, and output the sum of the two to obtain the difficult negative sample mining loss. The calculation process is expressed as: Where S + It is the output of the last deconvolution block. The subscript pos represents positive samples, neg represents negative samples, ratio represents the ratio of negative samples, Crossentropy represents cross entropy loss, top-k represents the highest k loss operations, and mean represents the average value.

6. A method for analyzing genomic regulatory elements based on deep learning according to claim 3, characterized in that: The output of the fully convolutional neural network coupled with the global average pooling model is a nucleotide-level prediction, which is post-processed to locate transcription factor binding sites and predict motifs: Step 1) Locate transcription factor binding sites A threshold is preset. If the predicted value of the nucleotide level output by the fully convolutional neural network coupled with the global average pooling model is greater than the threshold, the probability of all nucleotides is converted to 1, otherwise it is converted to 0; the number of label 1s is calculated using a sliding window of a preset short length, where the motif length < preset short length << sequence length; the region containing the largest number of label 1s is selected from the entire sequence as the location of the transcription factor binding site; Step 2) Predict motifs Repeat the above step 1) to locate the transcription factor binding site and obtain the localization region from the DNA sequence; score each sub-region of the localization region based on the training weight of the first convolutional layer, and select the sub-region with the highest score; align the selected sub-regions to calculate the position frequency matrix; match the position frequency matrix and the motifs in the database based on TOMTOM to achieve motif prediction.

7. The method for analyzing genomic regulatory elements based on deep learning according to claim 3, characterized in that: The FCN and GRU coupling model includes three convolution blocks, one GRU layer and three deconvolution blocks. It takes the normalized DNA sequence and the corresponding position weight matrix as input and includes two stages: encoding stage and decoding stage. The encoding stage is used to extract basic features of binding specificity from the input sequence and realize high-level encoding of features. It is based on three convolutional blocks and a GRU layer. Each convolutional block includes a convolutional layer, a ReLU layer, a maximum pooling layer, and a discard layer. The last convolutional block is connected to a GRU layer to capture the long-term dependencies between patterns and the spatial distances between DNA sequences. The decoding stage is used to learn the position of transcription factor binding sites in the input sequence and restore the size of the feature map. It is implemented based on three deconvolution blocks and three skip lines, where each deconvolution block consists of an upsampling layer, a batch normalization layer, a ReLU layer and a convolution layer. After the skip operation and upsampling, the low-level position information is directly transmitted from the encoder to the high-level semantic information of the same height in the decoder, providing more refined features for locating transcription factor binding sites. A sigmoid layer is connected after the convolution layer to infer the probability of each nucleotide belonging to each category.

8. The method for analyzing genomic regulatory elements based on deep learning according to claim 3, characterized in that: The output of the FCN and GRU coupled model is nucleotide-level predictions, which are post-processed to locate transcription factor binding sites and predict motifs: A threshold is preset. If the predicted value of the nucleotide level output by the FCN and GRU coupling model is greater than the threshold, the probability of all nucleotides is converted to 1, otherwise it is converted to 0; Use a sliding window of 20 bp to select the region containing the maximum number of labels 1 from the FCN and GRU coupling model output as the localization region; use the trained weights of the first convolutional layer to score each subregion of the localization region and select the subregion with the highest score; align and calculate the selected subregions to obtain the position frequency matrix; convert the position frequency matrix into the position weight matrix according to the following equation: Among them, f b,i represents the frequency of base b at position i, p b is the background frequency of base b in the genome; The learned position weight matrix representing the predicted motif is compared with the known motifs in the standard database through TOMTOM to achieve motif prediction.

9. The method for analyzing genomic regulatory elements based on deep learning according to claim 1, characterized in that: The input of the deep learning model includes Hi-C matrix, interactions involving target regions or proteins of interest, and corresponding accessible chromatin landscape data, and a three-layer CNN model is applied to extract features from the Hi-C matrix and the accessible chromatin landscape matrix, and output the identified whole-genome chromatin loops; the deep learning model includes three convolutional blocks, each of which consists of a convolutional layer, a ReLU layer, a dropout layer, and a global average pooling layer, the convolutional layer is used to directly capture local features, the global average pooling layer is used to capture global structural information from the Hi-C matrix and the chromatin landscape matrix, and the dropout layer is used to avoid falling into overfitting and reduce the complex adaptation relationship between neurons; Two fully connected layers of 64 neurons are used to fuse the features of the Hi-C matrix and the chromatin landscape matrix output by the last convolutional block. The output of the fully connected layer is connected to a sigmoid layer to output the probability of the candidate chromatin loop.

10. The method for analyzing genomic regulatory elements based on deep learning according to claim 1, characterized in that: The best deep learning model obtained through training is used to predict all potential chromatin loops in the corresponding chromosomes, which includes two stages. In the first stage, the best deep learning model is used to score all non-zero pixels, and the non-zero pixels indicate that there are potential chromatin loops on each chromosome. In the second stage, the candidate chromatin loops are merged, and only the candidate chromatin loops with a contact frequency greater than the average of all candidate chromatin loops are retained. A greedy algorithm is used to cluster all candidate chromatin loops, and the most representative pixels are selected as the identified chromatin loops.

Images