Methods, devices, media and programs for chromatin openness data mining based on long-read single-molecule sequencing data

Abstract

The present invention relates to a method, device, medium, and program for chromatin accessibility data mining based on long-read single-molecule sequencing data. The method comprises: extracting 6mA position information from the offline data to identify the positions of nucleosomes and MSPs; calculating the integrated FIRE signal value based on the estimated precision value assigned to each MSP to determine the FIRE peaks for each sample; merging the FIRE peaks between replicate samples within a group based on the physical position information of the FIRE peaks to obtain consistent FIRE peaks within each sample group; and merging the FIRE peaks of samples within a group based on the physical position information of the FIRE peaks to determine FIRE peaks that are unique between sample groups. In this way, chromatin accessibility information is mined within and between sample groups.

Description

Technical Field

[0001] The present invention generally relates to the processing of biological information, and in particular, to methods, computing devices, computer storage media, and computer program products for chromatin openness data mining based on long-read single-molecule sequencing data. Background Art

[0002] Long-read single-molecule sequencing data (such as, but not limited to, Fiber-seq data) can simultaneously provide information on three dimensions: DNA sequence, chromatin accessibility, and 5mC methylation. The following uses Fiber-seq data as an example to illustrate traditional methods for chromatin accessibility data mining based on long-read single-molecule sequencing data. Traditional chromatin accessibility data mining solutions for Fiber-seq data primarily include: Fiberseq-qc software, FiberTools software, and the Inferred Regulatory Elements from Fiber-seq Data (FIRE) software. Fiberseq-qc software primarily assesses the quality of 6mA methylation in Fiber-seq data. FiberTools software primarily extracts 6mA sites and MSP regions from sequencer-generated Fiber-seq data. FIRE software primarily identifies FIRE peaks (Fiber-seq Inferred Regulatory Element peaks, or FIRE peaks). The above traditional regulatory elements inferred based on Fiber-seq data can only complete specific chromatin accessibility analysis, have relatively simple functions, and can only be analyzed for a single sample.

[0003] In summary, the shortcomings of traditional chromatin accessibility data mining methods based on long-read single-molecule sequencing data are that they can only complete specific chromatin accessibility analysis, have relatively simple functions, and can only mine chromatin accessibility information of a single sample. Summary of the Invention

[0004] The present invention provides a method, computing device, computer storage medium and computer program product for chromatin accessibility data mining based on long-read single-molecule sequencing data, which can mine chromatin accessibility information within and between sample groups.

[0005] According to a first aspect of the present invention, a method for chromatin openness data mining based on long-read single-molecule sequencing data is provided. The method includes: extracting position information of methylation modification of the nitrogen atom at position 6 of adenine from the downstream data of long-read single-molecule sequencing, so as to identify the position information of nucleosomes and methyltransferase-sensitive regions based on the position information; calculating the integrated signal value of the regulatory element inferred based on Fiber-seq data according to the estimated accuracy value assigned to each methyltransferase-sensitive region, so as to determine the regulatory element peak value inferred based on Fiber-seq data for each sample; based on the determined physical position information of the regulatory element peak value inferred based on Fiber-seq data, merging the regulatory element peak value inferred based on Fiber-seq data between repeated samples within each sample group, so as to obtain the regulatory element peak value inferred based on Fiber-seq data with consistency within each sample group, each sample group including multiple samples; based on the physical position information of the regulatory element peak value inferred based on Fiber-seq data, merging the regulatory element peak value inferred based on Fiber-seq data of the samples within each sample group, so as to determine the regulatory element peak value inferred based on Fiber-seq data that is unique between multiple sample groups.

[0006] In some embodiments, the method for chromatin openness data mining based on long-read single-molecule sequencing data also includes: accumulating the supporting read lengths and unsupported read lengths of the regulatory element peaks inferred based on Fiber-seq data of the samples within each sample group, so as to determine the differences in the regulatory element peaks inferred based on Fiber-seq data between multiple sample groups based on the accumulated results of each of the multiple sample groups.

[0007] In some embodiments, obtaining the regulatory element peaks inferred based on Fiber-seq data for each sample group with intra-group consistency includes: merging the regulatory element peaks inferred based on Fiber-seq data between repeated samples in each sample group based on the physical location information of the regulatory element peaks inferred based on Fiber-seq data; retaining the regulatory element peaks inferred based on Fiber-seq data between samples whose overlapping length ratio is greater than a predetermined proportion as the regulatory element peaks inferred based on Fiber-seq data with intra-group consistency; in response to determining that the number of repeated samples in the group is greater than or equal to three, determining the regulatory element peaks inferred based on Fiber-seq data that appear in at least two samples among the retained regulatory element peaks inferred based on Fiber-seq data with intra-group consistency as the regulatory element peaks inferred based on Fiber-seq data with intra-group consistency; and annotating the regulatory element peaks inferred based on Fiber-seq data that are determined to be intra-group consistent.

[0008] In some embodiments, determining a group-specific regulatory element peak inferred based on Fiber-seq data between multiple sample groups includes: merging the regulatory element peak inferred based on Fiber-seq data within each sample group based on the physical location information of the regulatory element peak inferred based on Fiber-seq data, using the same method as the intra-group consistency regulatory element peak inferred based on Fiber-seq data (FIRE peaks); then determining whether the proportion of the intersection area of ​​the two sample groups in the current regulatory element peak inferred based on Fiber-seq data is greater than a predetermined intersection proportion threshold; in response to determining that the proportion of the two sample groups in the intersection area of ​​the current regulatory element peak inferred based on Fiber-seq data is greater than the predetermined intersection proportion threshold, determining that the current regulatory element peak inferred based on Fiber-seq data is a group-specific regulatory element peak inferred based on Fiber-seq data between the multiple sample groups; and performing functional and pathway enrichment analysis on the genes in the regulatory element peak inferred based on Fiber-seq data between the multiple sample groups.

[0009] In some embodiments, determining the differences in regulatory element peak values ​​inferred based on Fiber-seq data between multiple sample groups based on the cumulative results of each of the multiple sample groups includes: accumulating the supporting read lengths and unsupported read lengths of the regulatory element peak values ​​inferred based on Fiber-seq data of the samples within each sample group, so as to construct a quadruple table based on the cumulative results of the two sample groups; based on the constructed quadruple table, using a two-sided Fisher's exact test to determine the differences in regulatory element peak values ​​inferred based on Fiber-seq data between the two sample groups; and performing a multiple test via FDR on the differences in regulatory element peak values ​​inferred based on Fiber-seq data between the determined two sample groups, so as to select regulatory element peak values ​​inferred based on Fiber-seq data with an FDR less than or equal to 0.05 to determine as the differences in regulatory element peak values ​​inferred based on Fiber-seq data between the multiple sample groups.

