Structural variation detection algorithm, system and equipment based on third-generation sequencing data and generic genome and medium

Abstract

The invention discloses a structure variation detection algorithm, system and equipment based on three-generation sequencing data and a generic genome and a medium, and the detection algorithm comprises the following steps: detecting a snarl structure in a generic genome graph, and extracting reads corresponding to each snarl path from a gam comparison file; calculating the average coverage degree of all edges in each snarl path and the number of the edges with the coverage degree being 0; counting a path and a path direction which are possibly contained in each snarl, and reading information and path coverage information which are compared to the path; screening an optimal path and a second path according to the reads information and the path coverage information; and comparing the optimized optimal path and the second path with the reference path to obtain variation information. According to the method, the reads information corresponding to the path, the base coverage information of the path and the coverage information of the edge of the path are fused to serve as a basis for potential variation path selection, and the accuracy of third-generation sequencing data detection is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of gene sequence structural variation detection, and specifically relates to a structural variation detection algorithm, system, device and medium based on third-generation sequencing data and a pan-genome. Background Art

[0002] Variations in organisms generally include small variations and structural variations. Small variations include single nucleotide polymorphisms (SNPs) and small sequence insertions (INS) and sequence deletion variations (DEL). The latter two are generally collectively referred to as INDEL variations. SNPs mainly refer to DNA sequence diversity caused by single nucleotide variations at the genomic level. They are numerous and are important bases for studying genetic variations in human families and animal and plant strains. INDEL refers to the insertion or deletion of small fragment sequences at a certain position in the genome, and its length is usually below 50 bp. A structural variant (SV) is a genomic mutation involving 50 or more base pairs. SVs can generally be divided into 5 types, such as deletions, insertions, duplications, inversions, translocations. In addition, there are complex structural variations composed of these types. Although structural variations only account for a small part of all variations, they are usually relatively large. Existing analyses show that they have a high impact on gene expression. Structural variations can affect regulatory elements such as gene promoters and enhancers, or affect the structure of chromosomes, thereby changing the gene expression level. In some cases, structural variations can lead to the fusion of two genes, and this gene fusion may produce new fusion proteins with new functions or lose their original functions. Certain structural variations are also related to the occurrence of genetic diseases, such as certain types of cancer, hereditary heart diseases, etc. In addition, structural variations also play an important role in the process of species evolution. By changing the genomic structure, they promote the adaptability and diversity of species. Therefore, the study of structural variations is of great significance to human medicine and genetics.

[0003] The research on gene mutations depends on the technological progress of molecular biology. Undoubtedly, the DNA sequencing technology that has dominated genomics research in the past is the Illumina platform technology, which is a technology for generating short-read sequencing data. It can produce highly accurate (over 99.9%) sequencing data at low cost. However, the short-read sequencing technology has a major drawback in dealing with larger structural variations and genome assembly, that is, the limited length of reads, which makes it impossible for short reads to detect many types of structural variations. Therefore, long-read sequencing technologies have developed rapidly in recent years, mainly including ONT sequencing and PacBio sequencing, which can produce longer reads of thousands or tens of thousands of bases. Analyzing the genome with these reads can reveal more previously undiscovered genetic information.

[0004] A common method for detecting mutations is to first align the sequencing data to a reference genome and then search for mutations based on the alignment inconsistencies, such as the genotyping models of SVTyper and Delly. However, traditional reference genomes are not sensitive to mutations and have many drawbacks. For example, only one version of the mutation is selected for embedding at each mutation position, so they cannot well represent the genetic diversity. This will cause some sequencing data containing mutant genes to not align to the reference genome, and the subsequent mutation signals cannot be detected. To address the deficiencies of linear reference genomes, the concept of pan-genome has been proposed recently. The pan-genome can integrate the gene information of multiple samples and retain multiple mutant sequences at each mutation site. Compared with linear reference genes, the pan-genome graph can encode more mutation information to alleviate the drawbacks of linear reference genomes and improve the detection and genotyping of mutations. The pan-genome is usually represented by a directed graph, where nodes represent gene sequences and edges represent the connection relationships between sequences. Therefore, the alignment of sequencing data to the linear reference genome has also been extended to the alignment of sequencing data to the paths in the pan-genome graph.

