A whole genome haplotyping method for long read sequencing data

Abstract

This invention discloses a whole-genome haplotype typing method for long-read sequencing data. Specifically, it relates to a whole-genome haplotype typing method for long-read sequencing data. This invention addresses the technical challenge of existing methods requiring complete sequencing before typing, resulting in long analysis cycles and inability to meet rapid response requirements. It also overcomes the shortcomings of existing methods in terms of computational efficiency, noise resistance, and continuity. The process is as follows: 1. Determine the alleles supported by the read; 2. Obtain the final haplotype blocks; 3. Process the haplotype blocks sequentially to obtain a ligation queue, and determine the haplotype blocks based on the ligation queue; in adjacent haplotype block pairs, the preceding block is denoted as Block A, and the following block as Block B; determine whether Block A and Block B should be ligated; continue until all haplotype blocks are ligated, and output the typing results.

Description

Technical Field

[0001] This invention relates to a whole-genome haplotype typing method for long-read sequencing data. Background Technology

[0002] Haplotype typing is a key technology in genomics research and precision medicine, aiming to correctly assign genetic variations (such as single nucleotide polymorphisms, insertions, and deletions) in diploid organisms to paternal or maternal chromosomes. This technology is of great significance for elucidating allele-specific expression, the genetic mechanisms of complex diseases, and the evolutionary history of populations.

[0003] Traditional haplotype genotyping methods are primarily based on second-generation short-read sequencing technology. However, limited by the short read length, they struggle to traverse long repetitive sequences or highly homozygous regions, resulting in short and discontinuous haplotype blocks. In recent years, third-generation long-read sequencing technologies, represented by Pacific Biosciences and Oxford Nanopore Technologies, have rapidly developed. Their read lengths can reach tens or even hundreds of kb, directly traversing multiple heterozygous sites and providing abundant physical linkage information for constructing high-quality, continuous haplotypes. Currently, mainstream genotyping tools based on long-read sequencing mainly employ sequencing fragment-based genotyping strategies. They infer the phase relationships of variant sites by constructing mathematical models (such as the maximum cut problem and minimum error correction models), reducing over-reliance on population reference data or pedigree information.

[0004] However, the aforementioned existing technologies generally adopt an offline batch processing mode of "sequencing first, then analysis", which means that the genotyping program can only be started after the sequencing is completely completed and the full alignment file is generated. This makes the analysis process and data output completely disconnected in time.

[0005] Haplotype typing, as an important foundation of genomics and precision medicine, has made significant breakthroughs in recent years thanks to advancements in long-read sequencing technology, but it still faces many challenges in practical applications.

[0006] The complexity of the data itself remains a core factor affecting genotyping quality. While long-read sequencing data possesses the ability to traverse a wide range of variations, its raw error rate is still higher than that of short-read data, and data quality fluctuates significantly in regions with high GC content, complex repetitive sequences, and abundant structural variations. These factors lead to incomplete variant identification and allelic linkage information in some regions, thus affecting the accuracy and continuity of haplotypes. Furthermore, existing genotyping algorithms struggle to achieve an ideal balance between accuracy, efficiency, and scalability. As the scale of sequencing data continues to expand, mainstream algorithms experience a sharp increase in computational resource consumption when processing high-coverage or large-scale samples, making memory and CPU requirements bottlenecks in practical applications. More critically, current tools lack real-time capabilities and incur significant I / O overhead. Existing algorithms are all based on post-processing of the full dataset and cannot output results incrementally during sequencing, resulting in excessively long cycles from sample collection to result feedback, which fails to meet the timeliness requirements of clinical critical care diagnosis. Furthermore, these tools require repeated reading of massive amounts of alignment files (BAM / CRAM) during runtime. Frequent disk I / O operations and non-reusable intermediate computational states significantly limit analysis speed and waste computational resources. Finally, haplotype genotyping results generally exhibit block discontinuities and fragmentation. Limited by sequencing coverage, variant density, read distribution, and algorithm parameters, genotyping results often present as multiple short blocks, making it difficult to achieve complete haplotype reconstruction at the chromosome level. Cross-block connections rely on high-quality cross-reads, but the number of such reads is limited in real-world data, further exacerbating the haplotype block fragmentation problem. This not only affects the comprehensive analysis of allele linkages but also limits the in-depth application of haplotype information in the study of complex disease mechanisms and personalized medicine. Summary of the Invention

[0007] This invention aims to solve the technical problem that existing technologies require waiting for the complete sequencing process before genotyping can be performed, resulting in long analysis cycles and an inability to meet the needs of rapid response. At the same time, it overcomes the shortcomings of existing methods in terms of computational efficiency, noise resistance, and continuity, and proposes a whole-genome haplotype genotyping method for long-read sequencing data.

[0008] The specific process of a whole-genome haplotype typing method for long-read sequencing data is as follows:

[0009] Step 1: Receive sequencing reads and variant sites, and perform real-time multiple data filtering on the received sequencing reads and variant sites to obtain filtered sequencing reads and variant sites.

[0010] Based on the filtered sequencing reads and variant sites, effective reads are obtained;

[0011] Based on the effective reads, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site.

[0012] Each valid read segment is compared with the template A sequence to obtain the corresponding template A alignment result;

[0013] Each valid read segment is compared with the template B sequence to obtain the corresponding template B alignment result;

[0014] Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each valid read segment and the template A sequence;

[0015] Based on the alignment results of the corresponding template B, the edit distance 2 between each valid read segment and the template B sequence is calculated using the bit vector dynamic programming method;

[0016] By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined;

[0017] Step 2: Receive the selected SNP sites from Step 1 one by one. For newly added SNP sites... Perform dynamic grouping;

[0018] Construct an SNP graph;

[0019] Based on edge weights, the connection patterns between adjacent SNP sites within each encapsulated SNP group are determined, resulting in local monolithic blocks.

[0020] Further optimization of the local single-unit blocks yields the final single-unit block;

[0021] Step 3: Reorder and flatten the single-unit blocks output in Step 2 to obtain a queue to be connected, and determine the single-unit blocks based on the queue to be connected;

[0022] For adjacent single-block pairs, the first block in the adjacent single-block pair is denoted as Block A, and the second block is denoted as Block B;

[0023] Scan all sequencing reads aligned to the regions covered by Block A and Block B to identify all potential cross-block reads that span the end of Block A and the beginning of Block B.

[0024] A strict screening strategy is applied to the identified potential cross-block read segments to obtain the cross-block read segments;

[0025] For any two adjacent single-type blocks, Block A and Block B, there are only two logical possibilities for the single-type connection relationship between Block A and Block B; construct hypothetical models for each of the two possibilities.

[0026] Calculate the minimum error correction scores for Mode 1 and Mode 2;

[0027] Based on the minimum error correction MEC score of Mode 1 and the minimum error correction MEC score of Mode 2, determine whether Block A and Block B are connected.

[0028] Continue until all single-unit blocks are connected, then output the fractal result.

[0029] The beneficial effects of this invention are as follows:

[0030] This invention aims to provide a real-time whole-genome haplotype typing method for long-read sequencing data, in order to solve the technical problem that existing technologies must wait for the sequencing to be completely completed before typing can be performed, resulting in long analysis cycles and failure to meet the requirements for rapid response. At the same time, it overcomes the shortcomings of existing methods in terms of computational efficiency, noise resistance and continuity.

[0031] Specifically, this invention constructs a streaming processing architecture that receives mutation and sequencing data streams in real time and adopts an incremental graph model update strategy, breaking the limitations of traditional offline analysis and realizing "sequencing and genotyping on the same time". It can output high-precision haplotype results in real time during the sequencing process.

[0032] At the core algorithm level, this invention introduces a bit vector dynamic programming algorithm for local sequence re-alignment, effectively eliminating reference bias and correcting read support information. It utilizes an optimized SNP grouping strategy and a dynamic graph model to construct local haplotypes, decomposing the globally complex problem into manageable local problems, significantly reducing computational complexity while maintaining accuracy. Finally, it uses cross-block long read information and the minimum error correction (MEC) principle to connect local haplotype blocks in real time. This invention can significantly shorten the analysis cycle while ensuring high genotyping accuracy and constructing longer, more continuous haplotype blocks, making it suitable for various long-read sequencing datasets such as PacBio HiFi and Oxford Nanopore. Attached Figure Description

[0033] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0034] Specific Implementation Method 1: The specific process of this implementation method for whole-genome haplotype typing of long-read sequencing data is as follows:

[0035] This invention proposes a real-time whole-genome haplotype typing method for long-read sequencing data. Addressing the limitation of existing technologies that cannot output results in real-time during sequencing, this invention designs a streaming computing and incremental update architecture. The architecture is configured to interface with an upstream real-time variant detection module, capable of continuously receiving and processing its output variant site streams and sequencing read streams.

