A structural variation detection algorithm, system, device and medium based on third-generation sequencing data and pan-genome

Abstract

This invention discloses an algorithm, system, device, and medium for detecting structural variations based on third-generation sequencing data and a pan-genome. The detection algorithm includes: detecting SNARL structures in the pan-genome map; extracting reads corresponding to each SNARL path from the GAM alignment file; calculating the average coverage of all edges in each SNARL path and the number of edges with zero coverage; statistically analyzing the possible paths, path directions, reads aligned to the path, and path coverage information contained in each SNARL path; selecting the optimal path and the second path based on the read information and path coverage information; and comparing the optimized optimal path and the second path with the reference path to obtain variation information. This invention integrates the read information corresponding to the path, the base coverage information of the path, and the edge coverage information of the path as the basis for selecting potential variation paths, thus improving the accuracy of third-generation sequencing data detection.

Description

Technical Field

[0001] This invention belongs to the field of gene sequence structural variation detection technology, specifically involving a structural variation detection algorithm, system, device and medium based on third-generation sequencing data and pan-genome. Background Technology

[0002] Variations in organisms generally include minor variations and structural variations. Minor variations include single nucleotide polymorphisms (SNPs), small sequence insertions (INS), and deletions (DELs), the latter two often collectively referred to as INDEL variations. SNPs primarily refer to DNA sequence diversity caused by variations in a single nucleotide at the genomic level; they are numerous and are important evidence for studying genetic variation in human families and animal and plant strains. INDELs refer to the insertion or deletion of a small sequence fragment at a specific location in the genome, typically less than 50 bp in length. Structural variants (SVs) are genomic mutations involving 50 or more base pairs. SVs can generally be classified into five types: deletion, insertion, duplication, inversion, and translocation, as well as complex structural variations formed by combinations of these types. Although structural variations account for only a small portion of all variations, they are usually quite significant, and current analyses indicate they have a high impact on gene expression. Structural variations can affect regulatory elements such as gene promoters and enhancers, or influence chromosome structure, thereby altering gene expression levels. In some cases, structural variations can lead to the fusion of two genes. This fusion may produce new fusion proteins with novel functions or lose their original functions. Certain structural variations are also associated with the development of genetic diseases, such as certain types of cancer and hereditary heart disease. Furthermore, structural variations have played a crucial role in species evolution, promoting adaptability and diversity by altering genome structure. Therefore, the study of structural variations is of great significance to both human medicine and genetics.

[0003] The study of gene variation depends on advancements in molecular biology techniques. In the past, the dominant DNA sequencing technology in genomics research was undoubtedly the Illumina platform, a technology that generates short-read sequencing data. This technology produces highly accurate (over 99.9%) sequencing data at a low cost. However, short-read sequencing technology has a major drawback in handling larger structural variations and genome assembly: the limited read length prevents it from detecting many structural variations. Therefore, in recent years, long-read sequencing technologies have developed rapidly, primarily ONT sequencing and PacBio sequencing. These technologies can generate longer reads of thousands or even tens of thousands of bases. Analyzing these reads in the genome can reveal much previously undiscovered genetic information.

[0004] A common method for detecting variants is to first align sequencing data to a reference genome and then identify variants based on alignment discrepancies, as seen in SVTyper and Delly genotyping models. However, traditional reference genomes are not sensitive to variants and have several drawbacks. For example, they only select one version of the variant for embedding at each variant location, thus failing to accurately represent gene diversity. This can lead to sequencing data containing variant genes failing to align to the reference genome, resulting in the inability to detect variant signals. To address the shortcomings of linear reference genomes, the concept of pangenomes has recently been proposed. Pangenomes can integrate genetic information from multiple samples and retain multiple variant sequences at each variant site. Compared to linear reference genomes, pangenome graphs can encode more variant information, mitigating the drawbacks of linear reference genomes and improving variant detection and genotyping. Pangenomes are typically represented using directed graphs, where nodes represent gene sequences and edges represent connections between sequences. Therefore, the alignment of sequencing data to a linear reference genome is extended to the alignment of paths from sequencing data to the pangenome graph.

[0005] Currently, there are few tools that directly perform structural variant detection and genotyping based on pan-genome maps, with typical examples including VG and GraphTyper. However, VG directly selects the two pathways with the highest average coverage as potential variant pathways (the coverage of the corresponding pathway can be obtained based on the alignment information from sequencing data to the pathway), and then analyzes different cases for these two pathways. The method for detecting and genotyping structural variants is relatively simple, but if the selected pathway is incorrect, the subsequent analysis error will be significant. GraphTyper can only build pan-genome maps for analysis in small regions of the genome and cannot be directly applied to pan-genome maps of the entire chromosome. Furthermore, it is only suitable for second-generation sequencing data with relatively short sequencing lengths and is not applicable to third-generation sequencing data. As the application of third-generation sequencing data becomes more widespread, the shortcomings of GraphTyper are becoming increasingly apparent. Therefore, it is necessary to develop a structural variant detection algorithm based on third-generation sequencing data and pan-genome map structure to improve the accuracy of structural variant detection and provide a reliable foundation for subsequent analysis. Summary of the Invention