[0005] At present, there are not many tools for directly detecting and typing structural variations based on a pan-genome graph. Typical representatives include VG and GraphTyper. However, in VG, two paths with the highest average coverage are directly selected as potential variant paths (the coverage of the corresponding path can be obtained according to the alignment information from the sequencing data to the path), and then the different situations of these two paths are analyzed. The method for detecting and typing structural variations is relatively simple. If the selected paths are incorrect, the subsequent analysis error will be very large. GraphTyper can only build a pan-genome graph for analysis in a small region of the genome and cannot be directly applied to the pan-genome graph of the entire chromosome. Moreover, it is only applicable to second-generation sequencing data with short sequencing lengths and is not applicable to third-generation sequencing data. With the increasingly wide application of third-generation sequencing data, the disadvantages of GraphTyper are becoming increasingly obvious. Therefore, it is necessary to develop a structural variation detection algorithm based on third-generation sequencing data and the structure of the pan-genome graph to improve the accuracy of structural variation detection and provide a reliable basis for subsequent analysis. Summary of the Invention

[0006] To overcome the problem of low accuracy in detecting structural variations in the prior art, the object of the present invention is to provide a structural variation detection algorithm, system, device, and medium based on third-generation sequencing data and a pan-genome. The detection results of this algorithm are more reliable and real, improving the accuracy of structural variation detection.

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

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

[0009] Obtain a pan-genome graph;

[0010] Detect the snarl structure in the pan-genome graph, and extract the reads corresponding to each snarl path in the snarl structure from the gam alignment file; calculate the average coverage size of all edges in each snarl path and the number of edges with a coverage of 0; count the paths, path directions, reads information aligned to the paths, and path coverage information included in each snarl;

[0011] According to the reads information and path coverage information, screen the optimal path and the second path;

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

[0013] Compare the optimized optimal path and the second path with the reference path to obtain variant information.

[0014] Further, according to the reads information and the path coverage information, screening the optimal path and the second path includes the following steps: Sort the paths according to the number of reads aligned to the snarl path, select the path with the highest number of reads as the optimal path. If the number of reads is less than 5, then sort according to the calculated average coverage of the path, and select the path with the highest average coverage as the optimal path; Traverse the remaining paths in turn, calculate the number of unique reads aligned to each path, and select the path with the highest number of unique reads as the second path.

[0015] Further, optimizing the optimal path and the second path includes the following steps:

[0016] Sort all paths in ascending order of path coverage. If the coverages are the same, sort them in ascending order of path length. If the coverages and lengths of the paths are the same, sort them in descending order of the extraction order of the paths in the snarl, and then optimize the optimal path. If the number of reads of the optimal path is less than the set value, select the last one of the sorted paths as the optimal path; Then, according to the base coverage information obtained by using the vg alignment tool, calculate the path path_m with the largest average base coverage. If the base coverage of the path path_m with the largest average base coverage is unique and more than twice that of all other paths, and the average edge coverage of the path path_m with the largest average base coverage is similar to that of the optimal path, then the optimal path is the path path_m with the largest average base coverage; According to the selected optimal path, re-evaluate the path coverage of other paths and update the coverage; Then sort in ascending order of average edge coverage, ascending order of path length, and descending order of path order. Traverse the paths from back to front after sorting. Skip if it is the optimal path, otherwise select it as the second path; Then traverse the snarl path again. If there is a path whose path length is greater than the set value and 2 times larger than the path length of the second path, and the average edge coverage of this path is equivalent to that of the second path, then discard the original second path and use this path as the second path.

[0017] Further, comparing the optimal path and the second path with the reference path to obtain mutation information includes 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 the 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 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 proportion of edges with a coverage of 0 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 a mutation, and the second path is retained;

[0019] When the second path is the reference path, the optimal path is the variant path, and there is a variant under normal circumstances;

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

[0021] Furthermore, if the read counts of the optimal path and the second path are 0, the coverages of the optimal path and the second path are not 0, and the difference is not significant, then it is considered that there is a variant in the snarl region, and the second path is retained;

[0022] If the optimal path has only 2 nodes and the number of nodes in the second path is not less than 4, and the path coverage of the second path is greater than half of the optimal path, then it is considered that there is a variant in this snarl region, and the second path is retained at this time.

[0023] Furthermore, when the second path is the reference path, the optimal path is the variant path, and there is a variant under normal circumstances, including the following steps: If the difference in the number of bases between the optimal path and the reference path is less than the 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 more than 4 times the optimal path, the read count of the optimal path is less than 3, and the path coverage is less than 5, otherwise there is a potential variant in the optimal path, and the optimal path is retained; Then compare the second path and the reference path. If the difference in the number of bases is not significant, then there is no SV;

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

[0026] If the length of the second path is more than 4 times the optimal path and their coverages are comparable, then it is considered that there is a second variant, and the second path is retained.

[0027] Furthermore, according to the retained 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 to find the node sequence in the reference path in the pan-genome graph file and splice it forward or backward in turn according to the node direction; The variant sequence is to find the node sequence in the variant path in the pan-genome graph file and splice it forward or backward in turn according to the node direction to obtain the structural variant VCF output file.

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

