Topological association structural domain division method based on chromatin contact map

Abstract

The invention discloses a chromatin contact atlas-based topological correlation domain division method, which comprises the following steps of: acquiring and preprocessing chromatin contact atlas data, and constructing a normalized contact matrix; according to the normalized contact matrix, calculating single-segment Infomap entropy values of all possible division intervals; solving an Infomap entropy value of a minimum part in each state through dynamic programming recursion, and calculating a complete Infomap entropy value under each step number according to the Infomap entropy value of the minimum part in each state; selecting the step number corresponding to the minimum integral Infomap entropy value, and obtaining the optimal topological correlation structure domain division through backtracking; an approximation error of the partition structure is calculated to assess partition accuracy. According to the method, a machine learning model does not need to be trained, whole genome chromatin contact information can be fully utilized, optimal division of TAD is automatically achieved, the division accuracy and repeatability are improved, meanwhile, the calculation complexity and cost are reduced, and an efficient and reliable tool is provided for chromatin three-dimensional structure and gene regulation and control research.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of molecular biology research, and particularly relates to a topologically associating domain division method based on chromatin contact maps. BACKGROUND

[0002] Topologically associating domains (TADs) are chromatin regions in the genome that exhibit enriched intra-domain interactions and suppressed inter-domain interactions, and are widely considered as the basic units of chromosome three-dimensional organization and its functional regulation (such as transcription and replication). Accurate identification of topologically associating domains helps to analyze the relationship between chromatin three-dimensional structure and biological function, thereby promoting the understanding of gene regulation mechanisms and related diseases.

[0003] With the development of high-throughput chromosome capture (Hi-C) and other technologies, large-scale chromatin contact data can be obtained. These data provide a basis for the identification of topologically associating domains, but how to effectively utilize these contact data to accurately measure intra-domain and inter-domain interactions is still a key technical problem. According to the definition, topologically associating domains have two inherent characteristics: one is the enrichment of intra-domain "self-interaction", and the other is the suppression of inter-domain interaction. Therefore, the core of identifying topologically associating domains lies in reasonably quantifying and distinguishing the differences between intra-domain and inter-domain interactions.

[0004] In the prior art, Chinese patent CN118038990A discloses a multi-level chromatin topological domain identification method and system based on community discovery, relating to the technical field of biological gene recognition, which comprises: obtaining a raw Hi-C contact matrix of chromosome sequencing, pre-processing the Hi-C contact matrix to obtain an undirected weighted sub-network; using a community discovery algorithm to divide the undirected weighted sub-network into independent communities, searching for candidate double boundaries in all communities and forming a set; extracting boundary local feature representation of the candidate double boundaries, and inputting the boundary local feature representation into an MLP model to obtain a predicted score of each candidate double boundary; according to the predicted score, screening out candidate double boundaries that need to be retained, and the retained candidate double boundaries form a reliable boundary set; searching the candidate double boundaries in the reliable boundary set to identify multi-level chromatin topological domains TADs.

[0005] However, the method has the following limitations: 1. The existing method relies on community discovery algorithm and local bilateral boundary prediction, and the division result depends more on local features and prediction threshold setting, which cannot guarantee that the identified TAD is the globally optimal division, and local optimization may ignore the global chromatin interaction mode; 2. It is necessary to train the MLP model and design a multi-dimensional feature vector for prediction, which is sensitive to training data and threshold selection; When the training sample is insufficient or the resolution is inconsistent, boundary recognition bias may occur, reducing the recognition accuracy and repeatability, in addition, the MLP model training and feature calculation process have high computational complexity and long time consumption, increasing the cost and resource consumption of the method implementation;

[0006] Therefore, how to design a topological associated domain recognition method with high recognition accuracy and low time complexity has become a technical problem that needs to be solved in the current field. SUMMARY

[0007] The purpose of the present application is to overcome the defects of the prior art and provide a topological associated domain division method based on chromatin contact map.

[0008] The purpose of the present application can be achieved by the following technical solutions:

[0009] In one aspect, the present application provides a topological associated domain division method based on chromatin contact map, comprising the following steps:

[0010] Obtain the chromatin contact map data and perform preprocessing to obtain the normalized contact matrix;

[0011] Solve the Infomap entropy value corresponding to all division domains using the normalized contact matrix obtained after preprocessing;