[0006] To overcome the problem of low accuracy in structural variation detection in existing technologies, the purpose of this invention is to provide a structural variation detection algorithm, system, device, and medium based on third-generation sequencing data and pan-genome. The detection results of this algorithm are more reliable and accurate, thus improving the accuracy of structural variation detection.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0008] A structural variation detection algorithm based on third-generation sequencing data and pan-genome includes the following steps:

[0009] Obtain the pangenome map;

[0010] Detect SNARL structures in the pan-genome map and extract reads corresponding to each SNARL path from the GAM alignment file; calculate the average coverage of all edges in each SNARL path and the number of edges with zero coverage; and count the path, path direction, reads aligned to the path, and path coverage information contained in each SNARL.

[0011] Based on reads information and path coverage information, the optimal path and the second path are selected;

[0012] Optimize the optimal path and the second path;

[0013] The optimized optimal path and the second path are both compared with the reference path to obtain mutation information.

[0014] Furthermore, based on read information and path coverage information, the optimal path and the second path are selected, including the following steps: sort the paths according to the number of reads matched to the SNARL path, and select the path with the highest number of reads as the optimal path. If the number of reads is less than 5, sort them according to the calculated average coverage of the paths, and select the path with the highest average coverage as the optimal path. Iterate through the remaining paths in turn, calculate the number of unique reads matched to each path, and select the path with the highest number of unique reads as the second path.

[0015] Furthermore, the optimal path and the second path are optimized, including the following steps:

[0016] All paths are sorted in ascending order of path coverage; if coverage is the same, they are sorted in ascending order of path length; if both coverage and length are the same, they are sorted in descending order of path extraction from SNARL. Then, the optimal path is optimized. If the number of reads for the optimal path is less than a set value, the last path after sorting is selected as the optimal path. Next, based on the base coverage information obtained using the VG alignment tool, the path with the highest average base coverage, `path_m`, is calculated. If the base coverage of `path_m` is unique and greater than twice that of all other paths, and the edges of `path_m` have the highest average base coverage... If the average coverage is similar to the optimal path, then the optimal path is the path path_m with the largest average base coverage. Based on the selected optimal path, the path coverage of other paths is re-evaluated and updated. Then, the paths are sorted in ascending order of average edge coverage, ascending order of path length, and descending order of path sequence. After sorting, the paths are traversed from back to front. If it is the optimal path, it is skipped; otherwise, it is selected as the second path. Then, the SNARL paths are traversed again. If there is a path whose path length is greater than a set value and twice as long as the second path, and the average edge coverage of this path is comparable to that of the second path, then the original second path is discarded and this path is selected as the second path.

[0017] Furthermore, the optimal path is compared with the second path and the reference path to obtain variation information, including the following steps:

[0018] When the optimal path is the reference path, if the difference in the number of bases between the reference path and the second path is less than a set threshold, then the mutation in the SNARL region is a small mutation. If the path coverage of the optimal path is less than or equal to 2, the path coverage of the optimal path is greater than 4 times that of the second path, the number of reads of the optimal path is greater than 3 times that of the second path, and the number of edges with 0 coverage in the second path is more than 0.4 times the total number of edges, then there is no mutation in the SNARL region; otherwise, there is mutation, and the second path is retained.

[0019] When the second path is the reference path, the optimal path is the mutated path, and mutations are normal.

[0020] When neither the optimal path nor the second path is the reference path, compare the optimal path and the reference path. If the number of bases is not significantly different, then there is no SV.

[0021] Furthermore, if the number of reads for the optimal path and the second path is 0, the coverage of the optimal path and the second path is not 0, and the difference is not significant, then the SNARL region is considered to have a mutation, and the second path is retained.

[0022] If the optimal path has only 2 nodes and the second path has at least 4 nodes, and the path coverage of the second path is greater than half that of the optimal path, then this SNARL region is considered to have a mutation, and the second path is retained.

[0023] Furthermore, when the second path is the reference path, the optimal path is the variant path. Under normal circumstances, there is a variant, which includes the following steps: if the difference in the number of bases between the optimal path and the reference path is less than a set threshold, then the variant in the SNARL region is a small variant; if the coverage of the second path is more than 4 times the coverage of the optimal path, then the variant is ignored; otherwise, it is considered that there is a variant here, and the optimal path is retained.

[0024] If the coverage of the reference path is greater than 4 times that of the optimal path, the number of reads in the optimal path is less than 3, and the path coverage is less than 5, then there is a potential mutation in the optimal path, and the optimal path is retained. Then the second path and the reference path are compared. If the number of bases is not much different, then there is no SV.

[0025] If the coverage of the optimal path is more than 4 times that of the second path, the number of reads of the optimal path is more than 3 times that of the second path, and the number of edges with 0 coverage in the second path is more than 0.4 times the total number of edges, then there is no potential mutation in the second path; otherwise, it is considered that there is a second mutation and the second path is retained.