[0036] The core idea of ​​this invention is to transform the traditional "full post-processing" mode into a "data stream driven" mode. The system maintains a processing buffer in memory. The steps one (incremental re-alignment), two (dynamic local construction), and three (real-time global connection) described below are not executed sequentially all at once, but are repeatedly called and iteratively executed by the system as the sequencing data stream continues to be input.

[0037] Step 1: Receive sequencing reads and variant sites, and perform real-time multiple data filtering on the received sequencing reads and variant sites to obtain filtered sequencing reads and variant sites.

[0038] Based on the filtered sequencing reads and variant sites, effective reads are obtained;

[0039] Based on the effective reads, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site.

[0040] Each valid read segment is compared with the template A sequence to obtain the corresponding template A alignment result;

[0041] Each valid read segment is compared with the template B sequence to obtain the corresponding template B alignment result;

[0042] Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each valid read segment and the template A sequence;

[0043] Based on the alignment results of the corresponding template B, the edit distance 2 between each valid read segment and the template B sequence is calculated using the bit vector dynamic programming method;

[0044] By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined;

[0045] Step 2: Receive the selected SNP sites from Step 1 one by one. For newly added SNP sites... Perform dynamic grouping;

[0046] Construct an SNP graph; the edges between SNP sites in the SNP graph use allele information supported by reads;

[0047] Based on edge weights, the connection patterns between adjacent SNP sites within each encapsulated SNP group are determined, resulting in local monolithic blocks.

[0048] Further optimization of the local single-unit blocks yields the final single-unit block;

[0049] Step 3: Reorder and flatten the single-unit blocks output in Step 2 to obtain a queue to be connected, and determine the single-unit blocks based on the queue to be connected;

[0050] For adjacent single-block pairs, the first block in the adjacent single-block pair is denoted as Block A, and the second block is denoted as Block B;

[0051] Scan all sequencing reads aligned to the regions covered by Block A and Block B (e.g., Block A is [a1, a2] and Block B is [b1, b2], this region is the two regions) to identify all potential cross-block reads that span the end position of Block A and the beginning position of Block B.

[0052] To ensure the reliability of the connection evidence, this invention implements a strict screening strategy for the identified potential cross-block reads to obtain cross-block reads;

[0053] For any two adjacent single-type blocks, Block A and Block B, there are only two logical possibilities for the single-type connection relationship between Block A and Block B; construct hypothetical models for each of the two possibilities.

[0054] Calculate the minimum error correction (MEC) scores for Mode 1 and Mode 2;

[0055] Based on the minimum error correction MEC score of Mode 1 and the minimum error correction MEC score of Mode 2, determine whether Block A and Block B are connected.

[0056] Continue until all single-unit blocks are connected, then output the fractal result.

[0057] Specific Implementation Method Two: This implementation method differs from Specific Implementation Method One in that, in step one, sequencing reads and variant sites are received, and real-time multiple data filtering is performed on the received sequencing reads and variant sites to obtain filtered sequencing reads and variant sites.

[0058] Based on the filtered sequencing reads and variant sites, effective reads are obtained;

[0059] Based on the effective reads, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site.

[0060] Each valid read segment is compared with the template A sequence to obtain the corresponding template A alignment result;

[0061] Each valid read segment is compared with the template B sequence to obtain the corresponding template B alignment result;

[0062] Based on the alignment results of the corresponding template A, the edit distance 1 between each valid read segment and the template A sequence is calculated using the bit vector dynamic programming method;

[0063] Based on the alignment results of the corresponding template B, the edit distance 2 between each valid read segment and the template B sequence is calculated using the bit vector dynamic programming method;

[0064] By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined;

[0065] The specific process is as follows:

[0066] Step 11:

[0067] This invention is built on a real-time sequencing and analysis architecture. It does not directly perform whole-genome alignment and variant detection, but instead connects to the upstream real-time variant detection module in real time through a streaming interface.

[0068] The real-time variant calling framework (RVC, also known as https: / / github.com / CuiMiao-HIT / ) software outputs variant detection results for sequencing reads.

[0069] Receive two types of concurrent data streams sorted by genome coordinates from the output of real-time variant detection software;

[0070] The first type of data stream is sequencing reads; this is the read data that has completed preliminary alignment and is output in real time by upstream tools.

[0071] The second type of data stream is variant sites; these are candidate variant sites detected in real time by upstream tools; the variant site stream includes the genomic coordinates of the variant, reference bases, and variant base information;

[0072] Steps 1 and 2: To ensure the efficiency and accuracy of real-time processing, real-time multiple data filtering is performed on the sequencing reads and variant sites received in step 1 to obtain filtered sequencing reads and variant sites. Data that does not meet the conditions will be directly discarded without occupying subsequent computing resources.

[0073] Step 13: Based on the filtered sequencing reads and variant sites obtained in Step 12, obtain the effective reads;

[0074] Step 14: Based on the effective reads from Step 13, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site;

[0075] Step 15

[0076] Each valid read obtained in steps one and three is compared with the template A sequence to obtain the corresponding template A alignment result;

[0077] Each valid read obtained in steps one and three is compared with the template B sequence to obtain the corresponding template B alignment result;

[0078] Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each effective read segment obtained in steps one and three and the template A sequence;

[0079] Based on the alignment results of the corresponding template B, the bit vector dynamic programming method is used to calculate the edit distance 2 between each effective read segment obtained in steps one and three and the template B sequence;

[0080] By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined.

[0081] The other steps and parameters are the same as in Specific Implementation Method 1.

[0082] Specific Implementation Method 3: This implementation method differs from Specific Implementation Method 1 or 2 in that, in order to ensure the efficiency and accuracy of real-time processing, in steps 1 and 2, real-time multiple data filtering is performed on the sequencing reads and variant sites received in step 1 to obtain filtered sequencing reads and variant sites. Data that does not meet the conditions will be directly discarded without occupying subsequent computing resources.

[0083] The specific process is as follows:

[0084] Step 121: Perform real-time multiple data filtering on the sequencing reads received in Step 11 to obtain the filtered sequencing reads; the specific process is as follows:

[0085] 1) For the sequencing reads received in real time in step one The sequencing reads are preserved, and step 2) is executed.

[0086] In the sequencing reads received in real time in step one Delete the sequencing reads (sequencing reads come with their own...) ,like (Delete the read entry);

[0087] In minimap2, this value represents approximately 90% confidence in the read segment mapping location. The value calculation is based on the comparison tool's internal algorithm, taking into account multiple factors such as matching quality, uniqueness, and read segment length. Essentially, it is:

[0088]

[0089] in It represents the probability of a misalignment. Choosing a threshold of 5 is a balance point derived from multiple sets of experimental data, effectively filtering out ambiguous reads with multiple location mappings while retaining most reliable reads.

[0090] 2) Use SAMtools markdup or Picard MarkDuplicates tools to determine whether PCR repetitive sequences exist in the retained sequencing reads;

[0091] If present, delete the PCR repeat sequences to ensure that only one repeat sequence is retained (PCR repeat sequences within each sequencing read are labeled, so the repeat portion is deleted directly), and proceed to step 3).

[0092] If no PCR repeat sequence is found, proceed to step 3).

[0093] PCR duplication: Although the problem of PCR duplication is not as prominent in long-read sequencing as it is in short-read sequencing, duplicative reads can still be generated in some library preparation methods;

[0094] 3): Use the -F 0x900 parameter in SAMtools to determine whether there are sequences marked FLAG0x100 or FLAG 0x800 in the retained sequencing reads;

[0095] If they exist, the sequences marked with FLAG 0x100 or FLAG 0x800 will be deleted to obtain the filtered sequencing reads;

[0096] If it does not exist, the filtered sequencing reads are obtained;

[0097] Use the -F 0x900 parameter in SAMtools to filter out minor alignments (FLAG 0x100) and supplementary alignments (FLAG 0x800). In the FLAG system of BAM files, 0x100 (decimal 256) indicates a minor alignment, i.e., a non-optimal alignment position; 0x800 (decimal 2048) indicates a supplementary alignment, which is usually part of the original alignment but has been split.

[0098] Step 122: Perform real-time multiple data filtering on the variant sites received in step 11 to obtain the filtered variant sites;

[0099] The specific process is as follows:

[0100] Read the QUAL and FILTER fields from the variant sites, retain the variant sites marked as PASS in the QUAL and FILTER fields and whose quality values ​​are higher than a preset threshold (such as Q20), and delete the remaining variant sites;

[0101] Simultaneously, the frequency (AF) information of variant sites is obtained, and extreme variant sites with a frequency less than 0.2 or a frequency greater than 0.8 are removed to eliminate sequencing errors or chimeric interference, and only high-quality heterozygous variant sites are retained for subsequent genotyping.