[0012] Select the optimal topological associated domain division according to the minimum value of the Infomap entropy value corresponding to all division domains.

[0013] Further, the chromatin contact map data includes large-scale chromatin interaction frequency or intensity under multi-resolution, corresponding chromosome information, index information of specific chromosomes or specific genomic regions.

[0014] Further, the obtaining of the chromatin contact map data and the preprocessing thereof specifically includes:

[0015] Extract or construct the original contact count matrix from the contact map data, wherein the value of each element in the original contact count matrix is determined according to the large-scale chromatin interaction frequency or intensity under multi-resolution;

[0016] applying quality filtering on the original contact count matrix, the quality filtering comprising filtering based on mapping quality value (MAPQ) and removing or marking low quality or abnormal contact records, filtering threshold being MAPQ≥30;

[0017] applying balanced normalization on the quality filtered original contact count matrix, the balanced normalization being a Knight-Ruiz algorithm, and the normalized result being a normalized contact matrix A.

[0018] Further, the normalized contact matrix obtained after the preprocessing is used to solve the Infomap entropy values corresponding to all partition domains, specifically comprising:

[0019] According to the given step size parameter ∈ and the total interaction amount of the normalized contact matrix A, the value range of the step number is determined, the step number being a non-negative integer, and the continuous external interaction amount is discretized according to the step size parameter ∈;

[0020] The single segment Infomap entropy value of any possible continuous interval in the normalized contact matrix A is calculated in sequence;

[0021] A dynamic programming table is constructed with the right endpoint index of the interval and the step number as the state, and the minimum partial Infomap entropy value in each state is determined by recursion according to the right endpoint index and the step number dimension using the obtained single segment Infomap entropy value, and the precursor boundary index that produces the minimum value is recorded in the recursion process for backtracking, the minimum partial Infomap entropy value representing the minimum cumulative single segment entropy that can be reached under the given step number budget covering to the right endpoint;

[0022] For each value in the value range of the step number, the corresponding complete Infomap entropy value is calculated according to the corresponding minimum partial Infomap entropy value and the global encoding amount mapped by the step number; the global encoding amount is determined by the step number and the step size ∈.

[0023] Further, the value range of the step number is determined according to the given step size parameter ∈ and the total interaction amount of the normalized contact matrix A, and the formula is:

[0024]

[0025] Wherein, t represents the step number, and the value is a non-negative integer; represents the total interaction amount of the graph G, which is equal to the sum of all elements in the normalized contact matrix A; G is an undirected weighted graph corresponding to the normalized contact matrix A, the nodes in the graph being genomic bins divided according to the selected resolution, A ijrepresents the normalized interaction strength or weight between the i-th genomic bin and the j-th genomic bin, and n is the number of genomic bins obtained by dividing the genome at the selected resolution.

[0026] Further, the continuous interval refers to a section of index interval [i l ,i r ] in the row and column indexes of the normalized contact matrix A, which is adjacent in the linear order of the genome. l ,i r where i l ≤ i r ≤ n, and each continuous interval corresponds to a candidate topological domain on the genome, which is used to evaluate the local interaction tightness between sites in the interval and the interaction relationship outside the interval.

[0027] Further, the single-segment Infomap entropy value of any possible continuous interval in the normalized contact matrix A is calculated in sequence, and the formula is:

[0028]

[0029] where Q(i l ; i r ) represents the single-segment Infomap entropy value of the continuous interval [i l ,i r ], which is used to measure the information encoding contribution of the candidate domain, vol(G) represents the total interaction amount of the graph G, s(i l :i r ) represents the total interaction amount of the vertices in the continuous interval [i l ,i r ] and other vertices in the graph, g(i l :i r ) = vol(i l :i r ) - s(i l :i r ) represents the total interaction amount of the vertices in the continuous interval [i l ,i r ] and the vertices outside the interval.

[0030] Further, the minimum partial Infomap entropy value is:

[0031]

[0032] where ir denotes the index of the right end point of the interval, i prev denotes the index of the right end point of the candidate predecessor interval; t denotes the step number, when i r =1,

[0033] The dynamic programming table is constructed with the index of the right end point of the interval and the step number as states, and specifically comprises the following steps:

[0034] The index of the right end point of the interval i r is taken as a one-dimensional state, and the termination position of the candidate topological association domain is represented;

