Hierarchical TADs structural domain identification method and system based on interpretable deep learning

Abstract

The invention relates to the technical field of biological computing, in particular to a hierarchical TADs structural domain identification method and system based on interpretable deep learning, and the method comprises the steps: obtaining a Hi-C contact matrix, extracting basic TADs structural domain boundary points in the Hi-C contact matrix, and generating each hierarchy of candidate hierarchical TADs structural domains according to the boundary points and the maximum TADs structural domain limiting condition. And then training a deep learning model by taking the position relationship between any two chromatin fragments in the Hi-C contact matrix and all candidate hierarchical TADs structural domains and the genome distance between the two fragments as input and taking the interaction frequency between the two fragments as output. And then explaining the deep learning model, and screening out the most significant hierarchical TADs structural domain according to the contribution degree of each candidate hierarchical TADs structural domain to the whole contact frequency. The process is helpful to quantify the importance of each TADs structural domain on the whole Hi-C contact matrix, and provides an important basis for genome three-dimensional structure research and the like.

Description

Technical Field

[0001] The present invention relates to the field of biocomputing technology, and in particular to a hierarchical TADs domain identification method and system based on interpretable deep learning. Background Art

[0002] Chromatin of eukaryotic organisms is composed of DNA and proteins. Chromatin has a certain three-dimensional structure. They are tightly folded according to certain rules in the cell nucleus space with a diameter of only a few microns. The three-dimensional spatial structure and folding pattern of chromatin are closely related to key biological functions such as gene regulation, chromosome morphological changes, and genome stability. High-throughput chromosome conformation capture (Hi-C) technology is a sequencing technology that indirectly detects the spatial structure of chromatin in the whole genome. It measures the interaction frequency between chromatin fragments that gather at a certain distance on the linear position of the genome through steps such as fixation, shearing, cross-linking, and sequencing. These interaction frequencies can be used to infer the spatial position relationship between the corresponding chromatin fragments. Hi-C sequencing data are organized into a Hi-C contact matrix. Each subscript in the matrix identifies a chromatin segment (bin), and the value of the matrix element represents the interaction frequency between the corresponding bins.

[0003] Chromatin spatial domains mainly include A / B compartments, subcompartments, topologically associated domains (TADs), chromatin stripes, and chromatin loops. Among them, the TADs domain is shown in Hi-C data as follows: the frequency of interactions between chromatin fragments inside the TADs domain is significantly higher than that outside it. TADs have important biological functions. They regulate gene expression by limiting enhancer-promoter interactions. The destruction of TADs boundaries is associated with a variety of diseases, such as cancer and various limb deformities. Further studies have found that the TADs domain presents a complex hierarchical structure, which is manifested in two aspects: first, there may be a relatively small sub-TAD domain inside the larger parent TADs domain, that is, the sub-TAD domain is completely located inside the parent TADs domain; second, the two TADs domains overlap with each other, but there is no inclusion relationship, which is called the local intersection of TADs domains. The hierarchy of TADs has significantly different biological significance. For example, the external parent TADs domain remains almost constant during cell differentiation, while the internal daughter TADs domain undergoes drastic changes; the regulatory isolation effect of the parent TADs domain is weaker than that of the daughter TADs domain. Therefore, as an important component of the three-dimensional genome topology, the identification of hierarchical TADs domains is conducive to the construction of the three-dimensional structure of the genome. At the same time, by analyzing its hierarchical structure, it can provide a biological basis for the treatment of complex genetic diseases such as cancer and limb deformities.

[0004] Among the existing methods for identifying hierarchical TADs domains, GMAP models the contact frequency in the Hi-C contact matrix based on a statistical model, and combines the Gaussian mixture model with the proportion test to identify hierarchical TADs domains; 3DNetMod transforms the problem of identifying hierarchical TADs domains into a community detection problem of different sizes in the network, and uses the Louvain community detection algorithm to identify TADs domains; OnTAD defines a scoring function to optimize the possible TADs domain structures identified by the linear score to identify hierarchical TADs domains; SpectralTAD identifies the hierarchical structure of TADs domains through recursive spectral clustering. Although these methods can identify TADs domains in Hi-C data to a certain extent, there are still two problems: (1) A nested assumption is made for the hierarchical organization of TADs domains, that is, the hierarchical form of TADs can only be that small TADs domains are completely contained within large TADs domains, while the actual TADs hierarchical structure has various possibilities, including non-intersecting, nested, and partially intersecting TADs. (2) Only the identification results of TADs domains can be output, and there is a lack of quantitative evaluation of the importance of each TADs domain, and a lack of explanation for the quantitative relationship between the frequency of interaction between TADs domains and chromatin fragments. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a hierarchical TADs domain identification method and system based on explainable deep learning in view of the deficiencies in the above-mentioned prior art, so as to solve the technical problems that the current TADs domain identification method is incomplete in identifying hierarchical TADs domains and lacks quantitative evaluation of the importance of TAD domains.

