Method for analyzing biological sample data

Abstract

The invention provides a biological sample data analysis method. The method is applied to the technical field of data processing and comprises the steps that a generic genome graph and sequencing read length data of a biological sample are obtained, the generic genome graph is dynamically divided into a plurality of sub-graphs through hierarchical graph segmentation, and the segmentation density is automatically adjusted according to variation complexity of a genome region; performing topology-aware dynamic programming comparison on sequencing read length data in each sub-graph, and simultaneously processing sequence positions, node positions and path selection through a three-dimensional scoring matrix to generate corresponding sub-graph comparison results; combining all subgraph comparison results based on a minimum cost flow algorithm, and constructing a globally consistent genome comparison path; and detecting structural variation according to the global comparison path, and outputting a variation data set containing accurate breakpoint coordinates. In this way, the accuracy of biological sample data analysis can be improved.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a method for analyzing biological sample data. Background Technology

[0002] With the rapid development of high-throughput sequencing technology, obtaining whole-genome sequencing data of individuals or populations has become increasingly convenient. However, accurately interpreting this massive amount of data, especially precisely detecting complex genomic structural variations (such as insertions, deletions, inversions, and translocations), remains a major challenge for biomedical research and clinical diagnosis. Structural variations play a crucial role in genome evolution, the occurrence of complex diseases, and individual phenotypic differences; their accurate detection is essential for understanding life mechanisms and achieving precision medicine. Current mainstream biological sample data analysis methods are mainly based on linear reference genomes. However, a single linear reference cannot fully represent the genetic diversity within a population, especially in regions with high-density single nucleotide polymorphisms (SNPs), insertions, deletions, and complex structural variations. This mismatch leads to inaccurate alignment of sequencing reads in variant regions, resulting in false positives or false negatives in the detection of structural variations, severely reducing the accuracy of data analysis. To address these issues, pangenome graph models have emerged. It integrates the genome sequences of multiple individuals, incorporating known variants (especially common SNPs, Indels, and SVs in the population) as branching paths into the graph, theoretically providing more comprehensive reference information and improving the fidelity of sequencing read alignment. However, existing technologies still face significant challenges in practical applications of pan-genome graphs. Pan-genome graphs typically contain millions or even tens of millions of nodes and edges, forming extremely complex topologies.Performing precise full-graph dynamic programming alignments on such complex graph structures incurs enormous computational and memory overhead, especially when processing massive sequencing reads, resulting in low efficiency and difficulty in scaling to large-scale population analyses. Some optimization schemes attempt to statically partition large graphs into several subgraphs for distributed computation; however, current methods often employ static, fixed-size partitioning strategies, lacking the ability to account for the variational complexity of different genomic regions. The perception of complexity leads to unnecessary fine-grained segmentation in sparsely variable regions, increasing computational and communication overhead. In highly complex, variable-dense regions (such as repeat regions, inversions, tandem repeats, copy number variations, etc.), the boundaries of static segmentation often precisely intersect these complex regions. This prevents individual subgraphs from independently modeling complete local paths and variation options, forcing analysis tools to perform abrupt truncation at subgraph boundaries. This disrupts the topological integrity of local regions, obliterating the continuous path representation that complex variation patterns should have within those regions, inevitably leading to the loss of fine-structure information and inaccurate alignment results. Efficiently and accurately integrating the alignment results of segmented subgraphs globally is another key challenge. Existing methods (such as simple splicing or local path selection) cannot guarantee the topological consistency and path optimality of the merged path at the whole-graph level. Locally optimal paths within subgraphs may conflict or become discontinuous after global splicing, failing to form a smooth, unified global genome alignment path. AlignmentPath directly undermines the overall logical structural foundation upon which subsequent structural variation analysis relies. In repetitive regions, highly polymorphic regions, or regions with complex nested SVs, existing graph-based alignment methods are still limited by factors such as computational complexity, difficulties in static segmentation and integration, etc. Read length alignment in these regions is prone to errors, and path ambiguity is severe, significantly reducing the sensitivity and resolution for detecting small fragment insertions and deletions and for accurate breakpoint localization. The variant breakpoint coordinates output by existing tools are often unclear or significantly biased, hindering in-depth research on the biological mechanisms of structural variation and its clinical application.