[0035] The step number t is taken as another one-dimensional state, and the discretization step number under a given coding budget is represented;

[0036] In the dynamic programming recursion process, the minimum part Infomap entropy value H(i r ,t) corresponding to each state (i r ,t) is updated in sequence, and the predecessor index i prev that produces the optimal result is recorded, so that the optimal topological association domain division is finally obtained through backtracking.

[0037] Further, the complete Infomap entropy value is formulaed as:

[0038]

[0039] wherein, denotes the complete Infomap entropy value of the graph G, t denotes the step number, H(n,t) denotes the minimum part Infomap entropy value under the state of the index of the right end point being n and the step number being t, n is the number of genome bins obtained by selecting the resolution division, and t∈ is the global coding amount.

[0040] Further, the optimal topological association domain division is selected according to the minimum value of the Infomap entropy values corresponding to all division domains, and specifically comprises the following steps:

[0041] Within the value range of the step number t , the complete Infomap entropy value H(n,t) corresponding to each step number t value is calculated in sequence , and the sizes thereof are compared;

[0042] The step number t that makes the complete Infomap entropy value H(n,t) obtain the minimum value is selected * ;

[0043] The index of the right end point of the interval in the dynamic programming table is i r =n, and the step number is t=t *the state as the backtracking starting point, wherein n is the number of genomic bins obtained by selecting the resolution division;

[0044] from the state (i r , t), read the predecessor right endpoint index i prev in the dynamic programming table, save the interval [i prev +1, i r ] as an optimal topological association domain;

[0045] update the backtracking state, i r ← i prev , the step number repeat the backtracking step until the start position of the genome is covered;

[0046] The continuous interval sequence obtained by backtracking is arranged in the linear order of the genome, that is, the final optimal topological association domain division result is obtained.

[0047] Compared with the prior art, the present application has the following advantages:

[0048] (1) In the prior art, it is difficult to accurately quantify the intra-domain and inter-domain interactions when identifying topological association domains using the Hi-C contact matrix, resulting in inaccurate TAD boundary identification. The present application constructs a normalized chromatin contact matrix into a graph structure, calculates the single-section Infomap entropy value for all possible continuous intervals, and then determines the minimum cumulative Infomap entropy value through dynamic programming to select the optimal TAD division, thereby achieving accurate measurement of intra-domain enriched interactions and inter-domain insulated interactions and improving the accuracy and stability of TAD boundary identification.

[0049] (2) In the prior art, the training of the MLP model and the multi-dimensional feature vector are dependent on the candidate boundary screening process, which is greatly affected by the training data and threshold selection, resulting in deviation when the sample is insufficient or the resolution is inconsistent, high computational complexity, and high cost. The present application uses a dynamic programming method based on Infomap entropy, does not require training any model, automatically selects the optimal TAD boundary through systematic search and global optimal backtracking of the step number and step length parameter ∈, realizes robust identification under different resolutions and data volumes, significantly reduces the computational complexity and method dependence, and improves the repeatability and applicability.

[0050] (3) In the prior art, the search efficiency of the candidate boundary is low in multi-level TAD division, and it is difficult to obtain high-resolution, multi-level topological domains. The present application uses step number discretization and right endpoint index state recursion in the dynamic programming process, and combines the complete Infomap entropy minimization strategy to realize efficient optimal TAD division of the whole genome, thereby significantly improving the multi-level TAD identification speed and resolution under the premise of ensuring accuracy.

[0051] (4) The method proposed in the present application balances the calculation time and result accuracy by adjusting the step size parameter, when the step size is small, the result accuracy is high, which is suitable for the accurate identification task of small area chromatin data, when the step size is large, the calculation time is shorter, which is suitable for the identification task of chromatin data with more sites. The algorithm proposed in the present application has higher accuracy and lower calculation time, which can meet the needs of topological associated domain identification.

[0052] (5) The method proposed in the present application uses Infomap entropy as the objective function, in which the entropy value of each divided region quantifies the information flow in the interaction, which is used to measure the internal connection tightness of each region and the sparseness of external association; compared with the traditional method using structure entropy and modularity definition as the objective function, the algorithm proposed in the present application can meet the needs of accurate quantification of topological associated domain quality.