[0029] A pan-genome graph acquisition module for acquiring a pan-genome graph;

[0030] A data preprocessing module for detecting snarl structures in the pan-genome graph and extracting reads corresponding to each snarl path in the snarl structure from the gam alignment file; calculating the average coverage size of all edges in each snarl path and the number of edges with a coverage of 0; counting the paths, path directions, reads information aligned to the paths, and path coverage information that may be included in each snarl;

[0031] A variant path selection module for screening the optimal path and the second path according to the reads information and the path coverage information;

[0032] A path selection optimization module for optimizing the optimal path and the second path;

[0033] A path rough alignment module for comparing both the optimized optimal path and the second path with a reference path to obtain variant information.

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

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

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

[0037] The present invention utilizes the reads information corresponding to the paths, the base coverage information of the paths, and the coverage information of the edges of the paths, and combines the three as the basis for selecting potential variant paths, which is more reliable and accurate as a whole, overcoming the problem in the prior art that the analysis of snarls in the pan-genome graph by VG is relatively simple, and directly selects the two paths with the highest base coverage as candidate paths when selecting paths in the snarl, but sometimes the base coverage is unreliable, resulting in these two paths being unreliable and subsequent analysis being unreliable.

[0038] After selecting the optimal path, the present invention eliminates the influence of the optimal path on other paths, including the influence on reads and the influence on the coverage of edges, and then selects the second path, making the result more reliable and true. It overcomes the problem in the prior art that neither VG nor GraphTyper eliminates the influence of the optimal path on other paths in path selection, which may lead to an unrealistically high overall coverage of other paths and affect the accuracy of subsequent results.

[0039] Furthermore, after selecting the optimal path and the second path in the present invention, if the reads alignment situation is not good, the logic of path selection optimization will continue. The optimal path is compared with the path with the highest base coverage, and the second path is compared with paths whose length is more than twice that of it and whose coverage is comparable. The final candidate variant paths to be retained are determined according to specific strategies. The variant paths retained in this way will be more reliable and the subsequent results will be more accurate. It overcomes the problem in the prior art that after VG and GraphTyper select the optimal path and the second path, there is no logic of path selection optimization. However, there may be some paths with slightly lower coverage but better overall effects, and these paths should be more selected as candidate variant paths, thus affecting the result accuracy.

[0040] Furthermore, the present invention comprehensively considers the difference in the number of reads of paths and the multiple of coverage, and also considers the proportion of edges with coverage of 0 in candidate variant paths, making full use of the available information of candidate variant paths to make 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 according to the multiple difference of the base coverage of the paths, with a simple method and without considering complex situations.

[0041] Furthermore, the present invention also takes into account some complex situations, such as many reads being aligned to snarl paths that are not counted, resulting in very few reads corresponding to the counted paths, and the reads of the second path being very few but the length of the second path being much longer than the optimal path and the coverage being comparable. There are corresponding strategies to solve these problems, making the result more accurate. Brief Description of the Drawings

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

[0043] Figure 2 It is a flowchart of the variant path selection module;

[0044] Figure 3 It is a schematic diagram of the influence of the optimal path in snarl on the reads of other paths;

[0045] Figure 4 It is a flowchart of the path selection optimization module;

[0046] Figure 5 It is a schematic diagram of the influence of the optimal path in the snarl on the coverage of other paths;

[0047] Figure 6 It is a general box plot of the F1 scores of three structure variation detection tools;

[0048] Figure 7 It is a general radar chart of the results of three tools;

[0049] Figure 8 It is a bar chart of the f1 scores of three tools on different chromosomes;

[0050] Figure 9 It is the F1 score results of different tools for insertional variations of different sizes;

[0051] Figure 10 It is a schematic diagram of the structure of a structure variation detection system based on third-generation sequencing data and a pan-genome. Detailed implementation manners

[0052] For the convenience of understanding the present invention, the present invention will be described more comprehensively below with reference to the relevant drawings. The preferred embodiments of the present invention are shown in the drawings. However, the present invention can be implemented in many different forms and is not limited to the embodiments described herein. On the contrary, these embodiments are provided to make the disclosure of the present invention more thorough and comprehensive.

[0053] See Figure 1 , Figure 1 It is a schematic flowchart of the structure variation detection algorithm of the present invention based on third-generation sequencing data and a pan-genome, mainly including four steps: data preprocessing, variant path selection, path selection optimization, and path rough alignment.

[0054] Obtain the pan-genome graph;