[0006] The objective of the present invention is achieved by the following technical solutions: In a first aspect, the present invention provides a hierarchical TADs domain identification method based on interpretable deep learning, comprising: Obtain the Hi-C contact matrix and extract the basic TADs domain boundary points in the Hi-C contact matrix; Generate candidate hierarchical TADs domains of all hierarchical forms according to the boundary points of the basic TADs domain and the constraints of the maximum TADs domain; Taking the positional relationship between any two chromatin fragments and all candidate hierarchical TADs domains in the Hi-C contact matrix and the genomic distance between the two fragments as input, and the interaction frequency between the two fragments as output, a deep learning model is trained, wherein the deep learning model is used to represent the relationship between the candidate hierarchical TADs domains and the contact frequency of the chromatin fragments; The trained deep learning model is interpreted to evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix; the most significant hierarchical TADs domain is screened according to the contribution to obtain the final identified hierarchical TADs domain.

[0007] As a further improvement of the present invention, the basic TADs domain boundary points in the Hi-C contact matrix are extracted, specifically including: The TopDom algorithm was used to extract the basic TADs domain boundary points on the diagonal of the Hi-C contact matrix.

[0008] As a further improvement of the present invention, according to the restriction conditions of the basic TADs domain boundary points and the maximum TADs domain, all possible hierarchical forms of candidate hierarchical TADs domains are generated, specifically including: The basic TADs domain boundary points are combined by random combination to obtain candidate hierarchical TADs domains; The maximum allowable size of the TADs domain is determined according to a biological priori formula, and the candidate hierarchical TADs domains are screened according to the maximum allowable size of the TADs domain.

[0009] As a further improvement of the present invention, any two chromatin fragments in the Hi-C contact matrix, the candidate hierarchical TADs domains, and the genomic distance between any two chromatin fragments are used as input features, and the interaction frequency between the two chromatin fragments is used as an output feature to train a deep learning model, specifically including: For each contact frequency in the Hi-C contact matrix, an n+1-dimensional feature vector is constructed, which is used to indicate the positional relationship between the contact frequency and all candidate hierarchical TADs domains, where n is the number of candidate hierarchical TADs domains, and n>1; In the n+1-dimensional feature vector, if the corresponding contact frequency of the first n dimensions of the feature vector falls into the candidate hierarchical TADs domain, the value of the corresponding position of the feature vector is set to 1, otherwise, it is set to 0; if the n+1th dimension of the feature vector records the genomic distance between the two chromatin fragments involved in the contact frequency, d , then set to log 10 ( d +1); For the contact frequency, if the initial contact frequency is y , then the contact frequency after logarithmic transformation is log 10 ( y +1).

[0010] As a further improvement of the present invention, the deep learning model is formed based on a network framework of a one-dimensional convolutional layer and a fully connected layer, wherein the one-dimensional convolutional layer is used to extract the spatial position characteristics of each hierarchical TADs structural domain in the feature vector; the fully connected layer is used to further reduce the dimension of the spatial position characteristics and predict the value of the contact frequency.

[0011] As a further improvement of the present invention, the activation function in the deep learning model adopts a linear rectification function or a Sigmoid function.

[0012] As a further improvement of the present invention, the trained deep learning model is interpreted, specifically including: The trained deep learning model is interpreted using the SHAP interpretable algorithm, and the prediction results in the deep learning model are interpreted as the sum of the contribution value of each input feature and the benchmark output value; the contribution value of each input feature is used to indicate the importance of the input feature to the output; the benchmark output value is the prediction value of the deep learning model under the background input; the background input selection method is: n All dimension features are set to 0, and the last dimension genome distance is set to the maximum genome distance, which is the specified maximum TADs domain size; Calculate the SHAP values ​​of the first n dimensions; obtain an n×m contribution matrix, the columns of which are used for the contribution of each candidate hierarchical TADs domain to the corresponding contact frequency, m is the total number of contact frequencies, and n is the number of candidate hierarchical TADs domains.

[0013] As a further improvement of the present invention, the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix is ​​evaluated, specifically including: The negative elements in the contribution matrix are filtered out, and the remaining non-zero elements are averaged row by row to comprehensively evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix, and an n-dimensional TADs contribution column vector is obtained.

[0014] As a further improvement of the present invention, for the n-dimensional TADs contribution column vector, a threshold is selected according to the fractional threshold method, and the selected threshold is used to screen out the hierarchical TADs domains that contribute most significantly to the Hi-C contact matrix as the finally identified hierarchical TADs domains.

[0015] In a second aspect, the present invention provides a hierarchical TADs domain identification system based on interpretable deep learning, which is used to implement the above-mentioned hierarchical TADs domain identification method based on interpretable deep learning, comprising: Data preparation module, obtains Hi-C contact matrix; extracts basic TADs domain boundary points in Hi-C contact matrix; The data processing module generates candidate hierarchical TADs domains of all hierarchical forms according to the boundary points of the basic TADs domain and the restriction conditions of the maximum TADs domain; A model building and training module, which takes any two chromatin fragments in the Hi-C contact matrix, the candidate hierarchical TADs domains, and the genomic distance between any two chromatin fragments as input features, and the interaction frequency between the two chromatin fragments as output features, to build a deep learning model, wherein the deep learning model is used to represent the relationship between the candidate hierarchical TADs domains and the contact frequency of the chromatin fragments; The hierarchical TADs domain contribution calculation module is used to interpret the trained deep learning model and evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix; the hierarchical TADs domain is identified based on the contribution to obtain the final identified hierarchical TADs domain.