[0053] (6) The method proposed in the present application uses the idea of dynamic programming, in which the minimum entropy value corresponding to each divided region and the starting and ending sites of the region can be determined by solving the recursive formula and backtracking process; compared with the traditional method using the greedy idea, the algorithm proposed in the present application can meet the accuracy requirement of topological associated domain division. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 The flow chart of the topological associated domain division method of the embodiment of the present application is shown in the figure.

[0055] Figure 2 The flow chart of the method of step S2 of the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0056] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor should be within the scope of protection of the present application.

[0057] Embodiment 1:

[0058] The application provides a topological associated domain division method based on chromatin contact map, constructs a recursive formula by using the dynamic programming idea, measures the quality of the division area by calculating the Infomap entropy, and balances the precision and calculation efficiency by using the step length. First, the symmetric contact matrix is obtained by processing the chromatin contact map data, and the normalization processing is performed; second, the step number range is determined according to the matrix condition and the step length parameter, the optimal division in each case is determined by using all possible site and step number values through the dynamic programming recursive formula; finally, the approximate optimal division is determined according to the Infomap entropy values corresponding to all possible divisions calculated by the dynamic programming.

[0059] The embodiment provides a topological associated domain division method based on chromatin contact map, as shown in the formula (I), comprising the following steps: Figure 1

[0060] Step S1: obtaining chromatin contact map data and performing preprocessing to obtain a normalized contact matrix;

[0061] The chromatin contact map data comprises large-scale chromatin interaction times or intensity under multiple resolutions, corresponding chromosome information, index information of specific chromosomes or specific genomic regions.

[0062] The preprocessing is that all interactions are combined in the main sample and the duplicate sample, and the interactions are screened according to the mapping quality (MAPQ) filtering condition, and the MAPQ is greater than or equal to 30. The contact matrix is obtained by using the Knight-Ruiz (KR) normalization to process the original matrix data.

[0063] Step S2: solving the Infomap entropy values corresponding to all division domains by using the normalized contact matrix obtained after the preprocessing;

[0064] Step S2, as shown in the formula (II), comprises the following steps: Figure 2

[0065] Step S201: determining the step number value range according to the given step length parameter and the total interaction amount of the normalized contact matrix A, the step number is a non-negative integer, and is used for discretizing the continuous external interaction amount according to the step length parameter;

[0066] The step number value range is determined according to the given step length parameter and the total interaction amount of the normalized contact matrix A, and the formula is as follows:

[0067]

[0068] The step number t is a non-negative integer. ​​vol(G) represents the total interaction of graph G, which is equal to the sum of all elements in the normalized contact matrix A; G is the undirected weighted graph corresponding to the normalized contact matrix A, the nodes of which are genomic bins divided according to the selected resolution, A ij represents the element in the i-th row and the j-th column of the normalized contact matrix A, representing the normalized interaction intensity or weight between the i-th genomic bin and the j-th genomic bin, and n is the number of genomic bins divided according to the selected resolution.

[0069] Step S202: sequentially calculating the single-section Infomap entropy value of each possible continuous interval in the normalized contact matrix A; wherein the continuous interval refers to a section of adjacent indexes [i l ,i r ] in the row and column indexes of the normalized contact matrix A in the linear order of the genome, wherein i l ,i r are the left and right end point indexes of the interval, respectively, and satisfy 1≤i l ≤i r ≤n; each continuous interval corresponds to a candidate topological correlation domain on the genome, which is used to evaluate the local interaction tightness between sites in the interval and the interaction relationship outside the interval.

[0070] The single-section Infomap entropy value of each possible continuous interval in the normalized contact matrix A is sequentially calculated, and the formula is:

[0071]

[0072] wherein Q(i l ;i r ) represents the single-section Infomap entropy value of the continuous interval [i l ,i r ], which is used to measure the information encoding contribution of the candidate domain, vol(G) represents the total interaction of graph G, represents the total interaction between the vertices in the continuous interval [i l ,i r ] and the other vertices in the graph, represents the total interaction between the vertices in the continuous interval [i l ,i r ], and g(i l :i r )=vol(i l :i r )-s(i l :i r ) represents the total interaction between the vertices in the continuous interval [i l ,i r ] and the vertices outside the interval.

