Haploid differential methylation region identification method based on long read length sequencing

Abstract

The invention discloses a haplotype differential methylation region recognition method based on long-read-long sequencing, and relates to the field of bioinformatics and epigenetics, in particular to the haplotype differential methylation region recognition method based on long-read-long sequencing. The invention aims to solve the problems of low accuracy, poor applicability, high resource consumption and high time cost of the existing method for identifying the haplotype differential methylation region of long reading and long sequencing. The method comprises the following steps: acquiring sorted and compared BAM files of each person; performing data cleaning on the BAM file to obtain a file with methylation state CpG site information; extracting the feature of each CpG locus in the file with the methylation state CpG locus information; inputting the extracted features of each CpG locus into the trained long-short-term memory network, and outputting a typing result of each CpG locus by the trained long-short-term memory network model; and identifying a haplotype differential methylation region of the typing result of each CpG site in the file output in the step 3.

Description

Technical Field

[0001] The present invention relates to the fields of bioinformatics and epigenetics, and in particular to a method for identifying haplotype differentially methylated regions based on long-read sequencing. Background Art

[0002] DNA methylation is an important modification information in epigenetics, among which differentially methylated regions (DMRs) are used to characterize the differences in methylation levels between different samples or populations. DMRs can be defined from multiple perspectives, including: (1) population DMRs, which reflect methylation differences between different populations; (2) structural variation DMRs, which reflect changes in methylation levels in structural variation regions; and (3) haplotype differentially methylated regions (hDMRs), which are used to characterize methylation differences between haplotypes at the individual level. hDMRs are of great significance in revealing allele-specific regulation, haplotype genetic effects, and their impact on gene expression, and contribute to a deeper understanding of epigenetic mechanisms at the individual level. However, since haplotype resolution relies on high-quality sequencing data and accurate variant positioning, the identification of hDMRs faces challenges such as low resolution and complex operations.

[0003] Traditional hDMR identification methods mainly rely on whole-genome bisulfite sequencing (WGBS) or high-throughput methylation arrays (such as 450K and 850K arrays). These sequencing methods have the following drawbacks:

[0004] Limitations of sequencing technology: Sulfite treatment can cause DNA sequence damage, and the sequencing read length is short (e.g., the average WGBS read length is approximately 150 bp), making it difficult to accurately identify the methylation typing of long fragments, thus affecting the detection accuracy of hDMR.

[0005] Difficulties in data collection: Traditional hDMR identification methods typically rely on family data (paternal, maternal, and offspring) for haplotype typing in order to analyze the methylation information carried by haplotypes. However, in actual research, most samples lack parental information, which greatly limits the applicability of this method.

[0006] hDMR is the haplotype differentially methylated region;

[0007] Complex analysis process: Traditional methods rely on variant typing information to genotype methylation data, typically requiring variant detection followed by analysis in conjunction with methylation data. However, variant typing is complex, computationally intensive, and time-consuming, making it challenging for non-professional engineers. Summary of the Invention

[0008] The purpose of the present invention is to solve the problems of low accuracy, poor applicability, high resource consumption and high time cost in the existing methods for identifying haplotype differentially methylated regions (hDMRs) by long-read sequencing, and to propose a method for identifying haplotype differentially methylated regions based on long-read sequencing.

[0009] A method for identifying haplotype differentially methylated regions based on long-read sequencing. The specific process is as follows:

[0010] Step 1: Obtain the sorted and aligned BAM files for each person;

[0011] Each person corresponds to 1 BAM file;

[0012] Each BAM file contains M SNP sites;

[0013] Step 2: Clean the sorted and aligned BAM files to obtain files with CpG site information of methylation status;

[0014] Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information;

[0015] Step 3: Input the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 into the trained long short-term memory network LSTM. The trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information;

[0016] Step 4: Use the maximum likelihood regression method to identify the haplotype differential methylation regions of the typing results of each CpG site in the file output in step 3.

[0017] Preferably, in step 1, the sorted and aligned BAM files of each person are obtained;

[0018] Each BAM file contains M SNP sites; 100000≤M≤10000000;

[0019] The specific process is:

[0020] Step 11, obtain the human whole genome Fast5 electrical signal file from the sequencing chip R9 or R10 in the third-generation sequencing platform Oxford Nanopore;

[0021] Step 12: Use Guppy software to convert the Fast5 electrical signal file into a Fastq format file;

[0022] Step 13. Set Nanopolish software parameters: index;

[0023] Use Nanopolish software to create an index connection between the Fast5 electrical signal file in step 11 and the Fastq format file in step 12;

[0024] Step 14. Set the minimap2 software parameters: -a, -x; the value of -x is map-ont;

[0025] The Fastq format file of the human reference genome GRch38 and the index connected in step 13 is input into the minimap2 software. The minimap2 software converts the Fastq format file generated in step 12 into an aligned BAM file;

[0026] Step 15. Set the Samtools software parameters: sort;

[0027] Use Samtools software to sort the BAM file in step 14 to generate a sorted and aligned BAM file;

[0028] Step 16. Set the Samtools software parameters: index;

[0029] Use Samtools software to index the BAM file from step 15.

[0030] Preferably, in step 2, data cleaning is performed on the sorted and aligned BAM files to obtain a file with methylation status CpG site information;

[0031] Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information;

[0032] The specific process is:

[0033] Step 21: Clean the sorted and aligned BAM files to obtain files with methylation status CpG site information;

[0034] The specific process is:

[0035] Step 211: Remove low-quality reads from the sorted and aligned BAM files; the process is as follows:

[0036] Set Samtools software parameters: -q;

[0037] Samtools software was used to remove reads with alignment quality lower than 30 in the sorted and aligned BAM files;

[0038] Step 212: Remove the PCR duplicate fragments in the sorted and aligned BAM files retained in step 211; the process is:

[0039] Set PicardTools software parameters: MarkDuplicates, REMOVE_DUPLICATES, and the value of REMOVE_DUPLICATES is true;

[0040] Use PicardTools software to identify and remove PCR duplicate fragments in the sorted and aligned BAM files retained in step 211;

[0041] Step 213: Set the Nanopolish software parameter: --methylation, where the value of --methylation is CpG;