[0003] In summary, existing pan-genome graph-based analysis methods suffer from insufficient overall accuracy when processing large-scale biological sample data (especially when it is necessary to accurately detect structural variations in complex regions) due to their high computational complexity (especially the difficulty of global dynamic programming of complex paths), the use of static and rigid segmentation strategies that lead to the fragmentation of information in complex variation regions, and the lack of effective mechanisms to ensure the consistency of the globally optimal path after the integration of subgraph results. Summary of the Invention

[0004] This invention provides a method for analyzing biological sample data. The method includes:

[0005] S1: Obtain pangenome map and sequencing read data of biological samples. Dynamically divide the pangenome map into multiple sub-maps through hierarchical map segmentation. The segmentation density is automatically adjusted according to the variation complexity of the genomic region.

[0006] S2 performs topology-aware dynamic programming alignment of sequencing read data within each subgraph, and simultaneously processes sequence position, node position and path selection through a three-dimensional scoring matrix to generate corresponding subgraph alignment results;

[0007] S3, based on the minimum cost flow algorithm, merges all subgraph alignment results to construct a globally consistent genome alignment path;

[0008] S4 detects structural variations based on global path comparison and outputs a variation dataset containing precise breakpoint coordinates.

[0009] Furthermore, after S4 detects structural variations, the process includes: if the detected variation region contains low-quality bases with base quality values ​​below the base quality value threshold, a secondary validation is performed. Specifically, this involves: locating all sequencing reads within a 50bp range on both sides of the variation breakpoint; calculating the number of reads supported by the breakpoint using a local re-alignment algorithm, retaining the variation when the number of supported reads is ≥3; predicting the functional impact of the retained variations and labeling their clinical significance level using the ClinVar database; the secondary validation process is executed in parallel using GPUs, with each variation assigned an independent computation thread; if the variation is located in a known pathogenic gene region, an additional read length consistency check is initiated, utilizing the PacBioHiFi read length crossing breakpoint characteristic to verify the variation's authenticity.

[0010] Furthermore, the S1 layered graph segmentation process includes: when the variation density is >0.3, the segmentation level is automatically increased: the top layer segmentation generates 5-10 macro sub-graphs for quickly locating variation regions; the middle layer segmentation generates micro sub-graphs within the top layer sub-graphs to focus on high variation density regions; and the bottom layer segmentation generates nano sub-graphs for complex structural variation regions.

[0011] Furthermore, the layered graph segmentation process of S1 satisfies the following constraints: the number of connecting edges between subgraphs is less than 15% of the internal edges of the subgraph; when a branch region with a node out-degree greater than the node out-degree threshold is detected, the integrity of the branch path is preserved first, and path cutting is avoided by adding virtual nodes.

[0012] Furthermore, the three-dimensional scoring matrix of S2 includes topological consistency verification: if the current node v b With the predecessor node v b'If the paths are not consecutively adjacent on the topological path, a target penalty score is applied; when the score difference between multiple candidate paths is less than 2 points, the high-frequency path is selected first; for reads that cross the subgraph boundary, if the score difference between the alignments in adjacent subgraphs is greater than 10 points, a 100bp sequence window is extracted on both sides of the boundary; an extended subgraph is constructed; and a second dynamic programming alignment is performed in the extended subgraph.

[0013] Furthermore, following S3, it also includes: real-time monitoring of the time T consumed in comparing each subgraph. i When T max / T avg >2,T max For the maximum time consumption, T avg To ensure that the average time taken is less than the GPU utilization threshold, the subgraph with the highest time taken is further divided into several new subgraphs, which are then allocated to idle GPU threads. The GPU memory resources of the subgraphs that have been compared are released.

[0014] Furthermore, it also includes: correcting the matching score for low-quality bases whose base quality value is below the base quality value threshold: σ'=σ×(1-10 (-Q / 10) ), where σ is the original matching score and Q is the base quality value; if the frequency of occurrence of the base is greater than the frequency threshold, the bonus score is increased; when Q < 20 for 3 consecutive bases D, the topological constraints are relaxed for re-alignment; k-mer seed expansion is used to re-anchor the position; after re-alignment, candidate results with an alignment score improvement of < 5 are deleted.