[0055] Data preprocessing: Detect the snarl structure in the pan-genome graph (chrom.vg file), and extract the reads corresponding to the snarl paths in each snarl structure from the gam alignment file. The snarl structure in the figure is such a structure that removing all the 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 into a connected graph containing the nodes at both ends. The schematic diagram of the snarl structure is as Figure 1As shown in the data preprocessing module, like "bubbles" in the genome, the position of each snarl is a potential position for variation. Then, the average coverage size of all edges in each snarl path and the number of edges with a coverage of 0 are recalculated. When calculating the path coverage here, only the average coverage of the edges in the path is calculated, without considering the coverage of the bases in the nodes. Because according to the characteristics of the pan-genome graph, the coverage of the edges is more important. A high edge coverage indicates a high coverage of the path where the edge is located. In addition, the base sequences in some nodes are very long and contain repetitive sequences. The alignment of the middle part of the base sequences may not be very good, resulting in a situation where the coverage of two edges related to the node may be very high but the coverage of the bases inside the node is low, leading to a low overall base coverage of the path and affecting the accuracy of the results. Here, two types of coverage information appear. One is the average coverage of all bases in the path calculated using the vg alignment tool, and the other is the average coverage of all edges in the path recalculated. After that, unless otherwise specified, the path coverage that appears refers to the average coverage of all edges in the path. Finally, the paths, path directions, reads information aligned to the paths, and path coverage information that may be included in each snarl are counted and used for the next analysis.

[0056] Variant path selection: Two possible variant paths are initially selected, namely the optimal path and the second path. Since there are a very large number of paths in the snarls involved in SV, if pairwise combinations are made and the likelihood is calculated like in GraphTyper, the number of path combinations will be extremely large. So, the strategy is changed here. First, two paths are selected, and then they are compared with the reference path to detect the variation. Specifically, as Figure 2As shown in the figure, first find the possible optimal path based on the reads information and path coverage information: first sort the paths according to the number of reads aligned 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, it means that the alignment of the sequencing reads in this snarl region is not good. At this time, sort the paths according to the calculated average coverage, and select the path with the highest average coverage as the optimal path. Since the reads information is more reliable than the average coverage of the path (how many reads are aligned to the path is how many reads, which is more accurate. If a small part of the read sequence is aligned to a small part of the base sequence in the snarl path, and the other parts are not aligned, the coverage of this small part will also be counted when the coverage is counted, resulting in inaccurate results), so first of all, the reads aligned to the snarl path are taken as the basis. After the optimal path is obtained, the remaining paths are traversed in turn, the number of unique reads aligned to each path is calculated, and the path with the highest number of unique reads is selected as the second path. What is calculated here is the number of unique reads for each path rather than the number of all reads. Unique reads refer to those reads that appear in the snarl path but not on the optimal path. The number of unique reads is calculated because there is a situation where there are multiple paths that are very similar to the optimal path, with only a very small number of nodes that are different, and the coverage of these very few nodes is very low. If you simply select based on the number of reads on the path, these similar paths will be retained, but 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 of the path is used for judgment to eliminate the influence of the optimal path on other paths as much as possible, so that the overall preference is to select paths that are less similar to the optimal path. For example Figure 3 As shown in the figure, path1 corresponds to 8 reads, path2 corresponds to 5 reads, and path3 corresponds to 3 reads, but 4 of the 5 reads of path2 are repeated with path1. After path1 is selected as the optimal path, the number of unique reads of path2 is only 1, while the number of unique reads of path3 is still 3, so path3 should be selected as the second path. After the second path is selected, if the number of reads of the optimal path or the second path is 0, or 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 means that the alignment of the reads in this snarl region is not good, and then the path optimization logic is entered to use the path coverage for further analysis.