[0016] In a third aspect, the present invention provides a computer-readable storage medium storing one or more programs, wherein the one or more programs include instructions, which, when executed by a computing device, enable the computing device to perform the above-mentioned hierarchical TADs domain identification method based on interpretable deep learning.

[0017] In a fourth aspect, the present invention provides a computing device, comprising: One or more processors, a memory and one or more programs, wherein the one or more programs are stored in the memory and configured to be executed by the one or more processors, and the one or more programs include steps for executing the above-mentioned hierarchical TADs domain identification method based on interpretable deep learning.

[0018] The beneficial effects of the present invention are as follows: the present invention provides a hierarchical TADs domain identification method based on interpretable deep learning, which does not make assumptions or restrictions on the hierarchical form of the hierarchical TADs domain, and comprehensively considers hierarchical TADs of all structural forms; a deep learning model is established to accurately fit the relationship between all candidate hierarchical TADs domains and the frequencies in the Hi-C contact matrix, and interpret the quantitative relationship between the two; the deep learning model is interpreted, and in addition to identifying the hierarchical TADs domain, the importance of each hierarchical TADs domain can also be evaluated according to the SHAP value, which is convenient for use by biological and medical workers.

[0019] Furthermore, the present invention considers TADs domains at different levels and uses model interpretation to evaluate the contribution of candidate domains, thereby improving the accuracy of TADs identification.

[0020] Furthermore, the TopDom algorithm can more accurately identify the boundary points of TADs by analyzing the diagonal area of ​​the contact matrix. It uses the high interaction frequency inside the chromatin region and the low interaction frequency at the boundary to identify obvious boundary signals in the data. The algorithm can process Hi-C data of different resolutions, has strong adaptability, and can provide reliable results in genome studies of different scales.

[0021] Furthermore, by converting complex genomic distance and contact frequency information into numerical features, feature vectors can be used more effectively in machine learning algorithms. This approach simplifies the data processing process and improves computational efficiency. Through feature vectors, the model can more easily interpret the relationship between genomic distance and contact frequency. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0023] Figure 1 It is a flowchart of a hierarchical TADs domain identification method based on interpretable deep learning provided in an embodiment of the present invention.

[0024] Figure 2 It is a schematic diagram of generating candidate hierarchical TADs domains provided by an embodiment of the present invention.

[0025] Figure 3 It is an architectural diagram of the deep learning model provided by an embodiment of the present invention.

[0026] Figure 4 This is an architectural diagram of a hierarchical TADs domain identification system based on interpretable deep learning provided in an embodiment of the present invention.

[0027] Figure 5 It is a heat map of experimental results on simulated Hi-C data provided in an embodiment of the present invention; (a) is a heat map of the original Hi-C contact matrix; (b) is a comparison chart of the fitted contact frequency and the original contact frequency; (c) is the final hierarchical TADs domain recognition result of the present invention.

[0028] Figure 6 It is a heat map of experimental results on the experimental Hi-C data provided in the embodiment of the present invention; (a) is a heat map of the Hi-C contact matrix of chromosome 1 of the GM12878 cell line; (b) is a comparison chart of the fitted contact frequency and the original contact frequency; (c) is the final hierarchical TADs domain recognition result of the present invention.

[0029] Figure 7 It is a schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0030] In order to make the purpose and technical solution of the present invention clearer and easier to understand, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. The specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0031] The technical solution of the present invention will be clearly and completely described below in conjunction with the accompanying drawings and specific embodiments, wherein the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments.

[0032] Example 1 like Figure 1-Figure 3 As shown, this embodiment provides a hierarchical TADs domain identification method based on interpretable deep learning. The following is a specific implementation method.

[0033] S1: Obtain the Hi-C contact matrix and extract the basic TADs domain boundary points in the Hi-C contact matrix.

[0034] In this implementation, in order to remove the systematic deviation in the Hi-C contact matrix while retaining the basic TADs domain boundary point information, the Hi-C contact matrix was normalized by matrix transformation. The normalization methods include ICE (iterative correction and eigenvector decomposition) and KR (Knight-Ruiz) normalization.

[0035] In this implementation, KR normalization calculates the total frequency of rows and columns of the contact frequency matrix, gradually adjusts each element to eliminate sequencing depth and technical bias, and iterates multiple times until the matrix converges, ultimately generating a standardized Hi-C contact matrix.

[0036] ICE standardization corrects background noise and improves data accuracy by constructing an expectation matrix and correction factors, using an iterative process to continuously update the original contact frequency matrix, while using eigenvalue decomposition to ensure statistical characteristics.