[0010] In some embodiments, the method for chromatin openness data mining based on long-read single-molecule sequencing data also includes: among the regulatory element peaks inferred based on Fiber-seq data of each sample, based on the results after typing, selecting the regulatory element peaks inferred based on Fiber-seq data with a coverage depth greater than or equal to a predetermined coverage depth threshold, wherein the results after typing include at least haplotypes; counting the number of reads that support the regulatory element peaks inferred based on Fiber-seq data and the number of reads that do not support the regulatory element peaks inferred based on Fiber-seq data in each haplotype of the two haplotypes for constructing a quadruple table; based on the constructed quadruple table, using a two-sided Fisher's exact test to determine the difference in the regulatory element peaks inferred based on Fiber-seq data between the two haplotypes in a single sample; and among the differences in the regulatory element peaks inferred based on Fiber-seq data between the two haplotypes in a single sample, selecting the differences in the regulatory element peaks inferred based on Fiber-seq data between the two haplotypes that meet predetermined conditions through multiple tests. In some embodiments, identifying the position information of nucleosomes and methyltransferase-sensitive regions based on the position information includes: extracting the position tag information of the methylation modification of the nitrogen atom at position 6 of adenine from the offline data of Fiber-seq; filtering the position tag information about the GC-rich region in the sequence from the extracted position tag information of the methylation modification of the nitrogen atom at position 6 of adenine; identifying the nucleosome region and the non-nucleosome region based on the filtered position tag information through a trained prediction model; and obtaining an unmethylated region having a length greater than a predetermined length threshold, so as to use the unmethylated region to optimize the position information of the identified nucleosome region and the methyltransferase-sensitive region.

[0011] In some embodiments, determining the peak value of the regulatory element inferred based on the Fiber-seq data for each sample includes: assigning an estimated accuracy value to each methyltransferase-sensitive region; calculating the integrated signal value of the regulatory element inferred based on the Fiber-seq data according to the estimated accuracy value assigned to each methyltransferase-sensitive region; and performing a genome-wide correction on the calculated integrated signal value of the regulatory element inferred based on the Fiber-seq data to determine the peak value of the regulatory element inferred based on the Fiber-seq data for each sample.

[0012] In some embodiments, the method for chromatin openness data mining based on long-read single-molecule sequencing data also includes: obtaining information on regulatory elements in a predetermined database, so as to perform database annotation on the regulatory element peak values ​​inferred based on Fiber-seq data for each determined sample according to the information on the regulatory elements; based on the position information of multiple gene elements in the reference genome, determining the overlapping relationship between the positions of the gene elements and the positions of the regulatory element peak values ​​inferred based on Fiber-seq data for each sample, thereby performing position annotation on the regulatory element peak values ​​inferred based on Fiber-seq data; and extracting genes on the regulatory element peak values ​​inferred based on Fiber-seq data for functional and pathway enrichment analysis, so as to perform annotation on the regulatory element peak values ​​inferred based on Fiber-seq data.

[0013] According to a second aspect of the present invention, a computing device is further provided, comprising: at least one processing unit; and at least one memory coupled to the at least one processing unit and storing instructions for execution by the at least one processing unit, wherein the instructions, when executed by the at least one processing unit, cause the device to perform the method of the first aspect of the present invention.

[0014] According to a third aspect of the present invention, a non-transitory computer-readable storage medium is provided, wherein machine-executable instructions are stored on the non-transitory computer-readable storage medium, and when the machine-executable instructions are executed, the machine executes the method according to the first aspect of the present invention.

[0015] According to a fourth aspect of the present invention, a computer program product is further provided, wherein the computer program product comprises instructions, and when the instructions are executed by a machine, the method according to the first aspect of the present invention is implemented.

[0016] This summary is provided to introduce a selection of concepts in a simplified form that are further described below in the detailed description. It is not intended to identify key features or essential features of the invention, nor is it intended to limit the scope of the invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 A schematic diagram of a system for implementing a method for chromatin openness data mining based on long-read single-molecule sequencing data according to an embodiment of the present invention is shown.

[0018] Figure 2 A flowchart of a method for chromatin openness data mining based on long-read single-molecule sequencing data according to an embodiment of the present invention is shown.

[0019] Figure 3 A flow chart of a method for determining differences in FIRE peaks among multiple sample groups according to an embodiment of the present invention is shown.

[0020] Figure 4 A flowchart of a method for annotating FIRE peaks according to an embodiment of the present invention is shown.

[0021] Figure 5 A schematic diagram showing the location and proportion of FIRE peaks in different genetic elements.

[0022] Figure 6 A flow chart showing a method for determining differences in FIRE peaks among haplotypes according to an embodiment of the present invention is shown.

[0023] Figure 7 A flowchart of a method for identifying position information of nucleosomes and MSPs according to an embodiment of the present invention is shown.

[0024] Figure 8 A schematic diagram showing detected MSP information according to an embodiment of the present invention is shown.

[0025] Figure 9 The block diagram schematically shows an electronic device suitable for implementing the embodiments of the present invention.

[0026] Figure 10 A schematic diagram of a computing device 1000 for implementing a method for chromatin openness data mining based on long-read single-molecule sequencing data according to other embodiments of the present invention is shown.

[0027] In the various drawings, the same or corresponding reference numerals denote the same or corresponding parts. DETAILED DESCRIPTION

[0028] The preferred embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although preferred embodiments of the present invention are shown in the accompanying drawings, it should be understood that the present invention can be implemented in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided to make the present invention more thorough and complete and to fully convey the scope of the present invention to those skilled in the art.

[0029] As used herein, the term "including" and its variations represent open inclusion, i.e., "including but not limited to." Unless otherwise stated, the term "or" means "and / or." The term "based on" means "based at least in part on." The terms "one example embodiment" and "an embodiment" mean "at least one example embodiment." The term "another embodiment" means "at least one additional embodiment." The terms "first," "second," etc. may refer to different or identical objects.

[0030] As described above, the shortcomings of traditional chromatin accessibility data mining methods based on long-read single-molecule sequencing data are that they can only complete specific chromatin accessibility analysis, have relatively simple functions, and can only mine chromatin accessibility information of a single sample.

[0031] To at least partially address one or more of the aforementioned and other potential issues, exemplary embodiments of the present invention provide a data mining solution for long-read single-molecule sequencing data. This solution involves calculating an integrated FIRE signal value based on an estimated precision value assigned to each MSP to determine FIRE peaks for each sample; merging FIRE peaks between replicate samples within each sample group based on the physical location information of the determined FIRE peaks to obtain consistent FIRE peaks within each sample group, each sample group comprising multiple samples; and merging FIRE peaks within each sample group based on the physical location information of the FIRE peaks to determine group-specific FIRE peaks across multiple sample groups. This solution not only enables FIRE detection of individual samples but also determines consistent and group-specific FIRE peaks within each sample group. Consequently, the present invention enables mining of chromatin accessibility information within and across sample groups.

[0032] In this plan, the relevant terms are defined as follows:

[0033] Definition of terms

[0034] DNA methylation is an epigenetic modification in which methyl groups are covalently added to DNA molecules. This process can occur at positions such as the N-6 position of adenine (6-mA), the N-4 position of cytosine (4-mC), the N-7 position of guanine (7-mG), and the C-5 position of cytosine (5-mC).

[0035] 6mA

[0036] 6mA refers to the methylation of the nitrogen atom at the 6th position of adenine.

[0037] open chromatin

[0038] The highly folded chromatin structure exposes DNA sequences during replication and transcription. This exposed region is known as the open chromatin region, which allows for the binding of transcription factors and other regulatory elements, and is therefore closely related to transcriptional regulation. Disruption of this dense nucleosome structure allows cis-regulatory elements and trans-acting factors, such as promoters, enhancers, insulators, and silencers, to approach. This property is called chromatin accessibility, and this region is called open chromatin.

[0039] nucleosomes

[0040] The nucleosome is the fundamental structural unit of eukaryotic chromatin, composed of histones and approximately 200 base pairs of DNA. It is a spherical structure approximately 10 nm in diameter. Each nucleosome comprises a histone octamer and one molecule of histone H1. The assembly of DNA into nucleosomal structures provides the structural basis for key properties such as genetic material compaction, stability protection, the introduction of negative supercoiling, and selective gene expression.

