A method of typing assembly and identifying variations of a homologous polyploid genome

Abstract

This invention discloses a method for genotyping, assembling, and identifying variants in a homologous polyploid genome, relating to the field of genome sequencing technology. The genotyping and assembly method provided by this invention is based on second-generation sequencing data, Hi-C library data, and third-generation PacBio HIFI sequencing data. First, the optimal preliminary assembly result is selected through initial assembly screening. Then, the Hi-C library data and the optimal preliminary assembly result are mounted in Hi-C. Finally, collapsed regions in the mounted genome are filled to obtain the final version of the homologous polyploid genome. The variant identification method involves isolating several haplotype genomes from the homologous polyploid genome, comparing them pairwise, and identifying the variants between haplotype genomes. Using the method of this invention, the genome integrity is higher, and it is applicable to homologous polyploids with higher ploidy levels; it employs a novel redundancy and collapse assessment method, without relying on a reference genome or manual sequence addition; and it achieves faster and more economical genotyping, assembling, and variant identification.

Description

Technical Field

[0001] This invention relates to the field of genome sequencing technology, and in particular to the typing, assembly, and identification of variants in homologous polyploid genomes. Background Technology

[0002] Autopolyploids are polyploids that have additional chromosome sets originating from the same species, typically produced by the direct doubling of chromosomes in a diploid organism. Polyploids formed from chromosome doubling within the same species are called autopolyploids. Autopolyploids are mostly found in the plant kingdom, such as wheat, corn, and rice, which play an important role in agricultural production; they are relatively rare in the animal kingdom, appearing more often in certain fish or amphibians.

[0003] Genome assembly refers to the process of generating sequence fragments (reads) from the genome of a species using sequencing methods, and then splicing these fragments based on overlapping regions. First, longer contigs are assembled, and then these contigs are further spliced ​​into longer scaffolds that allow for gaps. By eliminating errors and gaps in the scaffolds, they are located on chromosomes, resulting in a high-quality whole-genome sequence. Traditional genome assembly, whether assembling a chimeric haplotype genome through redundancy removal or creating several sets, is insufficient to represent the entire complete homologous polyploid genome. Using traditional genome assembly techniques for homologous polyploid genome typing and assembly is prone to errors and deletions, failing to effectively preserve the genetic information in the homologous polyploid genome, especially in highly heterozygous regions. Summary of the Invention

[0004] To address the problems arising from incomplete assembly during homologous polyploid assembly using traditional genome assembly techniques, this invention provides a method for genotyping and identifying variations in homologous polyploid genomes. Specifically, this is achieved through the following techniques.

[0005] A method for typing and assembling a homologous polyploid genome includes the following steps:

[0006] Sequencing was performed on the test samples to obtain second-generation library data, Hi-C library data, and third-generation PacBio HIFI sequencing data; the genome size was evaluated and estimated.

[0007] The third-generation PacBio HIFI sequencing data were used for preliminary assembly, and the degree of similarity between the preliminary assembly results and the estimated genome size was evaluated to select the optimal preliminary assembly results.

[0008] Using the Hi-C library data and the optimal preliminary assembly results, Hi-C mounting was performed to obtain a preliminary Hi-C mounted genome.

[0009] By filling in the collapsed regions in the preliminary Hi-C mounted genome, the final version of the homologous polyploid genome is obtained.

[0010] Furthermore, based on the second-generation sequencing data, a survey analysis (KMER analysis) was performed using the software GenomeScope 2.0 (https: / / github.com / tbenavi1 / genomescope2.0) to obtain a genome size estimation result of 2n = NX (N≥4), with the ploidy parameter -p ploidy provided during the estimation; where n is the number of chromosomes in the gametes, X is the number of chromosomes of haplotypes in the autopolyploids, and N is the ploidy of the autopolyploids.

[0011] Furthermore, the method for preliminary assembly using the aforementioned third-generation PacBio HIFI sequencing data is as follows:

[0012] The third-generation PacBio HIFI sequencing data were initially assembled using four parameter modes of the assembly software hifiasm (https: / / github.com / chhylp123 / hifiasm).

[0013] The four parameter modes (two hifiasm non-fractal modes and two hifiasm fractal modes) are as follows:

[0014] (1) Hifiasm non-fractal mode (-l 0-n 3), that is, filtering units composed of ≤3 reads during the assembly process;

[0015] (2) Hifiasm non-fractal mode (-l 0-n 4), that is, filtering units composed of ≤4 reads during the assembly process;

[0016] (3) Hifiasm fracturing mode (-l 2-n 3), that is, filtering units composed of ≤3 reads during the assembly process;

[0017] (4) Hifiasm fracturing mode (-l 2-n 4), which means filtering units composed of ≤4 reads during the assembly process.

[0018] When using two hifiasm non-genotyping modes, the assembled genomes p_ctg and a_ctg are merged as the first preliminary assembly result; when using two hifiasm genotyping modes, the assembled genomes hap1_ctg and hap2_ctg are merged as the second preliminary assembly result.

[0019] Compare the continuity and completeness of the first preliminary assembly result and the second preliminary assembly result, and select the one that is closer to 1.2 times the estimated genome size as the optimal preliminary assembly result, that is, assemble a complete contig-level homologous polyploid genome.

[0020] The above assembly results do not remove redundant sequences so that the initial assembly version is large and complete enough to have enough contig sequences for subsequent Hi-C mounting.

[0021] It should be further noted that in the preliminary assembly method described above, phenotyping mode-l2 and non-phenotyping mode-l0 were used to obtain hap1.p_ctg / hap2.p_ctg and p_ctg / a_ctg, respectively. The assembly software hifiasm can only separate a maximum of two genome sets for polyploids. For example, even if two sets can be separated in autotetraploids, it can only guarantee that hap1 and hap2 each contain two haplotypes of different sizes.

[0022] The subsequent genotyping and assembly operations of this invention can ensure the complete assembly of four haplotypes. For the genotyping assembly method, hap1.p_ctg and hap2.p_ctg will be obtained. For the non-genotyping assembly method, two parts of contigs will be obtained: p_ctg (primary contigs) and a_ctg (composed of all contigs discarded during the generation of primary contigs). Although these two methods cannot accurately separate each haplotype, merging them together forms a complete homologous polyploid genome, sufficient for subsequent genotyping operations.

[0023] This method further employs the -n 3 / -n 4 mode for assembly in both operating modes, removing units at the ends of the assembly connection diagram that are not supported by multiple reads to varying degrees, thereby obtaining contig-level homologous polyploid genomes in four assembly modes.

[0024] From these four preliminary genome assembly versions, after comparing their continuity (N50) and integrity (BUSCO), we selected the version that was closest to the estimated complete homologous polyploid genome size of 1.2 times for subsequent mounting and typing operations.

[0025] Furthermore, the specific method for mounting Hi-C is as follows:

[0026] Hi-C library data was aligned to the optimal preliminary assembly result using Hi-C data alignment software (e.g., Bwa). Then, Hi-C signal clustering software (using polyploid mode - m diploid) was used to cluster and mount the preliminary assembled contig-level homologous polyploid genome (i.e., the optimal preliminary assembly result) according to the alignment results. Mapping software (e.g., Juciebox) was used to adjust homologous chromosomes together based on their interaction signals and determine the position and orientation of contigs within chromosomes (i.e., manual sequencing and orientation). During mapping, more than N interchromosomal interaction signals were identified as redundant and removed (i.e., redundancy removal) to obtain the Hi-C preliminary mounted genome.

[0027] Generally, Hi-C signal clustering software has different software parameters for haploid and polyploid, namely -m haploid and -m diploid respectively. The parameter used in polyploid mode is -m diploid.