[0026] If the length of the second path is more than four times that of the optimal path, and their coverage is comparable, then a second mutation is considered to exist, and the second path is retained.

[0027] Furthermore, based on the preserved path and the reference path, specific variation information is extracted. The variation position is the last position of the first node in the reference path. The reference sequence is the node sequence found in the reference path in the pan-genome map file, which is then assembled sequentially in either forward or reverse order according to the node direction. The variation sequence is the node sequence found in the reference variation path in the pan-genome map file, which is then assembled sequentially in either forward or reverse order according to the node direction to obtain the structural variation VCF output file.

[0028] A structural variation detection system based on third-generation sequencing data and pan-genome includes:

[0029] The pan-genome map acquisition module is used to acquire pan-genome maps;

[0030] The data preprocessing module is used to detect SNARL structures in the pan-genome map and extract reads corresponding to each SNARL path in the SNARL structure from the GAM alignment file; calculate the average coverage of all edges in each SNARL path and the number of edges with zero coverage; and count the possible paths, path directions, reads aligned to the path, and path coverage information in each SNARL.

[0031] The mutation path selection module is used to filter the optimal path and the second path based on reads information and path coverage information.

[0032] The path selection optimization module is used to optimize the optimal path and the second path;

[0033] The coarse path comparison module is used to compare the optimized optimal path and the second path with the reference path to obtain variation information.

[0034] An electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that the processor, when executing the computer program, implements the structural variation detection method based on third-generation sequencing data and pan-genome algorithm.

[0035] A computer-readable storage medium storing a computer program that, when executed by a processor, implements the structural variation detection method based on third-generation sequencing data and a pan-genome algorithm.

[0036] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0037] This invention utilizes the reads information corresponding to the path, the base coverage information of the path, and the edge coverage information of the path, and integrates the three as the basis for selecting potential variant paths. The overall method is more reliable and accurate, overcoming the problem in the prior art that VG's SNARL analysis of the pan-genome map is relatively simple. When selecting paths in SNARL, it directly selects the two paths with the highest base coverage as candidate paths. However, sometimes the base coverage is unreliable, which makes these two paths unreliable, and subsequent analyses will also be unreliable.

[0038] After selecting the optimal path, this invention eliminates the influence of the optimal path on other paths, including the influence of reads and edge coverage, before selecting the second path. This makes the results more reliable and realistic, overcoming the problem in existing technologies where VG and GraphTyper do not eliminate the influence of the optimal path on other paths in path selection, which may lead to the overall coverage of other paths being unrealistically high and affecting the accuracy of subsequent results.

[0039] Furthermore, in this invention, after selecting the optimal and second paths, if the read alignment results are unsatisfactory, a path selection optimization logic is continued. The optimal path is compared with the path with the highest base coverage, and the second path is compared with paths more than twice its length and with comparable coverage. The final candidate mutation paths are determined based on a specific strategy. This ensures more reliable mutation paths and more accurate subsequent results. It overcomes the problem in existing technologies where VG and GraphTyper, after selecting the optimal and second paths, lack path selection optimization logic, potentially leading to paths with slightly lower coverage but better overall performance. These paths should be selected as candidate mutation paths, thus affecting the accuracy of the results.

[0040] Furthermore, this invention comprehensively considers the difference in the number of reads and the multiple of coverage of the path, and also considers the proportion of edges with zero coverage in the candidate mutation path. It makes full use of the available information of the candidate mutation path, making the result more reliable. It overcomes the problem in the prior art that after VG selects the optimal path and candidate paths, it directly determines whether there is an SV based on the multiple of the base coverage of the path. The method is simple, but does not consider complex cases.

[0041] Furthermore, this invention also considers some complex situations, such as many reads being matched to uncounted SNARL paths, resulting in very few reads corresponding to the counted paths, and the second path having few reads but a length much greater than the optimal path and comparable coverage. There are corresponding strategies to solve these problems, making the results more accurate. Attached Figure Description

[0042] Figure 1 This is a schematic diagram of the structural variation detection algorithm based on third-generation sequencing data and pan-genome in this invention;

[0043] Figure 2 Flowchart for selecting the mutation path module;

[0044] Figure 3 This is a diagram illustrating the impact of the optimal path in SNARL on the reads of other paths.

[0045] Figure 4 Flowchart for the path selection optimization module;

[0046] Figure 5 This is a diagram illustrating the impact of the optimal path on the coverage of other paths in SNARL.

[0047] Figure 6 The overall box plot shows the F1 scores of the three structural variation detection tools;

[0048] Figure 7 An overall radar chart showing the results from the three tools;

[0049] Figure 8 Bar charts showing the F1 scores of the three tools on different chromosomes;

[0050] Figure 9 F1 scores for different sizes of insertion variants using different tools;

[0051] Figure 10 This is a schematic diagram of a structural variation detection system based on third-generation sequencing data and pan-genome. Detailed Implementation

[0052] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Preferred embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a thorough and complete understanding of the disclosure of the invention.