[0073] Step S203: constructing a dynamic programming table with the right endpoint index of the interval and the step number as the state, and determining the minimum partial Infomap entropy value in each state according to the obtained single-segment Infomap entropy value and the right endpoint index and the step number, and recording the precursor boundary index that produces the minimum value in the recursive process so as to backtrack, the minimum partial Infomap entropy value representing the minimum cumulative single-segment entropy that can be reached under the given step budget and covering the right endpoint;

[0074] wherein,, the minimum partial Infomap entropy value, the formula being:

[0075]

[0076] wherein, i r represents the right endpoint index of the interval, i prev represents the right endpoint index of the candidate precursor interval; t represents the step number, when i r = 1,

[0077] The dynamic programming table is constructed with the right endpoint index of the interval and the step number as the state, specifically including:

[0078] The right endpoint index i r of the interval is taken as a one-dimensional state, representing the termination position of the candidate topological association domain;

[0079] The step number t is taken as another-dimensional state, representing the discretization step number under the given coding budget;

[0080] In the dynamic programming recursive process, the minimum partial Infomap entropy value H(i r , t) corresponding to each state (i r , t) is updated in sequence, and the precursor index i prev that produces the optimal result is recorded, so as to finally obtain the optimal topological association domain division through backtracking.

[0081] Step S204: for each value in the value range of the step number, the corresponding complete Infomap entropy value is calculated according to the corresponding minimum partial Infomap entropy value and the global coding amount mapped by the step number; the global coding amount is determined by the step number and the step length ∈.

[0082] wherein, the complete Infomap entropy value, the formula being:

[0083]

[0084] wherein, H (G) represents the complete Infomap entropy value of graph G, t represents the step number, H (n, t) represents the minimum partial Infomap entropy value in the state of the right endpoint index n and the step number t, n is the number of genomic bins obtained by the selected resolution division, and t∈ is the global encoding amount.

[0085] Step S3: selecting the optimal topological correlation domain division according to the minimum value of the Infomap entropy values corresponding to all division domains, specifically including:

[0086] In the value range of the step number , the complete Infomap entropy value corresponding to each step number t is calculated in sequence, and the sizes are compared.

[0087] The step number t that makes the complete Infomap entropy value reach the minimum value is selected. * ;

[0088] Taking the state of the right endpoint index i r =n and the step number t=t * in the dynamic programming table as the backtracking starting point, wherein n is the number of genomic bins obtained by the selected resolution division.

[0089] Starting from the state (i r , t), the predecessor right endpoint index i prev recorded in the dynamic programming table is read, and the interval [i prev +1, i r ] is saved as an optimal topological correlation domain.

[0090] The backtracking state is updated as i r ←i prev , and the step number is repeated until the starting position of the genome is covered.

[0091] The continuous interval sequence obtained by backtracking is arranged in the linear order of the genome, that is, the final optimal topological correlation domain division result is obtained.

[0092] The embodiment also includes calculating the division structure approximation error ΔI according to the optimal topological correlation domain division result, and the expression is:

[0093]

[0094] Wherein, d is the number of regions of the selected division structure, ΔI is the approximation error of the selected division structure, and is used to measure the deviation between the division result and the ideal global encoding amount.

[0095] By calculating ΔI, the approximation degree of the optimal partition result can be quantified, so as to evaluate the accuracy of the TAD partition in information coding, and ensure that the partition result can reflect the interaction enrichment in the domain and maintain the overall consistency of the global interaction structure.

[0096] Embodiment 2:

[0097] The embodiment also provides a topologically associated domain partition system based on chromatin contact maps, comprising:

[0098] The data acquisition and preprocessing module is used for acquiring chromatin contact map data, including large-scale chromatin interaction times or intensities under multiple resolutions, corresponding chromosome information and genome region index information; performing quality filtering (for example, based on MAPQ >= 30) and normalization processing (for example, Knight-Ruiz algorithm) on the original data to generate a normalized contact matrix.

[0099] The Infomap entropy calculation module is used for calculating single-section Infomap entropy values of any possible continuous interval in the normalized contact matrix, and discretizing the external interaction amount according to the step parameter ∈, to provide input for dynamic programming.

[0100] The dynamic programming and optimal partition module is used for constructing a dynamic programming table with right endpoint index and step number as states, calculating minimum partial Infomap entropy values by recursion, and recording precursor boundary indexes; calculating complete Infomap entropy values for all step number values, and selecting a step number corresponding to the minimum value as an optimal partition basis, to obtain an optimal topologically associated domain partition result by backtracking.