[0041] Figure 1 FIG. 1 is a schematic diagram of a system 100 for implementing a method for chromatin openness data mining based on long-read single-molecule sequencing data according to an embodiment of the present invention. Figure 1 As shown, the system 100 includes a computing device 110 and a sequencing device 130. In some embodiments, the computing device 110 and the sequencing device 130 exchange data via a network (not shown).

[0042] The sequencing device 130 is used, for example, to prepare sample libraries and sequence target libraries, generating long-length single-molecule sequencing data (e.g., Fiber-seq data). For example, the sequencing device 130 releases cell nuclei from the sample, treats them with 6mA-MTase, extracts DNA, and then performs a DNA quality test. Once the DNA sample passes the test, Fiber-seq libraries are constructed. After library construction is complete, the library quality is tested. Once the library passes the test, sequencing is performed on the machine to generate Fiber-seq data.

[0043] The computing device 110 can be used to obtain consistent FIRE peaks within each sample group and inter-group-specific FIRE peaks based on the offline data from long-read single-molecule sequencing. In some embodiments, the computing device 110 is configured to identify nucleosome and MSP positional information and determine FIRE peaks for each sample. The computing device 110 is further configured to obtain consistent FIRE peaks within each sample group based on the physical positional information of the determined FIRE peaks, each sample group comprising multiple samples, and to determine inter-group-specific FIRE peaks based on the physical positional information of the FIRE peaks. In some embodiments, the computing device 110 includes one or more processing units, including specialized processing units such as GPUs, FPGAs, and ASICs, as well as general-purpose processing units such as CPUs. One or more virtual machines can also run on each computing device. The computing device 110 includes, for example, a nucleosome and MSP positional information identification unit 112, a FIRE peak determination unit 114 for each sample, a unit 116 for obtaining consistent FIRE peaks within each sample group, and a unit 118 for determining inter-group-specific FIRE peaks. The nucleosome and MSP position information identification unit 112 , the FIRE peaks determination unit 114 for each sample, the intra-group consistent FIRE peaks acquisition unit 116 , and the inter-group specific FIRE peaks determination unit 118 can be configured on one or more computing devices 110 .

[0044] The nucleosome and MSP position information identification unit 112 is configured to extract 6mA position information from the downstream data of long-read single-molecule sequencing, so as to identify the position information of nucleosomes and MSP based on the position information.

[0045] The FIRE peaks determination unit 114 for each sample is configured to calculate an integrated signal value of FIRE according to the estimated accuracy value allocated to each MSP, so as to determine the FIRE peak of each sample.

[0046] The intra-group consistent FIRE peaks acquisition unit 116 is configured to merge the FIRE peaks between repeated samples in each sample group based on the determined physical location information of the FIRE peaks, so as to obtain intra-group consistent FIRE peaks for each sample group, where each sample group includes multiple samples.

[0047] The inter-group specific FIRE peaks determining unit 118 is configured to merge the FIRE peaks of the samples within each sample group based on the physical location information of the FIRE peaks, so as to determine the inter-group specific FIRE peaks of the multiple sample groups.

[0048] The following will be combined Figure 2 A method for chromatin openness data mining based on long-read single-molecule sequencing data according to an embodiment of the present invention is described. Figure 2 FIG2 shows a flow chart of a method 200 for chromatin openness data mining based on long-read single-molecule sequencing data according to an embodiment of the present invention. It should be understood that the method 200 can be used, for example, in Figure 9 The electronic device 900 is described. Figure 1 The method 200 is executed at the described computing device 110. It should be understood that the method 200 may also include additional actions not shown and / or may omit actions shown, and the scope of the present invention is not limited in this respect.

[0049] At step 202 , the computing device 110 extracts 6mA position information from the downstream data of long-read single-molecule sequencing, so as to identify the position information of nucleosomes and MSPs based on the position information.

[0050] Regarding the positional tag information of 6mA, for example, it is the site distribution characteristics of 6mA. It should be understood that by analyzing the site distribution characteristics of 6mA, information on all open chromatin regions across the entire genome in that time and space is obtained. The degree of chromatin openness is closely related to gene transcriptional regulation. Once chromatin is open, it allows certain regulatory proteins (such as transcription factors) to bind to it.

[0051] The nucleosome is the basic structural unit of eukaryotic chromatin. It is composed of histones and approximately 200 base pairs of DNA, forming a spherical body with a diameter of approximately 10 nm. The DNA on the nucleosome is located in the non-open region of chromatin.

[0052] Methyltransferase-sensitive patches (MSPs) are regions of Fiber-seq reads with a high density of 6mA sites and are typically unobstructed by nucleosomes. Nucleosomes are connected by linker DNA, and these fiber sequence segments with a high density of 6mA sites on the linker DNA are defined as MSPs.

[0053] In some embodiments, the computing device 110 uses Fiberseq-qc software to infer 6mA sites, nucleosome positions, and MSP positions in Fiber-seq reads in the offline Fiber-seq data, and counts their number and distribution, thereby evaluating the quality of the Fiber-seq reads.

[0054] Regarding the method for identifying the position information of nucleosomes and MSPs, in some embodiments, for example, it includes: the computing device 110 extracts the 6mA position tag information from the Fiber-seq offline data; in the extracted 6mA position tag information, the position tag information about the GC-rich region in the sequence is filtered; based on the filtered position tag information, the nucleosome region and the non-nucleosome region are identified through the trained prediction model; and unmethylated regions with a length greater than a predetermined length threshold are obtained, so as to use the unmethylated regions to optimize the position information of the identified nucleosome region and MSP. The following will be combined with Figure 7 The method 700 for identifying the position information of nucleosomes and MSPs is described in detail and will not be repeated here.

[0055] At step 204 , the computing device 110 calculates the integrated signal value of FIRE according to the estimated accuracy value assigned to each MSP, so as to determine the FIRE peaks of each sample, thereby accurately determining the potential functional elements.

[0056] Regarding FIRE (Fiber-seq Inferred Regulatory Element), it is a regulatory element with regulatory characteristics inferred through bioinformatics methods based on the characteristics of MSP. These regulatory elements play a key role in gene regulatory networks. They can bind to transcription factors or other regulatory proteins to influence the transcriptional activity of target genes.

[0057] Regarding the method for determining the FIRE peaks of each sample, it includes, for example: the computing device 110 assigns an estimated precision value to each MSP; based on the estimated precision value assigned to each MSP, the integrated signal value of FIRE is calculated; and the calculated integrated signal value of FIRE is corrected on a genome-wide scale to determine the FIRE peaks of each sample.

[0058] Specifically, for example, first, the computing device 110 assigns an estimated accuracy value to each MSP.

[0059] The following describes the algorithm for calculating the estimated accuracy value in conjunction with formula (1).

[0060]

[0061] In the above formula (1), represents the number of “true positive” identifications with XGBoost scores no lower than the current element score in mixed positive labels. represents the number of “false positive” identifications with XGBoost scores no lower than the current element score in negative labels. represents the estimated accuracy value of the FIRE component.

[0062] Next, the computing device 110 calculates the integrated signal value of FIRE according to the estimated accuracy value assigned to each MSP.

[0063] Regarding the method for calculating the integrated signal value of FIRE, for example, the method includes: the computing device 110 calculates the integrated signal value of FIRE based on the number of FIRE elements at position g, the number of Fiber-seq reads, and the estimated precision value of the FIRE element. The algorithm for calculating the integrated signal value of the estimated precision value of FIRE is described below in conjunction with formula (2).