[0015] Furthermore, it also includes: if the read length is >10kb, then before S2, extract the first and last 1kb of k-mers, k=15, for topology-guided seed localization; locate the topological path of the k-mer in the pan-genome map; when the distance between the first and last paths is < read length × k, 1 < k < 1.5, directly generate candidate sub-maps; for regions with an error rate greater than the error rate threshold, increase the k-mer step size to 20bp.

[0016] Furthermore, S4 includes: statistically analyzing the distribution of breakpoint locations supporting read lengths to calculate breakpoint confidence intervals; determining breakpoint intervals using the 95th percentile; extracting all read lengths within the breakpoint interval for local assembly refinement; performing DeBruijn plot assembly; aligning the assembled sequence to the reference plot to determine precise breakpoint coordinates; and additionally checking sequence micro-homology when the breakpoint is located in a repetitive region.

[0017] Furthermore, it also includes: technical noise filtering of the final variant dataset, specifically: removing variant datasets with chain bias > 80%.

[0018] This invention utilizes a pan-genome map and adjusts the segmentation density according to variation complexity, subdividing high-variance regions and coarsely segmenting low-variance regions to adapt to the alignment needs of different regions and avoid alignment bias caused by unreasonable region division. By combining a three-dimensional scoring matrix of sequence position, node position, and path selection, it comprehensively assesses the alignment probability, reduces local misjudgments based solely on sequence similarity, and improves the precision of alignment within sub-maps. The minimum cost flow algorithm optimizes the merging of sub-map alignment results, balancing local and global costs, avoiding splicing errors, and ensuring the coherence and consistency of the whole-genome alignment path. By locating structural variations based on globally consistent alignment paths and accurately determining breakpoint coordinates using path information, it reduces breakpoint offsets or missed detections caused by alignment inconsistencies, thereby improving the accuracy of biological sample data analysis.

[0019] It should be understood that the description in the Summary of the Invention is not intended to limit the key or essential features of the embodiments of the present invention, nor is it intended to restrict the scope of the invention. Other features of the invention will become readily apparent from the following description. Attached Figure Description

[0020] The above and other features, advantages, and aspects of the various embodiments of the present invention will become more apparent from the accompanying drawings and the following detailed description. The drawings are provided for a better understanding of the invention and are not intended to limit the scope of the invention. In the drawings, the same or similar reference numerals denote the same or similar elements, wherein:

[0021] Figure 1 A flowchart illustrating a method for analyzing biological sample data according to an embodiment of the present invention is shown. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. 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.

[0023] Furthermore, the term "and / or" in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0024] Figure 1 A flowchart of a method for analyzing biological sample data according to an embodiment of the present invention is shown, the method comprising:

[0025] S1: Obtain pangenome map and sequencing read data of biological samples. Dynamically divide the pangenome map into multiple sub-maps through hierarchical map segmentation. The segmentation density is automatically adjusted according to the variation complexity of the genomic region.

[0026] S2 performs topology-aware dynamic programming alignment of sequencing read data within each subgraph, and simultaneously processes sequence position, node position and path selection through a three-dimensional scoring matrix to generate corresponding subgraph alignment results;

[0027] S3, based on the minimum cost flow algorithm, merges all subgraph alignment results to construct a globally consistent genome alignment path;

[0028] S4 detects structural variations based on global path comparison and outputs a variation dataset containing precise breakpoint coordinates.

[0029] According to embodiments of the present invention, a pan-genome map is used, and the segmentation density is adjusted according to the complexity of variation. This allows for the subdivision of highly variable regions and the coarse division of low-variable regions, adapting to the alignment needs of different regions and avoiding alignment bias caused by unreasonable region division. By combining a three-dimensional scoring matrix of sequence position, node position, and path selection, the alignment probability is comprehensively evaluated, reducing local misjudgments based solely on sequence similarity and improving the precision of alignment within sub-maps. The minimum cost flow algorithm optimizes the merging of sub-map alignment results, balancing local and global costs, avoiding splicing errors, and ensuring the coherence and consistency of the whole-genome alignment path. By locating structural variations based on globally consistent alignment paths and accurately determining breakpoint coordinates using path information, the breakpoint offset or missed detection caused by alignment inconsistencies is reduced, thereby improving the accuracy of biological sample data analysis.