[0057] Path selection optimization: When the read alignment in the snarl region is poor, the path selection strategy is optimized. Figure 2 As shown, the optimal path and the second path are optimized respectively. Specifically, Figure 4 As shown, all paths are sorted in ascending order according to path coverage. If the coverage is the same, they are sorted in ascending order according to path length, because the longer the path, the more difficult it is to achieve the same coverage and the higher the reliability. If the coverage and length of the path are the same, they are sorted in descending order according to the order of path extraction in snarl. Because the coverage and length of the path are sorted in ascending order, the better the overall rightward direction is. When all conditions are the same, the original front path tends to be selected, so the path with a small index is on the right. Then optimize the optimal path. If the number of reads of the optimal path selected previously is very small (less than the set value 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 the two methods is the same, which is more reliable). Then find the average base coverage obtained by the vg alignment tool and calculate the path path_m with the largest average base coverage. If the base coverage of the path path_m with the largest average base coverage is unique and greater than twice that of all other paths (including the optimal path), and the average coverage of the edges of the path path_m with the largest average base coverage is similar to the optimal path (that is, the average coverage of the edges of the path path_m with the largest average base coverage is greater than 0.75 times that of the optimal path), then the optimal path is replaced with 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 to exclude the impact of the optimal path on them. The impact of the optimal path on the reads of other paths has been excluded before, but the alignment of the reads in the snarl region is very poor at this time, which is not suitable for continued analysis with reads, but needs to be analyzed with path coverage, so it is necessary to continue to exclude the impact of the optimal path on the path coverage of other paths. This can prevent the nodes of some paths from being completely included in the optimal path. For those edges that are different from the optimal path, their coverage is originally 0, but due to the existence of the optimal path, the coverage of other edges may be relatively high, which will cause the overall coverage of these paths to be unrealistically high. 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 relatively high coverage of path2 is largely due to the influence of path1. The coverages of the two edges that are different from path1 in the middle are very low. Therefore, the influence of the duplicate edges of path1 should be excluded, and its coverage should be updated to 2. Finally, path3 is selected as the second path. Therefore, it is necessary to adjust the coverages of these paths to ensure their authenticity and accuracy. The method of updating the coverage is also relatively simple: if the edges in other paths are exactly the same as those in the optimal path, this edge needs to be ignored at this time to avoid the influence of the optimal path. After the coverage is updated, the paths are sorted again in ascending order of the average coverage of the edges, ascending order of the path length, and descending order of the path sequence. After sorting, traverse the paths from back to front. If it is the optimal path, skip it; otherwise, select it as the second path. Then traverse the snarl paths again. If there is a path whose path length is relatively large (by default, greater than the set value of 15 nodes) and is twice as large as the path length of the second path, and the average coverage of the edges of this path is comparable to that of the second path, then discard the original second path and keep this path as the second path. Because for paths with a larger length, it is more difficult to achieve the same path coverage as paths with a shorter length. If the length difference is more than twice and the path coverages are comparable, the longer path is more accurate and reliable and is more likely to be the true potential variant path. Thus, a relatively reliable potential variant path is obtained.

[0058] Rough path comparison: Further analyze the selected optimal path and the second path, and compare them with the reference path to obtain the final variant information. First, calculate the total base lengths of the reference path and the two selected paths (the optimal path and the second path) so as to judge whether there is an SV and the type of SV according to the length difference. Then, as Figure 1 shown in the path rough comparison module, analyze in three cases:

[0059] 1. The optimal path is the reference path: If the difference in the number of bases between the reference path (i.e., the optimal path) and the second path is not significant, less than 50 bases (i.e., the set threshold), then the variation in this snarl region is considered a small variation and is ignored. One of the following situations will cause the second path to be ignored, and it is considered that there is no variation in this snarl region: 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 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 proportion of edges with a coverage of 0 in the second path is more than 0.4 times the total number of edges. Otherwise, it is considered that there is a variation and the second path is retained. In addition, there are two special situations, namely: (1) If the number of reads of 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 indicates that there is a snarl path that has not been counted. Most reads are aligned to the uncounted snarl path, affecting the coverage of the counted path. Because the snarl region may be large and complex, containing many nodes and potential paths, it may not be completely counted. Here, the reference path must be counted, so the uncounted path must be a variant path, and many reads are aligned to the uncounted variant path. Therefore, it is considered that there is a variation in this snarl region at this time, and the second path is retained. (2) Since insertion mutations are more difficult to detect than deletion mutations, the following optimization is carried out for insertion mutations: If the optimal path has only 2 nodes and the number of nodes in the second path is not less than 4 (indicating an obvious insertion mutation), and the path coverage of the second path is greater than half of the optimal path, it is considered that there is a variation in this snarl region, and the second path is retained at this time. The influence of reads is no longer considered here because in insertion mutations, if the inserted sequence is very long, it is difficult for reads to cover and 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 there must be a variation under normal circumstances. If the difference in the number of bases between the optimal path and the reference path (i.e., the second path) is not significant, that is, less than 50 bases, then the variation in this snarl region is considered a small variation and is ignored. If the coverage of the second path is more than 4 times the coverage of the optimal path, then the optimal path may be selected according to the number of reads, but the overall number of reads is very small, and the variation is unreliable. At this time, the variation is ignored and no longer considered. Otherwise, it is considered that there is a variation here, and the optimal path is retained.