[0064]

[0065] In the above formula (2), represents the integrated signal value of FIRE, represents the number of FIRE elements at position g, represents the number of Fiber-seq reads at position g, and represents the estimated accuracy value of the FIRE element at position g.

[0066] Furthermore, the computing device 110 performs genome-wide correction on the calculated integrated signal value of FIRE to determine the FIRE peaks detected in each sample.

[0067] For example, the computing device 110 performs genome-wide multiple hypothesis testing correction (Bonferroni correction) on the calculated integrated signal value of FIRE to obtain a corrected significance threshold, so that the corrected significance threshold determines the FIRE peaks signal.

[0068] In some embodiments, the corrected significance threshold α=0.01.

[0069] In some embodiments, regions with the highest local FIRE scores and FDR values ​​less than 0.05 are searched in all genomic windows to identify the regions as FIRE peaks. It should be understood that FIRE peaks are generally considered to be potential functional elements because they show a higher enrichment of functional features than surrounding regions.

[0070] Table 1 below schematically shows information on the determined FIRE peaks of the samples.

[0071]

[0072] At step 206 , the computing device 110 merges the FIRE peaks between repeated samples within each sample group based on the determined physical location information of the FIRE peaks to obtain consistent FIREpeaks within each sample group, where each sample group includes multiple samples.

[0073] Regarding the method for obtaining FIRE peaks with intra-group consistency for each sample group, it includes, for example: the computing device 110 merges the FIRE peaks between repeated samples in each sample group based on the physical location information of the FIRE peaks; retains the FIRE peaks between samples whose overlapping length ratio is greater than a predetermined proportion as FIRE peaks with intra-group consistency; in response to determining that the number of repeated samples in the group is greater than or equal to three, determines, among the retained FIRE peaks with intra-group consistency, the FIRE peaks that appear in at least two samples as FIRE peaks with intra-group consistency; and annotates the FIRE peaks determined to be intra-group consistent.

[0074] Regarding the predetermined ratio, it is, for example but not limited to, 50%.

[0075] Regarding the method for annotating FIRE peaks determined to be intra-group consistent, it includes, for example: performing database annotation and position annotation on the FIRE peaks with intra-group consistency, and performing enrichment analysis on genes within the FIRE peaks with intra-group consistency.

[0076] At step 208 , the computing device 110 merges the FIRE peaks of the samples within each sample group based on the physical location information of the FIRE peaks, so as to determine the FIRE peaks specific to the sample groups.

[0077] Regarding the method for determining group-specific FIRE peaks among multiple sample groups, in some embodiments, it includes, for example: the computing device 110 merges the FIRE peaks within each sample group based on the physical location information of the FIRE peaks; determines whether the proportion of the two sample groups in the intersection area of ​​the current FIRE peaks is greater than a predetermined intersection proportion threshold; in response to determining that the proportion of the two sample groups in the intersection area of ​​the current FIRE peaks is greater than the predetermined intersection proportion threshold, determines that the current FIRE peaks are group-specific FIRE peaks among the multiple sample groups; and performs functional and pathway enrichment analysis on the genes in the FIRE peaks determined to be group-specific among the multiple sample groups.

[0078] Regarding the predetermined intersection ratio threshold, it is, for example but not limited to, 50%.

[0079] In this scheme, the integrated FIRE signal value is calculated based on the estimated precision value assigned to each MSP to determine the FIRE peaks for each sample; based on the physical location information of the determined FIRE peaks, the FIRE peaks between repeated samples within each sample group are merged to obtain consistent FIRE peaks within each sample group, each sample group including multiple samples; and based on the physical location information of the FIRE peaks, the FIRE peaks of the samples within each sample group are merged to determine group-specific FIRE peaks across multiple sample groups. This method not only enables FIRE detection of a single sample, but also determines consistent and group-specific FIRE peaks within the group. Thus, the method enables mining of chromatin accessibility information within and between sample groups.

[0080] The following will be combined Figure 3 A method 300 for determining differences in FIREpeaks between a plurality of sample groups according to an embodiment of the present invention is described. Figure 3 FIG. 3 is a flow chart showing a method 300 for determining differences in FIRE peaks between multiple sample groups according to an embodiment of the present invention. It should be understood that the method 300 may be used, for example, in Figure 9 The electronic device 900 is described. Figure 1The method 300 is executed at the described computing device 110. It should be understood that the method 300 may also include additional actions not shown and / or may omit actions shown, and the scope of the present invention is not limited in this respect.

[0081] At step 302 , the computing device 110 accumulates the supported read lengths and unsupported read lengths of the FIRE peaks of the samples within each sample group, so as to construct a quadruple table based on the accumulated results of the two sample groups.

[0082] At step 304 , the computing device 110 uses a two-sided Fisher's exact test based on the constructed quadruple table to determine the difference in FIRE peaks between the two sample groups.

[0083] At step 306 , the computing device 110 performs a multiple test based on the FDR on the differences in the FIRE peaks of the two sample groups, so as to select FIRE peaks with an FDR less than or equal to 0.05 to determine as the differences in the FIRE peaks between the multiple sample groups.

[0084] In some embodiments, the computing device 110 further performs functional annotation, pathway enrichment analysis, and motif analysis on the genes of the determined differential FIRE peaks between multiple sample groups.

[0085] In the above solution, the present invention can accurately determine the difference in FIRE peaks between two sample groups.

[0086] The following will be combined Figure 4 A method 400 for annotating FIRE peaks according to an embodiment of the present invention is described. Figure 4 FIG. 4 is a flow chart of a method 400 for annotating FIRE peaks according to an embodiment of the present invention. It should be understood that the method 400 may be used, for example, in Figure 9 The electronic device 900 is described. Figure 1 The method 400 is executed at the described computing device 110. It should be understood that the method 400 may also include additional actions not shown and / or may omit actions shown, and the scope of the present invention is not limited in this respect.

[0087] At step 402 , the computing device 110 obtains information on regulatory elements in a predetermined database, so as to perform database annotation on the determined FIRE peaks of each sample according to the information on the regulatory elements.

[0088] Regarding the predetermined database, it is, for example, a known regulatory function database. In some embodiments, the predetermined database is, for example, the ENCODE database. It should be understood that the ENCODE database is a comprehensive database of genomic functional elements, which contains epigenetic regulatory element information such as transcription factor binding sites, enhancers and promoters identified through ChIP-seq experiments, chromatin accessibility, histone modifications, and DNA methylation sites identified through techniques such as whole-genome bisulfite sequencing (WGBS). It should be understood that the predetermined database may also be one or more of the CTCF database, the H3K4me3 database, and the Enhancer database.

[0089] Regarding the method for performing database annotation, for example, the method includes: the computing device 110 determines whether there is a regional overlap between the position of the regulatory element in the predetermined database and the position of the FIRE peaks detected by each sample; if it is confirmed that there is a regional overlap between the position of the regulatory element in the predetermined database and the position of the FIRE peaks detected by each sample, then determining that the FIRE peaks detected by each sample are annotated to the database.

[0090] The following expression (3) schematically illustrates a situation where it is determined that there is a region overlap between the position of the regulatory element in the predetermined database and the position of the FIRE peaks detected in each sample.

[0091]

[0092] In the above expression (3), A represents the position of the FIRE peak detected in each sample. B represents the position of the regulatory element in the predetermined database. start represents the starting position. End represents the ending position. That is, A.start represents the starting position of the FIRE peak detected in each sample. A.end represents the ending position of the FIRE peak detected in each sample. B.start represents the starting position of the regulatory element in the predetermined database. B.end represents the ending position of the regulatory element in the predetermined database.

