Interaction DNA fragment identification method based on DNA methylation correlation coefficient

Abstract

The invention relates to an interaction DNA fragment recognition method based on a DNA methylation correlation coefficient, which recognizes interaction DNA fragments based on high correlation of methylation rates of two interaction DNA fragments, and can recognize DNA fragments possibly having interaction relationship in all experimental samples only by using DNA methylation sequencing data, so that the recognition efficiency is improved, and the detection cost is reduced. The interaction way comprises space folding, DNA cyclization, transcription factor combination, RNA interference and the like, the technical shortages that Hi-C can only identify a chromatin space structure and Chip-seq can only detect a DNA fragment combined with a specific transcription factor are made up, and the experiment and calculation cost is saved.

Description

Technical Field

[0001] The present invention relates to the field of bioinformatics, and in particular to a method for identifying interacting DNA fragments based on DNA methylation correlation coefficients. Background Art

[0002] DNA fragments in the cell nucleus can interact with distant DNA fragments of different chromosomes or the same chromosome through processes such as three-dimensional folding, DNA cyclization, transcription factor binding, and RNA interference, thereby affecting the expression of genes that are not similar in DNA sequence.

[0003] Commonly used methods for identifying interacting DNA fragments include high-throughput chromosome conformation capture sequencing (Hi-C), chromatin immunoprecipitation sequencing (Chip-Seq), and hypergeometric optimization of motif enrichment (Homer). Among them, Hi-C technology is based on the DNA fragment contact matrix obtained by sequencing, and then uses a variety of bioinformatics algorithms to identify the three-dimensional spatial structure of cell chromosomes, such as A / B compartments, topologically associated domains (TADs), and DNA loops. However, when using this technology to identify interacting DNA loops, there are extremely high requirements for sequencing depth and computational complexity, and it is difficult to identify DNA fragments that interact through transcription factors or RNA. Chip-seq can only be used to explore the binding of specific transcription factors to DNA fragments. If you want to explore the binding of all transcription factors to DNA fragments in specific cells, the experimental cost is too high. Homer can only identify possible DNA binding sites of transcription factors through mathematical algorithms. Therefore, how to identify all possible interacting DNA fragments in a specific cell or tissue with less experimental and computational costs is a current technical difficulty.

[0004] DNA methylation is a common epigenetic modification, mediated by DNA methyltransferases (DNMTs), which covalently bind the methyl group in S-adenosylmethionine to the 5'-terminal cytosine in CpG dinucleotides. A large number of studies have shown that the binding of transcription factors, RNA and DNA fragments is affected by the DNA methylation level of the fragment. The sequencing results of Methyl-HiC technology, which combines DNA methylation sequencing and Hi-C, also show that DNA fragments that are close or intersecting in space have synergistic changes in DNA methylation status. In addition, changes in DNA methylation can also reflect changes in gene expression levels. Currently, there is a large amount of DNA methylation sequencing data in public databases. Summary of the invention

[0005] The technical problem to be solved by the present invention is how to identify all possible interacting DNA fragments existing in specific cells or tissues with less experimental and computational costs.

[0006] The present invention provides a method for identifying interacting DNA fragments based on DNA methylation correlation coefficient, comprising: Step 1. Perform DNA methylation sequencing on multiple DNA fragment samples to obtain multiple bed sample files; Step 2. Integrate multiple bed sample files into one bed sample file according to the common CpG sites; Step 3. Calculate the Pearson correlation coefficient of the methylation rate of each CpG site with other CpG sites from the bed sample files integrated in step 2; and preset the Pearson correlation coefficient threshold , the maximum distance between adjacent CpG sites in a DNA fragment and the minimum number of CpG sites that a DNA fragment needs to contain ; Step 4. Set the screening conditions, and determine whether there is an interaction relationship between all DNA fragments according to the screening conditions, wherein the screening conditions are: , There are interactions between the DNA fragments that meet the above screening conditions, and all DNA fragments with interactions are integrated to obtain a set of interaction bed files.

[0007] Compared with the prior art, the present invention has the following advantages: the present invention identifies interacting DNA fragments based on the high correlation between the methylation rates of two interacting DNA fragments. Only by using DNA methylation sequencing data, DNA fragments that may have interacting relationships in all experimental samples can be identified. The interaction pathways include spatial folding, DNA cyclization, transcription factor binding, RNA interference, etc., which makes up for the technical shortcomings that Hi-C can only identify the spatial structure of chromatin and Chip-seq can only measure DNA fragments bound to specific transcription factors, saving experimental and computing costs.