[0102] Other steps and parameters are the same as in specific implementation method one or two.

[0103] Specific Implementation Method Four: This implementation method differs from one of Specific Implementation Methods One to Three in that, in step one and three, effective reads are obtained based on the filtered sequencing reads and variant sites obtained in step one and two.

[0104] The specific process is as follows:

[0105] Based on the assumption of coordinate ordering of sequencing reads and variant sites, the filtered sequencing reads obtained in steps one and two are placed into the sequencing read buffer, and the filtered variant sites obtained in steps one and two are placed into the variant site buffer.

[0106] When a new SNP site is received in step one, the genomic coordinates Pos of the new SNP site are recorded;

[0107] Compare the termination coordinates of each sequencing read in the sequencing read buffer with the termination coordinates of the new SNP site;

[0108] If the termination coordinate of a sequencing read in the buffer is less than the termination coordinate of a new SNP site (it can be determined that no new variant sites will fall within the coverage of the read. At this time, the variant coverage associated with the read is considered to be completely determined), calculate the number of variant sites K in the corresponding sequencing read; if K is less than 2, the corresponding sequencing read is determined to be unable to provide valid haplotype linkage information, and the corresponding sequencing read is removed from the buffer and discarded; if K is greater than or equal to 2, the corresponding sequencing read is determined to be a valid read and is sent to the subsequent re-alignment and genotyping process;

[0109] If the termination coordinates of all sequencing reads in the buffer are greater than or equal to the termination coordinates of the new SNP site, no action is taken.

[0110] The other steps and parameters are the same as those in one of the specific implementation methods one to three.

[0111] Specific Implementation Method 5: This implementation method differs from one of the specific implementation methods one to four in that, in step one to four, based on the effective read segment in step one to three, a window corresponding to each SNP site in the effective read segment is obtained, and two target sequence templates are constructed based on the window corresponding to each SNP site.

[0112] The specific process is as follows:

[0113] Since the initial alignment is based on a reference genome, reads carrying variant alleles may have lower alignment scores or minor misalignments near the variant sites, a phenomenon known as reference bias. To address this issue, this invention implements local realignment.

[0114] Step 141:

[0115] For each SNP site in the effective reads obtained in steps one and three, a fixed length of base sequence is extracted upstream and downstream of each SNP site, with the coordinates of each SNP site on the reference genome as the center (the preferred window size in this invention is 15bp, i.e., a total length of 31bp), to obtain the window (31bp) corresponding to each SNP site.

[0116] The SNP locus records the reference allele (REF) and the variant allele (ALT).

[0117] If the distance (10bp) from a certain SNP site to the upstream endpoint of the reference chromosome genome is less than a fixed length (15bp), then the base sequence from the center to the endpoint is extracted upstream from the coordinates of the SNP site on the reference genome, and a fixed length of base sequence is extracted downstream from the center (the window size is preferably 15bp in this invention), to obtain the window (10+1+15=26bp) corresponding to each SNP site.

[0118] If the distance (10bp) from a certain SNP site to the downstream endpoint of the reference chromosome genome is less than a fixed length (15bp), then the base sequence from the center to the endpoint is extracted downstream from the coordinates of the SNP site on the reference genome, and the fixed length base sequence is extracted upstream from the center (the window size is preferably 15bp in this invention), to obtain the window (10+1+15=26bp) corresponding to each SNP site.

[0119] Step 142: Construct two target sequence templates based on the window corresponding to each SNP site:

[0120] Template A sequence: The original sequence of the reference genome within the window (31 bp) obtained in step 141;

[0121] The central base of the window obtained in step 141 is the reference allele (REF).

[0122] The template A sequence is a reference sequence template;

[0123] Template B sequence: Replace the central base of the window obtained in step 141 with the variant allele (ALT), and keep the remaining sequences (15bp above and below) in the window consistent with the reference genome.

[0124] The template B sequence is a variant sequence template.

[0125] The other steps and parameters are the same as those in one of the specific implementation methods one to four.

[0126] Specific Implementation Method Six: This implementation method differs from one of Specific Implementation Methods One to Five in that, in step one five, each valid read segment obtained in step one three is compared with the template A sequence to obtain the corresponding template A comparison result;

[0127] Each valid read obtained in steps one and three is compared with the template B sequence to obtain the corresponding template B alignment result;

[0128] Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each effective read segment obtained in steps one and three and the template A sequence;

[0129] Based on the alignment results of the corresponding template B, the bit vector dynamic programming method is used to calculate the edit distance 2 between each effective read segment obtained in steps one and three and the template B sequence;

[0130] By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined;

[0131] The specific process is as follows:

[0132] Step 151

[0133] Each valid read obtained in steps one and three is compared with the template A sequence to obtain the corresponding template A alignment result;

[0134] Each valid read obtained in steps one and three is compared with the template B sequence to obtain the corresponding template B alignment result;

[0135] To efficiently determine whether each read covering an SNP supports a reference allele or a variant allele, this invention employs an optimized bit-vector dynamic programming algorithm. This algorithm is a variant of the Myers algorithm. It first extracts fixed-length bases before and after the SNP site to construct two sequences: one containing the reference allele and the other containing the variant allele, and then aligns these sequences with the reads.

[0136] Step 152

[0137] Based on the alignment results of the corresponding template A obtained in step 1-51, the bit vector dynamic programming method is used to calculate the edit distance 1 between each effective read segment obtained in step 1-3 and the template A sequence;

[0138] Based on the alignment results of the corresponding template B obtained in step 151, the bit vector dynamic programming method is used to calculate the edit distance 2 between each effective read segment obtained in step 13 and this template B sequence;

[0139] The specific process is as follows:

[0140] 1) Initialize the current score ;

[0141] make , Represents the target sequence The One location;

[0142] make , This represents the query sequence obtained in steps one and three. The Middle One location;

[0143] target sequence It can be either template A sequence or template B sequence;

[0144] 2) Obtain the target sequence The Character matching information at each position ; indicates as:

[0145] in, express The character matching vector;

[0146] The query sequence obtained in step one or three Each valid read segment is ;

[0147] like ,but =1, otherwise =0 ( (For A, C, T, G).

[0148] For characters in the alphabet ,like ,but The Middle The bit is 1 if it is not 0 otherwise. (For A, C, T, G).

[0149] Represents the target sequence The One location; Represents the query sequence The Middle One effective reading segment;

[0150] express With query sequence Matching results;

[0151] 3) Based on the first Character matching information at each position Calculate the vertical auxiliary vector ; indicates as: ;

[0152] in, express and bitwise OR operation;

[0153] Represents the vertical negative difference vector (-1); during initialization, let... It is a vector consisting entirely of zeros;

[0154] 4) Calculate the horizontal auxiliary vector based on the vertical auxiliary vector; expressed as: ;

[0155] in, express and bitwise OR operation;

[0156] Represents the vertical positive difference vector (+1); during initialization, set... It is a vector consisting entirely of 1s;

[0157] Represents a horizontal auxiliary vector;

[0158] This represents the bitwise AND operator. This represents the bitwise XOR operator.

[0159] 5) Calculate the horizontal positive difference vector based on the horizontal auxiliary vector. ; indicates as: ;

[0160] ;

[0161] in, express and bitwise OR operation;

[0162] Indicates bitwise inversion;

[0163] 6) Update the edit distance in real time using the highest bit of the horizontal difference vector; expressed as:

[0164] If the horizontal positive difference vector If the highest bit is 1, then the edit distance is... ;

[0165] If the horizontal negative difference vector If the highest bit is 1, then the edit distance is... ;

[0166] The highest position is the first one on the left;

[0167] For example, in a binary representation of 1010001, the leftmost '1' is the most significant bit.

[0168] 7) Update the vertical difference vector based on the real-time updated edit distance; expressed as:

[0169]

[0170] in, express Shift one position to the left; express Shift one position to the left;

[0171] 8) Order Repeat steps 2)-7) until the target sequence has been traversed. After all positions, the final edit distance This is the query sequence obtained in steps one and three. Each valid read segment Edit distance between the template A sequence and the template B sequence;

[0172] And select The central base corresponding to a smaller template A sequence or template B sequence is used as the base type supported by the read at the SNP site;

[0173] 9) Order Repeat steps 2)-8) until the query sequence obtained in steps 1 and 3 has been traversed. All valid reading segments ;

[0174] Step 153

[0175] The alleles supported by the read were determined by comparing edit distance 1 and edit distance 2. The process is as follows:

[0176] Set the edit distance threshold.

[0177] Alleles corresponding to edit distances less than the edit distance threshold of 1 are considered to be supported by reads;

[0178] Alleles corresponding to edit distances greater than or equal to the edit distance threshold of 1 are considered as reads that are not supported.