[0030] In some embodiments, after detecting structural variations in S4, the process further includes: if the detected variation region contains low-quality bases (i.e., base quality value Q < 20) or a coverage depth < 5× (i.e., the variation region is detected less than 5 times), then a secondary verification is performed. Specifically, this involves: locating all sequencing reads within a 50bp range on both sides of the variation breakpoint; calculating the number of reads supported by the breakpoint using a local re-alignment algorithm; retaining the variation when the number of supported reads is ≥ 3; predicting the functional impact of the retained variation and labeling its clinical significance level using the ClinVar database; wherein, the secondary verification process is executed in parallel using GPU, with each variation allocated an independent computation thread; if the variation is located in a known pathogenic gene region (such as BRCA1, TP53), then an additional read length consistency check is initiated, utilizing the PacBioHiFi read length crossing breakpoint characteristic to verify the authenticity of the variation. According to embodiments of the present invention, by initiating secondary validation on unreliable regions with low-quality bases and low coverage depth, false positive variants caused by sequencing noise or insufficient data are eliminated; by focusing on the critical 50bp region on both sides of the variant breakpoint, the number of supporting reads is calculated using a local re-alignment algorithm to avoid interference from global alignment and ensure that the variant has sufficient experimental data support; by combining functional impact prediction and ClinVar database annotation, true variants are screened from the perspective of biological function and clinical significance, eliminating interference from non-functional or benign variants, thereby improving the accuracy of biological sample data analysis.

[0031] For example, low-quality bases are defined as Phred scores < 20 (error probability > 1%), and the coverage depth threshold is set to 5× (the number of effective sequencing times for this site in whole-genome sequencing data < 5). Using the coordinates of the mutation breakpoint ± 50 bp as a window (GRCh38 reference coordinate system), all short reads within this 100 bp interval are extracted. Smith-Waterman local re-alignment is used (match score +2, mismatch penalty -4, gap opening penalty -6, extension penalty -1), and the number of reads with alignment scores > 15 that cross the breakpoint is counted, with a threshold of ≥ 3. For mutations within chr17:41,196,312-41,277,500 (BRCA1 locus) or chr17:7,577,984-7,591,939 (TP53 locus), at least two PacBioHiFi reads (≥ 10 kb, accuracy ≥ 99.9%) are required to support the breakpoint.

[0032] In some embodiments, the hierarchical graph segmentation process of S1 includes: when the variation density is >0.3, automatically increasing the segmentation level: top-level segmentation generates 5-10 macro-level sub-graphs (size 1-5Mb) for quickly locating variation regions; middle-level segmentation generates micro-level sub-graphs (size 200-500kb) within the top-level sub-graphs to focus on high variation density regions; bottom-level segmentation generates nano-level sub-graphs (size 20-50kb) for complex structural variation regions (such as repeating sequences and inversions). According to embodiments of the present invention, by dynamically adjusting the segmentation level through "variation density," resources are saved by using coarse analysis in regions with sparse variations, while more refined segmentation is applied to complex / high variation density regions. By concentrating computational resources and complex algorithms on known high variation density regions and regions with complex structural variations, potential noise and errors caused by unnecessary deep analysis in simple regions are avoided. The resolution of local sequence and variation identification is greatly improved by generating smaller units in complex regions, enabling more accurate detection, localization, and characterization of variations that are closely adjacent or embedded in complex genomic backgrounds. By isolating complex regions that are prone to alignment errors within smaller subgraphs for analysis, the propagation and accumulation of errors such as long-range alignment misjudgments and chimeric assembly are significantly reduced. By segmenting the massive genome into subgraphs of appropriate size and varying problem complexity, the analysis algorithm can more accurately optimize and adjust parameters for the characteristics of each subgraph, reducing the global bias caused by a "one-size-fits-all" strategy, thereby improving the accuracy of biological sample data analysis.