[0008] In a possible implementation, the bed sample file obtained in step 1 includes the DNA fragment number, the CpG site start position, the CpG site end position and the CpG site methylation rate β.

[0009] In a possible implementation, the expression formula for calculating the Pearson correlation coefficient of the methylation rate of each CpG site and other CpG sites in step 3 is: , In the formula, is the number of bed sample files, and Respectively represent Methylation rates of two CpG sites in a bed sample file.

[0010] In a possible implementation, the Pearson correlation coefficient threshold in step 3 is By specifying the significance level of the t-test , and then according to the number of samples n Calculation, the calculation formula is: , , , In the formula, is the probability density function of the t distribution, based on the probability density function of the t distribution Find the cumulative probability of the t distribution The critical t value , then based on Calculate the correlation coefficient threshold .

[0011] In a possible implementation, the interaction bed file in step 4 includes the DNA fragment number, the starting position of the CpG site, the ending position of the CpG site, the number of CpG sites contained in the DNA fragment, and a set of other DNA fragments that have an interaction relationship with the DNA fragment.

[0012] In a possible implementation, the interaction relationship between the DNA fragments in step 4 includes a positive feedback relationship and a negative feedback relationship. If the Pearson correlation coefficient is a positive number, the DNA fragments are in a positive feedback relationship. If the Pearson correlation coefficient is a negative number, the DNA fragments are in a negative feedback relationship.

[0013] In a possible implementation, the method further includes step 5. Obtaining two groups of interaction bed files through steps 1 to 4, integrating the two groups of interaction bed files according to the common DNA fragments, and comparing the differences between other DNA fragment sets with the same DNA fragments in the two groups of interaction bed files, and marking the difference types between other DNA fragment sets with the same DNA fragments in the two groups of interaction bed files, integrating and outputting the difference interaction bed files.

[0014] In a possible implementation, the differential interaction bed file outputted in step 5 includes the DNA fragment codes of the two DNA fragments that interact with each other in the two sets of interaction bed files, the starting positions of the CpG sites of the two DNA fragments, the ending positions of the CpG sites of the two DNA fragments, and the difference types between other DNA fragment sets in which the same DNA fragments in the two sets of interaction bed files have an interaction relationship.

[0015] In a possible implementation, in step 5, by comparing the differences between the sets of other DNA fragments in which the same DNA fragments in the two sets of interaction bed files have an interaction relationship, marking the difference types between the sets of other DNA fragments in which the same DNA fragments in the two sets of interaction bed files have an interaction relationship specifically includes: Step 501. The two groups of interaction bed files are group A interaction bed files and group B interaction bed files; Step 502. Intersect the set of other DNA fragments in the interaction bed file of group A and the interaction bed file of group B, and obtain three sets. Set 1 is the intersection of the set of other DNA fragments in the interaction bed file of group A and the interaction bed file of group B, and the set 2 is the subset of the set of other DNA fragments in the interaction bed file of group A, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B. Step 503. Set 1 is the interaction relationship between the same two DNA fragments in the interaction bed files of group A and group B, and the difference type is marked as interaction relationship A-interaction relationship B; The difference type of set 2 is labeled as interaction relationship A-irrelevant; The difference type of set three is labeled as irrelevant-interaction relationship B. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 Flow chart of specific embodiment 1 of the present invention; Figure 2 This is a flow chart of the interaction bed file of group A obtained through steps 1 to 4 in specific embodiment 2 of the present invention; Figure 3 This is a flow chart of the interaction bed file of group B obtained through steps 1 to 4 in specific embodiment 2 of the present invention; Figure 4 This is a flowchart of a specific embodiment 2 of the present invention, in which the interaction bed file of group A and the interaction bed file of group B are subjected to step 5 to obtain a difference interaction bed file. DETAILED DESCRIPTION

[0017] First, those skilled in the art should understand that these implementations are only used to explain the technical principles of the embodiments of the present invention, and are not intended to limit the protection scope of the embodiments of the present invention. Those skilled in the art can make adjustments to them as needed to adapt to specific application scenarios.

[0018] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments. Specific embodiment 1 Figure 1 Part a of the present invention is a specific embodiment 1 step 2. n A set of bed sample files after integrating bed sample files; Figure 1 Part b is the site matrix of calculating the Pearson correlation coefficient of the methylation rate of each CpG site with other CpG sites in step 3 of specific embodiment 1 of the present invention; Figure 1 Part c is the interaction bed file of step 4 of specific embodiment 1 of the present invention, which integrates all DNA fragments with interacting relationships.