[0037] Furthermore, the TopDom algorithm was used to extract the basic TADs domain boundary points on the diagonal of the Hi-C contact matrix, including: Different from the single diamond sliding window in the original TopDom algorithm, this implementation uses multiple sliding windows of different sizes to capture TADs domain boundary points of different scales. The maximum, minimum, and number of sliding window sizes are specified by the user.

[0038] First, for a specific sliding window, this implementation calculates the Hi-C average signal in different sliding windows, and uses the local weighted average curve fitting method to smooth the Hi-C average signal. Then, the local minimum and local maximum extraction algorithm is called to use the local minimum of the Hi-C average signal curve as the preliminary basic TADs domain boundary point, and filter out non-significant basic TADs domain boundary points whose local minimum and local maximum are too small. Finally, the basic TADs domain boundary points obtained from all sliding windows are taken as the union, and the boundary points with a genomic distance of ±1bin are merged into the same boundary point to reduce the accidental deviation of boundary point identification.

[0039] S2: Generate candidate hierarchical TADs domains of all hierarchical forms based on the boundary points of the basic TADs domain and the constraints of the maximum TADs domain.

[0040] Specifically, the boundary points of basic TADs domains are combined by random combination to obtain candidate hierarchical TADs domains; the maximum size of the allowed TADs domains is determined according to a biological priori formula, and the candidate hierarchical TADs domains are screened according to the maximum size of the allowed TADs domains.

[0041] Specifically, all TADs domain boundary points extracted in step 1 are randomly combined in pairs. k A total of 100 Hi-C contact matrices of the boundary points of 10 TADs domains can be generated. All possible candidate hierarchical TADs domains (k>1).

[0042] Furthermore, in order to make the finally identified hierarchical TADs domain more biologically meaningful, the user specifies the maximum size of the TADs domain allowed according to biological priors, and filters out TADs domains larger than the specified maximum size. The number of candidate hierarchical TADs domains that meet the requirements after filtering is recorded as n .

[0043] This implementation Figure 2 Taking the Hi-C contact matrix in as an example, the method of generating candidate hierarchical TADs domains is explained. It can be seen that Figure 2 There are a total of 8 basic TADs domain boundary points in the , so all possible candidate hierarchical TADs domains have The dotted line indicates the size of the largest TADs domain. Therefore, after filtering, there are 18 candidate hierarchical TADs domains that meet the requirements, and their serial numbers are marked in the upper right corner of the corresponding candidate hierarchical TADs domains in the figure.

[0044] S3: Construct a deep learning model that is used to represent the relationship between candidate hierarchical TADs domains and the contact frequency of chromatin fragments.

[0045] Specifically, the positional relationship between any two chromatin fragments and all candidate hierarchical TADs domains in the Hi-C contact matrix, as well as the genomic distance between the fragments (measured by the bins spanned by the two fragments) are used as feature inputs, and the interaction frequency between the two fragments is used as the result output. A deep learning model is established to map the relationship between the contact frequencies of candidate hierarchical TADs domains and chromatin fragments.

[0046] Specifically, this implementation constructs a n +1-dimensional feature vector, used to indicate the positional relationship between the contact frequency and all candidate hierarchical TADs domains, where n is the number of candidate hierarchical TADs domains, n>1.

[0047] Specifically, this implementation is for the front of the feature vector n If the corresponding contact frequency falls within a candidate hierarchical TADs domain, the value of the corresponding position of the feature vector is set to 1, otherwise, it is set to 0. For the last dimension of the feature vector, if the genomic distance between the two chromatin fragments involved in the contact frequency is d , then set to log10 ( d +1). For the contact frequency, the logarithmic transformed contact frequency is also used for modeling. If the initial contact frequency is y , then the transformed contact frequency is log10 ( y +1). The purpose of the logarithmic transformation is to make it easier for the model to capture the exponential decay relationship between contact frequency and genomic distance.

[0048] In this embodiment, the deep learning model is formed based on a network framework of a one-dimensional convolutional layer (1D-convolutional neural network, 1D-CNN) and a fully connected layer (fully connected layer, FC). Among them, the one-dimensional convolutional layer is used to extract the spatial position features of each hierarchical TADs domain in the feature vector; the fully connected layer is used to further reduce the dimension of the spatial position features and predict the value of the contact frequency.

[0049] Specifically, the activation function in the deep learning model of this implementation uses a linear rectified linear unit (ReLu) or Sigmoid to perform nonlinear transformation on the output of each layer of the network; the Dropout layer prevents overfitting by randomly deleting some neurons, and the Dropout probability is set in the range of 0.3~0.5. The loss function is set to the mean square error between the contact frequency predicted by the deep learning model and the actual contact frequency, that is:

[0050] In the formula, y i represents the actual contact frequency, x i Representation and y i The corresponding eigenvector, f (·) represents the mapping relationship defined by the deep learning model.

[0051] In this implementation, the deep learning model is Figure 3 As shown in the figure, the input of the model is an n+1-dimensional feature vector, and the output is the interaction frequency between chromatin fragments. The model is divided into an input layer, a hidden layer, and an output layer. The input layer is used to receive the input feature vector, the hidden layer is used to extract the spatial position relationship of the candidate hierarchical TADs domain, and the output layer is used to predict the interaction frequency.