[0028] Optionally, the Hi-C data comparison software (e.g., software Bwa) selects the -SP5M mode, uses a segment alignment method to take the alignment of the smallest query coordinate as the primary alignment, and marks shorter segment hits as secondary hits.

[0029] Optionally, the Hi-C signal clustering software (e.g., software 3d-dna) adopts the `--early-exit--mapq0-r 2-i 15000-m diploid` mode, exiting early after the first round of mounting with the lowest quality threshold, obtaining only clustering and preliminary sorting mounting results, i.e., the Hi-C preliminary genome mounting. This operating mode can significantly shorten the software running time, preserve as many interaction signals between homologous chromosomes as possible, and cluster homologous chromosomes together, facilitating subsequent manual fine-tuning using mapping software.

[0030] The following example uses autotetraploidy:

[0031] If there is one homotetraploid (with 2n = 4X = 48 chromosomes), this method will cluster the contig-level genome into 12 clusters. Each cluster actually contains 4 homologous chromosomes, but they are clustered together because the homologous chromosomes are too similar.

[0032] Unlike conventional Hi-C mapping, and noteworthyly, autopolyploid typing requires attention not only to intrachromosomal interactions but also to interchromosomal interactions. This method utilizes as many interchromosomal interaction signals as possible to guide autopolyploid typing.

[0033] It is worth mentioning that, in order to ensure the integrity of the genome, our contig-level genome (i.e., the optimal preliminary assembly result) is close to 1.2 times the size of the real genome, so there must be redundant sequences in it.

[0034] Due to the complexity of the genome, including repetitive and heterozygous regions, and the limitations of genome assembly software, some regions may assemble fewer than four sets, while others may assemble more than four sets; both are normal phenomena. These regions can be observed in mapping software through interchromosomal interaction signals. Regions with more than four sets appearing in the genome are identified as redundant sequences and removed during Hi-C mounting; regions with fewer than four sets are identified as collapsed regions and require subsequent filling in of the collapsed regions.

[0035] It should be noted that the aforementioned regions refer to any region on the genome. For example, a carotetraploid can be considered as a mixture of four extremely similar haplotypes. Any segment of the genome sequence within a carotetraploid should contain three other similar sequence regions. Due to the complexity of the genome, including repetitive and heterozygous regions, and limitations of genome assembly software, some regions may only assemble two haplotypes, while others may assemble eight. Therefore, redundancy removal and gap filling operations are necessary for the assembly results.

[0036] Furthermore, the method for filling the collapsed regions in the preliminary Hi-C mounted genome is as follows:

[0037] In a group of non-redundant chromosomes that are homologous to each other, the longest chromosome is selected as the ref, and the other homologous chromosomes are selected as queries. (Using gene alignment software, all queries are aligned to the ref.)

[0038] Based on the comparison results, collinearity blocks in the comparison are extracted, and regions that exist in the ref but not in the query are filled into the collapsed regions of the query.

[0039] Specifically, the ref and query are divided into several collinear blocks. The gap lengths between adjacent collinear blocks on the ref and query are compared. Based on the comparison results, the gap regions are divided into haplotype-specific regions, haplotype-collapsed regions, and non-obvious collapse regions.

[0040] For the haplotype collapsed region, the gap region between adjacent collinear blocks on the ref is extended 20kb to both ends as a candidate sequence to be filled into the query; the end of the candidate sequence is aligned to the end 20kb region of the corresponding collapsed region on each query, and the alignment method with the highest alignment quality value and the fewest skip alignments is selected to automatically fill the collapsed region of the query with the sequence on the ref.

[0041] No treatment is applied to haplotype-specific regions and regions with indistinct collapse.

[0042] The final version of the homologous polyploid genome, with the collapse filled, was obtained.

[0043] Furthermore, gene alignment software (such as Minimap2) uses the -ax asm5 mode to align other chromosomes (queries) with collapsed regions to the ref, resulting in an alignment sam file.

[0044] The alignment information is extracted from the SAM file, and secondary alignments and overlapping alignment intervals are removed from the alignment intervals. The ref and query are then divided into several non-redundant and one-to-one collinear blocks (collinear regions) according to the alignment intervals. Therefore, the area between two adjacent collinear blocks is the area where the query is not aligned to the ref.

[0045] If the unaligned blank region between two collinear blocks in the query is larger than the unaligned blank region between two collinear blocks in the reference, then the blank region in the query is an additional region assembled compared to the reference, belonging to a specific region relative to the reference, and is retained in the genome. Conversely, if the unaligned blank region between two collinear blocks in the query is smaller than the unaligned blank region between two collinear blocks in the reference, then it is determined to be a collapsed region that needs to be filled.

[0046] Specifically, a threshold of 5 can be set. If the length of the blank region on the reference exceeds 5 times the threshold of the blank region on the query, the blank region on the reference is extended 20kb in both directions as a candidate sequence to fill the collapsed region. Then, it is aligned to the last 20kb region of the corresponding collapsed region on the query. The alignment information with the highest alignment quality value and the fewest skip alignments is selected, and the sequence on the reference is used to fill the collapsed region of the query. If the length of the blank region on the reference exceeds the length of the blank region on the query but is less than 5 times the threshold, it is judged as an insignificant collapsed region, which may be due to differences caused by the evolution of homologous chromosomes over a long period of time, and is not processed.

[0047] The method for filling the collapsed region in the preliminary Hi-C mounted genome provided by this invention also compares the homologous chromosomes before and after the collapse filling, obtains the alignment PAF file using the -cx asm5 mode of minimap2 software, and draws a collinearity map based on the alignment position to verify the filling effect.

[0048] If the filling effect is unsatisfactory, the threshold parameters in the aforementioned collapse filling software need to be adjusted to achieve the best filling effect. This can fill the collapse to compensate for the defects of incomplete assembly without introducing too many specific sequences from the reference, thus avoiding redundancy. Finally, the Hi-C data is aligned to the filled homologous polyploid genome, and Hi-C heatmap software (such as HiCExplorer) is used to plot the Hi-C interaction signals across the entire genome, further visualizing the reliability of the homologous polyploid genome typing and assembly.

[0049] Furthermore, the final version of the homologous polyploid genome is also assessed for continuity, consistency, integrity, and redundancy (i.e., quality assessment).

[0050] Furthermore, genome continuity is evaluated by assessing the N50 size of the genome;

[0051] Second-generation data alignment software (such as Bwa software) is used to align the clean reads after quality control of second-generation high-throughput sequencing to the final version of the homologous polyploid genome. The consistency of the assembled genome is evaluated by statistically analyzing the alignment rate and coverage of the second-generation library data.

[0052] Using third-generation data alignment software (such as minimap2), clean reads after quality control of third-generation high-throughput sequencing were aligned to the final version of the homologous polyploid genome. The alignment rate and coverage of the third-generation PacBio HIFI sequencing data were statistically analyzed to evaluate the consistency of the assembled genome.

[0053] Wherein, the alignment rate represents the proportion of clean reads that have been aligned to the final version of the homologous polyploid genome out of the total number of clean reads; the coverage represents the proportion of the length of the reference genome fragment that has been aligned to the clean reads out of the total length of the final version of the homologous polyploid genome;

[0054] The completeness of genome assembly is determined by assessing the percentage of Complete BUSCOs using integrity assessment software (such as BUSCO software). This software evaluates the completeness and accuracy of genome assembly or prediction by comparing the predicted set of genes from an input genome or gene with a known set of single-copy genes.

[0055] The third-generation sequencing data (third-generation PacBio HIFI sequencing data) were then aligned to the final version of the homologous polyploid genome (non-primary alignments were filtered out so that each read was aligned to at most one location). The coverage depth of each site was calculated, the frequency of each depth was counted, and it was determined whether it conformed to a Poisson distribution, i.e., redundancy was determined.