[0101] The partition error calculation module is used for calculating the approximation error of the partition structure according to the optimal partition result, to evaluate the accuracy and reliability of the partition result in global information coding.

[0102] Through the cooperative work of the above modules, the system can automatically complete the whole process from chromatin contact map data to optimal topologically associated domain partition and error evaluation, reduce the dependence on artificial threshold selection and trained models, and improve the partition accuracy, repeatability and calculation efficiency.

[0103] Embodiment 3:

[0104] The following is a specific embodiment of the present application:

[0105] The given original contact count matrix is as follows:

[0106] 0 80 74 70 11 12 16 11 80 0 76 76 16 11 18 11 74 76 0 78 11 19 10 11 70 76 78 0 17 17 16 17 11 16 11 17 0 79 80 78 12 11 19 17 79 0 70 71 16 18 10 16 80 70 0 71 11 11 11 17 78 71 71 0

[0107] The interaction between eight genomic bins with serial numbers 0-7 is represented, and the number represents the interaction times, and the step size ∈ is set to 2.

[0108] The contact matrix A obtained by preprocessing is given as follows:

[0109] 0 0.6207 0.5843 0.5187 0.0862 0.0913 0.1215 0.0866 0.6207 0 0.5685 0.5334 0.1188 0.0793 0.1295 0.0821 0.5843 0.5685 0 0.5684 0.0848 0.1422 0.0747 0.0852 0.5187 0.5334 0.5684 0 0.1235 0.1199 0.1126 0.1241 0.0862 0.1188 0.0848 0.1235 0 0.5645 0.5703 0.5768 0.0913 0.0793 0.1422 0.1199 0.5645 0 0.5177 0.5446 0.1215 0.1295 0.0747 0.1126 0.5703 0.5177 0 0.5405 0.0866 0.0821 0.0852 0.1241 0.5768 0.5446 0.5405 0

[0110] For ease of presentation, the data provided in this example are kept to four decimal places.

[0111] Thus, the calculation of t is as follows: The value range of t is: t∈{0, 1, …, 17};

[0112] For all i l and i r , the value of Q(i l ; i r ) is calculated as follows:

[0113]

[0114]

[0115] Since i l ≤ i r , the lower triangular matrix has no available values, which are represented as NaN.

[0116] In step S203, the value of all H(i r , t) is calculated as follows:

[0117] inf inf inf inf inf inf inf 4.0654 inf inf inf inf inf inf inf inf 0.252 inf inf 1.8707 inf inf 3.8147 inf inf 0.6675 1.2484 inf 2.4851 3.1354 inf inf inf 0.5067 inf inf 2.1246 inf inf 3.6962 inf inf 0.908 1.4825 inf 2.5176 3.0808 inf inf inf 0.7586 1.3188 1.952 2.3706 3.1132 3.6705 inf inf inf 1.158 1.7363 2.5688 2.7527 3.3079 inf inf inf 1.0095 1.5726 2.1294 2.6175 3.1442 inf inf inf inf 1.4119 1.9657 2.5288 2.9835 inf inf inf inf 1.2633 1.8049 2.3644 2.8349 inf inf inf inf inf 1.6564 2.2007 2.7579 inf inf inf inf inf 1.5094 2.04 2.5959 inf inf inf inf inf inf 1.8914 2.4351 inf inf inf inf inf inf 1.7563 2.2837 inf inf inf inf inf inf inf 2.1351 inf inf inf inf inf inf inf inf

[0118] For ease of presentation, the row and column of H(i r , t) are exchanged here. According to the actual situation of the input matrix, there may be a case that no matter how the division is made, the current t cannot be obtained. In order to facilitate comparison, the H corresponding to the case that cannot exist is set to inf (i.e., positive infinity).

[0119] The value corresponding to the index i prev is as follows:

[0120]

[0121] Wherein, ‘-’ represents that the corresponding index does not exist.

[0122] In step S203, according to The minimum case is marked as:

[0123] inf inf inf inf inf inf inf 4.0654 inf inf inf inf inf inf inf inf 0.252 inf inf 1.8707 inf inf 3.8147 inf inf 0.6675 1.2484 inf 2.4851 3.1354 inf inf inf 0.5067 inf inf 2.1246 inf inf 3.6962 inf inf 0.908 1.4825 inf 2.5176 3.0808 inf inf inf 0.7586 1.3188 1.952 2.3706 3.1132 3.6705 inf inf inf 1.158 1.7363 2.5688 2.7527 3.3079 inf inf inf 1.0095 1.5726 2.1294 2.6175 3.1442 inf inf inf inf 1.4119 1.9657 2.5288 2.9835 inf inf inf inf 1.2633 1.8049 2.3644 2.8349 inf inf inf inf inf 1.6564 2.2007 2.7579 inf inf inf inf inf 1.5094 2.04 2.5959 inf inf inf inf inf inf 1.8914 2.4351 inf inf inf inf inf inf 1.7563 2.2837 inf inf inf inf inf inf inf 2.1351 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf

[0124] The corresponding index case is as follows:

[0125] - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 0 - - 3 - - 3 - - 0 0 - 3 3 - - - 1 1 0 4 4 4 - - - 1 3 0 4 3 - - - 2 3 3 5 3 - - - - 3 3 3 3 - - - - 3 3 4 3 - - - - - 3 4 4 - - - - - 4 4 5 - - - - - - 4 5 - - - - - - 5 6 - - - - - - - 6 - - - - - - - - ​

[0126] Thus the division result is [0, 1, 2, 3], [4, 5, 6, 7], i.e. the points corresponding to 0-3 and 4-7 are in two divisions respectively.

[0127] If the above functions are realized in the form of software function units and sold or used as independent products, they can be stored in a computer readable storage medium. Based on this understanding, the technical solutions of the present application or the parts of the present application that essentially contribute to the prior art or the parts of the technical solutions can be embodied in the form of software products. The computer software product is stored in a storage medium and includes a plurality of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in the various embodiments of the present application. The aforementioned storage medium includes a U disk, a mobile hard disk, a read-only memory (ROM, Read-Only Memory), a random access memory (RAM, Random Access Memory), a magnetic disk or an optical disk, and various media that can store program codes.

[0128] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art can easily think of various equivalent modifications or replacements within the technical scope disclosed by the present application, and these modifications or replacements should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A method for partitioning topologically related structural domains based on chromatin contact maps, characterized in that, Includes the following steps: Acquire chromatin contact map data and preprocess it to obtain a normalized contact matrix; The Infomap entropy values ​​corresponding to all partitioned structural domains are solved using the normalized contact matrix obtained after preprocessing. The optimal topological association domain partition is selected based on the minimum Infomap entropy value corresponding to all partitioned domains.

2. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 1, characterized in that, The chromatin contact map data includes the number or intensity of large-scale chromatin interactions at multiple resolutions, corresponding chromosome information, and index information of specific chromosomes or specific genomic regions.

3. A method for partitioning topologically related structural domains based on chromatin contact maps according to claim 1 or 2, characterized in that, The acquisition and preprocessing of chromatin contact map data specifically includes: The original contact count matrix is ​​extracted or constructed from the contact map data, wherein the value of each element in the original contact count matrix is ​​determined based on the number or intensity of large-scale chromatin interactions under multi-resolution conditions. The original contact count matrix is ​​subjected to quality filtering, which includes filtering and removing or marking low-quality or abnormal contact records based on the comparison quality value MAPQ, with a filtering threshold of MAPQ≥30; The original contact count matrix after quality filtering is subjected to balanced normalization processing, which is performed using the Knight-Ruiz algorithm. The normalization result is the normalized contact matrix A.

4. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 1, characterized in that, The process of using the preprocessed normalized contact matrix to solve for the Infomap entropy values ​​corresponding to all partitioned structural domains specifically includes: The range of step size is determined based on the given step size parameter ∈ and the total interaction quantity of the normalized contact matrix A. The step size is a non-negative integer and is used to discretize the continuous external interaction quantity according to the step size parameter ∈. For any possible continuous interval in the normalized contact matrix A, calculate the single-segment Infomap entropy value of that interval in turn; A dynamic programming table is constructed with the right endpoint index and the number of steps as the state. The minimum part of the Infomap entropy value in each state is determined by recursively using the obtained single-segment Infomap entropy value according to the right endpoint index and the number of steps. During the recursion process, the predecessor boundary index that generates the minimum value is recorded for backtracking. The minimum part of the Infomap entropy value represents the minimum cumulative single-segment entropy that can be achieved by covering the right endpoint under a given number of steps. For each value within the range of the step number, the corresponding complete Infomap entropy value is calculated based on the minimum Infomap entropy value and the global encoding amount obtained by mapping from that step number; the global encoding amount is determined by the step number and step size ∈.

5. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 4, characterized in that, The range of step count values ​​is determined based on the given step size parameter ∈ and the total interaction quantity of the normalized contact matrix A, using the following formula: Where t represents the number of steps, and its value is a non-negative integer; The total interaction quantity of graph G is equal to the sum of all elements in the normalized contact matrix A; G is an undirected weighted graph corresponding to the normalized contact matrix A, and the nodes in the graph are the genome bins obtained according to the selected resolution. ij The element in the i-th row and j-th column of the normalized contact matrix A represents the normalized interaction strength or weight between the i-th genomic bin and the j-th genomic bin, and n is the number of genomic bins obtained by the selected resolution.

6. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 4, characterized in that, The continuous interval refers to an index interval [i] that is adjacent in the row and column indices of the normalized contact matrix A according to the linear order of the genome. l i r ], where i l i r Let i be the left and right endpoint indices of the interval, respectively, satisfying 1 ≤ i l ≤i r ≤n; Each continuous interval corresponds to a candidate topological association domain on the genome, used to evaluate the local interaction tightness between sites within the interval and the interaction relationship with the outside of the interval.

7. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 4, characterized in that, The formula for calculating the single-segment Infomap entropy value of any possible continuous interval in the normalized contact matrix A is as follows: Where, Q(i) l i r ) represents a continuous interval [i l i r The single-segment Infomap entropy value of graph G is used to measure the information encoding contribution of the candidate structural domain. vol(G) represents the total interaction amount of graph G. Represents a continuous interval [i] l i r The total number of interactions between interior vertices and all other vertices in the entire graph. Represents a continuous interval [i] l i r The total number of interactions between internal vertices, g(i) l :i r ) = vol(i l :i r )-s(i l :i r ) represents a continuous interval [i l i r The total number of interactions between vertices inside the interval and vertices outside the interval.

8. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 4, characterized in that, The minimum entropy value of the Infomap is calculated using the following formula: Among them, i r Indicates the index of the right endpoint of the interval, i prev Indicates the right endpoint index of the candidate predecessor interval; t represents the step number, when i r When = 1, The construction of the dynamic programming table using the right endpoint index of the interval and the number of steps as the state specifically includes: Index i by the right endpoint of the interval r As a one-dimensional state, it represents the termination position of the candidate topological association structure domain; The number of steps t is used as another dimension of the state, representing the number of discretization steps under a given coding budget; During the dynamic programming recursion process, each state (i) is updated sequentially. r The minimum entropy value H(i) of the Infomap corresponding to t) r ,t), and record the predecessor index i that produced the optimal result. prev This is so that the optimal topological association structure domain partitioning can be obtained through backtracking.

9. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 4, characterized in that, The formula for the complete Infomap entropy value is: in, Let H(n,t) represent the complete Infomap entropy value of graph G, where t represents the step number, H(n,t) represents the minimum Infomap entropy value in the state with right endpoint index n and step number t, where n is the number of genome bins obtained by the selected resolution, and t ∈ represents the global coding amount.

10. The method for partitioning topologically related structural domains based on chromatin contact maps according to claim 1, characterized in that, The step of selecting the optimal topological association structural domain partition based on the minimum Infomap entropy value corresponding to all partitioned structural domains specifically includes: Within the range of step values Within this process, the entropy value of the complete Infomap corresponding to each step number t is calculated sequentially. And compare their sizes; Choose to make the full Infomap entropy value Number of steps t to obtain the minimum value * ; Let the right endpoint index of the dynamic programming table be i r =n, number of steps is t=t * The state is used as the starting point for backtracking, where n is the number of genome bins obtained by the selected resolution. From state (i r Starting from t), read the index i of the right endpoint of the predecessor record in the dynamic programming table. prev , the interval [i prev +1,i r This is stored as an optimal topological association structure domain; Update the backtracking state, and set i r ←i prev Steps Repeat the backtracking steps until the start of the genome is covered; The continuous interval sequences obtained by backtracking are arranged in a linear order according to the genome, which constitutes the final optimal topological association domain partitioning result.

Images