[0052] In the process of training the deep learning model, in order to speed up the convergence of the loss function and make the model more biologically meaningful, only the contact frequencies that fall within at least one candidate hierarchical TADs domain are used for training. n The eigenvector whose dimension is not all zero, the number of contact frequencies that meet this condition is m All contact frequencies are divided into training sets and test sets in a ratio of 7:3 to 9:1. The optimizer can use optimization algorithms with adaptive learning rates such as Adam and RMSprop. The size of the training batch can be set to 8, 16, 32, etc. according to the computing power of the computer. The initial learning rate is set to 1×10 -3 ~1×10 -2 .

[0053] S4: Interpret the trained deep learning model and evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix; identify the hierarchical TADs domain based on the contribution to obtain the final identified hierarchical TADs domain.

[0054] The specific model interpretation method of this implementation is: the trained deep learning model is interpreted using the SHAP interpretable algorithm, and the prediction results in the deep learning model are interpreted as the sum of the contribution value of each input feature and the benchmark output value. Among them, the contribution value of each input feature indicates the importance of the feature to the output; the benchmark output value is the prediction value of the deep learning model under the background input.

[0055] The mathematical representation of the SHAP interpretable algorithm implemented in this paper is:

[0056] Where: f is the deep learning model to be explained, x is the input feature vector, is the reference output value, The attributed value of each feature, that is, the SHAP value.

[0057] Background input selection mode: n The first dimension features are all set to 0, and the last dimension genome distance is set to the maximum genome distance, that is, the maximum TADs domain size specified. The purpose is to make the background input contain no TADs domain information, so that the interpretation of the model is completely based on the subsequent valid feature input, so as to improve the reliability of the model interpretation.

[0058] Considering that the last dimension of the input feature is the genomic distance information, which has no direct relationship with the contribution of the hierarchical TADs domain, this implementation only calculates the first n SHAP value of the dimension. Get a n × m A contribution matrix is ​​shown in Figure 1, where a column of the matrix represents the contribution of each candidate hierarchical TADs domain to the corresponding contact frequency. Where m is the total number of contact frequencies, and n is the number of candidate hierarchical TADs domains.

[0059] Furthermore, from a biological perspective, the hierarchical TADs domain can only have a positive contribution or no contribution to the contact frequency. Therefore, this implementation filters out the negative elements in the contribution matrix and averages the remaining non-zero elements row by row to comprehensively evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix, and obtains a n dimensional TADs contribution column vector. Each element of the vector represents the contribution of the corresponding candidate hierarchical TADs domain to the entire contact matrix, that is, their significance or importance to the entire Hi-C contact matrix.

[0060] For the n-dimensional TADs contribution column vector, a threshold is selected according to the fractional threshold method, and the selected threshold is used to screen out the hierarchical TADs domains that contribute most significantly to the Hi-C contact matrix as the final identified hierarchical TADs domains. In this implementation, the threshold selection strategy is the quantile threshold method, and the 25% quantile, 50% quantile, 75% quantile, etc. of the n-dimensional TADs contribution column vector can be selected as the threshold to screen out the hierarchical TADs domains with contributions higher than the threshold. The higher the threshold, the more significant the identified hierarchical TADs domains.

[0061] Example 2 This embodiment provides a hierarchical TADs domain identification system based on interpretable deep learning. The structure diagram of the hierarchical TADs domain identification system is as follows: Figure 4 As shown, it includes the following modules Data preparation module, obtains Hi-C contact matrix; extracts basic TADs domain boundary points in Hi-C contact matrix; The data processing module generates candidate hierarchical TADs domains of all hierarchical forms according to the boundary points of the basic TADs domain and the restriction conditions of the maximum TADs domain; The model building and training module uses any two chromatin fragments in the Hi-C contact matrix, the candidate hierarchical TADs domain, and the genomic distance between any two chromatin fragments as input features, and the interaction frequency between the two chromatin fragments as output features to build a deep learning model, which is used to represent the relationship between the candidate hierarchical TADs domain and the contact frequency of the chromatin fragments. Specifically, a deep learning model based on a one-dimensional convolutional layer (1D-CNN) and a fully connected layer (FC) is built according to the network parameter configuration specified by the user, and the training set and the test set are divided according to the specified ratio. The training data is input into the deep learning model in batches of a certain size according to the training parameter configuration, and the model is trained using a specified optimization algorithm according to the loss function to optimize the weights.

[0062] The hierarchical TADs domain contribution calculation module is used to interpret the trained deep learning model, evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix, identify the hierarchical TADs domain according to the contribution, and obtain the final identified hierarchical TADs domain. Specifically, it is used to extract the contribution of each candidate hierarchical TADs domain to each contact frequency using the SHAP model interpretation algorithm for the model with trained weights, and generate a contribution matrix. And according to the biological prior, the negative values ​​of the contribution matrix are filtered out, and the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix is ​​comprehensively evaluated.