[0053] See Figure 1 , Figure 1 This is a flowchart illustrating the structural variation detection algorithm based on third-generation sequencing data and pan-genome of the present invention, which mainly includes four steps: data preprocessing, variation path selection, path selection optimization, and coarse path comparison.

[0054] Obtain the pangenome map;

[0055] Data preprocessing: SNARL structures in the pan-genome graph (chrom.vg file) were detected, and reads corresponding to the SNARL paths in each SNARL structure were extracted from the GAM alignment file. The SNARL structure in the graph is such that removing all edges between the left and right parts of the pan-genome graph and the nodes at both ends of this structure will split the graph, forming a connected graph containing the nodes at both ends. A schematic diagram of the SNARL structure is shown below. Figure 1As shown in the data preprocessing module, each SNARL position is a potential location for mutations, much like "bubbles" in the genome. The average coverage of all edges in each SNARL path and the number of edges with zero coverage are then recalculated. When calculating path coverage, only the average coverage of edges within the path is considered, not the coverage of bases within nodes. This is because, according to the characteristics of pan-genome graphs, edge coverage is more important; high edge coverage indicates high path coverage. Furthermore, some nodes contain very long base sequences with repetitive sequences. The alignment of the middle base sequences may be poor, resulting in high coverage for the two edges related to a node but low coverage for the bases within the node itself, leading to low overall path base coverage and affecting the accuracy of the results. Two types of coverage information are presented here: one is the average coverage of all bases in the path calculated using the VG alignment tool, and the other is the recalculated average coverage of all edges in the path. Unless otherwise specified, the path coverage mentioned here refers to the average coverage of all edges in the path. Finally, the possible paths contained in each SNARL, path directions, reads aligned to the path, and path coverage information are compiled for further analysis.

[0056] Mutation Path Selection: Two possible mutation paths are initially selected: the optimal path and the second path. Since SV involves numerous SNAL paths, combining them pairwise and calculating likelihood as in GraphTyper would result in a vast number of path combinations. Therefore, the strategy is changed: first, two paths are selected, and then compared with a reference path to detect mutations. Specifically, as... Figure 2As shown, the first step is to identify the optimal path based on read information and path coverage information: Paths are first sorted according to the number of reads aligned to the SNARL path, and the path with the highest number of reads is selected as the optimal path. If the number of reads is less than 5, it indicates poor read alignment in this SNARL region. In this case, the paths are sorted according to the calculated average coverage, and the path with the highest average coverage is selected as the optimal path. Since read information is more reliable than average path coverage (the exact number of reads aligned to the path is more precise; however, if a small portion of a read's sequence aligns to a small portion of the SNARL path's base sequence, and the rest is not aligned, the coverage statistics will still include this small portion, leading to inaccurate results), the number of reads aligned to the SNARL path is prioritized. After identifying the optimal path, the remaining paths are iterated through sequentially, and the number of unique reads aligned to each path is calculated. The path with the highest number of unique reads is selected as the second path. This calculation focuses on the number of unique reads for each path, not the total number of reads. Unique reads are those that appear in the SNARL path but not on the optimal path. The reason for calculating unique reads is that there are situations where multiple paths are very similar to the optimal path, differing only in a few nodes, and these few nodes have very low coverage. If selection were solely based on the number of reads on a path, these similar paths would be retained, even though the reads on these similar paths are actually on the optimal path, and the higher average coverage is also due to the influence of the optimal path. Therefore, the number of unique reads for each path is used to minimize the influence of the optimal path on other paths, making the overall selection more inclined towards paths with lower similarity to the optimal path. Figure 3 As shown, path1 corresponds to 8 reads, path2 to 5 reads, and path3 to 3 reads. However, 4 out of the 5 reads in path2 are duplicates of those in path1. After selecting path1 as the optimal path, path2 has only 1 unique read, while path3 still has 3 unique reads. Therefore, path3 should be selected as the second path. After selecting the second path, if the number of reads in the optimal path or the second path is 0, or if the average coverage of the second path is very low and there are other paths with a coverage greater than 4 times that of the second path, it indicates that the alignment of reads in this SNARL region is not good. The path optimization logic is then used to further analyze the path coverage.