[0093] Table 2 below schematically shows the overlap ratio between the FIRE peaks detected for each sample and the positions of known regulatory elements in the ENCODE database.

[0094]

[0095] At step 404 , the computing device 110 determines the overlapping relationship between the positions of the plurality of genetic elements and the positions of the FIRE peaks of each sample based on the position information of the genetic elements in the reference genome, thereby performing position annotation for the FIRE peaks.

[0096] Figure 5 A schematic diagram showing the location and proportion of FIRE peaks in different genetic elements. Figure 5 This figure shows the distribution of regulatory elements (FIRE peaks) inferred from Fiber-seq data in Sample 1 across different genomic regions. This figure categorizes FIRE peaks and calculates the proportion of each category, reflecting the correspondence between chromatin accessibility and gene functional regions. The "Type" field in the figure includes exonic regions (exonic), intronic regions (intronic), intergenic regions (intergenic), upstream regions (upstream), downstream regions (downstream), UTR3 regions (UTR3), UTR5 regions (UTR5), and other unclassified regions (others). As can be seen, FIRE peaks are primarily distributed in UTR3 (39%), UTR5 (36%), and exonic regions (10%), indicating significant enrichment of chromatin accessibility signals in these regions. Furthermore, the proportion of FIRE peaks located in downstream regions (6%), intergenic regions (4%), other regions (3%), and upstream regions (1%) is relatively low, indicating relatively weak chromatin activity.

[0097] The genetic elements include, for example, multiple elements of the following: exon regions, intron regions, gene intervals, upstream intervals of genes, downstream intervals of genes, UTR3 regions of genes, UTR5 regions of genes, and the like.

[0098] At step 406 , the computing device 110 extracts genes on FIRE peaks for function and pathway enrichment analysis, so as to annotate the FIRE peaks.

[0099] Pathway enrichment analysis methods, for example, include: using topGO software to perform GO enrichment analysis on genes with FIRE peaks. Specifically, for example, the number of genes significantly enriched in each GO term is counted to perform secondary classification statistics for the number of genes in each GO term; and the top 20 terms of each category are selected to indicate the distribution of candidate genes in the GO annotation.

[0100] It should be understood that FIRE peaks are regions where the detected FIRE score is significantly higher than the background level, and these regions may be component regions with regulatory functions.

[0101] In the above scheme, by comparing the detected FIRE peaks with a database of known regulatory functions, the present invention facilitates the identification of regulatory elements. Furthermore, by annotating the genomic locations of FIRE peaks and performing enrichment analysis on the genes they affect, the present invention can understand the genomic origins of FIRE peaks and the pathways and functions they may affect.

[0102] The following will be combined Figure 6 A method 600 for determining differences in FIRE peaks among haplotypes according to an embodiment of the present invention is described. Figure 6 FIG. 6 is a flow chart showing a method 600 for determining differences in FIRE peaks in haplotypes according to an embodiment of the present invention. It should be understood that the method 600 may be used, for example, in Figure 9 The electronic device 900 is described. Figure 1 The method 600 is executed at the described computing device 110. It should be understood that the method 600 may also include additional actions not shown and / or may omit actions shown, and the scope of the present invention is not limited in this respect.

[0103] At step 602 , the computing device 110 selects FIRE peaks having a coverage depth greater than or equal to a predetermined coverage depth threshold from among the FIRE peaks of each sample based on the typing results, wherein the typing results include at least haplotypes.

[0104] Regarding the predetermined coverage depth threshold, it is, for example but not limited to, 10x.

[0105] At step 604 , the computing device 110 counts the number of reads supporting the FIRE peak and the number of reads not supporting the FIRE peak in each of the two haplotypes for constructing a quadruple table.

[0106] Table 3 below schematically shows the constructed quadruple table for the differences in FIRE peaks between two haplotypes.

[0107]

[0108] At step 606 , the computing device 110 uses a two-sided Fisher's exact test based on the constructed quadruple table to determine the difference in FIRE peaks between the two haplotypes in a single sample.

[0109] At step 608 , the computing device 110 selects the differences in FIRE peaks between the two haplotypes in the single sample that meet a predetermined condition through multiple testing.

[0110] Regarding the predetermined condition, it is, for example, Q value < 0.05.

[0111] Regarding the multiple test, for example, it is the Benjamini-Hochberg multiple test.

[0112] The Q-value calculation method, for example, includes the following steps: assuming that m (m is a natural number) FIRE peaks are tested for differences, resulting in m P-values; reordering the P-values ​​from largest to smallest, calculating P = [P1, P2, …, Pm]; calculating the Qi-value for each sorted Pi-value; and for each calculated Q = [Q1, Q2, …, Qm], if a Qi-value is greater than the previous Qi-1 value, assign it to Qi-1; otherwise, retain the corresponding Q-value. The resulting Q-value is called the corrected FDR.

[0113] The following formula (4) schematically shows a method for calculating the Qi value.

[0114]

[0115] In the above formula (4), Pi represents the value of the i-th element in P = [P1, P2, …, Pm]. m is the number of P values, and r is m, m-1, …, 1.

[0116] In the above scheme, the present invention can accurately determine the differences in FIRE peaks between two haplotypes.

[0117] The following will be combined Figure 7 A method 700 for identifying position information of nucleosomes and MSPs according to an embodiment of the present invention is described. Figure 7 FIG. 7 is a flow chart showing a method 700 for identifying the position information of nucleosomes and MSPs according to an embodiment of the present invention. It should be understood that the method 700 can be used, for example, in Figure 9 The electronic device 900 is described. Figure 1 The method 700 is executed at the described computing device 110. It should be understood that the method 700 may also include additional actions not shown and / or may omit actions shown, and the scope of the present invention is not limited in this respect.

[0118] At step 702 , the computing device 110 extracts 6 mA location tag information from the Fiber-seq offline data.

[0119] At step 704 , the computing device 110 filters the position tag information about the GC-rich region in the sequence from the extracted 6 mA position tag information.

[0120] At step 706 , the computing device 110 identifies nucleosome regions and non-nucleosome regions via the trained prediction model based on the filtered position label information.

[0121] Regarding the trained prediction model, it is, for example but not limited to, a Hidden Markov model (HMM) model.

[0122] For example, the computing device 110 randomly extracts 5,000 fiber seq reads and trains a Hidden Markov model (HMM) using the Baum-Welch algorithm to distinguish nucleosomes from non-nucleosomes.

[0123] It should be understood that MSP refers to a region in the fiber that is not blocked by nucleosomes and has a high density of 6mA modification sites. Table 4 below illustrates the MSP information detected for the samples.

[0124]

[0125] At step 708 , the computing device 110 obtains unmethylated regions having a length greater than a predetermined length threshold, so as to use the unmethylated regions to perform optimization on the identified nucleosome regions and the position information of the MSP.

[0126] Regarding the predetermined length threshold, it is, for example, 85 bp.

[0127] For example, the computing device 110 obtains unmethylated regions greater than 85 bp in length and uses the result to optimize the nucleosome and MSP results detected by the HMM model.

[0128] In the above scheme, the present invention can accurately identify the nucleosome region and MSP.

[0129] Figure 8 A schematic diagram showing detected MSP information according to an embodiment of the present invention is shown. Figure 8 The horizontal axis is the length of the detected MSP. Figure 8 The number of times the corresponding length is detected is represented by .

