Diploid typing method, system and assembly method based on hi-c data

Abstract

The application discloses a diploid typing method based on Hi-C data, comprising obtaining an existing typing assembly result; mapping the assembly result and obtaining a pure region and a hybrid region of a diploid genome; filtering a comparison result of the Hi-C reading and reserving a Hi-C comparison result containing a hybrid variation site; identifying a large section switching existing in the assembly result; and counting the number of Hi-C connections between the sections after cutting, thereby completing the typing of the diploid. The application further discloses a system for implementing the diploid typing method based on the Hi-C data, and an assembly method comprising the diploid typing method based on the Hi-C data. The application filters the Hi-C comparison result through the hybrid variation in the diploid genome, detects whether a large parent switching exists by using the filtered Hi-C comparison result, thereby realizing the final diploid typing method and assembly, and therefore the application has higher reliability and better accuracy.

Description

Technical Field

[0001] This invention belongs to the field of molecular biology, specifically relating to a diploid typing method, system, and assembly method based on Hi-C data. Background Technology

[0002] Haplotype-resolved genome assembly holds immense potential for advancing various fields of basic research. Recent advancements in long-read sequencing technologies have yielded numerous high-quality de novo genome assemblies. The Human Pangenome Reference Consortium recently published an initial draft of the Human Pangenome Reference, which includes 47 well-generated diploid assemblies. Currently, many long-read genome assembly tools tend to merge different haplotypes into a single consensus assembly without requiring additional sequencing data. However, this merging results in different haplotypes being mosaicked together, containing erroneous linkages not present in any haplotype, thus compromising the integrity of biological inferences. Therefore, some assembly tools attempt to address this issue by utilizing heterozygous variations between haplotypes to maintain local genotyping. However, due to limitations in low-heterozygous samples, these tools often produce genotypic blocks of finite length. To obtain longer contigs, these tools link shorter genotypic blocks together to form a primary assembly. Nevertheless, haplotypes within each contig may be randomly switched. The remaining genotypic blocks, also known as haplotypes, are relegated to alternative assemblies.

[0003] To generate haplotype-resolved genome assemblies from diploid samples, various tools require additional sequencing data, such as parental short-read data, template-ligated strand sequencing data, or Hi-C data. Trio-binning methods, such as TrioCanu and hifiasm(trio), utilize short-read data from the parents to identify parent-specific k-mers of length k. Using these k-mers, long-read sequence data from offspring can be separated into paternal and maternal data. While the genotypic assemblies generated by trio-binning methods are accurate, their practical application is limited by parental sequencing data, which is often unavailable. To address this limitation, single-sample haplotype-resolved assembly methods using template-ligated strand sequencing data or Hi-C data have emerged. Strand sequencing-based methods preserve the structural continuity of individual homologs. However, when extended to multiple species or individuals within a single species, the requirements of living cells and a single cycle of cell division pose challenges. On the other hand, proximity-linked data types such as Hi-C provide extensive connectivity information for genotyping across longer genomes and can be sequenced from the same sample. Hi-C data has proven to be a more general genotyping solution.

[0004] Using Hi-C data to genotype primary and alternate assemblies and extend them into two homologous haplotypes is an effective solution, especially when dealing with partially genotyped diploid assemblies generated by numerous assembly tools. The FALCON-Phase tool is used for primary / alternate patterns, utilizing Hi-C data for genotyping and extension. However, it skips SNP identification, which can affect genotyping accuracy. This method uses alignment quality and edit distance as filters to exclude uninformative Hi-C alignments. However, this filtering process also discards many informative Hi-C alignments, resulting in a significant reduction in the number of Hi-C alignments and thus lowering genotyping accuracy. Furthermore, FALCON-Phase lacks a mechanism to handle parental switching within contigs, further weakening its genotyping performance. During genotyping, FALCON-Phase groups contigs according to the primary contig and processes each group separately. For a given block, FALCON-Phase relies only on the genotyping information preceding it within the group and Hi-C connections. This method leads to the loss of valuable fractal information and Hi-C connection information behind the blocks, ultimately affecting the overall fractal accuracy. Summary of the Invention

[0005] One of the objectives of this invention is to provide a highly reliable and accurate diploid typing method based on Hi-C data.

[0006] The second objective of this invention is to provide a system for implementing the diploid typing method based on Hi-C data.

[0007] A third objective of this invention is to provide an assembly method that includes the diploid typing method based on Hi-C data.

[0008] The diploid typing method based on Hi-C data provided by this invention includes the following steps:

[0009] S1. Obtain the assembly results of the existing subtypes;