[0056] When the genome assembly results contain redundant or collapsed regions, the coverage depth of the redundant parts is only half the normal depth. Therefore, if there is a peak before the main peak (determined by the estimated genome size), it indicates the presence of redundancy (another possibility is the presence of highly heterozygous sequences); if there is a peak after the main peak, it indicates that repetitive sequences have been excessively merged, i.e., collapsed regions exist. When the genome assembly results have no redundancy (i.e., good genotyping), they conform to a Poisson distribution.

[0057] Optionally, the second-generation data alignment software adopts the mem alignment mode to align the second-generation sequencing data onto the genome, obtains the alignment bam file, and uses the samtools (https: / / github.com / samtools / samtools) depth software to obtain the depth of the second-generation reads at each site on the genome, and calculates the alignment rate of the second-generation sequencing data and the coverage of the second-generation data on the genome; the consistency evaluation of third-generation data is similar.

[0058] When assessing redundancy, the third-generation PacBio HIFI sequencing data were aligned to each haplotype (i.e., each haplotype genome in the final chromosome-level homologous polyploid genome), all non-primary alignments were filtered out, ensuring that each read was aligned to at most one site, and the coverage depth of all sites was evaluated to ensure that it conformed to a Poisson distribution.

[0059] Generally, if both the alignment rate and coverage are greater than 90%, then the consistency between the genome and sequencing data is considered to be without problems.

[0060] The present invention also provides a method for identifying variations in a homologous polyploid genome. First, a homologous polyploid genome is obtained using the above-mentioned typing and assembly method. Then, several haplotype genomes are isolated from the homologous polyploid genome. Next, the haplotype genomes are compared pairwise to finally identify the variations between the haplotype genomes.

[0061] Furthermore, the isolated haplotype genomes are paired and compared using software minimap2 or mummer. Based on the comparison results, structural rearrangements in collinear and non-collinear regions are summarized. Structural and functional annotations are performed on the variations between haplotypes.

[0062] Taking autotetraploids as an example:

[0063] By isolating each haplotype from the whole, four haplotype genomes (A1 / A2 / A3 / A4) are obtained. Pairwise alignment of these four haplotype genomes reveals collinearity and variation among the haplotypes.

[0064] For example, alignments are performed using A2 to A1, A3 to A2, and A4 to A3 to find collinear blocks and identify variations between haplotypes. Each haplotype is aligned using the alignment software minimap2 in -ax asm5--eqx mode, or the alignment software mummer in --maxmatch-c 500-b 500-l 50-g 90 mode. Based on the alignment information, collinearity, structural rearrangements (inversions, translocations, and duplications), local variations within structural rearrangements (SNPs, INDELs, CNVs, etc.), and unaligned regions are identified. The software first determines the longest collinear path. Since all non-collinear regions correspond to genomic regions rearranged between two genomes, the identification of collinearity also identifies all structural rearrangements.

[0065] Subsequently, all collinear regions in the alignments were classified as inversions, translocations, or duplications based on the conformations constituting the alignments. This method also identifies local variations within all collinear regions and structural rearrangements. Local variations include short variations and small insertions / deletions, as well as structural variations (such as large insertions / deletions, copy number variations, and highly differentiated regions). In summary, short variations are resolved from the constitutive alignments, while structural variations are inferred by comparing the overlap and gaps between consecutive alignments of collinear or rearranged regions.

[0066] After identifying all collinear blocks and variation information between haplotypes, all collinear regions and structural variation regions are displayed on a collinearity map of haplotype comparisons to visually present the differences between haplotypes. Furthermore, combined with subsequent genome functional annotation, genes affected by structural variations can be functionally enriched to identify significantly enriched functional categories, pathways, or processes, thereby revealing the biological significance of the variation.

[0067] Compared with the prior art, the advantages of the present invention are:

[0068] 1. The method for genotyping and identifying variants of homopolyploid genomes proposed in this invention can obtain complete homopolyploid genome sequences at the contig level using PacBioHIFI sequencing technology and various high-quality assembly methods, overcoming the incomplete situation in the previous preliminary genome assembly, which could only assemble one or two haplotypes.

[0069] 2. The method proposed in this invention has strong adaptability and can be applied to the typing and assembly of tetraploid, hexaploid, octoploid, and even higher ploidy homologous polyploids, providing more reliable and complete basic data for subsequent downstream genome analysis.

[0070] 3. The method proposed in this invention significantly shortens the software running time during Hi-C mounting and obtains as many interchromosomal interaction signals as possible with less sequencing cost to guide the mounting and typing of homologous polyploids.

[0071] 4. The method proposed in this invention automatically searches for collinear blocks during genome collapse filling and uses self-assembled homologous chromosomes to fill in incomplete regions. This overcomes the difficulties of relying on reference genomes and manual sequence addition during previous genome collapse filling processes, achieving a faster and more economical typing and assembly process.

[0072] 5. The method proposed in this invention adopts a novel assessment method for haplotype redundancy and collapse when evaluating the quality of genotyping. It uses the original PacBio HIFI sequencing data to align to haplotypes and plots depth and frequency. The peak shape is used to determine whether the haplotype contains redundancy or collapse.

[0073] 6. The method proposed in this invention uses haplotype comparison to analyze haplotypes when identifying variations. It further explores the differences between haplotypes within homologous polyploid genomes by finding collinear regions and simultaneously identifying structural rearrangements between haplotypes. Attached Figure Description

[0074] Figure 1 This is a flowchart illustrating the overall process for typing, assembling, and identifying variations in a homologous polyploid genome in a specific implementation.

[0075] Figure 2 A flowchart detailing the initial assembly of third-generation genome sequencing data (third-generation PacBio HIFI sequencing data);

[0076] Figure 3 Detailed flowchart for mounting Hi-C;

[0077] Figure 4 A flowchart detailing the process of filling and identifying variations in collapsed areas;

[0078] Figure 5 and Figure 6 These are collinearity maps between homologous chromosomes before and after the collapse filling process;

[0079] Figure 7 and Figure 8 These are the full and partial plots of the Hi-C interaction signal, respectively.

[0080] Figure 9 Depth-frequency plot for genome assembly results conforming to a Poisson distribution;

[0081] Figure 10 This is a visualization of haplotype collinearity and variation. Detailed Implementation

[0082] The technical solution of the present invention will be clearly and completely described below. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0083] The method for typing, assembling, and identifying variants in a homopolyploid genome provided by this invention comprises three parts: obtaining sequencing data of the homopolyploid genome, typing and assembling the genome, and identifying variants. The sequencing data of the homopolyploid genome includes three parts: second-generation sequencing data (second-generation library), Hi-C library data (Hi-C sequencing), and third-generation PacBio HIFI sequencing data.

[0084] The overall flowchart of the method for typing, assembly, and identification of variants provided by this invention is as follows: Figure 1 As shown.

[0085] In some embodiments of the present invention, the part-type assembly process is mainly as shown in steps one through six below. In general, it includes the following steps:

[0086] Step 1: Sequencing

[0087] Second-generation whole-genome sequencing, Hi-C sequencing, and PacBio HIFI sequencing were performed using the same biological sample.

[0088] Step 2: Genome Size Estimation

[0089] like Figure 4 As shown, a survey analysis was performed using second-generation whole-genome sequencing (second-generation sequencing data), and the size of the complete homologous polyploid genome was assessed using the software GenomeScope 2.0 (https: / / github.com / tbenavi1 / genomescope2.0).

[0090] Step 3: Contig-level complete genome assembly

[0091] like Figure 2 As shown, hifiasm (https: / / github.com / chhylp123 / hifiasm) was used to perform preliminary assembly based on third-generation PacBio HIFI sequencing data, utilizing four modes. The assembled contigs were merged with the contigs discarded during the assembly process to complete the preliminary assembly of the homologous polyploids, yielding the optimal preliminary assembly result.