[0063] This embodiment also includes a Hi-C contact matrix normalization module, a training data preparation module, and a hierarchical TADs domain screening and visualization module. Each module specifically includes: The Hi-C contact matrix normalization module is used to normalize the contact matrix using matrix transformation methods such as ICE and KR to remove systematic deviations.

[0064] The training data preparation module is used to identify the boundary points of basic TADs domains, and combine them into candidate hierarchical TADs domains based on the boundary points of basic TADs domains, and filter out the parts that do not conform to biological priors based on their size; the input feature vector of the model is constructed based on the positional relationship between the contact frequency and the candidate hierarchical TADs domains, and the genomic distance between the two chromatin fragments.

[0065] The hierarchical TADs domain screening and visualization module is used to post-process the contribution vectors of each hierarchical TADs domain output by the model interpretation, determine the threshold according to the quantile method, screen out the most significant hierarchical TADs domains above the threshold, and use the heatmap function in python or matlab to visualize the Hi-C contact matrix and annotate the hierarchical structure of the identified TADs domains.

[0066] The hierarchical TADs domain identification system based on interpretable deep learning of this embodiment corresponds to the hierarchical TADs domain identification method based on interpretable deep learning of the first embodiment mentioned above; wherein the Hi-C data standardization module and the training data preparation module correspond to steps 1 and 2, the model construction and training module corresponds to step 3, the hierarchical TADs domain contribution calculation module, and the hierarchical TADs domain screening and visualization module correspond to step 4, so they will not be repeated here.

[0067] Example 3 This embodiment, according to the method in Example 1, uses simulated Hi-C contact matrices with locally overlapping TADs domains and without locally overlapping TADs domains to carry out experiments, wherein the local overlap rate is 0.15 and the noise level is 20%. The number of samples is 5. This implementation calculates the average Jaccrad index and F1 score of the identified hierarchical TADs domain and the real TADs domain, and uses it as the recognition accuracy of the TADs domain, and compares the recognition accuracy with four representative methods of GMAP, 3DNetMod, OnTAD, and SpectralTAD. The comparison results are shown in Table 1. It can be seen that the method of the present invention has comprehensively surpassed the comparison method in both the Jaccrad index index and the F1 score index, and has significant advantages.

[0068] Table 1. Comparison of the accuracy of various methods in identifying hierarchical TADs domains on simulated Hi-C data

[0069] Furthermore, in order to intuitively demonstrate the process of identifying hierarchical TADs domains on simulated Hi-C data with locally overlapping TADs domains, the results of model fitting and the final results of hierarchical TADs domain identification are visualized, as shown in FIG. Figure 5 shown. Figure 5 (a) is the original Hi-C contact matrix heat map, Figure 5 (b) Comparison between the fitted contact frequency and the original contact frequency. It can be seen that the fitted value is very close to the original value, and the hierarchical TADs domain near the diagonal is also well captured by the model. Figure 5 (c) shows the final hierarchical TADs domain recognition result of the present invention. It can be seen that the method of the present invention successfully identifies the most significant hierarchical TADs domain in the Hi-C contact matrix and can effectively capture locally overlapping hierarchical TADs domains.

[0070] In addition, as another example, this embodiment also uses the 25Kb experimental Hi-C contact matrix of chromosome 1 of the GM12878 cell line to carry out experiments, where the number of samples is 5. Since it is impossible to know the true hierarchical TADs domain in the experimental Hi-C data, the average Jaccard index and F1 score of the hierarchical TADs domain between two samples are used as indicators to compare the recognition accuracy of each method. The results of the comparison are shown in Table 2. It can be seen that the method of the present invention has achieved the best performance in both the Jaccradindex indicator and the F1 score indicator.

[0071] Table 2. Comparison of the accuracy of various methods in identifying hierarchical TADs domains on experimental Hi-C data

[0072] Furthermore, in order to intuitively demonstrate the process of identifying hierarchical TADs domains in the experimental Hi-C data of the present invention, the results of model fitting and the final results of hierarchical TADs domain identification are visualized. Figure 6 As shown, similar to the results on simulated Hi-C data, the method of the present invention also achieved excellent contact frequency fitting and hierarchical TADs domain recognition effects on experimental Hi-C data.

[0073] Example 4 In another embodiment of the present invention, the present invention also provides a storage medium, specifically a computer-readable storage medium (Memory), which is a memory device in a terminal device for storing programs and data. It is understandable that the computer-readable storage medium here can include both the built-in storage medium in the terminal device and the extended storage medium supported by the terminal device, and can be any tangible medium containing or storing a program, which can be used by an instruction execution system, device or device or used in combination with it. The computer-readable storage medium provides a storage space, which stores the operating system of the terminal. In addition, one or more instructions suitable for being loaded and executed by the processor are also stored in the storage space, and these instructions can be one or more computer programs (including program codes). It should be noted that more specific examples (non-exhaustive list) of the computer-readable storage medium here include: an electrical connection with one or more wires, a portable disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the above.