[0179] Alleles corresponding to edit distances less than the edit distance threshold of 2 are considered to be read-supported.

[0180] Alleles corresponding to edit distances greater than or equal to the edit distance threshold of 2 are considered as reads that are not supported.

[0181] The bit vector method utilizes the bit manipulation capabilities of computers to encode each column of a dynamic programming matrix as a binary vector.

[0182] The other steps and parameters are the same as those in one of the specific implementation methods one to five.

[0183] Specific Implementation Method Seven: This implementation method differs from Specific Implementation Methods One through Six in that, in step two, the screened SNP sites from step one are received one by one. For newly arriving SNP sites... Perform dynamic grouping;

[0184] Construct an SNP graph;

[0185] Based on edge weights, the connection patterns between adjacent SNP sites within each encapsulated SNP group are determined, resulting in local monolithic blocks.

[0186] Further optimization of the local single-unit blocks yields the final single-unit block;

[0187] The specific process is as follows:

[0188] Step 21: Receive the filtered SNP sites from Step 15 one by one. For newly added SNP sites... Dynamic grouping is performed; the specific process is as follows:

[0189] Haplotype typing across the entire genome is typically an NP-hard problem with extremely high computational complexity. To reduce the problem size, enable parallel computing, and adapt to the continuous input characteristics of data streams, this invention transforms the traditional static grouping strategy into a dynamic greedy grouping mechanism. A "current SNP group" is maintained in memory and dynamically determined based on the real-time inflow of variant sites.

[0190] Receive the selected SNP sites one by one from step one five, and for newly added SNP sites... Perform the following truncation criteria:

[0191] (1) Physical length limit: The current SNP group is added to the new inflow Afterwards, the total physical length of the current group exceeds a preset threshold. (In a preferred embodiment of the present invention, (Set to 2,250,000 bp); this limit is to prevent memory overflow caused by a single computational task being too large.

[0192] (2) Number of sites limited: New inflows into the current SNP group Subsequently, the total number of SNPs in the current group exceeds a preset threshold. (In a preferred embodiment of the present invention, (Set to 2,250); this limitation is to control the computational complexity of graph model construction and ensure the efficiency of algorithm operation.

[0193] (3) Adjacent spacing restrictions: newly flowing in Compared with the last SNP site in the current SNP group The physical distance between them exceeds a preset threshold (In a preferred embodiment of the present invention, (Set to 300,000 bp); When two variant sites are too far apart, long-read sequencing data often cannot provide enough span reads to establish reliable physical linkage. In this case, forcibly breaking the grouping can reduce invalid computation and reduce the risk of erroneous linkage.

[0194] If any of the above truncation conditions are met, the current SNP group is determined to be completed. The current SNP group is immediately encapsulated as a task unit to be processed and pushed into the graph construction engine. Simultaneously, a new group is initialized to accommodate newly incoming SNPs. This mechanism ensures that the local monomer construction task is atomic and temporally independent.

[0195] To prevent insufficient coverage due to segment read delays, the system does not immediately calculate the defined SNP groups, but instead pushes them into the queue to be calculated.

[0196] Real-time monitoring of the minimum starting mapping coordinates of all unprocessed valid read segments in the buffer of step one and three, denoted as... For the SNP group that has been encapsulated and is to be calculated, the termination physical coordinates are... Continuous judgment and Size, if Greater than or equal to Continue monitoring; if Less than This means that the starting positions of existing valid read segments and subsequent new inflow read segments in the buffer are all located after the SNP group, and there is no physical possibility of overwriting them. At this point, it is determined that the read segment collection work of the encapsulated SNP group to be calculated has been completed, and the subsequent graph construction and path search process begins.

[0197] Step 22: Construct the SNP diagram; the specific process is as follows:

[0198] For each encapsulated SNP group obtained in step 2, a Directed Acyclic Graph (DAG) model is constructed. The DAG model is called the "SNP graph". This step directly uses the read segment data output in step 1.

[0199] Node definition in a directed acyclic graph: For each encapsulated SNP group obtained in step two, each SNP site within that group. (in (For group indexing), create two nodes in the graph, labeled as follows: and ;

[0200] in, The reference allele representing the current locus.

[0201] The variant allele representing the current locus.

[0202] Edge definition in a directed acyclic graph: Directed edges are established only between adjacent SNP sites; for adjacent sites... and There are four possible connecting edges:

[0203] side :connect and This indicates that both loci are reference alleles;

[0204] side :connect and This indicates that both loci are variant alleles;

[0205] side :connect and , indicating the reference allele and the variant allele;

[0206] side :connect and This indicates that the variant allele is linked to the reference allele.

[0207] Weight calculation: The weight of an edge reflects the strength of the support that the sequencing data provides for the current connectivity relationship;

[0208] Based on the SNP sites within each encapsulated SNP group obtained in step 21 and The read segment data is used to count the number of read segments that support each type of connection relationship, and the count is assigned to the corresponding edge weight.

[0209] For example, if there is Reading segment at the location Supporting reference bases and at the site Supporting variant bases, then the edge weight ;

[0210] If Reading segment at the location Supporting reference bases and at the site Support reference base, then edge weight ;

[0211] If Reading segment at the location Supports variant bases and at the site Supporting variant bases, then the edge weight ;

[0212] If Reading segment at the location Supports variant bases and at the site Support reference base, then edge weight ;

[0213] Steps 2 and 3: Determine the connection pattern between adjacent SNP sites within each encapsulated SNP group based on edge weights to obtain the connected SNP groups.

[0214] The specific process is as follows:

[0215] After constructing the SNP graph, the problem is transformed into finding two non-intersecting optimal paths in the graph, which represent haplotype sequences of two homologous chromosomes respectively. To solve this problem, this invention proposes a heuristic path search algorithm based on edge weights, which analyzes the connection patterns between adjacent SNP sites one by one.

[0216] Step 231: For any adjacent SNP pairs within each encapsulated SNP group... This examines two mutually exclusive single-type phase combination schemes (Phasing Configuration); the specific process is as follows:

[0217] Option 1 (parallel connection): Assume adjacent SNP pairs The two monomer types are 0-0 and 1-1, where 0-0 represents... The two haplotypes are the reference allele and the reference allele, respectively, and 1-1 indicates that... The two haplotypes are the variant allele and the variant allele, respectively; the formula for calculating the total support of Scheme 1 is:

[0218] ;

[0219] Option 2 (Cross-connection): Assume adjacent SNP pairs The two monomer types are 0-1 and 1-0, where 0-1 represents... The two monomer types are the reference base 0 and the variant base 1; 1-0 indicates The two monomer types are the variant base 1 and the reference base 0, respectively; the formula for calculating the total support of Scheme 2 is:

[0220] ;

[0221] in, This indicates the total support for Option 1;

[0222] express and The weights between them;

[0223] express and The weights between them;

[0224] This indicates the total support for Option 2;

[0225] express and The weights between them;

[0226] express and The weights between them;

[0227] Indicates the reference allele at the current locus;

[0228] This indicates the alternative allele at the current locus.

[0229] Indicates site Reference allele;

[0230] Indicates site Alternative alleles;

[0231] To ensure the reliability of the typing results, this invention introduces a rigorous decision-making mechanism to evaluate the reliability of the connections:

[0232] Step 232: To ensure the reliability of the classification results, this invention introduces a rigorous decision-making mechanism to evaluate the reliability of the connection;

[0233] Set minimum total support threshold ;

[0234] The total support of Option 1 Total support for Option 2 The total support is obtained by adding them together. ;

[0235] Total support With minimum total support threshold Compare;

[0236] like If so, the evidence is considered insufficient, and adjacent SNP pairs are considered... and No connection is established;

[0237] like Then proceed to step two or three.

[0238] Step Two Three Three

[0239] Set support ratio threshold ;

[0240] Based on the total support of Scheme 1 Total support for Option 2 Calculate the support ratio :

[0241]

[0242] if If the support for the two schemes is very close, then it is considered that the correct connection method cannot be determined, and adjacent SNP pairs... and No connection is established;

[0243] if The total support for option one is selected. Total support for Option 2 Adjacent SNP pairs in schemes with high overall support and Establish a connection;

[0244] Steps 2 and 4: The connections obtained in Steps 2 and 3... The group is further optimized to obtain the final single-unit block;

[0245] The specific process is as follows:

[0246] The initial single-unit blocks generated through the above path search may still contain a small number of errors. This invention further implements an optimization step based on read segment consistency. For the connected blocks obtained in steps two and three... For each monolithic block in the group, rescan all read segments that cover the current monolithic block and calculate the consistency score between each read segment and the monolithic block;

[0247] Ideally, a read from a chromosome should carry alleles that perfectly match the haplotype sequence of the corresponding chromosome. For example, suppose a haplotype block has haplotype information at four SNP sites: H1: [0,0,1,0] and H2: [1,1,0,1] (where 0 represents a reference base and 1 represents a variant base). If a read's detection result at these four sites is [REF, REF, ALT, REF] (i.e., [0,0,1,0]), then it perfectly supports H1, covering all four SNP sites; if the result is [ALT, ALT, REF, ALT] (i.e., [1,1,0,1]), then it perfectly supports H2.

[0248] The process is as follows:

[0249] For each read segment, calculate the relationship between each read segment and... The consistent number of sites is obtained Consistency score ; Statistics on each reading segment and The consistent number of sites is obtained Consistency score ;

[0250] Calculate support confidence :

[0251]

[0252] Read segments with a support confidence level greater than the threshold (0.65) are selected as high-consistency read segments;

[0253] in, Indicates segment coverage The quantity is counted when accepting SNPs and reading segments;

[0254] After selecting a set of reads with high confidence, these reads are used to analyze each read within a single-unit block. The site undergoes secondary verification; the process is as follows:

[0255] If a certain single-type block The phase of the current site conflicts with the phase of a high-consistency read that covers 75% (e.g., a threshold of 0.75) of the current site (if it covers the current site). If the phase of the allele at a high-consistency read at that locus is inconsistent with the phase of the locus, then the phase of the current high-consistency read locus is flipped (1 becomes 2, 2 becomes 1); if the conflict cannot be resolved, then the current high-consistency read locus is removed from the haplotype block to obtain the final haplotype block, so as to ensure that the final output local haplotype block has a high accuracy.

[0256] If a certain single-type block If the phase of a site does not conflict with the phase of 75% (e.g., a threshold of 0.75) of highly consistent reads covering the current site, then the phase is retained.

[0257] The other steps and parameters are the same as those in one of the specific implementation methods one to six.

[0258] Specific Implementation Method Eight: This implementation method differs from any of the specific implementation methods one to seven in that, in step three, the single-unit blocks output in step two are sequentially reordered and flattened to obtain a queue to be connected, and the single-unit blocks are determined based on the queue to be connected.

[0259] For adjacent single-unit block pairs, the preceding block in the adjacent single-unit block pair is denoted as... The last block is denoted as ;

[0260] Scan all matches and The covered area (e.g.) The interval is The following yes The sequencing reads (which are the two intervals) were used to identify all reads that spanned these two intervals in location. End position and The potential cross-block read segment at the start position;

[0261] To ensure the reliability of the connection evidence, this invention implements a strict screening strategy for the identified potential cross-block reads to obtain cross-block reads;

[0262] For any two adjacent single-type blocks and , and The single-unit connection relationship between them has only two logical possibilities; hypothetical models are constructed for each of the two possibilities;

[0263] Calculate the minimum error correction (MEC) scores for Mode 1 and Mode 2;

[0264] Based on the minimum error correction MEC scores of Mode 1 and Mode 2, determine... and To establish a connection;

[0265] Continue until all single-unit blocks are connected, then output the fractal result;

[0266] The specific process is as follows:

[0267] After step two, the system will output a series of highly accurate "local single-unit blocks" in parallel. Although these blocks contain... The phase relationship is determined, but the relative phase relationship between blocks has not yet been established. Given the read length advantage of third-generation sequencing data, this invention utilizes long reads capable of simultaneously spanning two adjacent haplotype blocks (i.e., cross-block reads) to connect discrete local blocks into chromosome-level global haplotypes based on the Minimum Error Correction (MEC) principle. To achieve real-time output, this step first constructs a strictly ordered haplotype block stream through reordering, and then uses a greedy strategy to call standard connection algorithms for processing. The connection strategy between blocks specifically includes four sub-processes: identification and screening of cross-block reads, construction of connection patterns, calculation of MEC scores, and connection decision.

[0268] Step 31: Reorder and flatten the single-unit blocks output in Step 2 to obtain a queue to be connected, and determine the single-unit blocks based on the queue to be connected;

[0269] The specific process is as follows:

[0270] Since parallel computing may cause the output of step two to be out of order (e.g., the third group is completed before the second group), this step maintains a "reordering buffer" and a counter for the expected group number;

[0271] (1) Reorder the single-unit blocks output in step two: buffer the out-of-order received blocks. Group, until the received SNP group number equals the expected group number;

[0272] (2) Flatten the reordered single-unit blocks: Once the final single-unit blocks obtained in step two and four are obtained in sequence, immediately flatten the final single-unit blocks obtained in step two and four (such as...) within They are pushed into the ligation queue in genomic order;

[0273] (3) Greedy loop: The connection engine maintains a pointer to the "global tail block";

[0274] Continuously pop new blocks from the queue to be connected, and treat the global tail block and the new block as single-type blocks, and send them to the four sub-processes of step 32, step 33, step 34 and step 35 in sequence for connection determination.

[0275] The global tail block is popped from the queue; it doesn't exist initially. The first block popped from the queue is the first global tail block. Then, a new block is popped. These two blocks are processed in steps 32 to 35. After processing is complete, the latest block is used as the global tail block.

[0276] Step 3.2: For adjacent pairs of single-unit blocks, the previous block in the adjacent pair is denoted as... The last block is denoted as ;

[0277] Scan all matches and The covered area (e.g.) The interval is The following yes The sequencing reads (which are the two intervals) were used to identify all reads that spanned these two intervals in location. End position and The potential cross-block read segment at the start position;

[0278] To ensure the reliability of the connection evidence, this invention implements a strict screening strategy for the identified potential cross-block reads to obtain cross-block reads;

[0279] The screening strategy is as follows:

[0280] (1) Effective coverage verification:

[0281] Requires potential cross-block read segments in At least one already classified site, and in It also covers at least one already classified type. site;

[0282] Potential cross-block reads that do not meet the requirements cannot provide phase linkage information and are therefore eliminated.

[0283] (2) Activity level check:

[0284] Calculate the overlap length of potential cross-block reads within each monolithic block;

[0285] The overlap length between the potential cross-block read base sequence length and the monomer block base sequence length must be at least one-fifth (1 / 5) of the total length of the monomer block.

[0286] Potential cross-block read segments that do not meet the requirements are removed;

[0287] This strategy aims to exclude low-quality reads that only briefly touch at the edges, ensuring that reads have sufficient anchoring length within both blocks to provide robust evidence of connectivity.

[0288] Step 3: For any two adjacent single-unit blocks and , and The single-unit connection relationship between them has only two logical possibilities; hypothetical models are constructed for each of the two possibilities; the specific process is as follows:

[0289] Mode 1 (KepT):

[0290] Monotype 1 and Haplotype 1 belongs to the same chromosome and forms a haplotype. ;

[0291] Monotype 2 and Haplotype 2 belongs to another chromosome, forming a haplotype. ;

[0292] Mode 2 (Reverse Mode):

[0293] Monotype 1 and The monomers are connected to form monomers. ;

[0294] Monotype 2 and The monomers 1 are connected to form monomers. ;

[0295] Steps 3 and 4: Calculate the minimum error correction scores for Mode 1 and Mode 2;

[0296] Step 3.5: Based on the minimum error correction MEC score of Mode 1 and the minimum error correction MEC score of Mode 2, determine... and To establish a connection;

[0297] Step 36: Repeat steps 31 to 35 to complete the connection of all single-unit blocks and output the fractal result.

[0298] The other steps and parameters are the same as those in any of the specific implementation methods one to seven.

[0299] Specific Implementation Method Nine: This implementation method differs from Specific Implementation Methods One through Eight in that, in steps three and four, the minimum error correction (MEC) scores for Mode One and Mode Two are calculated; the specific process is as follows:

[0300] To determine which ligation mode better reflects actual biological facts, this invention introduces Minimum Error Correction (MEC) as an evaluation function. The core idea of ​​the MEC principle is that a correct haplotype assembly scheme should minimize the number of bases inconsistent with the sequencing read (i.e., the number of errors requiring correction). The MEC scores for the two modes are calculated separately:

[0301] 1) Calculate the Minimum Error Correction (MEC) score for Mode 1; the specific process is as follows:

[0302] 11) Calculate the reading segment Source mismatch number The process is as follows:

[0303] exist Region and ( The number of inconsistent sites plus exist Region and ( The number of inconsistent sites;

[0304] 12) Calculate the reading segment Source Mismatch 2-2 The process is as follows:

[0305] exist Region and ( The number of inconsistent sites plus exist Region and ( The number of inconsistent sites;

[0306] 13) Reading paragraphs The minimum cost under mode is ;

[0307] 14) Accumulate the cost of all cross-block reads. The minimum error correction MEC score for Mode 1 is obtained.

[0308] 2) Calculate the Minimum Error Correction (MEC) score for Mode 2; the specific process is as follows:

[0309] 21) Calculate the reading segment Source mismatch number The process is as follows:

[0310] exist Region and ( The number of inconsistent sites plus exist Region and ( The number of inconsistent sites;

[0311] 22) Calculate the reading segment Source mismatch number The process is as follows:

[0312] exist Region and ( The number of inconsistent sites plus exist Region and ( The number of inconsistent sites;

[0313] 23) Reading paragraphs The minimum cost under Mode 2 is ;

[0314] 24) Accumulate the cost of all cross-block reads. The minimum error correction (MEC) score for mode two is obtained.

[0315] The other steps and parameters are the same as those in one of the specific implementation methods one to eight.

[0316] Specific Implementation Method Ten: This implementation method differs from Specific Implementation Methods One through Nine in that, in step three-five, the judgment is based on the minimum error correction MEC score of Mode One and the minimum error correction MEC score of Mode Two. and To establish a connection;

[0317] The specific process is as follows:

[0318] After obtaining the MEC scores for both modes, this invention makes a ligation decision by comparing the numerical differences between the two; to prevent erroneous ligations caused by sequencing noise, this invention introduces a ratio-based ( ) confidence assessment mechanism;

[0319] 1) Calculate the MEC ratio :

[0320]

[0321] The range of the ratio is [0, 1];

[0322] The score of the minimum error correction MEC for Mode 1;

[0323] The score of the Minimum Error Correction (MEC) for Mode 2;

[0324] 2) If MEC ratio A connection threshold greater than the preset threshold (set to 0.8 in a preferred embodiment of the invention) indicates that the error costs of the two connection modes are very similar, and the data provides insufficient evidence to distinguish which mode is correct (i.e., there is ambiguity). In this case, to ensure accuracy, the present invention selects... and Do not establish a connection, maintain and It is an independent disconnected state;