[0092] Step 4: Hi-C Mount

[0093] 1. For example Figure 3 As shown, the Hi-C data alignment software Bwa (https: / / github.com / lh3 / bwa) was used to align the raw Hi-C data to the complete contig-level homologous polyploid genome; 3d-dna (https: / / github.com / aidenlab / 3d-dna) or ALLHiC (https: / / github.com / tangerzhang / ALLHiC) was used to cluster the contigs and preliminarily anchor them.

[0094] Optionally, special MapQ quality values ​​and operating modes can be used to save sequencing costs and software running time.

[0095] 2. For example Figure 3 As shown, the Hi-C interaction signals were visualized using the mapping software Juciebox (https: / / github.com / aidenlab / juicebox) and further refined. Based on the interaction signals between homologous chromosomes, homologous chromosomes were aligned together, and the position and orientation of contigs within chromosomes were determined (i.e., manual sequencing and orientation). During mapping, interaction signals between chromosomes exceeding N were identified as redundant and removed (i.e., redundancy removal).

[0096] The initial Hi-C-mounted genome was finally obtained.

[0097] Step 5: Fill the collapsed area

[0098] 1. For example Figure 4 As shown, the longest homologous chromosome in the initial Hi-C-mounted genome is taken as the reference; other homologous chromosomes are used as queries. Each query is aligned to the reference, and collinear regions are extracted based on the alignment information.

[0099] 2. For example Figure 4 As shown, the ref and query are divided into several collinear blocks. Based on the comparison of the gap length between adjacent collinear blocks on the ref and query, the gap region is divided into: haplotype-specific region, haplotype-collapsed region, and insignificant collapse region.

[0100] (1) For haplotype collapse areas

[0101] Extend the gap region between adjacent collinear blocks on the ref by 20kb at both ends as a candidate sequence to be added to the query.

[0102] Align the end of the candidate sequence to the end 20kb region of the corresponding collapsed region on the query; select the alignment method with the highest alignment quality value and the fewest skip alignments, and automatically fill the sequence on the ref into the collapsed region of the query.

[0103] (2) No processing is done for haplotype-specific regions and regions with indistinct collapses; these are considered as true differences between haplotypes.

[0104] 3. For example Figure 4 As shown, after the collapse filling was completed, the original Hi-C data was re-aligned to the new genome (i.e., the final version of the homopolyploid genome with the collapse filled). The software HiCExplorer (https: / / github.com / deeptools / HiCExplorer) was used to display the interaction signals within and between chromosomes on the new genome, further verifying the effect of homopolyploid assembly typing.

[0105] Step Six: Quality Assessment

[0106] 1. Continuity

[0107] like Figure 1 , 4 As shown, the N50 size of the genome is generally used to evaluate the continuity of the genome; the larger the N50, the better the continuity of the genome.

[0108] 2. Consistency