[0010] S2. Map the assembly results obtained in step S1, and obtain the homozygous and heterozygous regions of the diploid genome according to the mapping relationship;

[0011] S3. Based on the heterozygous variations in the diploid genome, filter the Hi-C reading alignment results and retain the Hi-C alignment results containing heterozygous variation sites;

[0012] S4. Using the filtered Hi-C alignment results obtained in step S3, identify large-segment switching in the assembly results;

[0013] S5. Based on the data obtained in step S4, count the number of Hi-C connections between the cut fragments to complete the typing of diploids.

[0014] Step S2 involves mapping the assembly results obtained in step S1 and obtaining the homozygous and heterozygous regions of the diploid genome based on the mapping relationship. Specifically, this includes the following steps:

[0015] The alternative assemblies in the diploid assembly results are mapped to the primary assembly, and based on the mapping relationship, the primary assembly is divided into homozygous and heterozygous regions:

[0016] Map each contig in the alternative assembly to the contig in the primary assembly, and obtain the mapping fragment of each alternative contig on the primary contig.

[0017] For alternative contigs that have at least two mapping segments on the primary contig: sort the mapping segments in ascending order of their position on the primary contig, and record the mapping length of the i-th mapping segment as |S|. i |, i∈1,2,...,n, where n is the number of mapping segments;

[0018] f(i) is calculated using the following formula:

[0019]

[0020] In the formula, f(i) is the maximum mapping length to the i-th mapping segment; [S i ] represents the mapping length; gap(j,i) represents the interval between the i-th mapping segment and the j-th mapping segment;

[0021] Finally, the combination of segments that maximizes f(x) is found as the mapping region of the alternative contig on the primary contig, which is a heterozygous region. Regions on the primary contig that have no mapping relationship with the alternative contig are homozygous regions, thus dividing the primary assembly into homozygous and heterozygous regions.

[0022] Step S3, which involves filtering the Hi-C readings based on heterozygous variations in the diploid genome and retaining Hi-C alignment results containing heterozygous variation sites, specifically includes the following steps:

[0023] Based on heterozygous variations in the diploid genome, determine whether each pair of Hi-C readings contains a variant site:

[0024] The current Hi-C readings are mapped to the alternative assembly and the primary assembly, respectively, to identify heterozygous variants, and heterozygous variants that correspond to each other are retained according to the mapping relationship;

[0025] The Hi-C readings are aligned to the homozygous and heterozygous regions obtained in step S2, with only one optimal alignment result retained for each Hi-C reading. The optimal alignment result is defined as the one with the highest alignment score and the lowest error rate, denoted as M. b ;

[0026] In the retained comparison results, filtering is performed again: Hi-C readings that do not meet both conditions a and b are deleted, and the remaining comparison results are retained;

[0027] Condition a: The alignment quality is greater than q and the edit distance is less than e; q is the quality threshold and e is the distance threshold.

[0028] Condition b: Covers at least one heterozygous variant site;

[0029] For a pair of Hi-C readings, if the alignment result of one reading is filtered out, then the alignment result of the other reading will also be filtered out; the final alignment result is the Hi-C alignment result that includes heterozygous variant sites.

[0030] The quality threshold q is set to 1; the distance threshold e is set to 5.

[0031] Step S4, which uses the filtered Hi-C alignment results obtained in step S3 to identify large-scale switching in the assembly results, specifically includes the following steps:

[0032] Based on the filtered Hi-C alignment results obtained in step S3, plot the coverage curve of each candidate alternate contig and its corresponding heterozygous region for each site; the Hi-C pair coverage is defined as the area covered from the leftmost position of the Hi-C reading pair to the rightmost position of the Hi-C reading pair.

[0033] The following rules are used to determine large-scale handovers:

[0034] If there is a notch on the coverage curve, it is determined that there is a large segment of switching on the alternative contig;

[0035] The notch is defined as: the minimum point of the curve, and the coverage of the minimum point is less than the average coverage set ratio c.

[0036] The coverage ratio c is set to 1 / 15.

[0037] Step S5, which involves counting the number of Hi-C connections between the segments after slicing based on the data obtained in step S4, to complete the genotyping of diploid individuals, specifically includes the following steps:

[0038] For any two alternative contig fragments a and b, and their corresponding heterozygous regions a' and b', the following formula is used as the objective function for fractal analysis:

[0039] f(θ a ,θ b )=θa θ b (w ab +w a'b' -(w a'b +w ab' ))

[0040] In the formula θ a It is the haplotype of fragment a; θ b For the haplotype of segment b, θ a and θ b The value can be 1 or -1, representing the parent or mother lineage respectively; w ab w represents the Hi-C connection number of a and b. a'b' w represents the Hi-C connection number of a' and b'. a'b w represents the Hi-C connection number of a' and b. ab' Let a and b' be the Hi-C connection numbers;