[0033] For example, when the mutation density of the chromosomal region chr2:215,500,000-215,800,000 (300kb in length) is detected to be 0.35 (i.e., an average of 0.35 mutations per kb), hierarchical map segmentation is performed. Specifically, during top-level segmentation, the entire pan-genome map (e.g., 3.2Gb in size) is divided into 8 macroscopic submaps, each submap ranging in size from 400-500Mb (e.g., sub-submaps). Figure 1 :chr1:1-500,000,000; Subgraph 2:chr2:1-480,000,000, etc.); During middle-level segmentation, three micrographs are generated within the region of high variation density in the top-level subgraph (i.e., chr2:215,500,000-215,800,000). Figure 1:chr2:215,500,000-215,650,000 (150kb), micrograph 2:chr2:215,650,000-215,750,000 (100kb), micrograph 3:chr2:215,750,000-215,800,000 (50kb); During bottom-level segmentation, complex data was detected within micrograph 2. Structural variations (such as the inverted region chr2:215,700,000-215,720,000) are further segmented into two nanomaps: nanomap A: chr2:215,700,000-215,710,000 (10kb) and nanomap B: chr2:215,710,000-215,720,000 (10kb).

[0034] In some embodiments, the hierarchical graph segmentation process of S1 satisfies the following constraints: the number of connecting edges between subgraphs is less than 15% of the internal edges of the subgraph; when a high branch region (node ​​out-degree ≥ 4) is detected, the integrity of the branch path is preserved first, and path cutting is avoided by adding virtual nodes.

[0035] In some embodiments, the three-dimensional scoring matrix of S2 includes a topology consistency check: if the current node v b With the predecessor node v b' If the elements are not consecutively adjacent on the topological path, a penalty score δ is applied. path =-5; When the score difference among multiple candidate paths is <2 points, the high-frequency path (path with an occurrence frequency >60%) is selected first; For reads that cross subgraph boundaries, if the score difference in adjacent subgraphs is >10 points, a 100bp sequence window is extracted on both sides of the boundary; an extended subgraph (containing the original adjacent subgraphs and connecting paths) is constructed; secondary dynamic programming alignment is performed in the extended subgraph. According to the embodiments of the present invention, path continuity is ensured by penalizing non-contiguous adjacent nodes, reducing false alignments; the selection of reliable paths is enhanced by prioritizing high-frequency paths using consensus frequency, reducing the impact of random errors; secondary dynamic programming alignment by constructing an extended subgraph provides a more comprehensive sequence context, improves boundary region alignment, solves the bias caused by subgraph partitioning, and thus improves the accuracy of biological sample data analysis.

[0036] For example, when the comparison read length reaches node v b When (coordinates chr3:50, 100, 002), if its predecessor node v b' (chr3:50,099,997) requires passing through two intermediate nodes (not directly adjacent) in the pan-genome topological pathway, so δ is applied. pathThe penalty of -5 reduces the original alignment score from 32 to 27. In the chr5:180,200,500-180,200,800 region, candidate path A scores 45, and path B scores 43.5 (difference 1.5 < 2). If path A appears 620 times in 1000 reads (frequency 62% > 60%), then path A is preferred. Read R1 has an alignment score of 38 in subgraph G1 (chr7:55,000,000-55,200,000) and in the adjacent subgraph G2... If the score for (chr7:55, 200, 001-55, 400, 000) is 25 (difference 13 > 10), then the sequences of boundary region G1 (chr7:55, 199, 901-55, 200, 000) and boundary region G2 (chr7:55, 200, 001-55, 200, 100) are extracted, and an extended subgraph G1' (range chr7:55, 199, 901-55, 200, 100) containing two subgraphs and connecting paths is constructed. After re-alignment, the score is corrected to 41.