[0020] See also Figure 1 As shown, the embodiment of the present invention discloses a method for identifying interacting DNA fragments based on DNA methylation correlation coefficient, comprising: Step 1. n DNA fragment samples are subjected to DNA methylation sequencing to obtain multiple bed sample files. n=14; Each bed sample file includes the DNA fragment number (chrom), the starting position of the CpG site (start), the ending position of the CpG site (End), and the methylation rate of the CpG site β; Here, β is the percentage of the number of methylated reads at the current CpG site to all the reads covering the current CpG site.

[0021] Step 2. n The bed sample files are integrated according to the common CpG sites into a bed sample file with n+3 columns; the integrated bed sample file contains the DNA fragment number, the starting position of the CpG site, the ending position of the CpG site, the methylation rate β of DNA fragment sample 1 at the CpG site, the methylation rate β of DNA fragment sample 2 at the CpG site...DNA fragment sample n The methylation rate of CpG sites is β; the integrated bed sample file is as follows Figure 1 As shown in part a of the figure, Chrom in the first row is the number of the chromosome in the DNA fragment, Start is the starting position of the CpG site, End is the ending position of the CpG site, SA1, SA2, SA3...SA14 are the numbers of the DNA fragment samples, and 1, 2, 3...m in the first column are the numbers of the CpG sites; Step 3. Calculate the Pearson correlation coefficient of the methylation rate of each CpG site with other CpG sites from the bed sample files integrated in step 2; and preset the Pearson correlation coefficient threshold , the maximum distance between adjacent CpG sites in a DNA fragment and the minimum number of CpG sites that a DNA fragment needs to contain ; In this specific embodiment, the Pearson correlation coefficient threshold , the maximum distance between adjacent CpG sites in a DNA fragment and the minimum number of CpG sites that a DNA fragment needs to contain Specified based on staff experience, ,in: The expression formula of the Pearson correlation coefficient of the methylation rate of each CpG site and other CpG sites is: , In the formula, is the number of bed sample files, and Respectively represent Methylation rates of two CpG sites in a bed sample file.

[0022] The interaction relationship between DNA fragments includes positive feedback relationship and negative feedback relationship. If the Pearson correlation coefficient is a positive number, the DNA fragments are in a positive feedback relationship. If the Pearson correlation coefficient is a negative number, the DNA fragments are in a negative feedback relationship. At the same time, the correlation coefficients of each CpG site with other CpG sites form an m×m CpG site matrix, such as Figure 1 As shown in part b of , due to the symmetry of the correlation coefficient matrix, only half of the diagonal of the m×m CpG site matrix needs to be calculated to reduce the amount of calculation and avoid repeated calculations.

[0023] Step 4. Set the screening conditions, and determine whether there is an interaction relationship between all DNA fragments according to the screening conditions, wherein the screening conditions are: , There is an interaction relationship between the DNA fragments that meet the above screening conditions, and all the DNA fragments that have an interaction relationship are integrated to obtain a set of interaction bed files, wherein the interaction bed file includes the DNA fragment number, the starting position of the CpG site, the ending position of the CpG site, the number of CpG sites contained in the DNA fragment, and a set of other DNA fragments that have an interaction relationship with the DNA fragment. In this specific embodiment, the set of other DNA fragments that have an interaction relationship with the DNA fragment is expressed as "P / n-kCpG-DNA fragment position: start-end", P represents a positive feedback relationship, n represents a negative feedback relationship, and kCpG represents the number of CpG sites on the DNA fragment; Figure 1 As shown in part c, chr1, 10, 302, 3, p-3-chr2: 30-342, n-3-chr4: 210-342 indicates that there are 3 CpG sites at positions 10-30 of DNA fragment 1, and there is a positive feedback relationship with the 3 CpG sites at positions 30-342 of DNA fragment 2; there is a negative feedback relationship with the 3 CpG sites at positions 210-342 of DNA fragment 4; chr2, 30, 292, 4, n-3-chr4: 210-342, indicates that there are 4 CpG sites at positions 30-292 of DNA fragment 2, and there is a negative feedback relationship with the 3 CpG sites at positions 210-342 of DNA fragment 4.