[0057] Path selection optimization: When the alignment of reads in the SNARL region is poor, optimize the path selection strategy. For example... Figure 2 As shown, optimizations are performed on the optimal path and the second path, respectively. Specifically, as... Figure 4 As shown, all paths are sorted in ascending order of path coverage. If the coverage is the same, they are sorted in ascending order of path length, because the longer the path, the more difficult it is to achieve the same coverage, and the higher the reliability. If the path coverage and length are the same, they are sorted in descending order according to the order of path extraction in SNARL. Since both path coverage and length are sorted in ascending order, the further to the right, the better. When all conditions are the same, the path that was earlier in the original path is preferred, so the path with the smaller index is on the right. Then, the optimal path is optimized. If the number of reads for the previously selected optimal path is very small (less than the set value of 2), the last path after the above sorting is selected as the optimal path (the two may be the same; if they are the same, it means that the optimal path selected by both methods is the same, which is relatively reliable). Then, the path with the highest average base coverage (path_m) obtained using the VG alignment tool is found. If the base coverage of path_m with the highest average base coverage is unique and greater than twice that of all other paths (including the optimal path), and the average edge coverage of path_m with the highest average base coverage is similar to that of the optimal path (i.e., the average edge coverage of path_m with the highest average base coverage is greater than 0.75 times that of the optimal path), then the optimal path is replaced by path_m with the highest average base coverage. Based on the selected optimal path, the path coverage of other paths is re-evaluated to eliminate the influence of the optimal path on them. The influence of the optimal path on the reads of other paths has already been eliminated, but at this point, the SNARL region read alignment is very poor, making it unsuitable for further analysis using reads. Instead, path coverage analysis should be used, so it is necessary to continue eliminating the influence of the optimal path on the path coverage of other paths. This approach prevents certain paths from being completely contained within the optimal path. Edges that differ from the optimal path, which would originally have zero coverage, may have relatively high coverage due to the presence of the optimal path, leading to an unrealistically high overall coverage for these paths. For example... Figure 5As shown, path1 is the optimal path with an average coverage of 9.3, path2 has an average coverage of 7.3, and path3 has an average coverage of 5. However, the higher coverage of path2 is largely due to the influence of path1. The two edges in the middle that are different from those of path1 have very low coverage. Therefore, the influence of duplicate edges in path1 should be eliminated, and its coverage should be updated to 2. Finally, path3 is selected as the second path. Therefore, the coverage of these paths needs to be adjusted to ensure their authenticity and accuracy. The method for updating the coverage is relatively simple: if an edge in another path is exactly the same as one in the optimal path, this edge should be ignored to avoid the influence of the optimal path. After updating the coverage, the paths are sorted again in ascending order of average edge coverage, ascending order of path length, and descending order of path sequence. After sorting, the paths are traversed from back to front. If it is the optimal path, it is skipped; otherwise, it is selected as the second path. The SNARL paths are then traversed again. If a path exists that is significantly longer (by default, more than 15 nodes) than the second path and is twice as long, and its average edge coverage is comparable to the second path, then the original second path is discarded, and this new path is retained as the second path. This is because it is more difficult for longer paths to achieve the same overall path coverage compared to shorter paths. If the length difference is more than double and the path coverage is comparable, the longer path is more accurate and reliable, and is more likely to be a true potential mutation path. Thus, a relatively reliable potential mutation path is obtained.

[0058] Coarse path comparison: Further analysis is performed on the selected optimal and second paths, comparing them with the reference path to obtain the final mutation information. First, the total base length of the reference path and the two selected paths (optimal and second paths) is calculated to determine the presence and type of SVs based on the length difference. Then, as follows... Figure 1 The coarse path comparison module shows three scenarios for analysis:

[0059] 1. The optimal path is the reference path: If the number of bases in the reference path (i.e. the optimal path) and the second path is not much different, less than 50 bases (i.e., the set threshold), then the mutation in this SNARL region is considered to be a small mutation and is ignored. The second path is ignored and the SNARL region is considered to have no mutation if any of the following situations occur: the path coverage of the optimal path is very small (i.e. less than or equal to 2), the path coverage of the optimal path is greater than 4 times that of the second path, the number of reads in the optimal path is greater than 3 times that of the second path, or the number of edges with 0 coverage in the second path is more than 0.4 times the total number of edges. Otherwise, the mutation is considered to exist and the second path is retained. In addition, there are two special cases: (1) If the number of reads in the optimal path and the second path is 0, but the coverage of the optimal path and the second path is greater than 10, and the coverage of the optimal path is less than twice the coverage of the second path, it means that there is an uncounted SNARL path, and most reads are matched to the uncounted SNARL path, which affects the coverage of the counted path. Because the SNARL region may be large and complex, containing many nodes and potential paths, the statistics may not be complete. The reference path can definitely be counted here, so the path that is not counted must be the mutation path. Many reads are aligned to the mutation path that is not counted, so at this time, it is considered that this SNARL region has mutation, and the second path is retained. (2) Since insertion mutation is more difficult to detect than deletion mutation, optimization is performed for insertion mutation: if the optimal path has only 2 nodes and the number of nodes in the second path is not less than 4 (indicating that it is an obvious insertion mutation), and the path coverage of the second path is greater than half of that of the optimal path, then it is considered that this SNARL region has mutation, and the second path is retained. The influence of reads is no longer considered here, because in insertion mutation, if the insertion sequence is very long, it is difficult for reads to cover it and difficult to align to the correct position.

[0060] 2. The second path is the reference path: In this case, the optimal path is the variant path, and under normal circumstances, there will definitely be a variant. If the number of bases in the optimal path and the reference path (i.e., the second path) is not significantly different (less than 50 bases), then the variant in this SNARL region is considered a small variant and is ignored. If the coverage of the second path is more than four times that of the optimal path, then the optimal path may have been selected based on the number of reads, but the overall number of reads is very small, and the variant is unreliable. In this case, the variant is ignored. Otherwise, a variant is considered to exist here, and the optimal path is retained.