[0037] In some embodiments, after S3, the method further includes: real-time monitoring of the comparison time T of each subgraph. i When T appears max / T avg If the time taken is greater than 2 (maximum time taken > twice the average time taken) and the GPU utilization is less than 60%, then the most time-consuming subgraph is further divided into 2-4 nano-subgraphs (i.e., new subgraphs). These new subgraphs are then allocated to idle GPU threads, and the GPU memory resources of the subgraphs that have already been compared are released. According to this embodiment of the invention, by dividing the most time-consuming task when GPU utilization is low, idle GPU computing resources can be fully utilized. The newly generated subgraphs are allocated to idle threads, greatly increasing the system's parallel processing capability and avoiding the idleness of valuable computing resources. By promptly releasing the GPU memory resources occupied by the subgraphs that have already been compared, space is freed up for tasks that need to process newly divided subgraphs, reducing the possibility of program crashes, data loss, or incorrect results due to GPU memory exhaustion during high-load computing, thereby improving the accuracy of biological sample data analysis.

[0038] For example, during the execution of step S3, the comparison time (in seconds) of the 8 subgraphs is monitored in real time and recorded as T = [120, 95, 110, 105, 480, 115, 102, 98], and the maximum time T is calculated. max = 480 seconds (sub-figure 5), average time T avg = 153.125 seconds, ratio T max / T avg=3.13>2, the current GPU utilization monitoring value is 55%<60%; the time-consuming sub-image 5 (original size 450kb) is divided into 3 nano-sub-images: new sub-image 5a (150kb), new sub-image 5b (150kb), and new sub-image 5c (150kb). 3 idle GPU threads were detected (out of a total of 8 threads, 5 currently in use). The new sub-images are allocated to the idle threads, and the sub-images that have completed the comparison are released. Figure 1 The video memory resources of subgraph 3 (each occupying 1.2GB of video memory); the time taken for new subgraphs 5a / b / c are [155, 160, 158] seconds respectively, and the updated T max = 160 seconds, T avg =124.7 seconds, the ratio dropped to 1.28 < 2, and the GPU utilization increased to 82%.

[0039] In some embodiments, the method further includes: correcting the matching score for low-quality bases (base quality value Q < 30): σ' = σ × (1 - 10) (-Q / 10) ), where σ is the original matching score; if the base is located in a high-frequency path (occurrence frequency > 70%), the bonus score η = +2 is increased; when Q < 20 for three consecutive bases D, the topological constraints are temporarily relaxed (τ penalty is halved) for re-alignment; k-mer seed expansion is used to re-anchor the position; after re-alignment, candidate results with an alignment score improvement < 5 are deleted. According to embodiments of the present invention, by correcting the matching score of low-quality bases, the negative impact of low-quality bases on the matching score is reduced, and noise introduced by mismatches is avoided, thereby reducing false positives; by increasing the bonus score for high-frequency paths, common variants or sequence paths are preferentially selected, and prior knowledge is used to strengthen the real signal, reducing false negatives caused by rare noise, thereby improving the accuracy of biological sample data analysis.

[0040] In some embodiments, the method further includes: if the read length is >10kb, extracting the first and last 1kb k-mers (k=15) before S2 for topology-guided seed localization; locating the topological path of the k-mers in the pan-genome map; directly generating a candidate sub-map when the distance between the first and last paths is <read length × 1.2; increasing the k-mer step size to 20bp for high error rate regions (e.g., GC > 70%); and enabling fuzzy matching (allowing 1-2 mismatches) for low complexity regions. According to embodiments of the present invention, by extracting the first and last 1kb k-mers, error-prone regions in the middle of the read are avoided, ensuring that the localization seed is generated in a highly reliable region; by matching the topological path of the k-mers in the pan-genome map, structural information about the connection relationships between sequences is introduced, significantly reducing erroneous localization caused by repetitive sequences; and by generating a sub-map when the distance between the first and last paths is <read length × 1.2, unreasonable mapping results with contradictory first and last localizations are filtered out, thereby improving the accuracy of biological sample data analysis.