[0074] Computer readable storage media also include data signals propagated in baseband or as part of a carrier wave, which carry readable program codes. Such propagated data signals can take a variety of forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination of the above. The readable storage medium can also be any readable medium other than a readable storage medium, which can send, propagate, or transmit a program for use by or in conjunction with an instruction execution system, device, or device. The program code contained on the readable storage medium can be transmitted using any appropriate medium, including but not limited to wireless, wired, optical cable, RF, etc., or any suitable combination of the above.

[0075] Program code for performing the operations of the present invention may be written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Java, C++, etc., and conventional procedural programming languages ​​such as "C" or similar programming languages. The program code may be executed entirely on the user computing device, partially on the user device, as a separate software package, partially on the user computing device and partially on a remote computing device, or entirely on a remote computing device or server. In the case of a remote computing device, the remote computing device may be connected to the user computing device through any type of network, including a local area network (LAN) or a wide area network (WAN), or may be connected to an external computing device (e.g., through the Internet using an Internet service provider).

[0076] One or more instructions stored in a computer-readable storage medium may be loaded and executed by a processor to implement the corresponding steps of the hierarchical TADs domain identification method based on interpretable deep learning in the above embodiment; one or more instructions in the computer-readable storage medium are loaded and executed by the processor as follows: Obtain Hi-C contact matrix; extract basic TADs domain boundary points in Hi-C contact matrix; Generate candidate hierarchical TADs domains of all hierarchical forms according to the boundary points of the basic TADs domain and the constraints of the maximum TADs domain; Taking any two chromatin fragments, candidate hierarchical TADs domains and the genomic distance between any two chromatin fragments in the Hi-C contact matrix as input features, the input features are input into the trained hierarchical TADs domain recognition model, and the recognition model is used to represent the relationship between the contact frequency between the candidate hierarchical TADs domains and the chromatin fragments; The trained recognition model is interpreted to evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix; the hierarchical TADs domain is identified based on the interpreted recognition model to obtain the final identified hierarchical TADs domain.

[0077] Example 5 See also Figure 7 , the terminal device is a computer device, and the computer device 60 of this embodiment includes: a processor 61, a memory 62, and a computer program 63 stored in the memory 62 and executable on the processor 61. When the computer program 63 is executed by the processor 61, the hierarchical TADs domain identification method based on interpretable deep learning in the embodiment is implemented. To avoid repetition, it is not described one by one here. Alternatively, when the computer program 63 is executed by the processor 61, the functions of each model / unit in the computing system of the hierarchical TADs domain identification method based on interpretable deep learning in Example 1 are implemented. To avoid repetition, it is not described one by one here.

[0078] The computer device 60 may be a computing device such as a desktop computer, a notebook, a PDA, or a cloud server. The computer device 60 may include, but is not limited to, a processor 61 and a memory 62. Those skilled in the art will appreciate that Figure 7 This is only an example of the computer device 60 and does not constitute a limitation of the computer device 60. It may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, the computer device may also include input and output devices, network access devices, buses, etc.

[0079] The processor 61 may be a central processing unit (CPU), or other general-purpose processors, central processing units, graphics processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, data processing logic based on quantum computing, discrete hardware components, etc. A general-purpose processor may be a microprocessor or the processor may be any conventional processor, etc.

[0080] The memory 62 may be an internal storage unit of the computer device 60, such as a hard disk or memory of the computer device 60. The memory 62 may also be an external storage device of the computer device 60, such as a plug-in hard disk, a smart memory card (SmartMediaCard, SMC), a secure digital (SecureDigital, SD) card, a flash card (FlashCard), etc. equipped on the computer device 60.

[0081] Furthermore, the memory 62 may include both an internal storage unit of the computer device 60 and an external storage device. The memory 62 is used to store computer programs and other programs and data required by the computer device. The memory 62 may also be used to temporarily store data that has been output or is to be output.

[0082] Any reference to a memory, database or other medium used in the embodiments provided in this application may include at least one of a non-volatile and a volatile memory. Non-volatile memory may include a read-only memory (ROM), a magnetic tape, a floppy disk, a flash memory, an optical memory, a high-density embedded non-volatile memory, a resistive random access memory (ReRAM), a magnetic random access memory (MRAM), a ferroelectric random access memory (FRAM), a phase change memory (PCM), a graphene memory, etc. Volatile memory may include a random access memory (RAM) or an external cache memory, etc. As an illustration and not limitation, RAM may be in various forms, such as a static random access memory (SRAM) or a dynamic random access memory (DRAM), etc.

[0083] The database involved in each embodiment provided in this application may include at least one of a relational database and a non-relational database. The non-relational database may include a distributed database based on blockchain, etc., but is not limited thereto. The processor involved in each embodiment provided in this application may be a general-purpose processor, a central processing unit, a graphics processor, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., but is not limited thereto.

Claims