[0130] Figure 10 FIG2 shows a schematic diagram of a computing device 1000 for implementing a method for chromatin openness data mining based on long-read single-molecule sequencing data according to some other embodiments of the present invention. Figure 10As shown, the computing device 1000 includes: a sequencing data acquisition unit 1002, a data quality control unit 1004, a reference genome alignment unit 1006, a variation detection and phasing unit 1008, a 6mA site, nucleosome and MSP prediction and statistics unit 1010, a FIRE peaks determination unit 1012, a FIRE peaks difference determination unit 1014 between haplotypes, a FIRE peaks annotation unit 1016, an intra-group consistency FIRE peaks detection and statistics unit 1018, an inter-group specific FIRE peaks detection and statistics unit 1020, and an inter-group difference FIRE peaks detection and statistics unit 1022.

[0131] The sequencing data acquisition unit 1002 is used to acquire long-read single-molecule sequencing data of the sample to be tested. The long-read single-molecule sequencing data is, for example, data sequenced by PacBio. The sequences acquired from the sequencing data are, for example, HiFi reads.

[0132] Regarding the data quality control unit 1004, it is used to perform data quality control on the long read length single molecule sequencing data obtained. Regarding the reference genome comparison unit 1006, it is used to compare the sequencing data that has been quality controlled with the reference genome. For example, the obtained HiFi Reads are compared with the reference genome. The comparison step, for example, includes: using the pbmm2 module of the SMRTLink software to compare the HiFi Reads data with the reference genome. In some embodiments, the reference genome comparison unit 1006 is further used to perform statistics on the data comparison rate, sequencing depth and coverage based on the comparison results. For example, if the sequencing depth distribution graph is close to a Poisson distribution around the average sequencing depth, it is determined that the uniformity of the data is better.

[0133] The variation detection and phasing unit 1008 is used to detect and annotate variant types such as SNVs, InDels, CNVs, and SVs based on the alignment results of the sequencing data of the test sample with the reference genome; and to phase the SNV and InDel variant results. In some embodiments, the variation detection and phasing unit 1008 is used to detect SNVs and InDels based on the alignment results of the sequencing data of the test sample with the reference genome using DeepVariant software (DeepVariant software is an analysis software that uses deep neural networks to identify genetic variants from alignment files). Then, based on the SNV and InDel detection results, the SNV and InDel sites are annotated. For example, SNV annotations include basic variant site information annotations, gene and region information annotations, normal human database (frequency) annotations, disease database annotations, conservative (harmful) prediction annotations, and gene function and pathway annotations. The methods for detecting and annotating variant types such as InDels, CNVs, and SVs are not further described here. Regarding phasing (or haplotyping) of SNV and InDel variant results, for example, this includes phasing of SNV and InDel variant types using WhatsApp software. It should be understood that variant information output by variant detection software is generally isolated and cannot identify the parental origin of two or more alleles in the genotype information of a single individual. Therefore, completing variant phasing and haplotype construction is crucial for inferring the correct relationship between multiple alleles. In population genetics studies, phased variants can provide data supporting population structure, migration, and environmental pressures. At the individual level, phased variant data can aid clinical decision-making and play a crucial role in studying compound heterozygous genotypes, gene-specific expression, and disease pathogenesis.

[0134] Regarding the 6mA site, nucleosome and MSP prediction and statistics unit 1010, in some embodiments, it includes Figure 1 The nucleosome and MSP position information identification unit 1010 is used to predict 6mA sites, nucleosome positions, and MSP positions in Fiber-seq reads in the sequencing data, and to count and distribute these 6mA sites, nucleosome positions, and MSP positions. In some embodiments, the counted number and distribution can be used to assess the quality of Fiber-seq reads. In some embodiments, the 6mA site, nucleosome, and MSP prediction and statistics unit 1010 is also used to count correlations between samples.

[0135] Regarding the FIRE peaks determination unit 1012, it is used to determine the FIRE of each sample; determine the FIRE peaks with the highest local FIRE score and an FDR value less than 0.05 in all genomic windows; and perform summary statistics on the FIRE Peak. In some embodiments, the FIRE peaks determination unit 1012 includes, for example, Figure 1 The FIREpeaks determination unit 1012 is used for the FIREpeaks determination of each sample in the FIREpeaks determination unit. The FIREpeaks determination unit 1012 is also used for inter-sample correlation statistics. Regarding the summary statistics of FIREPeak, it includes, for example, statistics on the total number of FIRE peaks, the median FIRE peak length, the median FIRE coverage in the Peak, and the median of the integrated signal value of FIRE (e.g., FIRE score). Regarding the inter-sample correlation statistics, it is used, for example, to calculate the correlation between Fiber-seq samples. It should be understood that calculating the correlation between Fiber-seq samples has important biological and experimental significance, which can help evaluate the consistency and similarity between samples. For example, by counting the density of 6mA modifications in the genomic interval, the Pearson correlation coefficient between two samples can be calculated.

[0136] The unit 1014 for determining FIRE peak differences between haplotypes is used to analyze the differential FIRE peaks between two haplotypes in a single sample. For example, the unit 1014 uses FIRE software to analyze and count the differential FIRE peaks between the two haplotypes. The differential FIRE peaks between haplotype 1 and haplotype 2 in each sample are counted. In some embodiments, the unit 1014 for determining FIRE peak differences between haplotypes further includes: GO enrichment analysis of genes associated with the differential FIRE peaks between haplotypes, KEGG enrichment analysis of genes associated with the differential FIRE peaks between haplotypes, and motif analysis of the differential FIRE peaks between haplotypes.

[0137] The FIRE peaks annotation unit 1016 is used, for example, for analyzing the distribution of FIRE peaks within functional regions, GO enrichment analysis of FIRE peak-related genes, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of FIRE peak-related genes, and motif analysis of FIRE peaks. Taking the distribution analysis of FIRE peaks within functional regions as an example, the FIRE peaks annotation unit 1016 uses Annovar software to perform annotation based on the determined FIRE peak location distribution within the genome, obtaining the genomic location information of the FIRE peaks, thereby determining the locations of genetic elements that may have regulatory functions and performing statistical analysis. For another example, the FIRE peaks annotation unit 1016 is also used to analyze FIRE peaks within the ENCODE database. It should be understood that the ENCODE database is a comprehensive database of genomic functional elements, including transcription factor binding sites, enhancers and promoters, chromatin accessibility, and histone modifications identified through ChIP-seq experiments. In some embodiments, this includes information on epigenetic regulatory elements such as DNA methylation sites identified through techniques such as whole-genome bisulfite sequencing (WGBS). The overlap ratio of FIRE peaks identified in each sample with known regulatory elements in the ENCODE database was counted.

[0138] Regarding the intra-group consistency FIRE peaks detection and statistics unit 1018, in some embodiments, it includes Figure 1 The intra-group consistent FIRE peaks acquisition unit is included in the intra-group consistent FIRE peaks detection and statistics unit 1018 is used for merging consistent FIRE peaks of samples in the group; analyzing the intra-group consistent FIRE peaks in the ENCODE database; GO enrichment analysis of genes related to the intra-group consistent FIRE peaks; and KEGG enrichment analysis of genes related to the intra-group consistent FIRE peaks. Regarding the merging of consistent FIRE peaks of samples in the group, it includes, for example: merging FIRE peaks between repeated samples in the group, and retaining the common peaks with overlap of more than 50% between the two samples for subsequent analysis. If there are three or more repeats in the group, the FIRE peaks that appear in at least two samples must be analyzed subsequently.