[0041] For example, when processing read length R1 (12345bp), the first 1kb (1-1000bp) and the last 1kb (11346-12345bp) of 15-mer seeds are extracted, generating a total of 133 k-mers (66 at the beginning and 66 at the end, with 1 overlap). The topological paths of the first and last k-mers are located in the pan-genome map. The coordinates of the first k-mer path are chr3:50,100,000-50,100,015, and the coordinates of the last k-mer path are chr3:50,112,000-50,112,015. The path distance is 12015bp; the ratio of path distance to read length is calculated (12,015 / 12,345≈0.97<1.2), and a candidate subgraph set {subgraph A(chr3:50,090,000-50,120,000)} is directly generated; a GC=75% region is detected in subgraph A, and the k-mer step size is increased from 15bp to 20bp; for low-complexity regions in subgraph A (such as the AT repetition region chr3:50,105,000-50,105,050), fuzzy matching is enabled (allowing 2 mismatches).

[0042] In some embodiments, S4 includes: statistically analyzing the distribution of breakpoint locations supporting reads to calculate breakpoint confidence intervals; determining breakpoint intervals using the 95th percentile (e.g., chr1:100, 200-100, 500); extracting all reads within the breakpoint intervals for local assembly refinement; performing DeBruijn plot assembly (k=31); aligning the assembled sequences to a reference plot to determine precise breakpoint coordinates; and additionally checking sequence micro-homology (homologous sequences >5 bp) when the breakpoint is located in a repetitive region. According to embodiments of the present invention, by summarizing supporting data from multiple reads, random noise is reduced, providing statistically reliable breakpoint location estimates, thereby reducing the risk of false positives or false negatives; by defining breakpoint intervals based on high confidence, it is ensured that breakpoint locations are covered within a reasonable range, reducing mislocalization due to individual outliers; by utilizing all reads within the breakpoint region for local sequence reconstruction, incomplete information from individual reads is compensated for, improving sequence accuracy and variant detection precision, thereby enhancing the accuracy of biological sample data analysis.

[0043] In some embodiments, the method further includes: performing technical noise filtering on the final variant dataset, specifically: removing variants with chain bias > 80%.

[0044] In some embodiments, the method further includes: performing biological rationality verification on the final variant dataset, specifically: connecting gene annotation information on both sides of the breakpoint; if the breakpoint disrupts an exon splicing site, it is marked as a high-impact variant.

[0045] In some embodiments, the method further includes: performing population frequency correction on the final variant dataset, specifically: comparing with the 1,000 Genomes Database, marking common variants with a frequency >1%; and initiating phenotypic association analysis for rare variants (frequency <0.1%).

[0046] It should be noted that, for the sake of simplicity, the aforementioned method embodiments are all described as a series of actions. However, those skilled in the art should understand that the present invention is not limited to the described order of actions, as some steps can be performed in other orders or simultaneously according to the present invention. Secondly, those skilled in the art should also understand that the embodiments described in the specification are optional embodiments, and the actions and modules involved are not necessarily essential to the present invention. The acquisition, storage, and application of user personal information involved in the technical solution of the present invention comply with relevant laws and regulations and do not violate public order and good morals. It should be understood that the various forms of processes described above can be used to reorder, add, or delete steps. For example, the steps described in the present invention can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution of the present invention can be achieved, and this is not limited herein. The above specific embodiments do not constitute a limitation on the scope of protection of the present invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for analyzing biological sample data, characterized in that, include: S1: Obtain pangenome map and sequencing read data of biological samples. Dynamically divide the pangenome map into multiple sub-maps through hierarchical map segmentation. The segmentation density is automatically adjusted according to the variation complexity of the genomic region. S2 performs topology-aware dynamic programming alignment of sequencing read data within each subgraph, and simultaneously processes sequence position, node position and path selection through a three-dimensional scoring matrix to generate corresponding subgraph alignment results; S3, based on the minimum cost flow algorithm, merges all subgraph alignment results to construct a globally consistent genome alignment path; S4 detects structural variations based on global path comparison and outputs a variation dataset containing precise breakpoint coordinates.