[0109] like Figure 1 , 4 As shown, the software Bwa (https: / / github.com / lh3 / bwa) was used to process the second-generation library data ( Figure 4 The clean reads from the second-generation clean reads were compared back to the final chromosome-level homologous polyploid genome to calculate the alignment rate and coverage.

[0110] The alignment rate and coverage were calculated by back-aligning third-generation PacBio HIFI sequencing data (CCS data) to the homologous polyploid genome using minimap2 (https: / / github.com / lh3 / minimap2).

[0111] Generally, an alignment rate greater than 90% and a coverage greater than 90% indicate that the genome has a high degree of consistency with the original data.

[0112] 3. Completeness

[0113] like Figure 1 , 4As shown, the final chromosome-level homologous polyploid genome was compared with the conserved single-copy gene set protein library (i.e., the single-copy gene set included in the BUSCO software) using the software BUSCO (version: v5.4.6; parameter: –evalue 1e-05) to assess the integrity of the genome.

[0114] Generally, if more than 90% of single-copy genes are predicted, it indicates that the genome assembly is relatively complete.

[0115] 4. Redundancy.

[0116] like Figure 1 , 4 As shown, the software minimap2 (https: / / github.com / lh3 / minimap2) was used to re-align the third-generation PacBio HIFI sequencing data to each haplotype of the homologous polyploid genome (i.e., each haplotype genome in the final chromosome-level homologous polyploid genome); the depth of the original data coverage at each locus was calculated, the frequency of occurrence at each depth was counted, and a depth-frequency map was plotted. If the depth-frequency map conforms to a Poisson distribution, it indicates that each haplotype has no redundancy and no collapsed regions.

[0117] like Figure 1 As shown, in other embodiments of the present invention, a method for haplotype comparison and identification of variations is also provided based on the above-mentioned homologous polyploid genome typing and assembly results. The specific process is as follows:

[0118] Use minimap2 (https: / / github.com / lh3 / minimap2) or mummer (https: / / github.com / mummer4 / mummer) to perform pairwise alignments for each isolated haplotype.

[0119] By comparing information, we can summarize the structural rearrangements in collinear and non-collinear regions, and visualize the structural rearrangements in collinear and non-collinear regions between each group of homologous chromosomes.

[0120] We can also combine structural annotation to annotate variations between haplotypes and extract genes affected by these variations; we can also combine functional annotation to study which functional categories, pathways, or processes these affected genes may be enriched in. This step allows for a direct comparison of differences between haplotypes and extends this understanding to biological questions, contributing to a deeper understanding of the evolution and genetic mechanisms of plant and animal genomes.

[0121] Example

[0122] I. Sample Testing

[0123] 1. Check whether the sample (derived from fresh young blueberry leaves) contains foreign objects;

[0124] 2. High-quality genomic DNA was extracted using the modified CTAB / SDS method (Doyle and Doyle 1987) (Guobin et al. 2008);

[0125] 3. Genomic DNA was subjected to 0.75% agarose gel electrophoresis to detect the size of genomic DNA fragments and the degree of DNA degradation;

[0126] 4. Use a NanoDrop One spectrophotometer (Thermo Fisher Scientific) to determine DNA purity;

[0127] 5. Use a Qubit 3.0 fluorescence spectrometer (Life Technologies, Carlsbad, CA, USA) to detect the concentration of genomic DNA and perform precise quantification of DNA.

[0128] II. Sequencing

[0129] 1. Next-generation sequencing (e.g.) Figure 4 (As shown)

[0130] (1) Construction of the second-generation library

[0131] After the DNA samples passed the test, they were randomly broken using a Covaris ultrasonic disruptor.

[0132] Using the Nextera DNA Flex Library Prep Kit (Illumina, San Diego, CA, USA), Illumina sequencing libraries with insert sizes of 200-400 bp were constructed through steps such as end repair, addition of A tails, addition of sequencing adapters, purification, and PCR amplification.

[0133] (2) Document Quality Inspection

[0134] After the sequencing library is constructed, preliminary quantification is performed using the software Qubit 3.0.

[0135] After diluting the library, Agilent 2100 software was used for quality control. Once the insert size met expectations, the effective concentration of the library was accurately quantified using Q-PCR to ensure library quality.

[0136] (3) Sequencing

[0137] After the library passed the testing, high-throughput sequencing was performed using the Illumina NovaSeq 6000 platform (Illumina, San Diego, CA, USA) to obtain next-generation sequencing data.

[0138] (4) Statistics on the amount of second-generation sequencing data

[0139] The raw reads were filtered using the software Fastp (version: 0.21.0, parameter: --detect_adapter_for_pe), discarding low-quality reads. The filtering criteria are as follows:

[0140] Remove reads with an N base content exceeding 5%;

[0141] Remove reads with 50% or more low-quality (quality value less than or equal to 5) bases;

[0142] Remove the adapter sequence from the reads;

[0143] Remove repetitive sequences caused by PCR amplification.

[0144] The statistical results are shown in Table 1.

[0145] Table 1. Statistics on the amount of second-generation sequencing data.

[0146] Type Total_reads Total_bases GC_content Q20 Q30 rawdata 1,080,000,648 162,000,097,200 43.59％ 97.29％ 92.77％ cleandata 1,079,974,134 161,996,120,100 43.59％ 97.29％ 92.77％

[0147] 2. Hi-C library preparation and sequencing process (e.g.) Figure 3 (As shown)

[0148] Hi-C assisted assembly technology originates from chromosome conformation capture (3C) technology. Using the entire cell nucleus as the research object, it employs high-throughput sequencing technology combined with bioinformatics methods to study the spatial relationships of chromatin DNA across the entire genome. By capturing all DNA interaction patterns within chromatin, the genome draft is precisely located to chromatin, constructing complete chromatin sequence location structure information.

[0149] (1) Tissue sample testing: Fresh tissue samples are used.

[0150] (2) Hi-C Library Construction

[0151] ① Cell cross-linking: Using formaldehyde to fix the sample, cross-linking is performed between intracellular proteins and DNA, and between DNA molecules, to preserve their interaction relationships and maintain the 3D structure of the cell.

[0152] ② Restriction enzyme digestion: DNA is digested using restriction endonucleases, resulting in sticky ends on both sides of the cross-links. The commonly used restriction endonuclease is DPN II.

[0153] ③ End repair: Using the end repair mechanism, biotin-labeled bases are introduced to facilitate subsequent DNA purification and capture.

[0154] ④ Circularization: The repaired DNA is circularized, and the DNA fragments containing interactions are circularized to ensure that the location of interacting DNA can be determined during subsequent sequencing and analysis.

[0155] ⑤ DNA purification and capture: Decrosslink the DNA, purify the DNA, and break it into 300-700bp fragments; use streptavidin magnetic beads to capture DNA fragments containing interaction relationships for library construction.

[0156] ⑥ Library quality control: After library construction, the concentration and insert size of the library were tested using Qubit 3.0 and Agilent 2100 software, respectively. The effective concentration of the library was accurately quantified using Q-PCR to ensure library quality.

[0157] ⑦ Sequencing: After the library passes the inspection, high-throughput sequencing is performed using the Illumina platform with a read length of PE150 to obtain Hi-C library data.

[0158] (3) Statistics on Hi-C sequencing data volume

[0159] In addition to the second-generation data filtering mentioned above, the Clean Data was aligned to the reference genome and filtered using the software Jucier (version: 1.6). The filtering criteria are as follows:

[0160] ① Remove unmapped reads that failed to uniquely align to the reference genome at both ends;

[0161] ② Remove invalid pairs, including self circles and dangling ends;

[0162] The statistical results are shown in Table 2.

[0163] Table 2. Statistical results of Hi-C sequencing data volume.

[0164]

[0165] 3. Third-generation sequencing (e.g.) Figure 2 (As shown)

[0166] (1) Third-generation SMRT (Single Molecule Real Time) sequencing

[0167] ① DNA extraction: High-quality genomic DNA is extracted and subjected to strict quality control (Doyle and Doyle 1987) (Guobin et al. 2008). The genomic DNA is broken into appropriately sized fragments, and then damage repair, adapter ligation, fragment screening, and other operations are performed to construct a DNA library.

[0168] ② Library construction: After obtaining high-quality purified genomic DNA samples, PCR-free SMRT-Bell libraries were constructed and sequenced using the PacBio Sequel II sequencing platform.

[0169] ③ Sequencing: After fragment screening and quantification, the calculated concentration and volume of the library template and enzyme complex are transferred into the nanopores of the PacBio Sequel II sequencer for sequencing.

[0170] ④ Data Reading: The raw PacBio single-molecule sequencing data is recorded in the subreads.bam file, in which sequencing adapters and low-quality sequences have been removed; data reading uses PacBio's official software SMRTlink v11.0 to perform quality control and statistical analysis on the output data to obtain the final valid data.

[0171] (2) Statistics on the amount of third-generation sequencing data

[0172] Polymerase reads from sequencing were processed as follows to obtain subreads containing only the insert fragment. These subreads were then processed using SMRTLink v11.0 software to generate CCS reads (parameters: –min-passes = 3 –min-rq = 0.99), which are the third-generation PacBio HIFI sequencing data.

[0173] ① Remove Polymerase reads shorter than 50bp;

[0174] ②Remove Polymerase reads with a quality value below 0.8;

[0175] ③ Remove self-connected Polymerase reads;

[0176] ④ Remove the connector sequence from the Polymerase reads.

[0177] The statistical results are shown in Table 3.

[0178] Table 3. CCS Reads Statistical Results

[0179] num_seqs sum_len min_len avg_len max_len N50 3,723,167 68,225,527,142 93 18,324.6 50,622 19,347

[0180] III. Preliminary Genome Assembly

[0181] Hifiasm software (version: 0.19.5, parameters: -l 2 -n 4) was used, with only CCS Reads employed for classification and assembly. The assembly process is as follows: Figure 2 As shown.

[0182] 1. In the early stages, the software uses the OLC (Overlap Layout Consensus) algorithm for genome assembly. During the assembly process, all-versus-all alignments are first performed on the reads to identify overlaps between them, thus performing a clustering-like operation.

[0183] 2. Subsequent self-correction is performed on clustered reads (reads with overlap). The Hifiasm software performs three rounds of error correction by default.

[0184] PacBio sequencing errors are random, so reads can be used to correct each other.

[0185] 3. Perform an all-versus-all comparison on the corrected reads again to find overlaps between reads in order to construct a string graph.

[0186] To differentiate between different haplotypes, the software Hifiasm utilizes the differing sites (i.e., heterozygous sites) between haplotypes. At the read level, this means that most reads share the same bases, but some positions have different base types (at least three reads support each base type). These differing reads originate from different haplotypes, and haplotype construction can be performed using these reads.

[0187] Specifically, for homologous polyploid genomes:

[0188] (1) When the Hifiasm software is run in genotyping mode (-l 2), hap1_ctg and hap2_ctg are obtained (each containing two sets of genomes). Therefore, hap1_ctg + hap2_ctg needs to be used as the initial assembly version containing multiple sets of original haplotypes;

[0189] (2) When the software Hifiasm is run in non-fragmented mode (-l 0), it will get two parts of contigs: p_ctg (primary contigs) and a_ctg (composed of all contigs discarded when p_ctg is generated).

[0190] p_ctg stands for primary contigs, while a_ctg consists of all contigs discarded during the generation of primary contigs. Therefore, p_ctg + a_ctg is needed as the initial assembly version containing four sets of primitive haplotypes.

[0191] The preliminary assembly indicators are shown in Table 4.

[0192] Table 4 Preliminary Assembly Indicators Statistics

[0193] Item Value Total_length(bp) 2,758,649,574 total_length_withoutN(bp) 2,758,649,574 Total_number 1,685 GC_content 38.692％ N50(bp) 9,586,632 N90(bp) 2,704,003 Average (bp) 1,637,180.76 Median(bp) 54,848 Min(bp) 10,851 Max(bp) 29,666,133

[0194] IV. Hi-C Mount

[0195] Hi-C mounting process as follows: Figure 3 As shown. Specifically:

[0196] 1. Hi-C library data alignment at the contig level genome level

[0197] In the Hi-C library, DNA fragments containing interactions were sequenced using second-generation paired-end sequencing to obtain Reads1 and Reads2 (clean data). Reads1 and Reads2 were then aligned one-end to the contig sequence (i.e., the optimal preliminary assembly result) using Bwa alignment. Both reads uniquely aligned to the DNA fragments containing the reads in the genome draft, indicating the interaction between the two reads; these two read pairs are called UniqueMapped Paired-end Reads.

[0198] Because Hi-C reads are not all derived from the re-ligation of two spatially adjacent but linearly non-adjacent restriction enzyme fragments, the Hi-C data generated by sequencing has a high false positive rate. Therefore, it is necessary to separate valid Hi-C data (Valid Interaction Pairs) and invalid Hi-C data (Invalid Interaction Pairs) from UniqueMapped Paired-end Reads.

[0199] Valid Hi-C data (Valid Interaction Pairs) refers to sequenced reads derived from two spatially adjacent but linearly non-adjacent digested DNA fragments. It provides effective interaction information for subsequent genome assembly. Invalid Hi-C data (Invalid Interaction Pairs) does not provide inter-fragment interaction information and includes three types: self-circles, dangling ends, and dumped pairs.

[0200] When both ends of paired-end sequencing reads in a Hi-C library fall between adjacent restriction endonucleases, they can be categorized into three types based on the orientation of the paired reads aligned to the genome: Self Circle, Dangling End, and Dumped Pairs. These are all considered invalid interaction pairs. When both ends of paired-end sequencing reads in a Hi-C library fall between non-adjacent restriction endonucleases, they are identified as valid interaction pairs.

[0201] The final data retains only the valid interactions, i.e., the valid Hi-C data.

[0202] 2. Cluster mounting

[0203] Cluster mounting can be performed using the following two methods:

[0204] (1) Using software 3D-DNA clustering mounting

[0205] Using the software 3D-DNA (version: 180419, parameters: –early-exit –mapq 0-r 2-i 15000-mdiploid) and Jucier (version: 1.6), the pairwise interactions between Contigs (i.e., the pairwise interactions in the optimal preliminary assembly results) are converted into a specified binary file (i.e., a .hic file); then clustering is performed to generate an assembly file (.assembly).

[0206] (2) Using the ALLHiC software for clustering and mounting

[0207] First, using the software ALLHiC (version: 0.9.8, parameters: default parameters), the optimal preliminary assembly result of the contig sequences is initially clustered into X chromosome groups through agglomerative hierarchical clustering (bottom-up hierarchical clustering algorithm). Then, the contigs within the chromosome groups are ordered and oriented. Next, using the software 3D-DNA (version: 180419), the interaction relationships between each pair of contigs are converted into a specified binary file (.hic file), and the clustering generates an assembly file (.assembly).

[0208] 3. Manually sequence and phase the circuit to remove redundancy.

[0209] Manually reorder and reorient (.review.assembly) the pre-ordered and reoriented contigs (.assembly) using the Windows software Juciebox (version: 1.11.08). The following rules must be followed:

[0210] (1) If a sequence does not interact with a sequence of the same size at all, but interacts normally with other sequences, it is judged as a heterozygous sequence and needs to be removed manually.

[0211] If there are more than N sets of large sequence fragments (where N is the ploidy of the homologous polyploid) with obvious homologous signals, then the large sequence fragments beyond N sets are redundant sequences assembled from multiple sets and need to be removed manually.

[0212] (2) After sequencing and orientation, the redundant Contig sequences are removed and connected with 100 Ns to obtain the preliminary Hi-C mounted genome, which is the chromosome-level homologous polyploid genome before subsequent collapse filling.

[0213] V. Filling the collapsed area

[0214] After removing redundant regions and large, extra-assembled sequences from the chromosome-level genome, multiple haplotype genomes (Haplotype A1, Haplotype A2, Haplotype A3, Haplotype A4, etc.) are obtained. However, each haplotype will contain a collapsed region on a certain chromosome, which will be slightly shorter than other homologous chromosomes. This collapse phenomenon is normal. Since the assembly software can generate multiple sets of results for certain regions, it can naturally also generate fewer than N sets for some regions.

[0215] To fill this gap, we performed preliminary Hi-C-mounted genome mapping according to... Figure 4The process shown involves treating the collapsed area, as detailed below:

[0216] 1. Select the longest of the N homologous chromosomes of each chromosome as the ref, and you will get the longest and most complete haplotype genome, denoted as HaplotypeA1;

[0217] Other homologous chromosomes (Haplotype A2, Haplotype A3, and Haplotype A4, etc.) were used as queries. The queries were then mapped to the longest chromosome, Haplotype A1 (ref), using the software minimap2 (version: 2.24-r1122, parameters: -ax asm5–eqx).

[0218] 2. Extract collinear regions based on the comparison results.

[0219] The ref and query are divided into several collinear blocks. The gap lengths between adjacent collinear blocks on the ref and query are compared. Based on the comparison results, the gap regions are divided into haplotype-specific regions, haplotype-collapsed regions, and insignificant collapse regions.

[0220] 3. For haplotype-collinear regions, the gap between adjacent collinear blocks on the reference is extended 20kb in both directions to serve as candidate sequences to be filled into the query. The end of the candidate sequence is aligned to the end 20kb region of the corresponding collapsed region on the query. The alignment method with the highest alignment quality value and the fewest skip alignments is selected, and the sequence on the reference is automatically filled into the collapsed region of the query. Other regions (haplotype-specific regions and regions with insignificant collapse) are not processed.

[0221] 4. For chromosomes in the homologous chromosome sets before and after the collapse filling, pairwise alignments were performed using the software minimap2 (version: 2.24-r1122, parameter: -cx asm5). This involved comparing the initially Hi-C-mounted genome with the chromosomes in the homologous chromosome sets after collapse filling. Based on the positional coordinates of the alignment results, a collinearity map was plotted, which allowed observation of whether the collapsed region had been successfully filled.

[0222] Collinearity diagrams between homologous chromosomes before and after collapse filling are shown below. Figure 5 and Figure 6 As shown.

[0223] 5. Re-align the raw Hi-C data to the fully genotyped homologous polyploid genome (i.e., the new genome) after collapse filling using the software Bwa (version: 0.7.17-r1198-dirty, parameters: mem-A1-B4-t60-E50-L0-T 0). Use the software HiCExplorer (https: / / github.com / deeptools / HiCExplorer) to display the intrachromosomal and interchromosomal interaction signals on the final version of the homologous polyploid genome, generating interaction heatmaps (full and partial Hi-C interaction signal maps are shown below). Figure 7 and Figure 8 As shown in the figure, this further verifies the effect of homologous polyploid assembly typing.

[0224] VI. Quality Assessment of the Final Genome

[0225] The quality of genome assembly is generally assessed from several important aspects, such as continuity, consistency, completeness, and redundancy. Figure 1 , 4 As shown.

[0226] 1. Continuity assessment

[0227] N50 is the final N50 of the chromosome-level homologous polyploid genome. All chromosome sequences are sorted from longest to shortest and then added sequentially. When the total length of the added sequences reaches half the total genome length, the length of the last added sequence is the N50. This step is performed and evaluated using a Python script. The N50 size is generally used to evaluate genome continuity; a larger N50 indicates better genome continuity.

[0228] Tables 5-8 below show the statistical results of various indicators of the haplotype A1, haplotype A2, haplotype A3 and haplotype A4 genomes, respectively.

[0229] Table 5. Statistical table of haplotype A1 genomic indicators

[0230] Item Value Total_length(bp) 596,605,943 total_length_withoutN(bp) 596,572,443 Total_number 12 GC_content 38.231％ N50(bp) 48,493,429 N90(bp) 45,141,909 Average (bp) 49,717,161.92 Median(bp) 48,091,811.0 Min(bp) 42,063,948 Max(bp) 57,488,531

[0231] Table 6. Statistical table of haplotype A2 genomic indicators

[0232] Item Value Total_length(bp) 611,906,938 total_length_withoutN(bp) 611,875,438 Total_number 12 GC_content 38.324％ N50(bp) 50,146,325 N90(bp) 43,287,748 Average (bp) 50,992,244.83 Median(bp) 50,050,131.0 Min(bp) 43,054,162

[0233] Table 7. Statistical table of haplotype A3 genomic indicators

[0234] Item Value Total_length(bp) 610,355,175 total_length_withoutN(bp) 610,309,775 Total_number 12 GC_content 38.338％ N50(bp) 50,148,557 N90(bp) 46,260,919 Average (bp) 50,862,931.25 Median(bp) 49,679,923.0 Min(bp) 42,907,289 Max(bp) 59,622,604

[0235] Table 8. Statistical table of haplotype A4 genomic indicators

[0236]

[0237]

[0238] 2. Consistency assessment of second-generation sequencing data

[0239] The consistency of the assembled genome is assessed by statistically analyzing the second-generation (NGD) data alignment rate and coverage. The NDD data alignment rate (map_rate) represents the proportion of NDD clean reads aligned to the reference genome out of the total number of NDD clean reads. The NDD data coverage represents the proportion of the reference genome fragment length aligned to the total length of the reference genome out of the total reference genome length. When genome consistency is good, both the NDD alignment rate and coverage are generally greater than 90%.

[0240] The software Bwa (version: 0.7.17, parameter: mem) was used to map the clean reads after quality control of second-generation high-throughput sequencing onto the reference genome, and the alignment rate and coverage of the second-generation data were statistically analyzed.

[0241] The specific statistical results of the second-generation data of haplotypeA1, haplotypeA2, haplotypeA3, and haplotypeA4 genomes are shown in Table 9-12 below.

[0242] Table 9. Statistical table of haplotype A1 genome comparison using second-generation data.

[0243] Type Number map_rate 99.85 Average_depth 30.60 Coverage 90.93

[0244] Table 10. Statistical table of haplotype A2 genome comparison using second-generation data.

[0245] Type Number map_rate 99.86 Average_depth 30.83 Coverage 90.24

[0246] Table 11. Statistical table of haplotype A3 genome comparison using second-generation data.

[0247] Type Number map_rate 99.85 Average_depth 31.40 Coverage 90.78

[0248] Table 12. Statistical table of haplotype A4 genome comparison using second-generation data.

[0249] Type Number map_rate 99.85 Average_depth 31.61 Coverage 90.29

[0250] 3. Consistency assessment of third-generation sequencing data

[0251] The consistency of the assembled genome is assessed by statistically analyzing the alignment rate and coverage of third-generation (3G) data. The alignment rate (map_rate) represents the proportion of 3G clean reads aligned to the reference genome out of all 3G clean reads. The coverage represents the proportion of the reference genome fragment length aligned to the total length of the reference genome out of all 3G clean reads. When genome consistency is good, both the alignment rate and coverage are generally greater than 90%.

[0252] The third-generation PacBio HIFI sequencing data were mapped onto the reference genome using the software minimap2 (version: 2.24-r1122, parameter: -ax map-hifi), and the alignment rate and coverage of the three generations were statistically analyzed.

[0253] The specific statistical results of the three generations of haplotypeA1, haplotypeA2, haplotypeA3, and haplotypeA4 genomes are shown in Tables 13-16 below.

[0254] Table 13. Statistical table of haplotype A1 genome comparison across three generations.

[0255] Type Number map_rate 100.00 Average_depth 26.43 Coverage 99.80

[0256] Table 14. Statistical table of haplotype A2 genome comparison data from three generations.

[0257] Type Number map_rate 100.00 Average_depth 27.70 Coverage 99.81

[0258] Table 15. Statistical table of haplotype A3 genome comparison across three generations.

[0259] Type Number map_rate 100.00 Average_depth 27.70 Coverage 99.82

[0260] Table 16. Statistical table of haplotype A4 genome comparison across three generations.

[0261] Type Number map_rate 100.00 Average_depth 28.46 Coverage 99.79

[0262] 4. Genome integrity

[0263] The BUSCO software (version: v5.4.6; parameter: –evalue 1e-05) was used to assess the integrity of genome assembly. When genome integrity is good, the proportion of Complete BUSCOs in the BUSCO software assessment is generally greater than 90%.

[0264] The genomic integrity statistics of haplotypeA1, haplotypeA2, haplotypeA3, and haplotypeA4 are shown in Table 17-20 below.

[0265] Table 17 Haplotype A1 BUSCO Evaluation Statistics

[0266] Item Number Percent (%) Complete BUSCOs(C) 1549 96.0 Complete and single-copy BUSCOs(S) 1449 89.8 Complete and duplicated BUSCOs(D) 100 6.2 Fragmented BUSCOs(F) 0 0.0 Missing BUSCOs(M) 65 4.0 Total BUSCO groups searched 1614 100.0

[0267] Table 18 Haplotype A2 BUSCO Evaluation Statistics

[0268]

[0269]

[0270] Table 19 haplotype A3 BUSCO Evaluation Statistics

[0271] Item Number Percent(%) Complete BUSCOs(C) 1531 94.9 Complete and single-copy BUSCOs(S) 1372 85.0 Complete and duplicated BUSCOs(D) 159 9.9 Fragmented BUSCOs(F) 1 0.1 Missing BUSCOs(M) 82 5.1 Total BUSCO groups searched 1614 100.0

[0272] Table 20: haplotype A4 BUSCO Evaluation Statistics

[0273] Item Number Percent(%) Complete BUSCOs(C) 1517 94.0 Complete and single-copy BUSCOs(S) 1378 85.4 Complete and duplicated BUSCOs(D) 139 8.6 Fragmented BUSCOs(F) 0 0.0 Missing BUSCOs(M) 97 6.0 Total BUSCO groups searched 1614 100.0

[0274] 5. Redundancy assessment

[0275] Because there may be secondary alignments during the alignment process, meaning that a single read may be aligned to multiple locations on the genome, secondary alignments need to be filtered out in order to reflect the true genomic redundancy.

[0276] The third-generation sequencing data (third-generation PacBio HIFI sequencing data) were then aligned to the assembled genome (non-primary alignments were filtered out so that each read was aligned to at most one location), the coverage depth of each site was calculated, the frequency of each depth was counted, and it was determined whether it conformed to a Poisson distribution.

[0277] When genome assembly results contain redundancy, the coverage depth of the redundant portion is only half that of the normal depth. Therefore, if there is a peak before the main peak (determined by the estimated genome size), it indicates redundancy (another possibility is the presence of highly heterozygous sequences); if there is a peak after the main peak, it indicates that repetitive sequences have been excessively merged (collapsed). When the genome assembly results have no redundancy, they conform to a Poisson distribution, such as... Figure 9 As shown.

[0278] VII. Comparison of Haplotype Typing and Variation Identification

[0279] like Figure 1 , 4As shown, the isolated haplotype genomes were paired and compared using minimap2 (version: 2.24-r1122, parameters: -ax asm5–eqx) or mummer (version: 4.0.0rc1, parameters: --maxmatch-c 500-b500-l 50-g 90) software to identify genomic variations and plot collinearity maps to illustrate differences between haplotypes. The specific methods are as follows.

[0280] 1. Isolate each haplotype from the complete homologous haploid genome.

[0281] 2. Use minimap2 (https: / / github.com / lh3 / minimap2) or mummer (https: / / github.com / mummer4 / mummer) to perform pairwise comparisons on each isolated haplotype.

[0282] 3. By comparing the results, the structural rearrangements in collinear and non-collinear regions are summarized, and the structural rearrangements in collinear and non-collinear regions between each group of homologous chromosomes are visualized, such as... Figure 10 As shown.

[0283] 4. Combining structural annotation, variations between haplotypes were annotated, and genes affected by these variations were extracted. Structural annotation was achieved using existing software such as Stringtie (http: / / ccb.jhu.edu / software / stringtie) and TransDecoder (https: / / github.com / TransDecoder / TransDecoder) to predict coding genes based on transcriptome evidence, Exonerate (https: / / github.com / nathanweeks / exonerate) to predict coding genes based on homologous protein evidence, and Augustus (http: / / augustus.gobics.de / ) to predict coding genes de novo using a model trained on closely related species. Finally, the annotation results were integrated using the software EVidenceModeler (https: / / github.com / EVidenceModeler / EVidenceModeler / ?tab=readme-ov-file#evidencemodeler).

[0284] 5. Combine functional annotation to investigate which functional categories, pathways, or processes the genes affected by variations may be enriched in. Functional annotation utilizes existing software such as Diamond (…). https: / / github.com / bbuchfink / diamondThe protein sequences encoded by genes obtained from structural annotations were compared with existing protein databases Uniprot and NR for sequence similarity searches. Simultaneously, motif similarity searches were performed using InterProScan, the metabolic pathway databases KEGG (version: 2023-07-11), and Pfam (version: 2022-03-11).

[0285] The above detailed embodiments describe the implementation of the present invention; however, the present invention is not limited to the specific details described in the above embodiments. Within the scope of the claims and technical concept of the present invention, various simple modifications and changes can be made to the technical solution of the present invention, and these simple modifications all fall within the protection scope of the present invention.

Claims

1. A method for typing and assembling a homologous polyploid genome, characterized in that, Includes the following steps: Sequencing was performed on the test samples to obtain second-generation library data, Hi-C library data, and third-generation PacBio HIFI sequencing data; the genome size was evaluated and estimated. The third-generation PacBio HIFI sequencing data were used for preliminary assembly, and the degree of closeness between the preliminary assembly results and the estimated genome size was evaluated to select the optimal preliminary assembly results. Using the Hi-C library data and the optimal preliminary assembly results, Hi-C mounting was performed to obtain a preliminary Hi-C mounted genome. By filling in the collapsed regions in the preliminary Hi-C mounted genome, the final version of the homopolyploid genome is obtained; The specific method for preliminary assembly using the aforementioned third-generation PacBio HIFI sequencing data is as follows: The third-generation PacBio HIFI sequencing data were initially assembled using four parameter modes of the hifiasm software. The four parameter modes are: hifiasm non-fractal mode with filtering of cells consisting of ≤3 reads during assembly, hifiasm non-fractal mode with filtering of cells consisting of ≤4 reads during assembly, hifiasm fractal mode with filtering of cells consisting of ≤3 reads during assembly, and hifiasm fractal mode with filtering of cells consisting of ≤4 reads during assembly. When using two hifiasm non-genotyping modes, the assembled genomes p_ctg and a_ctg are merged as the first preliminary assembly result; when using two hifiasm genotyping modes, the assembled genomes hap1_ctg and hap2_ctg are merged as the second preliminary assembly result. Compare the continuity and completeness of the first preliminary assembly result and the second preliminary assembly result, and select the one that is closer to 1.2 times the estimated genome size as the optimal preliminary assembly result.

2. The method for typing and assembling a homologous polyploid genome according to claim 1, characterized in that, Based on the second-generation sequencing data, a survey analysis was performed using the software GenomeScope 2.0 to obtain genome size estimation results with 2n=NX and N≥4; where n is the number of chromosomes in the gametes, X is the number of chromosomes of haplotypes in the autopolyploids, and N is the ploidy of the autopolyploids.

3. The method for typing and assembling a homologous polyploid genome according to claim 1, characterized in that, The specific method for mounting Hi-C is as follows: Hi-C data alignment software was used to align the Hi-C library data to the optimal preliminary assembly result. Then, Hi-C signal clustering software was used to cluster and mount the optimal preliminary assembly result according to the alignment results. Map adjustment software was used to adjust homologous chromosomes together according to the interaction signals of homologous chromosomes and determine the position and orientation of contigs within chromosomes. During map adjustment, more than N chromosome interaction signals were identified as redundant and removed to obtain the Hi-C preliminary mounted genome.

4. The method for typing and assembling a homologous polyploid genome according to claim 1, characterized in that, The Hi-C data comparison software uses the -SP5M mode and adopts a segment alignment method to take the alignment of the smallest query coordinate as the primary alignment and mark the shorter segment hits as secondary hits. The Hi-C signal clustering software uses the `--early-exit --mapq 0 -r 2 -i 15000 -m diploid` mode to exit early after the first round of mounting with the lowest quality threshold, thus obtaining the preliminary Hi-C mounted genome.

5. The method for typing and assembling a homologous polyploid genome according to claim 1, characterized in that, The method for filling the collapsed regions in the preliminary Hi-C mounted genome is as follows: In a group of non-redundant chromosomes that are homologous to each other, the longest chromosome is selected as the ref, and the other homologous chromosomes are selected as queries. All queries are then matched to the ref. Based on the comparison results, the ref and query are divided into several collinear blocks. The gap lengths between adjacent collinear blocks on the ref and query are compared. Based on the comparison results, the gap regions are divided into haplotype-specific regions, haplotype-collapsed regions, and insignificant collapse regions. For the haplotype-collinear regions, the gap between adjacent collinear blocks on the ref is extended 20kb to both ends as candidate sequences to be filled into the query. The ends of the candidate sequences are aligned to the last 20kb region of the corresponding collapsed region on each query. The alignment method with the highest alignment quality value and the fewest skip alignments is selected, and the sequence on the ref is automatically filled into the collapsed region of the query. Haplotype-specific regions and inconspicuous collapsed regions are not processed. The final version of the homologous polyploid genome, with the collapse filled, was obtained.

6. The method for typing and assembling a homologous polyploid genome according to claim 1, characterized in that, The final version of the homologous polyploid genome was also evaluated for continuity, consistency, integrity, and redundancy.

7. The method for typing and assembling a homologous polyploid genome according to claim 6, characterized in that, Genomic continuity is assessed by evaluating the N50 size of the genome; The clean reads after quality control of second-generation high-throughput sequencing were aligned to the final version of the homologous polyploid genome using the Bwa software. The consistency of the assembled genome was evaluated by statistically analyzing the alignment rate and coverage of the second-generation library data. Using the software minimap2, clean reads after quality control of third-generation high-throughput sequencing were aligned to the final version of the homologous polyploid genome. The alignment rate and coverage of the third-generation PacBio HIFI sequencing data were statistically analyzed to evaluate the consistency of the assembled genome. Wherein, the alignment rate represents the proportion of clean reads that have been aligned to the final version of the homologous polyploid genome out of the total number of clean reads; the coverage represents the proportion of the length of the reference genome fragment that has been aligned to the clean reads out of the total length of the final version of the homologous polyploid genome. The percentage of Complete BUSCOs was assessed using the BUSCO software to determine the integrity of genome assembly. The third-generation PacBio HIFI sequencing data were then aligned to the final version of the homologous polyploid genome. The coverage depth of each site was calculated, the frequency of occurrence of each depth was counted, and it was determined whether it conformed to a Poisson distribution.

8. A method for identifying variations in a homologous polyploid genome, characterized in that, A homologous polyploid genome is obtained by using the typing and assembly method according to any one of claims 1-7. Several sets of haplotype genomes are isolated from the homologous polyploid genome. The haplotype genomes are compared pairwise to identify the variation between the haplotype genomes.

9. The method for identifying variations in a homologous polyploid genome according to claim 8, characterized in that, The isolated haplotype genomes were paired and compared using software minimap2 or mummer. Based on the comparison results, structural rearrangements in collinear and non-collinear regions were identified. Structural and functional annotations were performed on the variations between haplotypes.

Images