[0061] 3. Neither the optimal path nor the second path is a reference path: In this case, both the optimal and second paths need to be compared with the reference path. First, compare the optimal path and the reference path. Similarly, if the number of bases is not significantly different (less than 50 bases), it is considered that there is no potential mutation (SV). If one of the following conditions exists, it is considered that there is no potential mutation in the optimal path: the coverage of the reference path is greater than 4 times that of the optimal path; the number of reads (less than 3) and the path coverage (less than 5) of the optimal path are both very small. Otherwise, it is considered that there is a potential mutation in the optimal path, and the optimal path is retained. Then compare the second path and the reference path. Similarly, if the number of bases is not significantly different (less than 50 bases), it is considered that there is no SV. If one of the following conditions exists, it is considered that there is no potential mutation in the second path: the coverage of the optimal path is greater than 4 times that of the second path; the number of reads of the optimal path is greater than 3 times that of the second path; the number of edges with 0 coverage in the second path is more than 0.4 times the total number of edges. Otherwise, it is considered that there is a second mutation, and the second path is retained. There is a special case here: if the length of the second path is more than four times that of the optimal path, but their coverage is comparable (the coverage of the second path is 0.75 times that of the optimal path), then a second mutation is considered to exist, the second path is retained, and no further filtering using reads information is performed. This is because the number of nodes in the second path differs significantly from the optimal path, and the edge redundancy is also very low, indicating a weak correlation between the two paths. The coverage of the second path is not closely related to that of the optimal path, and a mutation is highly likely.

[0062] Finally, based on the preserved path and the reference path, specific variant information is extracted. If sequence alignment algorithms are used directly, the results may be very fine; a large SV might be broken down into multiple SVs and smaller variants. Therefore, a coarse alignment is performed using nodes in the path. The variant position is the last position of the first node in the reference path; the reference sequence is the node sequence found in the pan-genome map file, assembled sequentially in either forward (normal base sequence, such as ACGT) or reverse (complementary base sequence, such as TGCA) according to the node orientation; the variant sequence is the node sequence found in the pan-genome map file in the reference variant path, assembled sequentially in either forward or reverse according to the node orientation. This yields the final structural variant VCF output file.

[0063] To verify the accuracy of the algorithm in this invention, experiments were conducted on real samples. The sample used in this experiment was the HG002 sample. The Pan-SV and VG methods of this invention were run using the third-generation sequencing data of HG002, and GraphTyper was run using the second-generation sequencing data of HG002. The Truvari tool was used to detect the accuracy of structural variations.

[0064] The experimental results are shown in the attached figure. Figure 6The figure shows the overall box plot of the F1 scores of the three structural variation detection tools. As can be seen from the figure, whether in the high confidence region or in all regions, the overall F1 score of Pan-SV is better than that of VG and far superior to GraphTyper.

[0065] Figure 7 The overall radar chart of the results from the three tools shows that Pan-SV is closest to the upper right region and has the best performance.

[0066] Figure 8 The bar charts showing the f1 scores of the three tools on different chromosomes demonstrate that Pan-SV (the algorithm in this invention) yields the best results on all chromosomes.

[0067] Figure 9 The F1 scores of different tools for different sizes of insertion mutations show that Pan-SV performs best for different sizes of insertion mutations in the graph. Moreover, as the insertion mutation increases, the decrease in Pan-SV is smaller than that of VG. Although GraphTyper performs better as the insertion mutation increases, its overall F1 score is very low.

[0068] Advantages of this invention:

[0069] In existing technologies, VG's SNARL analysis of pan-genome maps is relatively simple. When selecting paths in SNARL, it directly selects the two paths with the highest base coverage as candidate paths. However, base coverage is sometimes unreliable, making these two paths unreliable, and subsequent analyses will also be unreliable. This invention utilizes the reads information corresponding to the path, the path's base coverage information, and the path's edge coverage information, integrating these three as the basis for selecting potential variant paths, resulting in a more reliable and accurate overall approach.

[0070] In existing technologies, neither VG nor GraphTyper excludes the influence of the optimal path on other paths during path selection, which may lead to an unrealistically high overall coverage of other paths, affecting the accuracy of subsequent results. This invention, however, after selecting the optimal path, excludes its influence on other paths, including the impact of reads and edge coverage, before selecting a second path, making the results more reliable and accurate.

[0071] In existing technologies, VG and GraphTyper do not perform path selection optimization after selecting the optimal and second paths. However, some paths with slightly lower coverage but better overall performance may exist, and these paths should be selected as candidate mutation paths. In this invention, after selecting the optimal and second paths, if the read alignment is unsatisfactory, path selection optimization logic continues. The optimal path is compared with the path with the highest base coverage, and the second path is compared with paths more than twice its length and with comparable coverage. The final candidate mutation path to be retained is determined based on a specific strategy. This results in more reliable retained mutation paths and more accurate subsequent results.