2. The method for analyzing biological sample data according to claim 1, characterized in that, After S4 detects structural variants, the process also includes: if the detected variant region contains low-quality bases with base quality values ​​below the base quality value threshold, a secondary validation is performed. Specifically, this involves: locating all sequencing reads within a 50bp range on both sides of the variant breakpoint; calculating the number of reads supported by the breakpoint using a local re-alignment algorithm, and retaining the variant when the number of supported reads is ≥3; predicting the functional impact of the retained variants and labeling their clinical significance level using the ClinVar database; the secondary validation process is executed in parallel using GPUs, with each variant assigned an independent computation thread; if the variant is located in a known pathogenic gene region, an additional read length consistency check is initiated, utilizing the PacBioHiFi read length crossing breakpoint characteristic to verify the variant's authenticity.

3. The method for analyzing biological sample data according to claim 2, characterized in that, The S1 layered graph segmentation process includes: when the variation density is >0.3, the segmentation level is automatically increased: the top layer segmentation generates 5-10 macro sub-graphs for quickly locating variation regions; the middle layer segmentation generates micro sub-graphs within the top layer sub-graphs to focus on high variation density regions; and the bottom layer segmentation generates nano sub-graphs for complex structural variation regions.

4. The method for analyzing biological sample data according to claim 3, characterized in that, The layered graph segmentation process of S1 satisfies the following constraints: the number of connecting edges between subgraphs is less than 15% of the internal edges of the subgraph; when a branch region with a node out-degree greater than the node out-degree threshold is detected, the integrity of the branch path is preserved first, and path cutting is avoided by adding virtual nodes.

5. The method for analyzing biological sample data according to claim 4, characterized in that, S2's three-dimensional scoring matrix includes topological consistency verification: if the current node v b With the predecessor node v b' If the paths are not consecutively adjacent on the topological path, a target penalty score is applied; when the score difference between multiple candidate paths is less than 2 points, the high-frequency path is selected first; for reads that cross the subgraph boundary, if the score difference between the alignments in adjacent subgraphs is greater than 10 points, a 100bp sequence window is extracted on both sides of the boundary; an extended subgraph is constructed; and a second dynamic programming alignment is performed in the extended subgraph.

6. The method for analyzing biological sample data according to claim 5, characterized in that, Following S3, it also includes: real-time monitoring of the time T spent comparing each subgraph. i When T max / T avg >2,T max For the maximum time consumption, T avg To ensure that the average time taken is less than the GPU utilization threshold, the subgraph with the highest time taken is further divided into several new subgraphs, which are then allocated to idle GPU threads. The GPU memory resources of the subgraphs that have been compared are released.

7. The method for analyzing biological sample data according to claim 6, characterized in that, Also includes: For low-quality bases whose base quality values ​​are below the base quality value threshold, a matching score correction is applied: σ'=σ×(1-10) (-Q / 10) ), where σ is the original matching score and Q is the base quality value; if the frequency of occurrence of the base is greater than the frequency threshold, the bonus score is increased; when Q < 20 for 3 consecutive bases D, the topological constraints are relaxed for re-alignment; k-mer seed expansion is used to re-anchor the position; after re-alignment, candidate results with an alignment score improvement of < 5 are deleted.

8. The method for analyzing biological sample data according to claim 7, characterized in that, Also includes: If the read length is >10kb, then before S2, extract the first and last 1kb k-mers of the read length, k=15, for topology-guided seed localization; locate the topological path of the k-mer in the pan-genome map; when the distance between the first and last paths is < read length × k, 1 < k < 1.5, directly generate candidate sub-maps; for regions with an error rate greater than the error rate threshold, increase the k-mer step size to 20bp.

9. The method for analyzing biological sample data according to claim 8, characterized in that, S4 includes: statistically analyzing the distribution of breakpoint locations for support reads to calculate breakpoint confidence intervals; determining breakpoint intervals using the 95th percentile; extracting all reads within the breakpoint interval for local assembly refinement; performing DeBruijn plot assembly; aligning the assembled sequences to the reference plot to determine precise breakpoint coordinates; and additionally checking sequence micro-homology when the breakpoint is located in a repetitive region.

10. The method for analyzing biological sample data according to claim 9, characterized in that, Also includes: Technical noise filtering is performed on the final mutated dataset, specifically by removing mutated datasets with chain bias > 80%.

Images