[0139] Regarding the inter-group specific FIRE peaks detection and statistics unit 1020, in some embodiments, it includes Figure 1The inter-group specific FIRE peaks determination unit is used for the determination and statistics of inter-group specific FIRE peaks in multiple samples; analysis of inter-group specific FIRE peaks in the ENCODE database; GO enrichment analysis of genes associated with inter-group specific FIRE peaks; and KEGG enrichment analysis of genes associated with inter-group specific FIRE peaks. For example, the inter-group specific FIRE peaks detection and statistics unit 1020 first merges the FIRE peaks of samples within a group, and then compares the unique and shared FIRE peaks between the two groups. If the FIRE peaks of the two groups have a greater than 50% overlap in their regions, they are considered shared FIRE peaks, and the rest are unique FIRE peaks for each group.

[0140] The Intergroup Difference FIRE Peaks Detection and Statistics Unit 1022 is used to determine and count intergroup difference FIREpeaks across multiple samples; analyze intergroup difference FIRE peaks in the ENCODE database; perform GO enrichment analysis on genes associated with intergroup difference FIREpeaks; and perform KEGG enrichment analysis on genes associated with intergroup difference FIRE peaks. For example, when performing intergroup difference analysis on multiple replicate samples, the Intergroup Difference FIRE Peaks Detection and Statistics Unit 1022 first merges the FIRE peak coverage of the samples, then uses the Fisher test to screen for differential FIRE peaks with an FDR <= 0.05, and annotates the genes associated with the differential FIRE peaks.

[0141] Figure 9 The block diagram of the electronic device 900 suitable for implementing the embodiment of the present invention is schematically shown. The electronic device 900 may be a device for implementing Figure 2 、 Figure 3 、 Figures 5 to 7 Methods 200-400, 600 to 700 are shown. Figure 9 As shown, electronic device 900 includes a central processing unit (i.e., CPU 901), which can perform various appropriate actions and processes according to computer program instructions stored in read-only memory (i.e., ROM 902) or computer program instructions loaded from storage unit 908 into random access memory (i.e., RAM 903). Various programs and data required for the operation of electronic device 900 can also be stored in RAM 903. CPU 901, ROM 902, and RAM 903 are connected to each other via bus 904. Input / output interface (i.e., I / O interface 905) is also connected to bus 904.

[0142] Multiple components in electronic device 900 are connected to I / O interface 905, including an input unit 906, an output unit 907, and a storage unit 908. CPU 901 executes the various methods and processes described above, such as methods 200-400, 600-700. For example, in some embodiments, methods 200-400, 600-700 may be implemented as a computer software program stored on a machine-readable medium, such as storage unit 908. In some embodiments, part or all of the computer program may be loaded and / or installed on electronic device 900 via ROM 902 and / or communication unit 909. When the computer program is loaded into RAM 903 and executed by CPU 901, one or more operations of methods 200-400, 600-700 described above may be performed. Alternatively, in other embodiments, CPU 901 may be configured to perform one or more actions of methods 200-400, 600-700 via any other suitable means (e.g., via firmware).

[0143] It should be further noted that the present invention may be a method, an apparatus, a system and / or a computer program product. The computer program product may include a computer-readable storage medium carrying computer-readable program instructions for executing various aspects of the present invention.

[0144] A computer-readable storage medium can be a tangible device that can hold and store instructions used by an instruction execution device. Computer-readable storage media can be, for example, but not limited to, an electrical storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination thereof. More specific examples (a non-exhaustive list) of computer-readable storage media include: a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanical encoding device, such as a punch card or a raised structure within a groove on which instructions are stored, and any suitable combination thereof. As used herein, a computer-readable storage medium is not to be construed as a transient signal per se, such as a radio wave or other freely propagating electromagnetic wave, an electromagnetic wave propagating through a waveguide or other transmission medium (e.g., a light pulse passing through a fiber optic cable), or an electrical signal transmitted through an electrical wire.

[0145] The computer-readable program instructions described herein can be downloaded from a computer-readable storage medium to each computing / processing device, or downloaded to an external computer or external storage device via a network, such as the Internet, a local area network, a wide area network, and / or a wireless network. The network can include copper transmission cables, fiber optic transmission, wireless transmission, routers, firewalls, switches, gateway computers, and / or edge servers. The network adapter card or network interface in each computing / processing device receives the computer-readable program instructions from the network and forwards the computer-readable program instructions to be stored in the computer-readable storage medium in each computing / processing device.

[0146] The computer program instructions for performing the operations of the present invention may be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, state setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Smalltalk, C++, and conventional procedural programming languages ​​such as "C" or similar programming languages. The computer-readable program instructions may be executed entirely on the user's computer, partially on the user's computer, as a stand-alone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In the case of a remote computer, the remote computer may be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or may be connected to an external computer (e.g., via the Internet using an Internet service provider). In some embodiments, the state information of the computer-readable program instructions is used to personalize an electronic circuit, such as a programmable logic circuit, a field programmable gate array (FPGA), or a programmable logic array (PLA), so that the electronic circuit can execute the computer-readable program instructions to implement various aspects of the present invention.

[0147] Various aspects of the present invention are described herein with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each block of the flowcharts and / or block diagrams, and combinations of blocks in the flowcharts and / or block diagrams, can be implemented by computer-readable program instructions.

[0148] These computer-readable program instructions can be provided to a processor in a voice interaction device, a general-purpose computer, a special-purpose computer, or a processing unit of another programmable data processing device, thereby producing a machine such that when these instructions are executed by the processing unit of the computer or other programmable data processing device, a device is generated that implements the functions / actions specified in one or more blocks in the flowchart and / or block diagram. These computer-readable program instructions can also be stored in a computer-readable storage medium, where these instructions cause the computer, programmable data processing device, and / or other device to operate in a specific manner. Thus, the computer-readable medium storing the instructions comprises an article of manufacture that includes instructions for implementing various aspects of the functions / actions specified in one or more blocks in the flowchart and / or block diagram.

[0149] Computer-readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device so that a series of operational steps are performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the instructions executed on the computer, other programmable data processing apparatus, or other device to implement the functions / actions specified in one or more blocks in the flowchart and / or block diagram.

[0150] The flow charts and block diagrams in the accompanying drawings show the possible architecture, functions and operations of the devices, methods and computer program products according to multiple embodiments of the present invention. In this regard, each box in the flow chart or block diagram can represent a part of a module, program segment or instruction, and the part of this module, program segment or instruction contains one or more executable instructions for realizing the logical function of the specification. In some alternative implementations, the functions marked in the box can also occur in a sequence different from that marked in the accompanying drawings. For example, two consecutive boxes can actually be executed substantially in parallel, and they can sometimes be executed in the opposite order, depending on the functions involved. It should also be noted that each box in the block diagram and / or flow chart, and the combination of the boxes in the block diagram and / or flow chart can be implemented with a special hardware-based system that performs the function or action of the specification, or can be implemented with a combination of special hardware and computer instructions.

[0151] While various embodiments of the present invention have been described above, the above descriptions are intended to be illustrative, non-exhaustive, and not limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is selected to best explain the principles of the embodiments, their practical applications, or technological improvements in the marketplace, or to enable others skilled in the art to understand the embodiments disclosed herein.

[0152] The above are merely optional embodiments of the present invention and are not intended to limit the present invention. Those skilled in the art will readily appreciate that the present invention is susceptible to various modifications and variations. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the present invention.

Claims