[0024] In addition, regarding the Pearson correlation coefficient threshold In this specific embodiment, the significance level of t-test is specified. , and then according to the number of DNA fragment samples n Calculation, the calculation formula is: , , , In the formula, is the probability density function of the t distribution, based on the probability density function of the t distribution Find the cumulative probability of the t distribution The critical t value , then based on Calculate the correlation coefficient threshold . Specific embodiment 2 Figure 2 Part a in it is a set of bed sample files obtained by performing step 2 in specific embodiment 2 of the present invention; Figure 2 Part b is a set of site matrices obtained by performing step 3 in specific embodiment 2 of the present invention to obtain the Pearson correlation coefficient; Figure 2 Part c is another set of interaction bed files obtained by performing step 4 in specific embodiment 2 of the present invention.

[0026] Figure 3 Part a in the figure is another set of bed sample files obtained by performing step 2 in specific embodiment 2 of the present invention; Figure 3 Part b is another set of site matrices obtained by performing step 3 in specific embodiment 2 of the present invention to obtain the Pearson correlation coefficient; Figure 3 Part c is another set of interaction bed files obtained by performing step 4 in specific embodiment 2 of the present invention.

[0027] This specific embodiment adds the following on the basis of Embodiment 1: Step 5. Obtain two groups of interaction bed files by implementing step 1 to step 4 in specific implementation 1, integrate the two groups of interaction bed files according to the common DNA fragments, and integrate and output the difference interaction bed files by comparing the differences between other DNA fragment sets with the same DNA fragments in the two groups of interaction bed files, and marking the difference types between other DNA fragment sets with the same DNA fragments in the two groups of interaction bed files; the difference interaction bed files contain the DNA fragment codes of the two DNA fragments that interact in the two groups of interaction bed files, the starting positions of the CpG sites of the two DNA fragments, the ending positions of the CpG sites of the two DNA fragments, and the difference types between other DNA fragment sets with the same DNA fragments in the two groups of interaction bed files, such as Figure 4 As shown; Regarding the differences between the sets of other DNA fragments that have an interaction relationship between the same DNA fragments in two sets of interaction bed files, the types of differences between the sets of other DNA fragments that have an interaction relationship between the same DNA fragments in two sets of interaction bed files are specifically marked as follows: Step 501. The two groups of interaction bed files are group A interaction bed files and group B interaction bed files; Step 502. Intersect the set of other DNA fragments in the interaction bed file of group A and the interaction bed file of group B, and obtain three sets. Set 1 is the intersection of the set of other DNA fragments in the interaction bed file of group A and the interaction bed file of group B, and the set 2 is the subset of the set of other DNA fragments in the interaction bed file of group A, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B. Step 503. Set 1 is the interaction relationship between the same two DNA fragments in the interaction bed files of group A and group B, and the difference type is marked as interaction relationship A-interaction relationship B; Specifically, if Figure 2 and Figure 3 As shown, first, two groups are used n DNA fragment samples were subjected to step 1 to obtain 2 groups n bed sample files; 2 groups n The bed sample files are integrated in step 2 to obtain group A bed sample files and group B bed sample files, such as Figure 2 Part a and Figure 3 As shown in part a of ; then, the bed sample files of group A and group B are subjected to step 3 to obtain the Pearson correlation coefficient matrix of the bed sample files of group A and the coefficient matrix of the bed sample files of group B, as shown in Figure 2 Part b and Figure 3 As shown in part b of ; then, the Pearson correlation coefficient matrix of the bed sample file of group A and the coefficient matrix of the bed sample file of group B are obtained through step 4 to obtain the interaction bed file of group A and the interaction bed file of group B, as shown in Figure 2 Part c and Figure 3 As shown in part c of the figure, group A: chr1, 10, 302, 3, p-3-chr2: 30-342, n-3-chr4: 210-342; group B: chr1, 10, 302, 3, p-3-chr3: 110-302; after integration in step 5, as shown in Figure 4As shown, set 1 is an empty set, set 2 is chr1, 10, 302, 3, p-3-chr2: 30-342; chr1, 10, 302, 3, n-3-chr4: 210-342, and set 3 is chr1, 10, 302, 3, p-3-chr3: 110-302; among them, chr1, 10, 302, chr2, 30, 342, p-no, represents the 10-302 position of chromosome 1 and the 30-342 position of chromosome 2 The difference type between the interaction bed file of group A and the interaction bed file of group B is positive feedback relationship-irrelevance; no-p indicates a positive feedback relationship from irrelevance of group A to group B, no-n indicates a negative feedback relationship from irrelevance of group A to group B, pn indicates a positive feedback relationship from group A to negative feedback relationship of group B, p-no indicates a positive feedback relationship from group A to irrelevance of group B, np indicates a negative feedback relationship from group A to positive feedback relationship of group B, and n-no indicates a negative feedback relationship from group A to irrelevance of group B.