[0325] If the MEC ratio is less than or equal to the connectivity threshold, it indicates that one mode is significantly superior to the other; comparison and Size:

[0326] like If the value is smaller, then accept the hold mode; and merge, All The monomer type remains unchanged;

[0327] like If smaller, then accept the inverted mode; and Merge, but for All The single-type markers are flipped (1 becomes 2, 2 becomes 1) to achieve the correct physical connection.

[0328] The other steps and parameters are the same as those in any of the specific implementation methods one to nine.

[0329] Through the above process, this invention can progressively assemble and merge tens of thousands of local haplotype blocks across the entire genome based on long-read evidence, ultimately outputting chromosome-level haplotype typing results with high continuity and excellent N50 index. Experimental results show that, using the ligation strategy of this invention, the haplotype N50 can be increased to over 2.78 Mb on Oxford Nanopore sequencing data, which is significantly better than the results without global ligation.

[0330] This invention may have other embodiments. Without departing from the spirit and essence of this invention, those skilled in the art can make various corresponding changes and modifications according to this invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A whole-genome haplotype typing method for long-read sequencing data, characterized in that: The specific process of the method is as follows: Step 1: Receive sequencing reads and variant sites, and perform real-time multiple data filtering on the received sequencing reads and variant sites to obtain filtered sequencing reads and variant sites. Based on the filtered sequencing reads and variant sites, effective reads are obtained; Based on the effective reads, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site. Each valid read segment is compared with the template A sequence to obtain the corresponding template A alignment result; Each valid read segment is compared with the template B sequence to obtain the corresponding template B alignment result; Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each valid read segment and the template A sequence; Based on the alignment results of the corresponding template B, the edit distance 2 between each valid read segment and the template B sequence is calculated using the bit vector dynamic programming method; By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined; Step 2: Receive the selected SNP sites from Step 1 one by one. For newly added SNP sites... Perform dynamic grouping; Construct an SNP graph; Based on edge weights, the connection patterns between adjacent SNP sites within each encapsulated SNP group are determined, resulting in local monolithic blocks. Further optimization of the local single-unit blocks yields the final single-unit block; Step 3: Reorder and flatten the single-unit blocks output in Step 2 to obtain a queue to be connected, and determine the single-unit blocks based on the queue to be connected; For adjacent single-block pairs, the preceding block in the adjacent single-block pair is denoted as Block A, and the following block is denoted as Block B; Scan all sequencing reads aligned to the regions covered by Block A and Block B to identify all potential cross-block reads that span the end of Block A and the beginning of Block B. A strict screening strategy is applied to the identified potential cross-block read segments to obtain the cross-block read segments; For any two adjacent single-type blocks, Block A and Block B, there are only two logical possibilities for the single-type connection relationship between Block A and Block B; construct hypothetical models for each of the two possibilities. Calculate the minimum error correction scores for Mode 1 and Mode 2; Based on the minimum error correction MEC score of Mode 1 and the minimum error correction MEC score of Mode 2, determine whether Block A and Block B are connected. Continue until all single-unit blocks are connected, then output the fractal result.

2. The whole-genome haplotype typing method for long-read sequencing data according to claim 1, characterized in that: In step one, sequencing reads and variant sites are received, and real-time multiple data filtering is performed on the received sequencing reads and variant sites to obtain filtered sequencing reads and variant sites. Based on the filtered sequencing reads and variant sites, effective reads are obtained; Based on the effective reads, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site. Each valid read segment is compared with the template A sequence to obtain the corresponding template A alignment result; Each valid read segment is compared with the template B sequence to obtain the corresponding template B alignment result; Based on the alignment results of the corresponding template A, the edit distance 1 between each valid read segment and the template A sequence is calculated using the bit vector dynamic programming method; Based on the alignment results of the corresponding template B, the edit distance 2 between each valid read segment and the template B sequence is calculated using the bit vector dynamic programming method; By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined; The specific process is as follows: Step 11: The real-time variant detection software RVC outputs variant detection results for sequencing reads; Receive two types of concurrent data streams sorted by genome coordinates from the output of real-time variant detection software; The first type of data stream consists of sequencing reads; The second type of data stream consists of variant sites; The variant site stream includes the genomic coordinates of the variant, reference bases, and information on the variant bases; Step 1 and 2: Perform real-time multiple data filtering on the sequencing reads and variant sites received in real time in Step 1 to obtain the filtered sequencing reads and variant sites. Step 13: Based on the filtered sequencing reads and variant sites obtained in Step 12, obtain the effective reads; Step 14: Based on the effective reads from Step 13, obtain the window corresponding to each SNP site in the effective reads, and construct two target sequence templates based on the window corresponding to each SNP site; Step 15 Each valid read obtained in steps one and three is compared with the template A sequence to obtain the corresponding template A alignment result; Each valid read obtained in steps one and three is compared with the template B sequence to obtain the corresponding template B alignment result; Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each effective read segment obtained in steps one and three and the template A sequence; Based on the alignment results of the corresponding template B, the bit vector dynamic programming method is used to calculate the edit distance 2 between each effective read segment obtained in steps one and three and the template B sequence; By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined.

3. The whole-genome haplotype typing method for long-read sequencing data according to claim 2, characterized in that: In steps one and two, real-time multiple data filtering is performed on the sequencing reads and variant sites received in real time in step one to obtain filtered sequencing reads and variant sites. The specific process is as follows: Step 121: Perform real-time multiple data filtering on the sequencing reads received in Step 11 to obtain the filtered sequencing reads; the specific process is as follows: 1) For the sequencing reads received in real time in step one) The sequencing reads are preserved, and step 2) is executed. In the sequencing reads received in real time in step one Delete the sequencing reads; 2) Use SAMtools markdup or Picard MarkDuplicates tools to determine whether PCR repetitive sequences exist in the retained sequencing reads; If present, delete the PCR repeat sequences to ensure that only one repeat sequence is retained, and proceed to step 3). If no PCR repeat sequence is found, proceed to step 3). 3): Use the -F 0x900 parameter in SAMtools to determine whether there are sequences marked FLAG0x100 or FLAG 0x800 in the retained sequencing reads; If they exist, the sequences marked with FLAG 0x100 or FLAG 0x800 will be deleted to obtain the filtered sequencing reads; If it does not exist, the filtered sequencing reads are obtained; Step 122: Perform real-time multiple data filtering on the variant sites received in step 11 to obtain the filtered variant sites; The specific process is as follows: Read the QUAL and FILTER fields from the variant sites, retain the variant sites marked as PASS in the QUAL and FILTER fields that have a quality value higher than the preset threshold, and delete the remaining variant sites. At the same time, frequency information of variant sites is obtained, and extreme variant sites with a frequency of less than 0.2 or greater than 0.8 are removed.