[0042] Input the sorted and aligned BAM files retained in step 212, the Fastq format files indexed in step 13, and the human reference genome GRch38 into Nanopolish software, and Nanopolish software outputs a file with methylation status CpG site information;

[0043] Step 22: Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information.

[0044] Preferably, in step 22, based on the file with the methylation status CpG site information, the features of each CpG site in the file with the methylation status CpG site information are extracted;

[0045] The specific process is:

[0046] Step 221: Count the number of methylation states of 0 and 1 for each CpG site in the file with methylation status CpG site information obtained in step 21;

[0047] The methylation level β of each CpG site is calculated based on the methylation status corresponding to each CpG site. j ;

[0048] Step 222: Find the number of the first allele G1 and the number of the second allele G2 corresponding to the position information of each CpG in the file with the methylation status CpG site information saved in step 213 in the sorted and aligned BAM file saved in step 212;

[0049] Step 223: defining each CpG in the file with methylation status CpG site information saved in step 213 as a central CpG;

[0050] The number of CpGs in the 2kb region upstream and downstream of the central CpG and the methylation level information of each CpG in the 2kb region upstream and downstream of the central CpG are saved as a 2kb information file corresponding to the CpG site;

[0051] Step 224: Calculate the mean and variance of the methylation level β values of all CpG sites in the 2 kb information file corresponding to the CpG sites saved in step 223;

[0052] Step 225: Calculate the relative position Pos of each CpG site and the central CpG site in the 2kb information file corresponding to the CpG site saved in step 223 i ;

[0053] Step 226: Compare the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223 with the human reference genome GRch38 to determine whether the central CpG site is located in the exon, intron, promoter, or enhancer region of the human reference genome GRch38, and perform annotation.

[0054] Step 227: The methylation level β of each CpG site obtained in step 221 is converted to j , the number of the first allele G1 and the number of the second allele G2 found in step 222, the 2kb information file corresponding to the CpG site saved in step 223, the mean and variance of the methylation level β values of all CpG sites in the 2kb information file corresponding to the CpG site calculated in step 224, and the relative position Pos of each CpG site to the central CpG site in the 2kb information file corresponding to the CpG site calculated in step 225 i The annotation result of step 226 is used as the feature of each CpG site in the file with methylation status CpG site information.

[0055] Preferably, in step 221, the number of methylation states of 0 and 1 for each CpG site in the file with methylation status CpG site information obtained in step 21 is counted;

[0056] The methylation level β of each CpG site is calculated based on the methylation status corresponding to each CpG site. j ; expressed as:

[0057]

[0058] Among them, β j is the methylation level of the j-th CpG site:

[0059] M c is the number of methylation states of 1 corresponding to the j-th CpG site, M Tis the number of methylation states of 0 corresponding to the j-th CpG site;

[0060] The methylation status of a CpG site is 0, indicating unmethylation;

[0061] The methylation status of a CpG site is 1, indicating methylation.

[0062] Preferably, the relative position Pos of each CpG site and the central CpG site in the 2kb information file corresponding to the CpG site saved in step 223 is calculated in step 225. i ; The specific process is:

[0063]

[0064] in,

[0065] CpG center The position of the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223;

[0066] CpG i The position of the i-th CpG site excluding the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223;

[0067] Preferably, in step 3, the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 are input into a trained long short-term memory network LSTM, and the trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information;

[0068] The specific process is:

[0069] Step 31: Get the training set;

[0070] Step 32: Using the features of each CpG site in the file with methylation status CpG site information in the training set as the input of the long short-term memory network model, and using the typing results of each CpG site in the file with methylation status CpG site information in the training set as the output of the long short-term memory network model;

[0071] The cross-entropy loss function is used to optimize the parameters of the long short-term memory network model, and the Adam optimizer is combined to perform gradient updates until the cross-entropy loss function converges to obtain a trained long short-term memory network model.

[0072] Step 33: Input the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 into the trained long short-term memory network LSTM, and the trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information.

[0073] Preferably, the training set is obtained in step 31; the specific process is:

[0074] Step 311: Obtain the sorted and aligned BAM files of each person;

[0075] Each person corresponds to 1 BAM file;

[0076] Each BAM file contains M SNP sites;

[0077] Step 312: performing data cleaning on the sorted and aligned BAM files to obtain a file with methylation status CpG site information;

[0078] Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information;

[0079] Step 313: Using the sequencing tool NanoMethPhase, each CpG site in the file containing the methylation status CpG site information is typed into haplotype 1 and haplotype 2, thereby obtaining a CpG typing file.

[0080] Preferably, in step 4, the maximum likelihood regression method is used to identify the haplotype differential methylation region hDMR of the typing result of each central CpG site in the methylation status CpG site information file output in step 3; the specific process is:

[0081] Step 41: Find candidate hDMR blocks; the process is:

[0082] Step 411: Divide all CpG sites in the CpG typing file output in step 3 into windows, and obtain seed sites in each window; the process is as follows:

[0083] Use a 4 kb window to slide all CpG sites in the CpG typing file output in step 3, and select the CpG site with the largest Δβ value in each window as the seed site;

[0084] The calculation process of the Δβ value of the CpG site is:

[0085] Methylation level β of haplotype 1 based on each CpG site hap1 and haplotype 2 methylation levels β hap2 Calculate Δβ = β for each CpG sitehap1 -β hap2 ;

[0086] Step 412: Set i=1;

[0087] Step 413: Expand the hDMR region toward both sides with the seed site as the center within the i-th window;

[0088] The specific process is:

[0089] 1) Copy the seed site in the i-th window and place it in the empty hDMR region, and delete the seed site from the i-th window. If the absolute value of the Δβ value of the deleted seed site is greater than or equal to 0.1, retain the seed site in the hDMR region and execute 2);

[0090] Copy the seed site in the i-th window and place it in the empty hDMR region, and delete the seed site from the i-th window. If the absolute value of the Δβ value of the deleted seed site is less than 0.1, the seed site is not retained in the hDMR region, and the hDMR region is not expanded. Execute step 414.