[0028] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed by the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims

1. A method for identifying interacting DNA fragments based on DNA methylation correlation coefficient, comprising: Step 1. Perform DNA methylation sequencing on multiple DNA fragment samples to obtain multiple bed sample files; Step 2. Integrate multiple bed sample files into one bed sample file according to the common CpG sites; Step 3. Calculate the Pearson correlation coefficient of the methylation rate of each CpG site with other CpG sites from the bed sample files integrated in step 2; And preset the Pearson correlation coefficient threshold , the maximum distance between adjacent CpG sites in a DNA fragment and the minimum number of CpG sites that a DNA fragment needs to contain ; Step 4. Set the screening conditions, and determine whether there is an interaction relationship between all DNA fragments according to the screening conditions, wherein the screening conditions are: ， There are interactions between the DNA fragments that meet the above screening conditions, and all DNA fragments with interactions are integrated to obtain a set of interaction bed files.

2. A method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 1, characterized in that: The bed sample file obtained in step 1 includes the DNA fragment number, the starting position of the CpG site, the ending position of the CpG site and the methylation rate β of the CpG site.

3. A method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 1, characterized in that: The expression formula for calculating the Pearson correlation coefficient of the methylation rate of each CpG site and other CpG sites in step 3 is: ， In the formula, is the number of bed sample files, and Respectively represent Methylation rates of two CpG sites in a bed sample file.

4. The method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 1, characterized in that: The Pearson correlation coefficient threshold in step 3 By specifying the significance level of the t-test , and then according to the number of DNA fragment samples n Calculation, the calculation formula is: ， ， ， In the formula, is the probability density function of the t distribution, based on the probability density function of the t distribution Find the cumulative probability of the t distribution The critical t value , then based on Calculate the correlation coefficient threshold .

5. The method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 1, characterized in that: The interaction bed file in step 4 includes the DNA fragment number, the starting position of the CpG site, the ending position of the CpG site, the number of CpG sites contained in the DNA fragment, and a set of other DNA fragments that have an interaction relationship with the DNA fragment.

6. A method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 5, characterized in that: The interaction relationship between the DNA fragments in step 4 includes a positive feedback relationship and a negative feedback relationship. If the Pearson correlation coefficient is a positive number, the DNA fragments are in a positive feedback relationship. If the Pearson correlation coefficient is a negative number, the DNA fragments are in a negative feedback relationship.

7. The method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 1, characterized in that: The method also includes step 5. obtaining two groups of interaction bed files through steps 1 to 4, integrating the two groups of interaction bed files according to the common DNA fragments, and comparing the differences between other DNA fragment sets with interaction relationships between the same DNA fragments in the two groups of interaction bed files, and marking the difference types between other DNA fragment sets with interaction relationships between the same DNA fragments in the two groups of interaction bed files, integrating and outputting the difference interaction bed files.

8. The method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 7, characterized in that: The differential interaction bed files output in step 5 include the DNA fragment codes of the two DNA fragments that interact with each other in the two sets of interaction bed files, the starting positions of the CpG sites of the two DNA fragments, the ending positions of the CpG sites of the two DNA fragments, and the difference types between other DNA fragment sets that have an interaction relationship with the same DNA fragments in the two sets of interaction bed files.

9. The method for identifying interacting DNA fragments based on DNA methylation correlation coefficient according to claim 7, characterized in that: In step 5, by comparing the differences between the sets of other DNA fragments in which the same DNA fragments in the two sets of interaction bed files have an interaction relationship, marking the difference types between the sets of other DNA fragments in which the same DNA fragments in the two sets of interaction bed files have an interaction relationship specifically includes: Step 501. The two groups of interaction bed files are group A interaction bed files and group B interaction bed files; Step 502. Intersect the set of other DNA fragments in the interaction bed file of group A and the interaction bed file of group B, and obtain three sets. Set 1 is the intersection of the set of other DNA fragments in the interaction bed file of group A and the interaction bed file of group B, and the set 2 is the subset of the set of other DNA fragments in the interaction bed file of group A, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B, and the set 3 is the subset of the set of other DNA fragments in the interaction bed file of group B. Step 503. Set 1 is the interaction relationship between the same two DNA fragments in the interaction bed files of group A and group B, and the difference type is marked as interaction relationship A-interaction relationship B; The difference type of set 2 is labeled as interaction relationship A-irrelevant; The difference type of set three is labeled as irrelevant-interaction relationship B.

Images