4. The whole-genome haplotype typing method for long-read sequencing data according to claim 3, characterized in that: In step one and three, effective reads are obtained based on the filtered sequencing reads and variant sites obtained in step one and two. The specific process is as follows: Place the filtered sequencing reads obtained in steps one and two into the sequencing read buffer, and place the filtered variant sites obtained in steps one and two into the variant site buffer. When a new SNP site is received in step one, the genomic coordinates Pos of the new SNP site are recorded; Compare the termination coordinates of each sequencing read in the sequencing read buffer with the termination coordinates of the new SNP site; If the termination coordinate of a sequencing read in the buffer is less than the termination coordinate of a new SNP site, calculate the number of variant sites K in the corresponding sequencing read; if K is less than 2, the corresponding sequencing read is determined to be unable to provide valid haplotype linkage information, and the corresponding sequencing read is removed from the buffer and discarded; if K is greater than or equal to 2, the corresponding sequencing read is determined to be a valid read. If the termination coordinates of all sequencing reads in the buffer are greater than or equal to the termination coordinates of the new SNP site, no action is taken.

5. The whole-genome haplotype typing method for long-read sequencing data according to claim 4, characterized in that: In step one four, based on the effective read segment from step one three, a window corresponding to each SNP site in the effective read segment is obtained, and two target sequence templates are constructed based on the window corresponding to each SNP site. The specific process is as follows: Step 141: For each SNP site in the effective reads obtained in Step 13, take the coordinates of each SNP site on the reference genome as the center, and extract a fixed length of base sequence upstream and downstream of each center to obtain the window corresponding to each SNP site. SNP loci contain reference alleles and variant alleles; If the distance between a certain SNP site and the upstream endpoint of the reference chromosome genome is less than a fixed length, then the base sequence from the center to the endpoint is extracted from the coordinates of the SNP site on the reference genome, and a fixed length of base sequence is extracted from the center to the downstream, thus obtaining the window corresponding to each SNP site. If the distance between a certain SNP site and the downstream endpoint of the reference chromosome genome is less than a fixed length, then the base sequence from the center to the endpoint is extracted from the coordinates of the SNP site on the reference genome, and the fixed-length base sequence is extracted from the center to the upstream, thus obtaining the window corresponding to each SNP site. Step 142: Construct two target sequence templates based on the window corresponding to each SNP site: Template A sequence: The original sequence of the reference genome within the window obtained in step 141; The central base of the window obtained in step 141 is the reference allele; The template A sequence is a reference sequence template; Template B sequence: Replace the central base of the window obtained in step 141 with the variant allele (ALT), and keep the remaining sequences in the window consistent with the reference genome except for the center; The template B sequence is a variant sequence template.

6. The whole-genome haplotype typing method for long-read sequencing data according to claim 5, characterized in that: In step one five, each valid read segment obtained in step one three is compared with the template A sequence to obtain the corresponding template A alignment result; Each valid read obtained in steps one and three is compared with the template B sequence to obtain the corresponding template B alignment result; Based on the alignment results of the corresponding template A, the bit vector dynamic programming method is used to calculate the edit distance 1 between each effective read segment obtained in steps one and three and the template A sequence; Based on the alignment results of the corresponding template B, the bit vector dynamic programming method is used to calculate the edit distance 2 between each effective read segment obtained in steps one and three and the template B sequence; By comparing edit distance 1 and edit distance 2, the alleles supported by the read were determined; The specific process is as follows: Step 151: Align each valid read segment obtained in Step 13 with the template A sequence to obtain the alignment result of the corresponding template A; Each valid read obtained in steps one and three is compared with the template B sequence to obtain the corresponding template B alignment result; Step 152: Based on the alignment results of the corresponding template A obtained in Step 151, use the bit vector dynamic programming method to calculate the edit distance 1 between each effective read segment obtained in Step 13 and this template A sequence; Based on the alignment results of the corresponding template B obtained in step 151, the bit vector dynamic programming method is used to calculate the edit distance 2 between each effective read segment obtained in step 13 and this template B sequence; The specific process is as follows: 1) Initialize the current score ; make , Represents the target sequence The One location; make , This represents the query sequence obtained in steps one and three. The Middle One location; target sequence It can be either template A sequence or template B sequence; 2) Obtain the target sequence The Character matching information at each position ; indicates as: in, express The character matching vector; The query sequence obtained in step one or three Each valid read segment is ; like ,but =1, otherwise =0; Represents the target sequence The One location; Represents the query sequence The Middle One effective reading segment; express With query sequence Matching results; 3) Based on the first Character matching information at each position Calculate the vertical auxiliary vector ; indicates as: ; in, express and Bitwise OR operation; Represents the vertical negative difference vector; during initialization, let... It is a vector consisting entirely of zeros; 4) Calculate the horizontal auxiliary vector based on the vertical auxiliary vector; expressed as: ; in, express and Bitwise OR operation; Represents the vertical positive difference vector; during initialization, let... It is a vector consisting entirely of 1s; Represents a horizontal auxiliary vector; This represents the bitwise AND operator. This represents the bitwise XOR operator. 5) Calculate the horizontal positive difference vector based on the horizontal auxiliary vector. ; indicates as: ; ； in, express and Bitwise OR operation; Indicates bitwise inversion; 6) Update the edit distance in real time using the highest bit of the horizontal difference vector; expressed as: If the horizontal positive difference vector If the highest bit is 1, then the edit distance is... ; If the horizontal negative difference vector If the highest bit is 1, then the edit distance is... ; The highest position is the first one on the left; 7) Update the vertical difference vector based on the real-time updated edit distance; expressed as: in, express Shift one position to the left; express Shift one position to the left; 8) Order Repeat steps 2)-7) until the target sequence has been traversed. After considering all positions, the final edit distance Score is the query sequence obtained in step one or three. Each valid read segment Edit distance between the template A sequence and the template B sequence; And select the central base corresponding to the template A sequence or template B sequence with a smaller score as the base type supported by the read at the SNP site; 9) Order Repeat steps 2)-8) until the query sequence obtained in steps 1 and 3 has been traversed. All valid reading segments ; Step 153: By comparing edit distance 1 and edit distance 2, determine the alleles supported by the read segment. The process is as follows: Set the edit distance threshold. Alleles corresponding to edit distances less than the edit distance threshold of 1 are considered to be supported by reads; Alleles corresponding to edit distances greater than or equal to the edit distance threshold of 1 are considered as reads that are not supported. Alleles corresponding to edit distances less than the edit distance threshold of 2 are considered to be read-supported. Alleles corresponding to edit distances greater than or equal to the edit distance threshold of 2 are considered as reads that are not supported.