[0091] 2) Calculate the Δβ value of the CpG site adjacent to the left of the deleted seed site and the Δβ value of the CpG site adjacent to the right of the deleted seed site in the i-th window;

[0092] If the Δβ value of the CpG site adjacent to the left of the deleted seed site is greater than the Δβ value of the CpG site adjacent to the right of the deleted seed site, the distance between the CpG site adjacent to the left of the seed site and the left edge of the hDMR region is calculated; if the distance is less than 500 bp, the CpG site adjacent to the left of the deleted seed site is placed in the leftmost region of the hDMR, and 3) is executed; if the distance is greater than or equal to 500 bp, the CpG site adjacent to the left is retained in the i-th window and is not placed in the leftmost region of the hDMR, and the left side of the hDMR region is no longer extended; the distance between the CpG site adjacent to the right of the seed site and the right edge of the hDMR region is calculated; if the distance is less than 500 bp, the CpG site adjacent to the right of the deleted seed site is placed in the rightmost region of the hDMR, and 3) is executed; if the distance is greater than or equal to 500 bp, the CpG site adjacent to the right is retained in the i-th window and is not placed in the rightmost region of the hDMR, and the right side of the hDMR region is no longer extended, and 4) is executed;

[0093] If the Δβ value of the CpG site adjacent to the left of the deleted seed site is less than or equal to the Δβ value of the CpG site adjacent to the right of the deleted seed site, the distance between the CpG site adjacent to the right of the seed site and the right boundary of the hDMR region is calculated. If the distance is less than 500 bp, the CpG site adjacent to the right of the deleted seed site is placed in the rightmost region of the hDMR, and step 3) is executed. If the distance is greater than or equal to 500 bp, the CpG site adjacent to the right is retained in the i-th window and is not placed in the rightmost region of the hDMR, and the right side of the hDMR region is no longer extended. The distance between the CpG site adjacent to the left of the seed site and the left boundary of the hDMR region is calculated. If the distance is less than 500 bp, the CpG site adjacent to the left of the deleted seed site is placed in the leftmost region of the hDMR, and step 3) is executed. If the distance is greater than or equal to 500 bp, the CpG site adjacent to the left is retained in the i-th window and is not placed in the leftmost region of the hDMR, and the left side of the hDMR region is no longer extended, and step 4) is executed.

[0094] 3) Calculate the mean of the absolute values of the Δβ values of all CpG sites in the hDMR region;

[0095] If the mean of the absolute values of the Δβ values of all CpG sites in the hDMR region is greater than or equal to 0.1, the newly inserted CpG site is retained in the hDMR region, and the retained CpG site is deleted from the i-th window, and step 2 is executed;

[0096] If the mean of the absolute values of the Δβ values is less than 0.1, the hDMR region does not retain the newly inserted CpG site, and the hDMR region is no longer expanded; step 414 is executed;

[0097] 4) Determine whether the number of CpG sites in the hDMR region is greater than or equal to 5;

[0098] If the number of CpG sites is greater than or equal to 5, retain the hDMR region; and execute step 414;

[0099] If the number of CpG sites is less than 5, delete the hDMR region; and execute step 414;

[0100] Step 414: In the i-th window, the point with the largest absolute value of the new Δβ is used as the seed site, and the hDMR region is expanded in both directions with the seed site as the center. The specific process is as follows:

[0101] Deleting the CpG sites placed in the hDMR region in step 413 within the i-th window;

[0102] The CpG site corresponding to the maximum absolute value of the Δβ value in the i-th window after deletion is selected as the seed site;

[0103] Repeat steps 413 to 414 until the number of CpG sites in the i-th window is less than 5, and then execute step 415;

[0104] Step 415: Let i=i+1, repeat steps 413 to 414, and then execute step 416;

[0105] Step 416: Repeat step 415 until all windows are calculated and all hDMR regions are obtained;

[0106] Each hDMR region is defined as a candidate hDMR region;

[0107] Step 42: Perform Wilcoxon test and multiple statistical test on all candidate hDMR regions obtained in step 41 to obtain final candidate hDMR regions.

[0108] Preferably, in step 42, Wilcoxon test and multiple statistical test are performed on all candidate hDMR regions obtained in step 41 to obtain final candidate hDMR regions; the specific process is:

[0109] Step 421: Perform a Wilcoxon test on all candidate hDMR regions to obtain p-values, and select candidate hDMR regions corresponding to p < 0.05. The specific process is as follows:

[0110] The methylation level of haplotype 1 for each CpG site in all candidate hDMR regions is β hap1 and haplotype 2 methylation levels β hap2 The Wilcoxon test was performed to obtain the p value, and the candidate hDMR regions corresponding to p < 0.05 were selected;

[0111] Step 422: Perform multiple verification and correction on the candidate hDMR region selected in step 421 to obtain the final candidate hDMR region. The specific process is as follows:

[0112] The Benjamini-Hochberg method is used to correct the p-values of the candidate hDMR regions selected in step 421, and candidate hDMR regions with p<0.05 are retained to obtain the final candidate hDMR regions.

[0113] The beneficial effects of the present invention are:

[0114] Improve detection accuracy: Utilize long-read sequencing data to overcome the impact of short-read sequencing on methylation typing accuracy and improve the accuracy of hDMR detection.

[0115] Reduced data requirements: No need to rely on family data, suitable for large-scale population studies, and improve the applicability of the method.

[0116] Simplified analysis process: The method of the present invention avoids the complex variant typing steps and directly analyzes the methylation data, reducing computing resource consumption and improving computing efficiency.

[0117] Application prospects: This invention can be widely used in the fields of genomic epigenetic analysis, imprinted gene research, and disease-related methylation marker screening. It is particularly suitable for DNA methylation research at haplotype resolution, providing a new research method for the epigenetic regulatory mechanism of complex diseases.

[0118] hDMR recognition model:

[0119] The long short-term memory (LSTM) network based on deep learning can directly identify hDMRs without relying on variant typing information, greatly simplifying the data analysis process;

[0120] Efficient computing:

[0121] The method of the present invention directly analyzes methylation status based on long-read sequencing data, thereby improving the recognition accuracy of hDMR regions; without the need for additional family data support, it reduces the limitations of sample data and improves generalization capabilities.