[0061] 3. Neither the optimal path nor the second path is the reference path: In this case, both the optimal path and the second path need to be compared with the reference path. First, compare the optimal path with the reference path. Similarly, if the difference in the number of bases is small (less than 50 bases), it is considered that there is no SV. If any of the following situations exist, it is considered that there is no potential variation in the optimal path: the coverage of the reference path is more than 4 times that of the optimal path, the number of reads in the optimal path (less than 3) and the path coverage (less than 5) are both very small. Otherwise, it is considered that there is a potential variation in the optimal path, and the optimal path is retained. Then compare the second path with the reference path. Similarly, if the difference in the number of bases is small (less than 50 bases), it is considered that there is no SV. If any of the following situations exist, it is considered that there is no potential variation in the second path: the coverage of the optimal path is more than 4 times that of the second path, the number of reads in the optimal path is more than 3 times that of the second path, and the proportion of edges with a coverage of 0 in the second path is more than 0.4 times the total number of edges. Otherwise, it is considered that there is a second variation, and the second path is retained. There is also a special case here. If the length of the second path is more than 4 times that of the optimal path, but their coverages are comparable (the coverage of the second path is more than 0.75 times the coverage of the optimal path), it is considered that there is a second variation, and the second path is retained without filtering with reads information. Because the number of nodes in the second path and the optimal path differ a lot, and the edge duplication degree is very small, indicating that the association between the two paths is very small, and the coverage of the second path has little association with the optimal path, there is likely a variation.

[0062] Finally, according to the retained paths and the reference path, specific variant information is extracted. If directly using sequence alignment algorithms for alignment, the alignment results may be very detailed, and a large SV may be split into multiple SVs and small variants. Therefore, here, a rough alignment is performed using the nodes in the path. Among them, the variant position is the last position of the first node in the reference path; the reference sequence is to find the node sequence in the reference path in the pan-genome graph file and splice it forward (normal base sequence, such as ACGT) or backward (reverse complement sequence of the base sequence, such as TGCA) according to the node direction in turn; the variant sequence is to find the node sequence in the variant path in the pan-genome graph file and splice it forward or backward according to the node direction in turn. Thus, the final structural variant VCF output file is obtained.

[0063] To verify the accuracy of the algorithm in the present invention, experimental verification was carried out on real samples. The sample selected for this experiment was the HG002 sample. The third-generation sequencing data of HG002 was used to run Pan-SV and VG of the present invention, the second-generation sequencing data of HG002 was used to run GraphTyper, and the Truvari tool was used for structural variant accuracy detection.

[0064] The experimental results are shown in the attached figure. Figure 6The overall box plot of the F1 scores of three structural variant detection tools. As can be seen from the figure, whether in high-confidence regions or 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 of the three tools. As can be seen, Pan-SV is closest to the upper right region and has the best effect.

[0066] Figure 8 The bar chart of the f1 scores of the three tools on different chromosomes. As can be seen, the results of Pan-SV (the algorithm in the present invention) are the best on all chromosomes.

[0067] Figure 9 The F1 score results of different tools for insertions of different sizes. As can be seen, for insertions of different size ranges in the figure, the results of Pan-SV are the best. Moreover, as the size of the insertion variant increases, the decline of Pan-SV is smaller than that of VG. Although the effect of GraphTyper becomes better as the size of the insertion variant increases, its overall F1 score is very low.

[0068] Advantages of the present invention:

[0069] In the prior art, the analysis of snarl in the pan-genome graph by VG 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, resulting in these two paths being unreliable and subsequent analysis being unreliable. In contrast, the present invention utilizes the read information corresponding to the paths, the base coverage information of the paths, and the coverage information of the edges of the paths, and combines the three as the basis for selecting potential variant paths, which is more reliable and accurate overall.

[0070] In the prior art, neither VG nor GraphTyper excludes the influence of the optimal path on other paths in path selection, which may lead to an unrealistically high overall coverage of other paths and affect the accuracy of subsequent results. In contrast, after selecting the optimal path, the present invention excludes the influence of the optimal path on other paths, including the influence of reads and the coverage of edges, and then selects the second path, making the results more reliable and realistic.

[0071] In the prior art, after VG and GraphTyper select the optimal path and the second path, there is no logic for optimizing path selection. However, there may be some paths with slightly lower coverage but better overall effects, and these paths should be more preferably selected as candidate variant paths. In the present invention, after the optimal path and the second path are selected, if the reads alignment situation is not good, the logic for optimizing path selection will continue. The optimal path is compared with the path having the highest base coverage, and the second path is compared with the path whose length is more than twice that of the second path and whose coverage is comparable. The final candidate variant path to be retained is determined according to specific strategies. The variant paths retained in this way will be more reliable, and the subsequent results will be more accurate.

[0072] In the prior art, after VG selects the optimal path and the candidate path, it directly determines whether there is an SV according to the multiple difference in the base coverage of the paths. The method is simple and does not consider complex situations. The present invention comprehensively considers the number of reads and the multiple difference in coverage of the paths, and also considers the proportion of the edges with coverage of 0 in the candidate variant paths, making full use of the available information of the candidate variant paths to make the results more reliable. In addition, the present invention also considers some complex situations, such as many reads being aligned to the snarl paths that are not counted, resulting in very few reads corresponding to the counted paths, and the second path having very few reads but the length of the second path being much longer than that of the optimal path and the coverage being comparable. There are corresponding strategies to solve these problems, making the results more accurate.