7. The whole-genome haplotype typing method for long-read sequencing data according to claim 6, characterized in that: In step two, the filtered SNP sites from step one are received one by one, and for newly incoming SNP sites... Perform dynamic grouping; construct SNP graph; determine the connection pattern between adjacent SNP sites within each encapsulated SNP group based on edge weights to obtain local monolithic blocks; further optimize the local monolithic blocks to obtain the final monolithic blocks; The specific process is as follows: Step 21: Receive the filtered SNP sites from Step 15 one by one. For newly added SNP sites... Dynamic grouping is performed; the specific process is as follows: Receive the selected SNP sites one by one from step one five, and for newly added SNP sites... Perform the following truncation criteria: (1) Physical length limit: The current SNP group is added to the new inflow Afterwards, the total physical length of the current group exceeds a preset threshold. ; (2) Number of sites limited: New inflows into the current SNP group Subsequently, the total number of SNPs in the current group exceeds a preset threshold. ; (3) Adjacent spacing restrictions: newly flowing in Compared with the last SNP site in the current SNP group The physical distance between them exceeds a preset threshold ; If any of the above truncation conditions are met, the "current SNP group" is determined to be completed. The current SNP group is immediately encapsulated as a task unit to be processed and pushed into the graph construction engine. Simultaneously, a new group is initialized to accommodate newly incoming SNPs. ; Real-time monitoring of the minimum starting mapping coordinates of all unprocessed valid read segments in the buffer of step one and three, denoted as... For the SNP group that has been encapsulated and is to be calculated, the termination physical coordinates are... Continuous judgment and Size, if Greater than or equal to Continue monitoring; if Less than This means that the starting positions of existing valid read segments and subsequent new inflow read segments in the buffer are all located after the SNP group, and there is no physical possibility of overwriting them; at this point, it is determined that the collection of read segments of the encapsulated SNP group to be calculated has been completed. Step 22: Construct the SNP diagram; the specific process is as follows: Node definition in a directed acyclic graph: For each encapsulated SNP group obtained in step two, each SNP site within that group. Create two nodes in the graph and label them as follows: and ; in, The reference allele representing the current locus; The variant alleles representing the current locus; Edge definition in a directed acyclic graph: Directed edges are established only between adjacent SNP sites; for adjacent sites... and There are four possible connecting edges: side :connect and This indicates that both loci are reference alleles; side :connect and This indicates that both loci are variant alleles; side :connect and , indicating the reference allele and the variant allele; side :connect and This indicates that the variant allele is linked to the reference allele. Weight calculation: Based on the SNP sites within each encapsulated SNP group obtained in step two-one. and The read segment data is used to count the number of read segments that support each type of connection relationship, and the count is assigned to the corresponding edge weight. Steps two and three: Based on edge weights, determine the connection patterns between adjacent SNP sites within each encapsulated SNP group to obtain the connected SNP groups; the specific process is as follows: Step 231: For any adjacent SNP pairs within each encapsulated SNP group... This examines two mutually exclusive single-type phase combination schemes; the specific process is as follows: Option 1: Assume adjacent SNP pairs The two monomer types are 0-0 and 1-1, where 0-0 represents... The two haplotypes are the reference allele and the reference allele, respectively, and 1-1 indicates that... The two haplotypes are the variant allele and the variant allele, respectively; the formula for calculating the total support of Scheme 1 is: ； Option 2: Assume adjacent SNP pairs The two monomer types are 0-1 and 1-0, where 0-1 represents... The two monomer types are the reference base 0 and the variant base 1; 1-0 indicates The two monomer types are the variant base 1 and the reference base 0, respectively; the formula for calculating the total support of Scheme 2 is: ； in, This indicates the total support for Option 1; express and Weights between them; express and Weights between them; This indicates the total support for Option 2; express and Weights between them; express and Weights between them; Indicates the reference allele at the current locus; Indicates the variant allele at the current locus; Indicates site Reference allele; Indicates site Variable alleles; Step Two Three Two, Set minimum total support threshold ; The total support of Option 1 Total support for Option 2 The total support is obtained by adding them together. ; Total support With minimum total support threshold Compare; like If so, the evidence is considered insufficient, and adjacent SNP pairs are considered... and No connection is established; like Then proceed to step two or three. Step Two Three Three Set support ratio threshold ; Based on the total support of Scheme 1 Total support for Option 2 Calculate the support ratio : if If the correct connection method cannot be determined, then adjacent SNP pairs are considered to be... and No connection is established; if The total support for option one is selected. Total support for Option 2 Adjacent SNP pairs in schemes with high overall support and Establish a connection; Steps 2 and 4: The connections obtained in Steps 2 and 3... The group is further optimized to obtain the final single-unit block; the specific process is as follows: After the connection obtained in steps two and three For each monolithic block in the group, rescan all read segments that cover the current monolithic block, and calculate the consistency score between each read segment and the monolithic block; the process is as follows: For each read segment, calculate the relationship between each read segment and... The consistent number of sites is obtained Consistency score ; Statistics on each reading segment and The consistent number of sites is obtained Consistency score ; Calculate support confidence : Read segments with a confidence level greater than a threshold are selected as high-consistency read segments; in, Indicates segment coverage quantity; After selecting a set of highly consistent read segments, these read segments are used to process each segment within a single-unit block. The site undergoes secondary verification; the process is as follows: If a certain single-type block If the phase of a site conflicts with the phase of a 75% high-consistency read covering the current site, the phase of the current high-consistency read site is flipped and corrected; if the conflict cannot be resolved, the current high-consistency read site is removed from the haplogroup, resulting in the final haplogroup. If a certain single-type block If the phase of a site does not conflict with the phase of a 75% high-consistency read covering the current site, the phase is retained.

8. A whole-genome haplotype typing method for long-read sequencing data according to claim 7, characterized in that: In step three, the single-unit blocks output in step two are reordered and flattened sequentially to obtain a queue to be connected, and the single-unit blocks are determined based on the queue to be connected. For adjacent single-unit block pairs, the preceding block in the adjacent single-unit block pair is denoted as... The last block is denoted as ; Scan all matches and Sequencing reads covering the covered region identified all reads that spanned [location missing]. End position and Potential cross-block read segment at the starting position; A strict screening strategy is applied to the identified potential cross-block read segments to obtain the cross-block read segments; For any two adjacent single-type blocks and , and The single-unit connection relationship between them has only two logical possibilities; hypothetical models are constructed for each of the two possibilities; Calculate the minimum error correction scores for Mode 1 and Mode 2; Based on the minimum error correction MEC scores of Mode 1 and Mode 2, determine... and To establish a connection; Continue until all single-unit blocks are connected, then output the fractal result; The specific process is as follows: Step 31: Reorder and flatten the single-unit blocks output from Step 2 to obtain a queue to be connected, and determine the single-unit blocks based on the queue to be connected; the specific process is as follows: (1) Reorder the single-unit blocks output in step two: buffer the out-of-order received blocks. Group, until received The group number equals the expected group number; (2) Flatten the reordered haplotype blocks: Once the final haplotype blocks obtained in step two and four are obtained in sequence, the final haplotype blocks obtained in step two and four are immediately pushed into the queue to be joined in genomic position order. (3) Continuously pop new blocks from the queue to be connected, and treat the global tail block and the new block as single-type blocks; Step 3.2: For adjacent pairs of single-unit blocks, the previous block in the adjacent pair is denoted as... The last block is denoted as ; Scan all matches and Sequencing reads covering the covered region identified all reads that spanned [location missing]. End position and The potential cross-block read segment at the start position; A strict screening strategy is applied to the identified potential cross-block read segments to obtain the cross-block read segments; The screening strategy is as follows: (1) Effective coverage verification: Requires potential cross-block read segments in At least one already classified site, and in It also covers at least one already classified type. site; Potential cross-block reads that do not meet the requirements cannot provide phase linkage information and are therefore eliminated. (2) Activity level check: Calculate the overlap length of potential cross-block reads within each monolithic block; The overlap length between the potential cross-block read base sequence length and the monomer block base sequence length must be at least one-fifth of the total length of the monomer block; Potential cross-block read segments that do not meet the requirements are removed; Step 3: For any two adjacent single-unit blocks and , and The single-unit connection relationship between them has only two logical possibilities; hypothetical models are constructed for each of the two possibilities; the specific process is as follows: Mode 1: Monotype 1 and Haplotype 1 belongs to the same chromosome and forms a haplotype. ; Monotype 2 and Haplotype 2 belongs to another chromosome, forming a haplotype. ; Mode 2: Monotype 1 and The monomers are connected to form monomers. ; Monotype 2 and The monomers 1 are connected to form monomers. ; Steps 3 and 4: Calculate the minimum error correction scores for Mode 1 and Mode 2; Step 3.5: Based on the minimum error correction MEC score of Mode 1 and the minimum error correction MEC score of Mode 2, determine... and To establish a connection; Step 36: Repeat steps 31 to 35 to complete the connection of all single-unit blocks and output the fractal result.

9. A whole-genome haplotype typing method for long-read sequencing data according to claim 8, characterized in that: In steps three and four, the minimum error correction scores for Mode 1 and Mode 2 are calculated; the specific process is as follows: 1) Calculate the Minimum Error Correction (MEC) score for Mode 1; the specific process is as follows: 11) Calculate the reading segment Source mismatch number The process is as follows: exist Region and The number of inconsistent sites plus exist Region and The number of inconsistent sites; 12) Calculate the reading segment Source Mismatch 2-2 The process is as follows: exist Region and The number of inconsistent sites plus exist Region and The number of inconsistent sites; 13) Reading paragraphs The minimum cost under mode is ; 14) Accumulate the cost of all cross-block reads. The minimum error correction MEC score for Mode 1 is obtained. 2) Calculate the Minimum Error Correction (MEC) score for Mode 2; the specific process is as follows: 21) Calculate the reading segment Source mismatch number The process is as follows: exist Region and The number of inconsistent sites plus exist Region and The number of inconsistent sites; 22) Calculate the reading segment Source mismatch number The process is as follows: exist Region and The number of inconsistent sites plus exist Region and The number of inconsistent sites; 23) Reading paragraphs The minimum cost under Mode 2 is ; 24) Accumulate the cost of all cross-block reads. The minimum error correction (MEC) score for mode two is obtained.

10. A whole-genome haplotype typing method for long-read sequencing data according to claim 9, characterized in that: In step three, the minimum error correction MEC score of mode one and the minimum error correction MEC score of mode two are used to determine... and To establish a connection; The specific process is as follows: 1) Calculate the MEC ratio : The range of the ratio is [0, 1]; The score of the minimum error correction MEC for Mode 1; The score of the Minimum Error Correction (MEC) for Mode 2; 2) If MEC ratio If the connection threshold is greater than the preset threshold, select... and Do not establish a connection, maintain and It is an independent disconnected state; If the MEC ratio is less than or equal to the connection threshold, compare and Size: like If the value is smaller, then accept the hold mode; and merge, All The monomer type remains unchanged; like If smaller, then accept the inverted mode; and Merge, but for All The single-type markers are flipped.

Images