[0122] This paper integrates long-read sequencing technology, deep learning models, and weighted statistical testing to construct an efficient and accurate hDMR identification process. Compared with traditional methods, the overall analysis efficiency is improved by approximately 60%. In large-scale sample testing, the hDMR identification accuracy rate can reach 90%, fully demonstrating the practical value and broad prospects of this invention in epigenetic research and biomarker screening. BRIEF DESCRIPTION OF THE DRAWINGS

[0123] Figure 1 It is a flow chart of the present invention. DETAILED DESCRIPTION

[0124] Specific embodiment 1: This embodiment is a method for identifying haplotype differential methylation regions based on long-read sequencing. The specific process is as follows:

[0125] Step 1: Obtain the sorted and aligned BAM files for each person;

[0126] Each person corresponds to 1 BAM file;

[0127] Each BAM file contains M SNP sites;

[0128] Step 2: Clean the sorted and aligned BAM files to obtain files with CpG site information of methylation status;

[0129] Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information;

[0130] Step 3: Input the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 into the trained long short-term memory network LSTM. The trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information;

[0131] Step 4: Use the maximum likelihood regression method to identify the haplotype differentially methylated regions (hDMRs) of the typing results of each CpG site in the file output in step 3.

[0132] Specific embodiment 2: This embodiment differs from specific embodiment 1 in that: in step 1, the sorted and aligned BAM files of each person are obtained;

[0133] Each BAM file contains M SNP sites; 100000≤M≤10000000;

[0134] The specific process is:

[0135] Step 11, obtain the human whole genome Fast5 electrical signal file from the sequencing chip R9 or R10 in the third-generation sequencing platform Oxford Nanopore;

[0136] Step 12: Use Guppy software to convert the Fast5 electrical signal file into a Fastq format file;

[0137] Step 13. Set Nanopolish software parameters: index;

[0138] Use Nanopolish software to establish an index connection between the Fast5 electrical signal file in step 11 and the Fastq format file in step 12 (establish a connection relationship between the two files);

[0139] Step 14. Set the minimap2 software parameters: -a, -x; the value of -x is map-ont;

[0140] The Fastq format file of the human reference genome GRch38 and the index connected in step 13 is input into the minimap2 software. The minimap2 software converts the Fastq format file generated in step 12 into an aligned BAM file;

[0141] Step 15. Set the Samtools software parameters: sort;

[0142] Use Samtools software to sort the BAM file in step 14 to generate a sorted and aligned BAM file;

[0143] Step 16. Set the Samtools software parameters: index;

[0144] Use Samtools software to index the BAM file from step 15.

[0145] Other steps and parameters are the same as those in the first embodiment.

[0146] Specific embodiment three: This embodiment differs from specific embodiment one or two in that: in step 2, the sorted and aligned BAM files are cleaned to obtain a file with methylation status CpG site information;

[0147] Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information;

[0148] The specific process is:

[0149] Step 21: Clean the sorted and aligned BAM files to obtain files with methylation status CpG site information to ensure the accuracy of downstream analysis;

[0150] The specific process is:

[0151] Step 211: Remove low-quality reads from the sorted and aligned BAM files; the process is as follows:

[0152] Set Samtools software parameters: -q;

[0153] Samtools software was used to remove reads with alignment quality lower than 30 in the sorted and aligned BAM files;

[0154] BAM files contain reads, and reads have their own alignment quality;

[0155] Step 212: Remove the PCR duplicate fragments in the sorted and aligned BAM files retained in step 211; the process is:

[0156] Set Picard Tools software parameters: MarkDuplicates, REMOVE_DUPLICATES, and the value of REMOVE_DUPLICATES is true;

[0157] Picard Tools software was used to identify and remove PCR duplicate fragments in the sorted and aligned BAM files retained in step 211 to avoid duplication affecting the calculation of methylation levels;

[0158] Step 213: Set the Nanopolish software parameter: --methylation, where the value of --methylation is CpG;

[0159] Input the sorted and aligned BAM files retained in step 212, the Fastq format files indexed in step 13, and the human reference genome GRch38 into Nanopolish software, and Nanopolish software outputs a file with methylation status CpG site information;

[0160] Step 22: Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information.

[0161] Other steps and parameters are the same as those in the first or second embodiment.

[0162] Specific embodiment 4: This embodiment differs from any one of specific embodiments 1 to 3 in that: in step 22, based on the file with the methylation status CpG site information, the features of each CpG site in the file with the methylation status CpG site information are extracted;

[0163] The specific process is:

[0164] Step 221: Count the number of methylation states of 0 and 1 for each CpG site in the file with methylation status CpG site information obtained in step 21;

[0165] The methylation level β of each CpG site is calculated based on the methylation status corresponding to each CpG site. j ;

[0166] Step 222: Find the number of the first allele G1 and the number of the second allele G2 corresponding to the position information of each CpG in the file with the methylation status CpG site information saved in step 213 in the sorted and aligned BAM file saved in step 212;

[0167] Step 223: defining each CpG in the file with methylation status CpG site information saved in step 213 as a central CpG;

[0168] The number of CpGs in the 2kb region upstream and downstream of the central CpG and the methylation level information of each CpG in the 2kb region upstream and downstream of the central CpG are saved as a 2kb information file corresponding to the CpG site;

[0169] Step 224: Calculate the mean and variance of the methylation level β values of all CpG sites in the 2 kb information file corresponding to the CpG sites saved in step 223;

[0170] Step 225: Calculate the CpG site information file corresponding to the CpG site saved in step 223. i ) and central CpG sites (CpG center )'s relative position Posi;

[0171] Step 226: Compare the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223 with the human reference genome GRch38 to determine whether the central CpG site is located in the exon, intron, promoter, or enhancer region of the human reference genome GRch38, and perform annotation.

[0172] Step 227: The methylation level β of each CpG site obtained in step 221 is converted to j , the number of the first allele G1 and the number of the second allele G2 found in step 222, the 2kb information file corresponding to the CpG site saved in step 223, the mean and variance of the methylation level β values of all CpG sites in the 2kb information file corresponding to the CpG site calculated in step 224, and the methylation level of each CpG site (CpG i ) and central CpG sites (CpG center )'s relative position Pos i The annotation result of step 226 is used as the feature of each CpG site in the file with methylation status CpG site information.