1. A hierarchical TADs domain identification method based on interpretable deep learning, characterized in that: include: Obtain the Hi-C contact matrix and extract the basic TADs domain boundary points in the Hi-C contact matrix; Generate candidate hierarchical TADs domains of all hierarchical forms according to the boundary points of the basic TADs domain and the constraints of the maximum TADs domain; Taking the positional relationship between any two chromatin fragments and all candidate hierarchical TADs domains in the Hi-C contact matrix and the genomic distance between the two fragments as input, and the interaction frequency between the two fragments as output, a deep learning model is trained, wherein the deep learning model is used to represent the relationship between the candidate hierarchical TADs domains and the contact frequency of the chromatin fragments; The trained deep learning model is interpreted to evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix; the most significant hierarchical TADs domain is screened according to the contribution to obtain the final identified hierarchical TADs domain.

2. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 1 is characterized in that: Extract the basic TADs domain boundary points in the Hi-C contact matrix, including: The TopDom algorithm was used to extract the basic TADs domain boundary points on the diagonal of the Hi-C contact matrix.

3. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 1 is characterized in that: According to the constraints of the basic TADs domain boundary points and the maximum TADs domain, all possible hierarchical forms of candidate hierarchical TADs domains are generated, including: The basic TADs domain boundary points are combined by random combination to obtain candidate hierarchical TADs domains; The maximum allowable size of the TADs domain is determined according to a biological priori formula, and the candidate hierarchical TADs domains are screened according to the maximum allowable size of the TADs domain.

4. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 1 is characterized in that The deep learning model is trained using the positional relationship between any two chromatin fragments and all candidate hierarchical TADs domains in the Hi-C contact matrix, as well as the genomic distance between the two fragments, and the interaction frequency between the two fragments as output. The specific optimization includes: For each contact frequency in the Hi-C contact matrix, an n+1-dimensional feature vector is constructed, which is used to indicate the positional relationship between the contact frequency and all candidate hierarchical TADs domains, where n is the number of candidate hierarchical TADs domains, and n>1; In the n+1-dimensional feature vector, if the corresponding contact frequency of the first n dimensions of the feature vector falls into the candidate hierarchical TADs domain, the value of the corresponding position of the feature vector is set to 1, otherwise, it is set to 0; if the n+1th dimension of the feature vector records the genomic distance between the two chromatin fragments involved in the contact frequency, d , then set to log 10 ( d +1); For the contact frequency, if the initial contact frequency is y , then the contact frequency after logarithmic transformation is log 10 ( y +1).

5. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 4 is characterized in that: The deep learning model is formed based on a network framework of a one-dimensional convolutional layer and a fully connected layer, wherein the one-dimensional convolutional layer is used to extract the spatial position features of each hierarchical TADs domain in the feature vector; The fully connected layer is used to further reduce the dimension of the spatial position feature and predict the value of the contact frequency.

6. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 5 is characterized in that: The activation function in the deep learning model uses a linear rectification function or a Sigmoid function.

7. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 1 is characterized in that: Explain the trained deep learning model, including: The trained deep learning model is interpreted using the SHAP interpretable algorithm, and the prediction results in the deep learning model are interpreted as the sum of the contribution value of each input feature and the benchmark output value; the contribution value of each input feature is used to indicate the importance of the input feature to the output; the benchmark output value is the prediction value of the deep learning model under background input; the background input selection method is: the first n-dimensional features are all set to 0, and the last-dimensional genome distance is set to the maximum genome distance, that is, the specified maximum TADs domain size; Calculate the SHAP values ​​of the first n dimensions; obtain an n×m contribution matrix, the columns of which are used for the contribution of each candidate hierarchical TADs domain to the corresponding contact frequency, m is the total number of contact frequencies, and n is the number of candidate hierarchical TADs domains.

8. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 7 is characterized in that: Evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix, including: The negative elements in the contribution matrix are filtered out, and the remaining non-zero elements are averaged row by row to comprehensively evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix, and an n-dimensional TADs contribution column vector is obtained.

9. The hierarchical TADs domain identification method based on interpretable deep learning according to claim 8 is characterized in that: For the n-dimensional TADs contribution column vector, a threshold is selected according to the fractional threshold method, and the selected threshold is used to screen out the hierarchical TADs domains that contribute most significantly to the Hi-C contact matrix as the final identified hierarchical TADs domains.

10. A hierarchical TADs domain identification system based on interpretable deep learning, used to implement the hierarchical TADs domain identification method based on interpretable deep learning according to any one of claims 1 to 9, characterized in that: include: Data preparation module, to obtain Hi-C contact matrix; Extract basic TADs domain boundary points in Hi-C contact matrix; The data processing module generates candidate hierarchical TADs domains of all hierarchical forms according to the boundary points of the basic TADs domain and the constraints of the maximum TADs domain; A model building and training module, which takes the positional relationship between any two chromatin fragments and all candidate hierarchical TADs domains in the Hi-C contact matrix and the genomic distance between the two fragments as input, and takes the interaction frequency between the two fragments as output to build a deep learning model, wherein the deep learning model is used to represent the relationship between the candidate hierarchical TADs domains and the contact frequency of the chromatin fragments; The hierarchical TADs domain contribution calculation module is used to interpret the trained deep learning model and evaluate the contribution of each candidate hierarchical TADs domain to the overall contact frequency in the Hi-C contact matrix; the hierarchical TADs domain is identified based on the contribution to obtain the final identified hierarchical TADs domain.

Images