[0041] For all segments to be segmented, the following formula is used as the objective function:

[0042]

[0043] In the formula For all segments to be classified, the haplotype is denoted as ;

[0044] Finally, the objective function is obtained using a random algorithm. The maximum value, the corresponding haplotype of each segment is its final haplotype, and the segments from the same parent are obtained;

[0045] The number of Hi-C connections between fragments originating from the same parent is greater than a first preset value, while the number of connections between fragments from different parents is less than a first preset value. Fragments belonging to the heterozygous region originate from either the father or mother parent. Based on the number of Hi-C connections, fragments belonging to the heterozygous region are classified into two categories. The number of Hi-C connections between fragments of the same category is greater than a second preset value, while the number of Hi-C connections between fragments of different categories is less than a second preset value. Fragments in the homozygous region do not carry parent information and are not classified.

[0046] The haplotype mentioned above is specifically obtained by using the following steps:

[0047] (1) Randomly set the haplotype of each segment;

[0048] (2) Randomly select a segment to be fractalized and change its haplotype: If the objective function If the value can be increased, the change is retained; otherwise, it is not retained.

[0049] (3) Repeat step (2) until the objective function is achieved. The value reaches a local maximum, and the following is determined:

[0050] If the local maximum value is greater than the current global maximum value, then the current global maximum value is set as the local maximum value.

[0051] (4) Randomly change the haplotype of a set of segments to be classified and the segments connected to them, and use step (2) to obtain a new local maximum.

[0052] (5) Repeat steps (2) to (4) until the set number of times to obtain the optimal local maximum value;

[0053] (6) Make the objective function When the global maximum value is reached, the haplotype of each segment is taken as the final haplotype of that segment, thus obtaining the final haplotype of the segment.

[0054] This invention also provides a system for implementing the diploid genotyping method based on Hi-C data, comprising a data acquisition module, a region acquisition module, a readout filtering module, a switch identification module, and a genotyping module; the data acquisition module, region acquisition module, readout filtering module, switch identification module, and genotyping module are connected in series; the data acquisition module is used to acquire the assembly results of existing genotypes and upload the data information to the region acquisition module; the region acquisition module is used to map the acquired assembly results according to the received data information, and obtain the homozygous and heterozygous regions of the diploid genome according to the mapping relationship, and upload the data information to the readout filtering module; the readout filtering module is used to filter the Hi-C readout alignment results according to the received data information and the heterozygous variations in the diploid genome, retain the Hi-C alignment results containing heterozygous variation sites, and upload the data information to the switch identification module; the switch identification module is used to identify large-segment switches in the assembly results according to the received data information and the obtained filtered Hi-C alignment results, and upload the data information to the genotyping module; the genotyping module is used to count the number of Hi-C connections between the cut fragments according to the received data information, thereby completing the diploid genotyping.

[0055] The present invention also provides an assembly method including the diploid typing method based on Hi-C data, further comprising the following steps:

[0056] S6. Based on the typing results obtained in step S5, the fragments from the same parent are clustered and connected into an assembly result to obtain the final diploid assembly result.

[0057] Step S6, based on the typing results obtained in step S5, involves clustering fragments from the same parent and connecting them into a single assembly result to obtain the final diploid assembly result. This specifically includes the following steps:

[0058] Based on the genotyping results of each fragment obtained in step S5, fragments with the same haplotype are fragments from the same parent; fragments with the same haplotype are grouped into one category, resulting in two categories of fragments, from the father and mother respectively;

[0059] For each class of fragments, heterozygous and homozygous fragments are sequentially connected according to the position of the fragment on the primary contig or the mapping position of the fragment on the primary contig, to obtain a new contig, which is the final diploid assembly result.

[0060] The diploid typing method, system, and assembly method based on Hi-C data provided by this invention filters Hi-C alignment results through heterozygous variations in the diploid genome, and then uses the filtered Hi-C alignment results to detect whether there is a large segment of parental switching, thereby realizing the final diploid typing method and assembly. Therefore, this invention has higher reliability and better accuracy. Attached Figure Description

[0061] Figure 1 This is a schematic diagram of the method flow for the classification method of the present invention.

[0062] Figure 2 This is a schematic diagram of the functional modules of the system of the present invention.

[0063] Figure 3 This is a schematic diagram of the assembly method of the present invention. Detailed Implementation

[0064] First, the terminology used in this invention is explained as follows:

[0065] Long-read sequencing, also known as third-generation sequencing, is a next-generation sequencing technology that contrasts with traditional Sanger sequencing and next-generation sequencing. It is characterized by longer read lengths and higher error rates. Representative long-read sequencing platforms include the PacBio Sequel II, PacBio Revio, Nanopore MinION, and Nanopore GridION sequencers. Long-read sequencing technology is currently widely used in various bioinformatics fields, including de novo assembly, variant detection, gene expression analysis, and disease research.

[0066] Hi-C (Hi-Cell) is a molecular biology technique used to study the three-dimensional structure of chromosomes and the spatial organization of the genome. Hi-C is a high-throughput sequencing technology capable of detecting and interpreting the interactions and spatial locations between different regions within chromosomes. Hi-C technology utilizes steps such as cross-linking, cutting, ligation, and sequencing to explore the spatial structure of chromosomes in the cell nucleus. First, chromosomes within the cell are cross-linked, connecting closely spaced DNA fragments together. Next, the DNA is cut into smaller fragments by enzymes, and the cross-linked DNA fragments are ligated to form a new DNA molecule containing sequences from different sites. Subsequently, these ligated DNA fragments are sequenced to obtain the sequences of interacting DNA fragments.

[0067] Haplotype analysis: Haplotypes typically refer to the combinations of identical alleles on homologous chromosomes at multiple consecutive loci on a chromosome. Haplotype analysis involves analyzing genomic data to determine the combinations of alleles at specific loci in an individual or population. Haplotype analysis provides accurate and reliable haplotype information and has wide applications in disease research, drug development, and personalized medicine.

[0068] Heterozygous variation: At a specific gene locus, an individual possesses two different alleles. Every human cell has two sets of chromosomes, and unless it's on the sex chromosomes, each gene locus has two alleles inherited from both parents. When these two alleles are different, meaning the genes from the parents are not entirely identical, this situation is called heterozygosity.

[0069] Coverage: This typically refers to the number of times or the depth to which a specific region or the entire genome in sequencing data has been sequenced. Coverage can be used to measure the depth or completeness of genome sequencing.

[0070] Parental switching: During diploid genome assembly, due to a lack of sufficient parental information, fragments from the father and mother are incorrectly spliced ​​together, forming a contig of parental mixing.

[0071] Diploid typing refers to the process in diploid organisms of identifying the combination of alleles from multiple consecutive genomes on homologous chromosomes to generate two haploid genomes, one from the father and one from the mother.

[0072] like Figure 1 The diagram shown is a flowchart of the typing method of the present invention: This diploid typing method based on Hi-C data disclosed in the present invention includes the following steps:

[0073] S1. Obtain the assembly results of the existing subtypes;

[0074] S2. Map the assembly results obtained in step S1, and obtain the homozygous and heterozygous regions of the diploid genome based on the mapping relationship; specifically including the following steps:

[0075] The alternative assemblies in the diploid assembly results are mapped to the primary assembly, and based on the mapping relationship, the primary assembly is divided into homozygous and heterozygous regions:

[0076] Map each contig in the alternative assembly to the contig in the primary assembly, and obtain the mapping fragment of each alternative contig on the primary contig.

[0077] For a candidate alternate contig with at least two mapping segments on the primary contig: sort the mapping segments in ascending order of their position on the primary contig, and record the mapping length of the i-th mapping segment as |S|. i |, i∈1,2,...,n, where n is the number of mapping segments;

[0078] f(i) is calculated using the following formula:

[0079]

[0080] In the formula, f(i) is the maximum mapping length to the i-th mapping segment; [S i ] represents the mapping length; gap(j,i) represents the interval between the i-th mapping segment and the j-th mapping segment;

[0081] Finally, the combination of segments that maximizes f(x) is found, which is the mapping region of the alternative contig on the primary contig, i.e., the heterozygous region. The region on the primary contig that has no mapping relationship with the alternative contig is the homozygous region. Thus, the primary assembly is divided into homozygous and heterozygous regions.

[0082] S3. Based on heterozygous variations in the diploid genome, filter the Hi-C reading alignment results and retain Hi-C alignment results containing heterozygous variation sites; specifically including the following steps:

[0083] Based on heterozygous variations in the diploid genome, determine whether each pair of Hi-C readings contains a variant site:

[0084] The current Hi-C readings are mapped to the alternative assembly and the primary assembly, respectively, to identify heterozygous variants, and heterozygous variants that correspond to each other are retained according to the mapping relationship;

[0085] The Hi-C readings are aligned to the homozygous and heterozygous regions obtained in step S2, with only one optimal alignment result retained for each Hi-C reading. The optimal alignment result is defined as the one with the highest alignment score and the lowest error rate, denoted as M. b ;