[0073] In the present invention, multiple types of information of the paths in the pan-genome graph are jointly used for selecting candidate variant paths: using the reads information corresponding to the paths, the base coverage information of the paths, and the coverage information of the edges of the paths, and using the characteristics of the three types of information together as the basis for selecting potential variant paths. When selecting paths, the influence of the optimal path on other paths is excluded, improving the reliability of path selection: after the optimal path is selected, the influence of the optimal path on other paths, including the influence of reads and the influence of the edge coverage, will be excluded, and then the second path is selected, making the results more reliable and real. After the two candidate variant paths, namely the optimal path and the second path, are selected, the logic for optimizing path selection continues to retain more reliable paths: after the optimal path and the second path are selected, if the reads alignment situation is not good, the optimization of path selection will continue. The optimal path is compared with the path having the highest base coverage, and the second path is compared with the path whose length is more than twice that of the second path and whose coverage is comparable. The final candidate variant path to be retained is determined according to specific strategies. The present invention makes full use of the available information of the candidate variant paths to determine whether there is an SV in the paths and considers some complex situations.

[0074] The present invention comprehensively considers the number of reads of a path and the multiple difference in coverage, and also considers the proportion of edges with a coverage of 0 in the candidate variant path, making full use of the available information of the candidate variant path to make the result more reliable. In addition, the present invention also takes into account some complex situations, such as a large number of reads being aligned to a snarl path that is not counted, resulting in very few reads corresponding to the counted path, and the second path having very few reads but the length of the second path being much greater than that of the optimal path and the coverage being comparable. There are corresponding strategies to solve these problems, making the result more accurate.

[0075] Example 2

[0076] See Figure 10 , a structural variant detection system based on third-generation sequencing data and a pan-genome, including: a pan-genome graph acquisition module, a data preprocessing module, a variant path selection module, a path selection optimization module, and a path rough alignment module.

[0077] The pan-genome graph acquisition module is used to acquire a pan-genome graph;

[0078] The data preprocessing module is used to detect the snarl structure in the pan-genome graph and extract the reads corresponding to each snarl path from the gam alignment file; calculate the average coverage size of all edges in each snarl path and the number of edges with a coverage of 0; count the paths, path directions, read information aligned to the paths, and path coverage information that may be included in each snarl;

[0079] The variant path selection module is used to screen the optimal path and the second path according to the read information and the path coverage information;

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

[0081] The path rough alignment module is used to compare the optimized optimal path and the second path with the reference path respectively to obtain variant 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. When the processor executes the computer program, the structural variant detection method based on third-generation sequencing data and a pan-genome algorithm is implemented.

[0084] Example 4

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

[0086] The above description is only for the best embodiments of the present invention, but it should not be construed as a limitation of the claims. The present invention is not limited to the above embodiments, and its specific structure is allowed to vary. Any changes made within the protection scope of the independent claims of the present invention are within the protection scope of the present invention.

[0087] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the technical field to which the present invention belongs. The terms used in the description of the present invention herein are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The term "and / or" used herein includes any and all combinations of one or more of the related listed items.

Claims

1. A structural variation detection algorithm based on third-generation sequencing data and pan-genome, characterized in that: The following steps are involved: Obtain pan-genome maps; Detect the snarl structure in the pan-genome graph and extract the reads corresponding to each snarl path in the snarl structure from the gam alignment file; Calculate the average coverage size of all edges in each snarl path and the number of edges with coverage of 0; Statistics of the paths, path directions, reads mapped to the paths, and path coverage information contained in each snarl; Based on the reads information and pathway coverage information, the optimal pathway and the second pathway are selected; Optimizing the optimal path and the second path; The optimized best path and the second path are compared with the reference path to obtain variation information.

2. The structural variation detection algorithm based on third-generation sequencing data and pan-genome according to claim 1, characterized in that: According to the reads information and path coverage information, the optimal path and the second path are screened, including the following steps: sorting the paths according to the number of reads mapped to the snarl path, selecting the path with the highest number of reads as the optimal path, if the number of reads is less than 5, sorting the paths according to the calculated average coverage, and selecting the path with the highest average coverage as the optimal path; traversing the remaining paths in turn, calculating the number of unique reads mapped to each path, and selecting the path with the highest number of unique reads as the second path.