[0173] The other steps and parameters are the same as those in the first to third embodiments.

[0174] Specific embodiment 5: This embodiment differs from any one of specific embodiments 1 to 4 in that: in step 221 , the number of methylation states of 0 and 1 for each CpG site in the file with methylation state CpG site information obtained in step 21 is counted;

[0175] The methylation level β of each CpG site is calculated based on the methylation status corresponding to each CpG site. j ; expressed as:

[0176]

[0177] Among them, β j is the methylation level of the j-th CpG site:

[0178] Mc is the number of methylation states of 1 corresponding to the j-th CpG site, M T is the number of methylation states of 0 corresponding to the j-th CpG site;

[0179] The methylation status of a CpG site is 0, indicating unmethylation;

[0180] The methylation status of a CpG site is 1, indicating methylation.

[0181] The other steps and parameters are the same as those in the first to fourth embodiments.

[0182] Specific embodiment 6: This embodiment differs from any one of the specific embodiments 1 to 5 in that: in step 225, each CpG site (CpG i ) and central CpG sites (CpG center )'s relative position Pos i ; The specific process is:

[0183]

[0184] in,

[0185] CpG center The position of the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223;

[0186] CpG i The position of the i-th CpG site excluding the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223;

[0187] The other steps and parameters are the same as those in the first to fifth embodiments.

[0188] Specific embodiment seven: This embodiment differs from any one of specific embodiments one to six in that: in step 3, the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 are input into a trained long short-term memory network LSTM, and the trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information;

[0189] The specific process is:

[0190] Step 31: Get the training set;

[0191] Step 32: Using the features of each CpG site in the file with methylation status CpG site information in the training set as the input of the long short-term memory network model, and using the typing results of each CpG site in the file with methylation status CpG site information in the training set as the output of the long short-term memory network model;

[0192] The cross-entropy loss function is used to optimize the parameters of the long short-term memory network model, and the Adam optimizer is combined to perform gradient updates until the cross-entropy loss function converges to obtain a trained long short-term memory network model.

[0193] Step 33: Input the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 into the trained long short-term memory network LSTM, and the trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information.

[0194] Haplotype 1 of each CpG site contains 0, 1 or more numbers with a methylation state of 1 and 0, 1 or more numbers with a methylation state of 0;

[0195] A bidirectional LSTM (BiLSTM) structure is used to simultaneously learn contextual information and improve the model's ability to analyze methylation patterns.

[0196] The first two layers of the bidirectional LSTM contain 128 neurons each, and the final layer contains 64 neurons, ensuring the model's effectiveness in complex data patterns. An attention mechanism is used to further enhance focus on key CpG sites and improve the accuracy of DMC recognition.

[0197] The cross entropy loss function is used to optimize the model parameters, and the Adam optimizer is used for gradient update.

[0198] The other steps and parameters are the same as those in the first to sixth embodiments.

[0199] Specific embodiment eight: This embodiment differs from any one of specific embodiments one to seven in that: the training set is obtained in step 31; the specific process is:

[0200] Step 311: Obtain the sorted and aligned BAM files of each person;

[0201] Each person corresponds to 1 BAM file;

[0202] Each BAM file contains M SNP sites;

[0203] The process is the same as step 1;

[0204] Step 312: performing data cleaning on the sorted and aligned BAM files to obtain a file with methylation status CpG site information;

[0205] Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information;

[0206] The process is the same as step 2;

[0207] Step 313: Using the sequencing tool NanoMethPhase, each CpG site in the file containing the methylation status CpG site information is typed into haplotype 1 and haplotype 2, thereby obtaining a CpG typing file.

[0208] The other steps and parameters are the same as those in the first to seventh embodiments.

[0209] Specific embodiment 9: This embodiment differs from any one of specific embodiments 1 to 8 in that: in step 4, the maximum likelihood regression method is used to identify the haplotype differential methylation region hDMR of the typing result of each central CpG site in the methylation status CpG site information file output in step 3; the specific process is:

[0210] Step 41: Find candidate hDMR blocks; the process is:

[0211] Step 411: Divide all CpG sites in the CpG typing file output in step 3 into windows, and obtain seed sites in each window; the process is as follows:

[0212] Use a 4 kb window to slide all CpG sites in the CpG typing file output in step 3, and select the CpG site with the largest Δβ value in each window as the seed site;

[0213] The calculation process of the Δβ value of the CpG site is:

[0214] Methylation level β of haplotype 1 based on each CpG site hap1 and haplotype 2 methylation levels β hap2 Calculate Δβ = β for each CpG site hap1 -β hap2 ;

[0215] Step 412: let i=1;

[0216] Step 413: Expand the hDMR region toward both sides with the seed site as the center within the i-th window;

[0217] The specific process is:

[0218] 1) Copy the seed site in the i-th window and place it in the empty hDMR region, and delete the seed site from the i-th window. If the absolute value of the Δβ value of the deleted seed site is greater than or equal to 0.1, retain the seed site in the hDMR region and execute 2);

[0219] Copy the seed site in the i-th window and place it in the empty hDMR region, and delete the seed site from the i-th window. If the absolute value of the Δβ value of the deleted seed site is less than 0.1, the seed site is not retained in the hDMR region, and the hDMR region is not expanded. Execute step 414.

[0220] 2) Calculate the Δβ value of the CpG site adjacent to the left of the deleted seed site and the Δβ value of the CpG site adjacent to the right of the deleted seed site in the i-th window;