[0086] In the retained comparison results, filtering is performed again: Hi-C readings that do not meet both conditions a and b are deleted, and the remaining comparison results are retained;

[0087] Condition a: The alignment quality is greater than q and the edit distance is less than e; q is the quality threshold (preferably 1), and e is the distance threshold (preferably 5);

[0088] Condition b: Covers at least one heterozygous variant site;

[0089] For a pair of Hi-C readings, if the alignment result of one reading is filtered out, then the alignment result of the other reading will also be filtered out; the final alignment result is the Hi-C alignment result that includes heterozygous variant sites.

[0090] S4. Using the filtered Hi-C alignment results obtained in step S3, identify large-segment switching in the assembly results; specifically, this includes the following steps:

[0091] Based on the filtered Hi-C alignment results obtained in step S3, plot the coverage curve of each candidate alternate contig and its corresponding heterozygous region for each site; the Hi-C pair coverage is defined as the area covered from the leftmost position of the Hi-C reading pair to the rightmost position of the Hi-C reading pair.

[0092] The following rules are used to determine large-scale handovers:

[0093] If there is a notch on the coverage curve, it is determined that there is a large segment of switching on the candidate alternate contig.

[0094] The notch is defined as: the minimum point of the curve, and the coverage of the minimum point is less than the average coverage set ratio c; the preferred value of the coverage ratio c is 1 / 15.

[0095] S5. Based on the data obtained in step S4, count the number of Hi-C connections between the segments after cutting, thereby completing the genotyping of diploids; specifically including the following steps:

[0096] For any two alternative contig fragments a and b, and their corresponding heterozygous regions a' and b', the following formula is used as the objective function for fractal analysis:

[0097] f(θ a ,θ b )=θ a θ b (w ab +w a'b' -(w a'b +w ab' ))

[0098] In the formula θ a It is the haplotype of fragment a; θ b For the haplotype of segment b, θ a and θ b The value can be 1 or -1, representing the parent or mother lineage respectively; w ab w represents the Hi-C connection number of a and b. a'b' w represents the Hi-C connection number of a' and b'. a'b w represents the Hi-C connection number of a' and b. ab' Let a and b' be the Hi-C connection numbers;

[0099] For all segments to be segmented, the following formula is used as the objective function:

[0100]

[0101] In the formula For all segments to be classified, the haplotype is denoted as ;

[0102] Finally, the objective function is obtained using a random algorithm. The maximum value, the haplotype of each corresponding segment is its final haplotype, and the segments from the same parent are obtained;

[0103] The number of Hi-C connections between fragments originating from the same parent is greater than a first preset value, while the number of connections between fragments from different parents is less than a first preset value. In practice, the number of Hi-C connections between fragments originating from the same parent will be much greater than the number of connections between fragments from different parents. Fragments belonging to the heterozygous region originate from either the father or mother parent. Based on the number of Hi-C connections, fragments belonging to the heterozygous region are classified into two categories. The number of Hi-C connections between fragments of the same category is greater than a second preset value, while the number of Hi-C connections between fragments of different categories is less than a second preset value. Fragments in the homozygous region do not carry parent information and are not classified.

[0104] The haplotype mentioned above is specifically obtained by using the following steps:

[0105] (1) Randomly set the haplotype of each segment;

[0106] (2) Randomly select a segment to be fractalized and change its haplotype: If the objective function If the value can be increased, the change is retained; otherwise, it is not retained.

[0107] (3) Repeat step (2) until the objective function is achieved. The value reaches a local maximum, and the following is determined:

[0108] If the local maximum value is greater than the current global maximum value, then the current global maximum value is set as the local maximum value.

[0109] (4) Randomly change the haplotype of a set of segments to be classified and the segments connected to them, and use step (2) to obtain a new local maximum.

[0110] (5) Repeat steps (2) to (4) until the set number of times to obtain the optimal local maximum value;

[0111] (6) Make the objective function When the global maximum value is reached, the haplotype value of each segment is the final haplotype of that segment, thus obtaining the final haplotype of the segment.

[0112] The present invention will be further described below with reference to an embodiment:

[0113] The input assembly results are the human (HG002) genome, whole genome sequencing reads, and corresponding Hi-C sequencing reads assembled using the PECAT assembly tool;

[0114] Determining homozygous and heterozygous regions: In this embodiment, the primary assembly has a total of 571 heterozygous regions and 774 homozygous regions;

[0115] Filtering Hi-C alignment results: In this embodiment, a total of 4,399,346 heterozygous variant sites were identified, and 32,332,595 Hi-C pairs were retained;

[0116] Identifying large-scale parent-gene switching: In this embodiment, parent-gene switching was identified in 6 heterozygous segments;