[0072] In existing technologies, after VG selects the optimal and candidate paths, it directly determines whether an SV exists based on the multiple difference in base coverage between the paths. This method is simple but does not consider complex cases. This invention, however, comprehensively considers both the number of reads and the multiple difference in coverage between paths, and also takes into account the proportion of edges with zero coverage in the candidate mutation paths, fully utilizing the available information of the candidate mutation paths to make the results more reliable. Furthermore, this invention also considers some complex cases, such as many reads aligning to uncounted SNARL paths, resulting in very few reads corresponding to the counted paths, and second paths having few reads but a length much greater than the optimal path with comparable coverage. These issues are addressed with corresponding strategies, making the results more accurate.

[0073] This invention combines multiple pieces of information from the pan-genome graph to select candidate variant paths: it utilizes read information, base coverage information, and edge coverage information corresponding to each path, using these three types of information together as the basis for selecting potential variant paths. During path selection, the influence of the optimal path on other paths is eliminated, improving the reliability of path selection: after selecting the optimal path, its influence on other paths, including read and edge coverage, is eliminated before selecting the second path, making the result more reliable and accurate. After selecting the optimal and second candidate variant paths, a path selection optimization logic continues to retain more reliable paths: if read alignment is unfavorable after selecting the optimal and second paths, path selection optimization continues, comparing the optimal path with the path with the highest base coverage, and comparing the second path with paths more than twice the length and with comparable coverage. The final candidate variant path to be retained is determined based on a specific strategy. This invention fully utilizes the available information of candidate variant paths to determine whether SVs exist within a path and considers some complex situations.

[0074] This invention comprehensively considers the difference in read count and coverage multiples between paths, and also takes into account the proportion of edges with 0 coverage in candidate mutation paths, making full use of the available information of candidate mutation paths and thus making the results more reliable. Furthermore, this invention also considers some complex situations, such as many reads being matched to uncounted SNARL paths, resulting in very few reads corresponding to the counted paths, and second paths having few reads but a length much greater than the optimal path with comparable coverage. These situations are addressed with corresponding strategies, making the results more accurate.

[0075] Example 2

[0076] See Figure 10 A structural variation detection system based on third-generation sequencing data and pan-genome includes: a pan-genome map acquisition module, a data preprocessing module, a variation path selection module, a path selection optimization module, and a path coarse comparison module.

[0077] The pan-genome map acquisition module is used to acquire pan-genome maps;

[0078] The data preprocessing module is used to detect SNARL structures in the pan-genome map and extract reads corresponding to each SNARL path from the GAM alignment file; calculate the average coverage of all edges in each SNARL path and the number of edges with zero coverage; and count the paths that may be contained in each SNARL, the path direction, the reads aligned to the path, and the path coverage information.

[0079] The mutation path selection module is used to filter the optimal path and the second path based on reads information and path coverage information.

[0080] The path selection optimization module is used to optimize the optimal path and the second path;

[0081] The coarse path comparison module compares the optimized optimal path, the second path, and the reference path to obtain variation information.

[0082] Example 3

[0083] An electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the structural variation detection method based on third-generation sequencing data and pan-genome algorithm.

[0084] Example 4

[0085] A computer-readable storage medium storing a computer program, characterized in that, when executed by a processor, the computer program implements the structural variation detection method based on third-generation sequencing data and pan-genome algorithm.

[0086] The above description is only of the preferred embodiment of the present invention and should not be construed as limiting the scope of the claims. The present invention is not limited to the above embodiments, and variations in its specific structure are permitted. All variations made within the scope of the independent claims of the present invention are also within the scope of protection of the present invention.

[0087] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

Claims

1. A structural variation detection algorithm based on third-generation sequencing data and pan-genome, characterized in that, Includes the following steps: Obtain the pangenome map; Detect snarl structures in the pangenome map and extract reads corresponding to each snarl path in the snarl structure from the gam alignment file; Calculate the average coverage of all edges in each SNARL path, and the number of edges with zero coverage; The statistics include the path, path direction, reads aligned to the path, and path coverage information contained in each SNARL; Based on read information and path coverage information, the optimal path and the second path are selected: the paths are sorted according to the number of reads matched to the SNARL path, and the path with the highest number of reads is selected as the optimal path. If the number of reads is less than 5, the paths are sorted according to the calculated average coverage, and the path with the highest average coverage is selected as the optimal path. The remaining paths are traversed in turn, the number of unique reads matched to each path is calculated, and the path with the highest number of unique reads is selected as the second path. Optimize the optimal and second paths: Sort all paths in ascending order of path coverage; if coverage is the same, sort in ascending order of path length; if both coverage and length are the same, sort in descending order of path extraction from SNARL. Then optimize the optimal path. If the number of reads for the optimal path is less than a set value, select the last path after sorting as the optimal path. Then, based on the base coverage information obtained using the VG alignment tool, calculate the path with the highest average base coverage, path_m. If the path with the highest average base coverage, path_m, is unique and has a base coverage greater than twice that of all other paths, and the path with the highest average base coverage... If the average edge coverage of path_m is similar to that of the optimal path, then the optimal path is path_m with the largest average base coverage. Based on the selected optimal path, the path coverage of other paths is re-evaluated and updated. Then, the paths are sorted in ascending order of average edge coverage, ascending order of path length, and descending order of path sequence. After sorting, the paths are traversed from back to front. If it is the optimal path, it is skipped; otherwise, it is selected as the second path. Then, the SNARL paths are traversed again. If there is a path whose path length is greater than a set value and twice as long as the second path, and the average edge coverage of this path is comparable to that of the second path, then the original second path is discarded and this path is selected as the second path. The optimized optimal path and the second path are both compared with the reference path to obtain mutation information.