[0221] If the Δβ value of the CpG site adjacent to the left of the deleted seed site is greater than the Δβ value of the CpG site adjacent to the right of the deleted seed site, the distance between the CpG site adjacent to the left of the seed site and the left edge of the hDMR region is calculated; if the distance is less than 500 bp, the CpG site adjacent to the left of the deleted seed site is placed in the leftmost region of the hDMR, and 3) is executed; if the distance is greater than or equal to 500 bp, the CpG site adjacent to the left is retained in the i-th window and is not placed in the leftmost region of the hDMR, and the left side of the hDMR region is no longer extended; the distance between the CpG site adjacent to the right of the seed site and the right edge of the hDMR region is calculated; if the distance is less than 500 bp, the CpG site adjacent to the right of the deleted seed site is placed in the rightmost region of the hDMR, and 3) is executed; if the distance is greater than or equal to 500 bp, the CpG site adjacent to the right is retained in the i-th window and is not placed in the rightmost region of the hDMR, and the right side of the hDMR region is no longer extended, and 4) is executed;

[0222] If the Δβ value of the CpG site adjacent to the left of the deleted seed site is less than or equal to the Δβ value of the CpG site adjacent to the right of the deleted seed site, the distance between the CpG site adjacent to the right of the seed site and the right boundary of the hDMR region is calculated. If the distance is less than 500 bp, the CpG site adjacent to the right of the deleted seed site is placed in the rightmost region of the hDMR, and step 3) is executed. If the distance is greater than or equal to 500 bp, the CpG site adjacent to the right is retained in the i-th window and is not placed in the rightmost region of the hDMR, and the right side of the hDMR region is no longer extended. The distance between the CpG site adjacent to the left of the seed site and the left boundary of the hDMR region is calculated. If the distance is less than 500 bp, the CpG site adjacent to the left of the deleted seed site is placed in the leftmost region of the hDMR, and step 3) is executed. If the distance is greater than or equal to 500 bp, the CpG site adjacent to the left is retained in the i-th window and is not placed in the leftmost region of the hDMR, and the left side of the hDMR region is no longer extended, and step 4) is executed.

[0223] 3) Calculate the mean of the absolute values of the Δβ values of all CpG sites in the hDMR region;

[0224] If the mean of the absolute values of the Δβ values of all CpG sites in the hDMR region is greater than or equal to 0.1, the newly inserted CpG site is retained in the hDMR region, and the retained CpG site is deleted from the i-th window, and step 2 is executed;

[0225] If the mean of the absolute values of the Δβ values is less than 0.1, the hDMR region does not retain the newly inserted CpG site, and the hDMR region is no longer expanded; step 414 is executed;

[0226] 4) Determine whether the number of CpG sites in the hDMR region is greater than or equal to 5;

[0227] If the number of CpG sites is greater than or equal to 5, retain the hDMR region; and execute step 414;

[0228] If the number of CpG sites is less than 5, delete the hDMR region; and execute step 414;

[0229] Step 414: In the i-th window, the point with the largest absolute value of the new Δβ is used as the seed site, and the hDMR region is expanded in both directions with the seed site as the center. The specific process is as follows:

[0230] In the i-th window, the CpG site placed in the hDMR region in step 413 is deleted (regardless of whether the hDMR region is ultimately retained, as long as the CpG site corresponding to the hDMR region is placed, it is deleted);

[0231] The CpG site corresponding to the maximum absolute value of the Δβ value in the i-th window after deletion is selected as the seed site;

[0232] Repeat steps 413 to 414 until the number of CpG sites in the i-th window is less than 5, and then execute step 415;

[0233] One window may not obtain an hDMR region, may obtain one hDMR region, or may obtain multiple hDMR regions;

[0234] Step 415: Let i=i+1, repeat steps 413 to 414, and then execute step 416;

[0235] Step 416: Repeat step 415 until all windows are calculated and all hDMR regions are obtained;

[0236] Each hDMR region is defined as a candidate hDMR region;

[0237] Step 42: Perform Wilcoxon test and multiple statistical test on all candidate hDMR regions obtained in step 41 to obtain final candidate hDMR regions.

[0238] The other steps and parameters are the same as those in the specific implementation modes 1 to 8-1.

[0239] Specific embodiment ten: This embodiment differs from any one of specific embodiments one to nine in that: in step 42, Wilcoxon test and multiple statistical test are performed on all candidate hDMR regions obtained in step 41 to obtain the final candidate hDMR regions; the specific process is as follows:

[0240] Step 421: Perform a Wilcoxon test on all candidate hDMR regions to obtain p-values, and select candidate hDMR regions corresponding to p<0.05 (one candidate hDMR region corresponds to one p-value). The specific process is as follows:

[0241] The methylation level of haplotype 1 for each CpG site in all candidate hDMR regions is β hap1 and haplotype 2 methylation levels β hap2 The Wilcoxon test was performed to obtain the p value, and the candidate hDMR regions corresponding to p < 0.05 were selected (one candidate hDMR region corresponds to one p value);

[0242] Step 422: Perform multiple verification and correction on the candidate hDMR region selected in step 421 to obtain the final candidate hDMR region. The specific process is as follows:

[0243] The Benjamini-Hochberg method is used to correct the p-values of the candidate hDMR regions selected in step 421, and candidate hDMR regions with p<0.05 are retained to obtain the final candidate hDMR regions.

[0244] Ensure that the hDMRs finally identified are statistically significant and biologically reliable.

[0245] The other steps and parameters are the same as those in the specific implementation modes 1 to 9-1.

[0246] The present invention may have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A method for identifying haplotype differentially methylated regions based on long-read sequencing, characterized by: The specific process of the method is: Step 1: Obtain the sorted and aligned BAM files for each person; Each person corresponds to 1 BAM file; Each BAM file contains M SNP sites; Step 2: Clean the sorted and aligned BAM files to obtain files with CpG site information of methylation status; Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information; Step 3: Input the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 into the trained long short-term memory network LSTM. The trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information; Step 4: Use the maximum likelihood regression method to identify the haplotype differential methylation regions of the typing results of each CpG site in the file output in step 3.

2. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 1, characterized in that: In step 1, the sorted and aligned BAM files of each person are obtained; Each BAM file contains M SNP sites; 100000≤M≤10000000; The specific process is: Step 11, obtaining the human whole genome Fast5 electrical signal file from the sequencing chip R9 or R10 in the third-generation sequencing platform Oxford Nanopore; Step 12: Use Guppy software to convert the Fast5 electrical signal file into a Fastq format file; Step 13. Set Nanopolish software parameters: index; Use Nanopolish software to create an index connection between the Fast5 electrical signal file in step 11 and the Fastq format file in step 12; Step 14. Set the minimap2 software parameters: -a, -x; the value of -x is map-ont; The Fastq format file of the human reference genome GRch38 and the index connected in step 13 is input into the minimap2 software. The minimap2 software converts the Fastq format file generated in step 12 into an aligned BAM file; Step 15. Set the Samtools software parameters: sort; Use Samtools software to sort the BAM file in step 14 to generate a sorted and aligned BAM file; Step 16. Set the Samtools software parameters: index; Use Samtools software to index the BAM file from step 15.

3. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 2, characterized in that: In step 2, the sorted and aligned BAM files are cleaned to obtain files with methylation status CpG site information; Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information; The specific process is: Step 21: Clean the sorted and aligned BAM files to obtain files with methylation status CpG site information; The specific process is: Step 211: Remove low-quality reads from the sorted and aligned BAM files; the process is as follows: Set Samtools software parameters: -q; Samtools software was used to remove reads with alignment quality lower than 30 in the sorted and aligned BAM files; Step 212: Remove the PCR duplicate fragments in the sorted and aligned BAM files retained in step 211; the process is: Set Picard Tools software parameters: MarkDuplicates, REMOVE_DUPLICATES, The value of REMOVE_DUPLICATES is true; Use Picard Tools software to identify and remove PCR duplicate fragments in the sorted and aligned BAM files retained in step 211; Step 213: Set the Nanopolish software parameter: --methylation, where the value of --methylation is CpG; Input the sorted and aligned BAM files retained in step 212, the Fastq format files indexed in step 13, and the human reference genome GRch38 into Nanopolish software, and Nanopolish software outputs a file with methylation status CpG site information; Step 22: Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information.

4. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 3, characterized in that: In step 22, based on the file with the methylation status CpG site information, the features of each CpG site in the file with the methylation status CpG site information are extracted; The specific process is: Step 221: Count the number of methylation states of 0 and 1 for each CpG site in the file with methylation status CpG site information obtained in step 21; The methylation level β of each CpG site is calculated based on the methylation status corresponding to each CpG site. j ; Step 222: Find the number of the first allele G1 and the number of the second allele G2 corresponding to the position information of each CpG in the file with the methylation status CpG site information saved in step 213 in the sorted and aligned BAM file saved in step 212; Step 223: defining each CpG in the file with methylation status CpG site information saved in step 213 as a central CpG; The number of CpGs in the 2kb region upstream and downstream of the central CpG and the methylation level information of each CpG in the 2kb region upstream and downstream of the central CpG are saved as a 2kb information file corresponding to the CpG site; Step 224: Calculate the mean and variance of the methylation level β values of all CpG sites in the 2 kb information file corresponding to the CpG sites saved in step 223; Step 225: Calculate the relative position Pos of each CpG site and the central CpG site in the 2kb information file corresponding to the CpG site saved in step 223 i ; Step 226: Compare the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223 with the human reference genome GRch38 to determine whether the central CpG site is located in the exon, intron, promoter, or enhancer region of the human reference genome GRch38, and perform annotation. Step 227: The methylation level β of each CpG site obtained in step 221 is converted to j , the number of the first allele G1 and the number of the second allele G2 found in step 222, the 2kb information file corresponding to the CpG site saved in step 223, the mean and variance of the methylation level β values of all CpG sites in the 2kb information file corresponding to the CpG site calculated in step 224, and the relative position Pos of each CpG site to the central CpG site in the 2kb information file corresponding to the CpG site calculated in step 225 i The annotation result of step 226 is used as the feature of each CpG site in the file with methylation status CpG site information.

5. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 4, characterized in that: In step 221, the number of methylation states of each CpG site in the file with methylation status CpG site information obtained in step 21 is counted as 0 and 1; The methylation level β of each CpG site is calculated based on the methylation status corresponding to each CpG site. j ; Expressed as: in, β j is the methylation level of the j-th CpG site: M c is the number of methylation states of 1 corresponding to the j-th CpG site; M T is the number of methylation states of 0 corresponding to the j-th CpG site; The methylation status of a CpG site is 0, indicating unmethylation; The methylation status of a CpG site is 1, indicating methylation.

6. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 5, characterized in that: In step 225, the relative position Pos of each CpG site and the central CpG site in the 2kb information file corresponding to the CpG site saved in step 223 is calculated. i ; The specific process is: in, CpG center The position of the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223; CpG i The position of the i-th CpG site excluding the central CpG site in the 2 kb information file corresponding to the CpG site saved in step 223 .

7. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 6, characterized in that: In step 3, the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 are input into the trained long short-term memory network LSTM, and the trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information; The specific process is: Step 31: Get the training set; Step 32: Using the features of each CpG site in the file with methylation status CpG site information in the training set as the input of the long short-term memory network model, and using the typing results of each CpG site in the file with methylation status CpG site information in the training set as the output of the long short-term memory network model; The cross-entropy loss function is used to optimize the parameters of the long short-term memory network model, and the Adam optimizer is combined to perform gradient updates until the cross-entropy loss function converges to obtain a trained long short-term memory network model. Step 33: Input the features of each CpG site in the file containing the methylation status CpG site information extracted in step 2 into the trained long short-term memory network LSTM, and the trained long short-term memory network model outputs the typing results of each CpG site in the file containing the methylation status CpG site information.

8. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 7, characterized in that: In step 31, a training set is obtained; The specific process is: Step 311: Obtain the sorted and aligned BAM files of each person; Each person corresponds to 1 BAM file; Each BAM file contains M SNP sites; Step 312: performing data cleaning on the sorted and aligned BAM files to obtain a file with methylation status CpG site information; Based on the file with the methylation status CpG site information, extract the features of each CpG site in the file with the methylation status CpG site information; Step 313: Using the sequencing tool NanoMethPhase, each CpG site in the file containing the methylation status CpG site information is typed into haplotype 1 and haplotype 2, thereby obtaining a CpG typing file.

9. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 8, characterized in that: In step 4, the maximum likelihood regression method is used to identify the haplotype differential methylation region hDMR of the typing result of each central CpG site in the methylation status CpG site information file output in step 3; The specific process is: Step 41: Find candidate hDMR blocks; The specific process is: Step 411: Divide all CpG sites in the CpG typing file output in step 3 into windows, and obtain a seed site in each window; The specific process is: Use a 4 kb window to slide all CpG sites in the CpG typing file output in step 3, and select the CpG site with the largest Δβ value in each window as the seed site; The calculation process of the Δβ value of the CpG site is: Methylation level β of haplotype 1 based on each CpG site hap1 and haplotype 2 methylation levels β hap2 Calculate Δβ = β for each CpG site hap1 -β hap2 ; Step 412: Set i=1; Step 413: Expand the hDMR region toward both sides with the seed site as the center within the i-th window; The specific process is: 1) Copy the seed site in the i-th window and place it in the empty hDMR region, and delete the seed site from the i-th window. If the absolute value of the Δβ value of the deleted seed site is greater than or equal to 0.1, retain the seed site in the hDMR region and execute 2); Copy the seed site in the i-th window and place it in the empty hDMR region, and delete the seed site from the i-th window. If the absolute value of the Δβ value of the deleted seed site is less than 0.1, the seed site is not retained in the hDMR region, and the hDMR region is not expanded. Execute step 414. 2) Calculate the Δβ value of the CpG site adjacent to the left of the deleted seed site and the Δβ value of the CpG site adjacent to the right of the deleted seed site in the i-th window; If the Δβ value of the CpG site adjacent to the left of the deleted seed site is greater than the Δβ value of the CpG site adjacent to the right of the deleted seed site, the distance between the CpG site adjacent to the left of the seed site and the left boundary of the hDMR region is calculated; if the distance is less than 500 bp, the CpG site adjacent to the left of the deleted seed site is placed in the leftmost region of the hDMR, and 3) is executed; if the distance is greater than or equal to 500 bp, the CpG site adjacent to the left is retained in the i-th window and is not placed in the leftmost region of the hDMR, and the left side of the hDMR region is no longer extended; Calculate the distance between the CpG site adjacent to the right of the seed site and the right edge of the hDMR region. If the distance is less than 500 bp, place the CpG site adjacent to the right of the deleted seed site in the rightmost region of the hDMR and proceed to step 3). If the distance is greater than or equal to 500 bp, retain the CpG site adjacent to the right in the i-th window and do not place it in the rightmost region of the hDMR. The right side of the hDMR region is no longer extended, and proceed to step 4). If the Δβ value of the CpG site adjacent to the left of the deleted seed site is less than or equal to the Δβ value of the CpG site adjacent to the right of the deleted seed site, the distance between the CpG site adjacent to the right of the seed site and the right boundary of the hDMR region is calculated. If the distance is less than 500 bp, the CpG site adjacent to the right of the deleted seed site is placed in the rightmost region of the hDMR, and step 3) is executed. If the distance is greater than or equal to 500 bp, the CpG site adjacent to the right is retained in the i-th window and is not placed in the rightmost region of the hDMR. In this case, the right side of the hDMR region is no longer extended. Calculate the distance between the CpG site adjacent to the left of the seed site and the left edge of the hDMR region. If the distance is less than 500 bp, place the CpG site adjacent to the left of the deleted seed site in the leftmost region of the hDMR and proceed to step 3). If the distance is greater than or equal to 500 bp, retain the CpG site adjacent to the left in the i-th window and do not place it in the leftmost region of the hDMR. The left side of the hDMR region is no longer extended, and proceed to step 4). 3) Calculate the mean of the absolute values of the Δβ values of all CpG sites in the hDMR region; If the mean of the absolute values of the Δβ values of all CpG sites in the hDMR region is greater than or equal to 0.1, the newly inserted CpG site is retained in the hDMR region, and the retained CpG site is deleted from the i-th window, and step 2 is executed; If the mean of the absolute values of the Δβ values is less than 0.1, the hDMR region does not retain the newly inserted CpG site, and the hDMR region is no longer expanded; step 414 is executed; 4) Determine whether the number of CpG sites in the hDMR region is greater than or equal to 5; If the number of CpG sites is greater than or equal to 5, retain the hDMR region; and execute step 414; If the number of CpG sites is less than 5, delete the hDMR region; and execute step 414; Step 414: In the i-th window, the point with the largest absolute value of the new Δβ is used as the seed site, and the hDMR region is expanded toward both sides with the seed site as the center; The specific process is: Deleting the CpG sites placed in the hDMR region in step 413 within the i-th window; The CpG site corresponding to the maximum absolute value of the Δβ value in the i-th window after deletion is selected as the seed site; Repeat steps 413 to 414 until the number of CpG sites in the i-th window is less than 5, and then execute step 415; Step 415: Let i=i+1, repeat steps 413 to 414, and then execute step 416; Step 416: Repeat step 415 until all windows are calculated and all hDMR regions are obtained; Each hDMR region is defined as a candidate hDMR region; Step 42: Perform Wilcoxon test and multiple statistical test on all candidate hDMR regions obtained in step 41 to obtain final candidate hDMR regions.

10. The method for identifying haplotype differentially methylated regions based on long-read sequencing according to claim 9, characterized in that: In step 42, Wilcoxon test and multiple statistical test are performed on all candidate hDMR regions obtained in step 41 to obtain final candidate hDMR regions; The specific process is: Step 421: Perform a Wilcoxon test on all candidate hDMR regions to obtain p-values, and select candidate hDMR regions corresponding to p < 0.05; The specific process is: The methylation level of haplotype 1 for each CpG site in all candidate hDMR regions is β hap1 and haplotype 2 methylation levels β hap2 The Wilcoxon test was performed to obtain the p value, and the candidate hDMR regions corresponding to p < 0.05 were selected; Step 422: Perform multiple verification and correction on the candidate hDMR region selected in step 421 to obtain a final candidate hDMR region; The specific process is: The Benjamini-Hochberg method is used to correct the p-values of the candidate hDMR regions selected in step 421, and candidate hDMR regions with p<0.05 are retained to obtain the final candidate hDMR regions.

Images