[0117] Finally, genotyping was performed: the number of Hi-C connections between the cut fragments was counted based on the retained Hi-C alignment results. Since the number of connections between fragments from the same parent is much greater than the number of connections between fragments from different parents, heterozygous fragments were divided into two categories based on this characteristic: one from the paternal parent and one from the maternal parent. Homozygous regions should appear in both the paternal and maternal parents.

[0118] In this embodiment, the obtained genotyping results are compared with the FALCON-Phase method. Specifically, using the genotyping method of this invention and the FALCON-Phase method with the same input data, the sequence attributes of the genotyping sequences, genotyping switching errors, Hamming errors, genotyping errors, and gene integrity are compared. The quality of the diploid genome genotyped in this invention is evaluated using these methods.

[0119] The following compares the human genotyping results generated by eight examples of this invention with the genotyping results obtained by the FALCON-Phase method. The sequencing data reads used in the five examples of this invention were obtained from the Vertebrate Genome Project (https: / / www.genomeark.org / vgp-all / ), the Human Pangenome Project (https: / / humanpangenome.org / ), the NCBI SRA database (https: / / www.ncbi.nlm.nih.gov / sra), and the official Nanopore sequencing website (https: / / nanoporetech.com / ), specifically mHomSap3CLR, HG00733 CLR, HG002 HiFi, NA12878 Nanopore, HG00733 Nanopore, HG002 Nanopore, HG002 simplx, and HG002 duplex. HG002 simplx and HG002 duplex used R10 reagents, while the other three Nanopore data used R9 reagents. Specific information for the eight sequencing datasets is shown in Table 1.

[0120] Table 1. Illustration of sequencing dataset information

[0121]

[0122] Tables 2 and 3, containing five Nanopore datasets and three PacBio datasets respectively, compare the genotyping tools in terms of genome size, N50, genotyping block N50, switching error rate, Hamming error rate, genotyping error rate, and gene integrity. N50 refers to the genome N50 and genotyping block N50 obtained by sorting DNA sequences in haplotypes from longest to shortest, summing the sequences from longest to shortest, and adding the sum of the last summed sequences when the sum reaches 50% of the total length of all sequences. Switching error rate and Hamming error rate were calculated by Merqury using second-generation parental data. Genotyping error rate refers to the percentage of paternal (maternal) fragments that are genotyped as maternal (paternal) fragments. Gene integrity was calculated by BUSCO.

[0123] Table 2. Schematic diagram of the typing results of five Nanopore datasets.

[0124]

[0125] #R10 dataset. The top of each row represents haplotype 1, and the bottom row represents haplotype 2.

[0126] Table 3. A schematic diagram of the typing results for the three PacBio datasets.

[0127]

[0128] *Since the Hi-C alignment results of FALCON-Phase on mHomSap3 CLR and HG002 HiFi data are unavailable, the classification error rate cannot be calculated. The top of each row indicates haplotype 1, and the bottom row indicates haplotype 2.

[0129] As can be seen from the results in Tables 2 and 3, the classification results obtained in this embodiment of the invention have several advantages over FALCON-Phase: the classification blocks are longer, as reflected in classification block N50; the accuracy is higher, as reflected in the Hamming error rate and classification error rate. The Hamming error rate and classification error rate are important indicators of the accuracy of classification tools, and the Hamming error rate and classification error rate of this embodiment of the invention are much lower than those of FALCON-Phase.

[0130] like Figure 2The diagram shows the functional modules of the system of this invention: The system for implementing the diploid genotyping method based on Hi-C data disclosed in this invention includes a data acquisition module, a region acquisition module, a reading filtering module, a switching identification module, and a genotyping module; these modules are connected in series. The data acquisition module acquires the assembly results of existing genotypes and uploads the data information to the region acquisition module. The region acquisition module maps the acquired assembly results according to the received data information and obtains the homozygous and heterozygous regions of the diploid genome based on the mapping relationship. The data is then uploaded to the readout filtering module. The readout filtering module filters the Hi-C readout alignment results based on the received data and heterozygous variations in the diploid genome, retaining Hi-C alignment results containing heterozygous variation sites, and uploads the data to the switching identification module. The switching identification module identifies large-segment switching in the assembly results based on the received data and the filtered Hi-C alignment results, and uploads the data to the typing module. The typing module counts the number of Hi-C connections between the cut fragments based on the received data, thus completing the typing of the diploid.

[0131] like Figure 3 The diagram shown is a flowchart of the assembly method of the present invention: The assembly method disclosed in this invention, which includes the diploid typing method based on Hi-C data, includes the following steps:

[0132] S1. Obtain the assembly results of the existing subtypes;

[0133] S2. Map the assembly results obtained in step S1, and obtain the homozygous and heterozygous regions of the diploid genome according to the mapping relationship;

[0134] S3. Based on the heterozygous variations in the diploid genome, filter the Hi-C reading alignment results and retain the Hi-C alignment results containing heterozygous variation sites;

[0135] S4. Using the filtered Hi-C alignment results obtained in step S3, identify large-segment switching in the assembly results;

[0136] S5. Based on the data obtained in step S4, count the number of Hi-C connections between the segments after cutting, thereby completing the typing of diploids;

[0137] S6. Based on the typing results obtained in step S5, cluster the fragments from the same parent and connect them into an assembly result to obtain the final diploid assembly result; specifically including the following steps:

[0138] Based on the genotyping results of each fragment obtained in step S5, fragments with the same haplotype are fragments from the same parent; fragments with the same haplotype are grouped into one category, resulting in two categories of fragments, from the father and mother respectively;

[0139] For each class of fragments, heterozygous and homozygous fragments are sequentially connected according to the position of the fragment on the primary contig or the mapping position of the fragment on the primary contig, to obtain a new contig, which is the final diploid assembly result.

Claims

1. A diploid typing method based on Hi-C data, comprising the following steps: S1. Obtain the assembly results of the existing subtypes; the assembly results include the main assembly and alternative assemblies; S2. Map the assembly results obtained in step S1, and obtain the homozygous and heterozygous regions of the diploid genome according to the mapping relationship; S3. Based on heterozygous variations in the diploid genome, filter the Hi-C reading alignment results and retain Hi-C alignment results containing heterozygous variation sites; specifically including the following steps: Based on heterozygous variations in the diploid genome, determine whether each pair of Hi-C readings contains a variant site: The current Hi-C readings are mapped to the alternative assembly and the primary assembly, respectively, to identify heterozygous variants, and heterozygous variants that correspond to each other are retained according to the mapping relationship; The Hi-C readings are aligned to the homozygous and heterozygous regions obtained in step S2, with only one optimal alignment result retained for each Hi-C reading. The optimal alignment result is defined as the one with the highest alignment score and the lowest error rate. ; In the retained comparison results, filtering is performed again: Hi-C readings that do not meet both conditions a and b are deleted, and the remaining comparison results are retained; Condition a: The alignment quality is greater than q and the edit distance is less than e; q is the quality threshold and e is the distance threshold. Condition b: Covers at least one heterozygous variant site; For a pair of Hi-C readings, if the alignment result of one reading is filtered out, then the alignment result of the other reading will also be filtered out; the final alignment result is the Hi-C alignment result that includes heterozygous variant sites. S4. Using the filtered Hi-C alignment results obtained in step S3, identify large-scale switching in the assembly results; S5. Based on the data obtained in step S4, count the number of Hi-C connections between the cut fragments to complete the typing of diploids.

2. The diploid typing method based on Hi-C data according to claim 1, characterized in that... Step S2 involves mapping the assembly results obtained in step S1 and obtaining the homozygous and heterozygous regions of the diploid genome based on the mapping relationship. Specifically, this includes the following steps: The alternative assemblies in the diploid assembly results are mapped to the primary assembly, and based on the mapping relationship, the primary assembly is divided into homozygous and heterozygous regions: Map each contig in the alternative assembly to the contig in the primary assembly, and obtain the mapping fragment of each alternative contig on the primary contig. For alternative contigs with at least two mapping segments on the primary contig: sort the mapping segments in ascending order of their position on the primary contig, and record the mapping length of the i-th mapping segment as follows: , , The number of mapped segments; The following formula is used to calculate the result. : In the formula This represents the maximum mapping length up to the i-th mapping segment; The interval between the i-th and j-th mapping segments; Finally, find the one that makes The largest fragment combination serves as the mapping region of the alternative contig on the primary contig, which is a heterozygous region. Regions on the primary contig that have no mapping relationship with the alternative contig are homozygous regions, thus dividing the primary assembly into homozygous and heterozygous regions.

3. The diploid typing method based on Hi-C data according to claim 2, characterized in that... The quality threshold q is set to 1; the distance threshold e is set to 5.

4. The diploid typing method based on Hi-C data according to claim 3, characterized in that... Step S4, which uses the filtered Hi-C alignment results obtained in step S3 to identify large-scale switching in the assembly results, specifically includes the following steps: Based on the filtered Hi-C alignment results obtained in step S3, plot the coverage curve of each candidate alternate contig and its corresponding heterozygous region for each site; the Hi-C pair coverage is defined as the area covered from the leftmost position of the Hi-C reading pair to the rightmost position of the Hi-C reading pair. The following rules are used to determine large-scale handovers: If there is a notch on the coverage curve, it is determined that there is a large segment of switching on the alternative contig; The notch is defined as: the minimum point of the curve, and the coverage of the minimum point is less than the average coverage set ratio c.