3. The structural variation detection algorithm based on third-generation sequencing data and pan-genome according to claim 1, characterized in that: Optimizing the optimal path and the second path includes the following steps: All paths are sorted in ascending order according to path coverage. If the coverage is the same, they are sorted in ascending order according to path length. If the coverage and length of the paths are the same, they are sorted in descending order according to the order of path extraction in snarl, and then the optimal path is optimized. If the number of reads of the optimal path is less than the set value, the last path after sorting is selected as the optimal path; then, based on the base coverage information obtained by the vg alignment tool, the path path_m with the largest average base coverage is calculated. If the base coverage of the path path_m with the largest average base coverage is unique and greater than twice that of all other paths, and the edges of the path path_m with the largest average base coverage are If the average coverage of is close to the optimal path, the optimal path is the path path_m with the largest average base coverage; based on the selected optimal path, re-evaluate the path coverage of other paths and update the coverage; then sort them in ascending order of average edge coverage, ascending order of path length, and descending order of path order, and traverse the path from back to front after sorting. If it is the optimal path, skip it, otherwise use it as the second path; then traverse the snarl path again, if there is a path whose path length is greater than the set value and is twice as long as the path length of the second path, and the average edge coverage of this path is equivalent to that of the second path, then discard the original second path and use this path as the second path.

4. The structural variation detection algorithm based on third-generation sequencing data and pan-genome according to claim 1, characterized in that: Comparing the optimal path with the second path and the reference path to obtain variation information includes the following steps: When the optimal path is the reference path, if the difference between the number of bases of the reference path and the second path is less than the set threshold, the variation in the snarl region is a small variation. 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 edges with a coverage of 0 in the second path account for more than 0.4 times the total number of edges, then there is no variation in the snarl region. Otherwise, there is variation and the second path is retained. When the second path is the reference path, the optimal path is the variant path, and there is a variation under normal circumstances; When neither the optimal path nor the second path is the reference path, the optimal path and the reference path are compared. If the number of bases is not much different, there is no SV.

5. The structural variation detection algorithm based on third-generation sequencing data and pan-genome according to claim 4, characterized in that: If the number of reads of 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 large, it is considered that the snarl region has a mutation and the second path is retained; If the optimal path has only 2 nodes and the number of nodes in the second path is not less than 4, and the path coverage of the second path is greater than half of the optimal path, it is considered that there is a variation in this snarl region, and the second path is retained.

6. The structural variation detection algorithm based on third-generation sequencing data and pan-genome according to claim 5, characterized in that: When the second path is the reference path, the optimal path is the variant path, and normally there is a variant, including the following steps: if the difference between the number of bases in the optimal path and the reference path is less than the set threshold, the variant in the snarl region is a small variant; if the coverage of the second path is greater 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 pathway is greater than 4 times that of the optimal pathway, the number of reads in the optimal pathway is less than 3, and the pathway coverage is less than 5, otherwise there is a potential variation in the optimal pathway, and the optimal pathway is retained; then the second pathway is compared with the reference pathway, and if the number of bases is not much different, there is no SV; If the 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, and the number of edges with 0 coverage in the second path accounts for more than 0.4 times the total number of edges, then there is no potential variation in the second path. Otherwise, it is considered that there is a second variation and the second path is retained; If the length of the second path is greater than 4 times that of the optimal path and the coverage of the two is equivalent, it is considered that a second variation exists and the second path is retained.

7. The structural variation detection algorithm based on third-generation sequencing data and pan-genome according to claim 1, characterized in that: According to the retained 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 in the reference path found in the pan-genome map file, and is spliced ​​forward or backward in sequence according to the node direction. The variant sequence is the node sequence found in the reference variant path in the pan-genome map file, which is spliced ​​forward or backward in sequence according to the node direction to obtain the structural variant VCF output file.

8. A structural variation detection system based on third-generation sequencing data and pan-genome, characterized in that: include: A pan-genome map acquisition module, used to acquire a pan-genome map; Data preprocessing module, used to detect the snarl structure in the pan-genome graph and extract the reads corresponding to each snarl path in the snarl structure from the gam alignment file; Calculate the average coverage size of all edges in each snarl path and the number of edges with coverage of 0; Count the paths, path directions, reads information aligned to the paths, and path coverage information that may be contained in each snarl; The variation pathway selection module is used to select the optimal pathway and the second pathway based on the reads information and pathway coverage information; A path selection optimization module, used to optimize the best path and the second path; The path rough comparison module is used to compare the optimized best path and the second path with the reference path to obtain variation information.

9. 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, the method for detecting structural variations based on third-generation sequencing data and a pan-genome algorithm as described in any one of claims 1 to 7 is implemented.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the method for detecting structural variations based on third-generation sequencing data and a pan-genome algorithm as described in any one of claims 1 to 7 is implemented.

Images