2. The structural variation detection algorithm based on third-generation sequencing data and pan-genome as described in claim 1, characterized in that, The optimal path is compared with the second path and the reference path to obtain mutation information, including the following steps: When the optimal path is the reference path, if the difference in the number of bases between the reference path and the second path is less than a set threshold, then the mutation in the SNARL region is a small mutation. If the path coverage of the optimal path is less than or equal to 2, the path coverage of the optimal path is greater than 4 times that of the second path, the number of reads of the optimal path is greater than 3 times that of the second path, and the number of edges with 0 coverage in the second path is more than 0.4 times the total number of edges, then there is no mutation in the SNARL region; otherwise, there is mutation, and the second path is retained. When the second path is the reference path, the optimal path is the mutated path, and mutations are normal. When neither the optimal path nor the second path is the reference path, compare the optimal path and the reference path. If the number of bases is not significantly different, then there is no structural variation.

3. The structural variation detection algorithm based on third-generation sequencing data and pan-genome as described in claim 2, characterized in that, If the number of reads for the optimal path and the second path is 0, the coverage of the optimal path and the second path is not 0, and the difference is not significant, then the SNARL region is considered to have a mutation, and the second path is retained. If the optimal path has only 2 nodes and the second path has at least 4 nodes, and the path coverage of the second path is greater than half that of the optimal path, then this SNARL region is considered to have a mutation, and the second path is retained.

4. The structural variation detection algorithm based on third-generation sequencing data and pan-genome as described in claim 3, characterized in that, When the second path is the reference path, the optimal path is the variant path. Under normal circumstances, there is a variant, which includes the following steps: if the difference in the number of bases between the optimal path and the reference path is less than a set threshold, then the variant in the SNARL region is a small variant; if the coverage of the second path is more than 4 times the coverage of the optimal path, then the variant is ignored; otherwise, it is considered that there is a variant here, and the optimal path is retained. If the coverage of the reference path is greater than 4 times that of the optimal path, the number of reads in the optimal path is less than 3, and the path coverage is less than 5, then there is a potential mutation in the optimal path, and the optimal path is retained. Then the second path and the reference path are compared. If the number of bases is not much different, then there is no structural mutation. If the coverage of the optimal path is more than 4 times that of the second path, the number of reads of the optimal path is more than 3 times that of the second path, and the number of edges with 0 coverage in the second path is more than 0.4 times the total number of edges, then there is no potential mutation in the second path; otherwise, it is considered that there is a second mutation and the second path is retained. If the length of the second path is more than four times that of the optimal path, and their coverage is comparable, then a second mutation is considered to exist, and the second path is retained.

5. The structural variation detection algorithm based on third-generation sequencing data and pan-genome as described in claim 1, characterized in that, Based on the preserved path and the reference path, specific variant information is extracted. The variant position is the last position of the first node in the reference path. The reference sequence is the node sequence found in the pangenome map file and assembled sequentially in either the forward or reverse direction according to the node orientation. The variant sequence is the node sequence found in the parameter variant path in the pangenome map file. It is then assembled sequentially in either the forward or reverse direction according to the node direction to obtain the structural variant VCF output file.

6. A structural variation detection system based on third-generation sequencing data and pan-genome, characterized in that, Includes a module for executing the structural variation detection algorithm of claim 1: The pan-genome map acquisition module is used to acquire pan-genome maps; The data preprocessing module is used to detect snarl structures in the pangenome map and extract reads corresponding to each snarl path in the snarl structure from the gam alignment file. Calculate the average coverage of all edges in each SNARL path, and the number of edges with zero coverage; Analyze the possible paths contained in each SNARL, the path direction, the reads aligned to the path, and the path coverage information; The mutation path selection module is used to filter the optimal path and the second path based on reads information and path coverage information. The path selection optimization module is used to optimize the optimal path and the second path; The coarse path comparison module is used to compare the optimized optimal path and the second path with the reference path to obtain variation information.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the structural variation detection algorithm based on third-generation sequencing data and pan-genome as described in any one of claims 1 to 5.

8. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the structural variation detection algorithm based on third-generation sequencing data and pan-genome as described in any one of claims 1 to 5.

Images