5. The diploid typing method based on Hi-C data according to claim 4, characterized in that... Step S5, which involves counting the number of Hi-C connections between the segments after slicing based on the data obtained in step S4, to complete the genotyping of diploid individuals, specifically includes the following steps: For any two alternative contig fragments a and b, and their corresponding heterozygous regions and The following formula is used as the objective function for fractal analysis: In the formula It is a haplotype of fragment a; It is the haplotype of fragment b. and The value can be 1 or -1, representing the parent or mother line respectively; Let a and b be the Hi-C connection numbers; for and Hi-C connection count; for The number of Hi-C connections between b and . For a and Hi-C connection count; For all segments to be segmented, the following formula is used as the objective function: In the formula For all segments to be classified, the haplotype is denoted as ; Finally, the objective function f( is obtained using a random algorithm.) The maximum value of ) is the haplotype of each corresponding segment, which is its final haplotype, and the segments from the same parent are obtained; The number of Hi-C connections between fragments originating from the same parent is greater than a first preset value, while the number of connections between fragments from different parents is less than a first preset value. Fragments belonging to the heterozygous region originate from either the father or mother parent. Based on the number of Hi-C connections, fragments belonging to the heterozygous region are classified into two categories. The number of Hi-C connections between fragments of the same category is greater than a second preset value, while the number of Hi-C connections between fragments of different categories is less than a second preset value. Fragments in the homozygous region do not carry parent information and are not classified.

6. The diploid typing method based on Hi-C data according to claim 5, characterized in that... The haplotype mentioned above is specifically obtained by using the following steps: (1) Randomly set the haplotype of each segment; (2) Randomly select a segment to be classified and change its haplotype: If the objective function If the value can be increased, the change is retained; otherwise, it is not retained. (3) Repeat step (2) until the objective function is achieved. The value reaches a local maximum, and the following is determined: If the local maximum value is greater than the current global maximum value, then the current global maximum value is set as the local maximum value. (4) Randomly change the haplotype of a set of segments to be classified and the segments connected to them, and use step (2) to obtain a new local maximum value; (5) Repeat steps (2) to (4) until the set number of times to obtain the optimal local maximum value; (6) Make the objective function f( When the global maximum value is reached, the haplotype of each segment is taken as the final haplotype of that segment, thus obtaining the final haplotype of the segment.

7. A system for implementing the diploid typing method based on Hi-C data as described in any one of claims 1 to 6, characterized in that... It includes a data acquisition module, a region acquisition module, a reading filtering module, a switching identification module, and a genotyping module; the data acquisition module, region acquisition module, reading filtering module, switching identification module, and genotyping module are connected in series; the data acquisition module is used to acquire the assembly results of existing genotypes and upload the data information to the region acquisition module; the region acquisition module is used to map the acquired assembly results according to the received data information, and obtain the homozygous and heterozygous regions of the diploid genome according to the mapping relationship, and upload the data information to the reading filtering module; The reading filtering module is used to filter the Hi-C reading alignment results based on the received data information and the heterozygous variations in the diploid genome, retain the Hi-C alignment results containing heterozygous variation sites, and upload the data information to the switching recognition module. The switching identification module is used to identify large-scale switching in the assembly result based on the received data information and the filtered Hi-C comparison results, and then upload the data information to the typing module. The typing module is used to count the number of Hi-C connections between the segments after cutting based on the received data information, thereby completing the typing of diploids.

8. An assembly method comprising the diploid typing method based on Hi-C data as described in any one of claims 1 to 6, characterized in that... It also includes the following steps: S6. Based on the typing results obtained in step S5, the fragments from the same parent are clustered and connected into an assembly result to obtain the final diploid assembly result.

9. The assembly method according to claim 8, characterized in that... Step S6, based on the typing results obtained in step S5, involves clustering fragments from the same parent and connecting them into a single assembly result to obtain the final diploid assembly result. This specifically includes the following steps: Based on the genotyping results of each fragment obtained in step S5, fragments with the same haplotype are fragments from the same parent; fragments with the same haplotype are grouped into one category, resulting in two categories of fragments, from the father and mother respectively; For each class of fragments, heterozygous and homozygous fragments are sequentially connected according to the position of the fragment on the primary contig or the mapping position of the fragment on the primary contig, to obtain a new contig, which serves as the final diploid assembly result.

Images