1. A method for chromatin openness data mining based on long-read single-molecule sequencing data, characterized in that: include: Extracting positional information of methylation modification of the nitrogen atom at position 6 of adenine from the downstream data of long-read single-molecule sequencing, so as to identify the positional information of nucleosomes and methyltransferase-sensitive regions based on the positional information; Based on the estimated precision value assigned to each methyltransferase-sensitive region, the integrated signal value of the regulatory elements inferred based on Fiber-seq data was calculated to determine the peak of regulatory elements inferred based on Fiber-seq data for each sample; Based on the determined physical position information of the regulatory element peak value inferred based on the Fiber-seq data, merging the regulatory element peak values ​​inferred based on the Fiber-seq data between repeated samples within each sample group to obtain consistent regulatory element peak values ​​inferred based on the Fiber-seq data within each sample group, wherein each sample group includes multiple samples; Based on the physical position information of the regulatory element peaks inferred based on Fiber-seq data, the regulatory element peaks inferred based on Fiber-seq data of the samples within each sample group are merged to determine the group-specific regulatory element peaks inferred based on Fiber-seq data among multiple sample groups.

2. The method according to claim 1, characterized in that Also includes: The supported read lengths and unsupported read lengths of the regulatory element peaks inferred based on Fiber-seq data of the samples within each sample group are accumulated, so as to determine the differences in the regulatory element peaks inferred based on Fiber-seq data between multiple sample groups based on the accumulated results of multiple sample groups.

3. The method according to claim 1, characterized in that The regulatory element peaks inferred from Fiber-seq data that were consistent within each sample group included: The physical position information of the regulatory element peaks inferred based on Fiber-seq data was merged among the repeated samples in each sample group; The regulatory element peaks inferred based on Fiber-seq data between samples with a coincidence length ratio greater than a predetermined ratio are retained as the regulatory element peaks inferred based on Fiber-seq data with intra-group consistency; In response to determining that the number of duplicate samples within the group is greater than or equal to three, among the retained regulatory element peaks inferred based on Fiber-seq data that are consistent within the group, determining the regulatory element peaks inferred based on Fiber-seq data that appear in at least two samples as regulatory element peaks inferred based on Fiber-seq data that are consistent within the group; and Regulatory element peaks inferred from Fiber-seq data that were determined to be consistent within the group were annotated.

4. The method according to claim 1, wherein Identify group-specific peaks of regulatory elements inferred from Fiber-seq data across multiple sample groups, including: Based on the physical position information of the regulatory element peaks inferred based on Fiber-seq data, the regulatory element peaks inferred based on Fiber-seq data within each sample group were merged; Determine whether the proportion of the two sample groups in the intersection region of the regulatory element peaks currently inferred based on the fiber-seq data is greater than a predetermined intersection proportion threshold; In response to determining that a proportion of an intersection region of the regulatory element peak currently inferred based on the Fiber-seq data for the two sample groups is greater than a predetermined intersection proportion threshold, determining that the regulatory element peak currently inferred based on the Fiber-seq data is a regulatory element peak inferred based on the Fiber-seq data that is unique between the multiple sample groups; and Functional and pathway enrichment analysis was performed on genes in the peaks of regulatory elements inferred from Fiber-seq data that were identified as group-specific across multiple sample groups.

5. The method according to claim 2, characterized in that The differences in the peak values ​​of regulatory elements inferred from Fiber-seq data between the multiple sample groups are determined based on the cumulative results of the multiple sample groups, including: The supported read lengths and unsupported read lengths of the regulatory element peaks inferred based on Fiber-seq data of samples within each sample group are accumulated to construct a quadruple table based on the accumulated results of the two sample groups; Based on the constructed quadruple table, a two-sided Fisher's exact test was used to determine the differences in the peak values ​​of regulatory elements inferred from Fiber-seq data between the two sample groups; and The differences in the peak values ​​of regulatory elements inferred based on Fiber-seq data between the two determined sample groups were subjected to multiple tests using FDR, so as to select the peak values ​​of regulatory elements inferred based on Fiber-seq data with FDR less than or equal to 0.05 to determine as the differences in the peak values ​​of regulatory elements inferred based on Fiber-seq data between the multiple sample groups.

6. The method according to claim 1, characterized in that Also includes: Among the regulatory element peaks inferred based on Fiber-seq data for each sample, selecting regulatory element peaks inferred based on Fiber-seq data with a coverage depth greater than or equal to a predetermined coverage depth threshold based on the results of typing, wherein the results of typing at least include haplotypes; The number of reads supporting the peak of the regulatory element inferred based on Fiber-seq data and the number of reads not supporting the peak of the regulatory element inferred based on Fiber-seq data in each of the two haplotypes were counted to construct a quadruple table; Based on the constructed quadruple table, a two-sided Fisher's exact test was used to determine the differences in the peak values ​​of regulatory elements inferred based on Fiber-seq data between two haplotypes in a single sample; as well as Among the differences in the peak values ​​of regulatory elements inferred based on Fiber-seq data between two haplotypes in a single sample, the differences in the peak values ​​of regulatory elements inferred based on Fiber-seq data between the two haplotypes that meet predetermined conditions are selected through multiple testing.

7. The method according to claim 1, characterized in that Identifying the position information of nucleosomes and methyltransferase-sensitive regions based on the position information includes: Extract the position tag information of the methylation modification of the nitrogen atom at position 6 of adenine from the Fiber-seq data; Among the extracted position tag information of methylation modification of the nitrogen atom at position 6 of adenine, the position tag information about the GC-rich region in the sequence is filtered; Based on the filtered position label information, identifying nucleosome regions and non-nucleosome regions through a trained prediction model; and Unmethylated regions having a length greater than a predetermined length threshold are obtained so as to optimize the position information of the identified nucleosome regions and methyltransferase-sensitive regions using the unmethylated regions.

8. The method according to claim 1, characterized in that The regulatory element peaks inferred from Fiber-seq data for each sample were determined to include: Each methyltransferase-sensitive region was assigned an estimated accuracy value; Calculate the integrated signal value of regulatory elements inferred from Fiber-seq data based on the estimated precision value assigned to each methyltransferase-sensitive region; and The calculated integrated signal values ​​of regulatory elements inferred from Fiber-seq data were corrected genome-wide to determine the peak number of regulatory elements inferred from Fiber-seq data for each sample.

9. The method according to claim 1, characterized in that Also includes: Acquiring information on regulatory elements in a predetermined database, so as to perform database annotation on the regulatory element peak values ​​inferred based on the Fiber-seq data for each determined sample according to the information on the regulatory elements; Based on the position information of multiple genetic elements in the reference genome, determining the overlapping relationship between the position of the genetic element and the position of the regulatory element peak inferred based on the Fiber-seq data of each sample, thereby annotating the position of the regulatory element peak inferred based on the Fiber-seq data; as well as Genes at the regulatory element peaks inferred based on Fiber-seq data were extracted for functional and pathway enrichment analysis, so as to annotate the regulatory element peaks inferred based on Fiber-seq data.

10. A computing device, characterized in that include: at least one processing unit; At least one memory, the at least one memory being coupled to the at least one processing unit and storing instructions for execution by the at least one processing unit, the instructions, when executed by the at least one processing unit, causing the apparatus to perform the steps of the method according to any one of claims 1 to 9.

11. A computer-readable storage medium, characterized in that A computer program is stored on a computer-readable storage medium, and when the computer program is executed by a machine, the method according to any one of claims 1 to 9 is implemented.

12. A computer program product, characterized in that The computer program product comprises instructions, which, when executed by a machine, implement the method according to any one of claims 1 to 9.

Images