A method for identifying whole genome structural variations based on third-generation sequencing

Abstract

The embodiment of the application discloses a whole genome structural variation identification method based on third-generation sequencing, comprising: aligning the data to be sequenced with a preset reference genome, and sorting the preset alignment file to construct an index; for each sequencing read, the base alignment condition is analyzed, SV identification is carried out, and whether the signal of the variable DNA fragment is included is identified; the error and overlapping variable DNA fragment signal is corrected; the corrected sequence is output, compared with the preset reference genome, SV identification is carried out, and the identification result is output; and the data to be sequenced is typed according to the identification result. The application can improve the detection accuracy and sensitivity, and the alignment boundary for high-error third-generation sequencing sequences; improve the identification accuracy of complex or longer variable DNA fragments; accurately distinguish adjacent variable DNA fragments with close distances; and accurately type the variable DNA fragments.

Description

Technical Field

[0001] This invention relates to the field of genome variation identification technology, and in particular to a method for identifying whole-genome structural variations based on third-generation sequencing. Background Technology

[0002] Structural variation (SV) generally refers to DNA fragment variations longer than 50 bp. Based on the type of variation, it can be classified as deletion, duplication, insertion, inversion, translocation, etc. Structural variations affect gene transcription and translation in various ways, thereby causing various hereditary diseases. When structural variations occur in the coding region of a gene, they alter gene transcription and translation; when structural variations occur in non-coding regions, they affect the regulatory functions of gene expression regulatory elements through positional effects. With the rapid development of sequencing technology and the continuous reduction in sequencing costs, more and more structural variations have been found to be associated with human hereditary diseases and even cancer. For example, Down syndrome caused by trisomy 21, Cri-du-chat syndrome caused by deletion on the short arm of chromosome 5, and learning disabilities caused by deletion variations at 17q21.31.

[0003] Currently, SV detection methods based on second-generation short-read sequencing can be categorized into read-pair methods, read-depth methods, split-read methods, and sequence assembly methods, or combinations of these methods. However, due to the short read lengths of second-generation data, many structural variants are much longer than the read length, thus failing to effectively detect all structural variants and resulting in a high rate of false positives.

[0004] The rapid development of third-generation long-read sequencing in the past two years (mainly including PacBio single-molecule real-time sequencing and Oxford Nanopore Sequencing) has provided a possibility for improving the detection rate and accuracy of structural variants. Currently, software used for identifying third-generation sequencing data includes pbSV, cuteSV, and Sniffles. The typical process involves aligning sequencing data to the genome, identifying signals containing variant DNA fragments (SVs), clustering all SV signals, and finally averaging the SV signals for each cluster. This method is significantly affected by sequencing and alignment errors, especially inaccurate identification of SV boundaries; it cannot identify complex SVs; it cannot genotype SVs of the same type but different lengths; it cannot accurately distinguish closely spaced adjacent SVs; and it cannot accurately identify SVs exceeding the length of the sequencing data. Therefore, this technology requires further development and improvement. Summary of the Invention

[0005] In view of this, the purpose of this invention is to provide a method for identifying whole-genome structural variations based on third-generation sequencing, which can improve the detection accuracy and sensitivity, as well as the alignment boundary, for third-generation sequencing sequences with high error rates; improve the identification accuracy of complex or long variant DNA fragments (SVs) (exceeding the length of sequencing data); accurately distinguish adjacent variant DNA fragments (SVs) that are close in distance; and accurately genotype variant DNA fragments (SVs).

[0006] The embodiments of the present invention are implemented as follows:

[0007] A method for identifying whole-genome structural variations based on third-generation sequencing, comprising:

[0008] S11: Align the sequencing data with the preset reference genome, sort the preset alignment files, and build an index.

[0009] S21 analyzes the base alignment (cigar string) of each sequencing read to identify SVs and determine whether they contain a signal of a variant DNA fragment (SV).

[0010] S31 corrects erroneous and overlapping variant DNA fragment signals.

[0011] S41 outputs the polished sequence, compares it with the preset reference genome, performs SV identification, and outputs the identification results.

[0012] S51, based on the identification results, the sequencing data is genotyped.

[0013] In a preferred embodiment of the present invention, in S21 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of parsing the base alignment information (cigar string) of each sequencing read to identify whether it contains a variant DNA fragment (SV) includes:

[0014] S211, If ​​the sequencing read contains an insertion or deletion fragment of 30 bp or more, then store the sequencing read and the interval position of the insertion or deletion fragment in a list;

[0015] S212, If the sequencing read contains a clip of 100 bp or more, store the interval positions of the sequencing read and the clip in a list.

[0016] S213, if its cigar string includes more than one different part aligned to different positions (Supplement), then store the sequencing read and multiple alignment positions in a list.

[0017] In a preferred embodiment of the present invention, in S31 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the correction of erroneous and overlapping variant DNA fragment signals includes:

[0018] S311 calculates the alignment depth of the sequencing data in 500bp windows. If the alignment depth of a certain window is much greater than the average depth, the variant DNA fragment (SV) signal caused by the secondary alignment of that window is filtered out to reduce false variant DNA fragment (SV) signals caused by alignment errors.

[0019] S312 clusters overlapping variant DNA fragment (SV) signals.

[0020] S313: For a given cluster, calculate the consistency between signals. If the consistency is high, take the median value and output it, then delete the cluster.

[0021] S314. Further merge clusters. Iterate over any two clusters. If the signals of the two clusters contain overlapping sequencing reads, and there are more than two overlapping sequencing reads, then merge the two clusters into one cluster.

[0022] S315 filters the clusters.

[0023] S316 outputs a variant DNA fragment (SV) of one type for each cluster.

[0024] In a preferred embodiment of the present invention, in S31 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the correction of erroneous and overlapping variant DNA fragment signals further includes:

[0025] S317, Repeat the filtering of clusters until there are no clusters.

[0026] In a preferred embodiment of the present invention, in S312 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of clustering overlapping variant DNA fragment (SV) signals includes:

[0027] S3121 sorts the list of variant DNA fragment (SV) signals in ascending order by chromosome, start coordinate, and end coordinate.

[0028] S3122: The current variant DNA fragment (SV) signal is compared cyclically with the previous variant DNA fragment (SV) signal. If the two signals are on the same chromosome and the distance is within 500 bp, the current signal is added to the cluster of the previous signal; otherwise, it is used as a new clustering starting point.

[0029] In a preferred embodiment of the present invention, in S313 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of calculating the consistency between signals for a cluster, and if the consistency is high, taking the median value and outputting the cluster, includes:

[0030] S3131 checks whether the types of variant DNA fragment (SV) signals are consistent; if they are inconsistent, the cluster is skipped.

[0031] S3132 sequentially detects the variance of the start coordinates, the variance of the end coordinates, and the variance of the length of the variant DNA fragment (SV) signal. If any variance is greater than 50 (parameter), the cluster is skipped.

[0032] S3133: Take the median value of the starting coordinate as the starting position of the mutant DNA fragment (SV), the median value of the ending coordinate as the ending position of the mutant DNA fragment (SV), and the median value of the length as the length of the mutant DNA fragment (SV), and output the mutant DNA fragment (SV).

[0033] S3134, Remove the variant DNA fragment (SV) cluster from the cluster list output above.

[0034] In a preferred embodiment of the present invention, in S315 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the filtering of clusters includes:

[0035] S3151, calculate the minimum start position and maximum end position of the cluster variant DNA fragment (SV) signal, and use them as the signal interval of the cluster.

[0036] S3152, based on the depth information obtained in S311, calculate the average depth of the cluster signal interval.

[0037] S3153, if the number of reads contained in the variant DNA fragment (SV) signal of the cluster is less than 30% of the average depth, or less than 3, it is considered that the variant DNA fragment (SV) signal is caused by an excessively high error rate or an alignment error, and the cluster is deleted.

[0038] In a preferred embodiment of the present invention, in S316 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of outputting a variant DNA fragment (SV) of one type for each cluster includes:

[0039] S3161, calculate the minimum start position and the maximum end position of the cluster variant DNA fragment (SV) signal, and use them as the signal interval of the cluster.

[0040] S3162 checks whether the signal range of the cluster can be crossed by any sequencing read in the cluster.

[0041] S3163: If there is coverage, the sequencing reads are used as the reference sequence for the cluster; if there is no coverage, the sequencing reads of the cluster are paired, the graph is assembled, and the longest contig sequence is output as the reference sequence for the cluster.

[0042] S3164 aligns all sequencing reads to the reference sequence of the cluster and filters out alignments containing variant DNA fragment (SV) signals. The reference sequence is then polished using the filtered consistency alignments.

[0043] S3165 outputs the polished sequence.

[0044] In a preferred embodiment of the present invention, in S316 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of outputting a variant DNA fragment (SV) of one type for each cluster further includes:

[0045] S3166, re-align all sequencing reads in the cluster to the above-output variant DNA fragments, filter out sequencing reads without variant DNA fragment (SV) signals, and delete the signals contained in the filtered sequencing reads from the cluster.

[0046] S3167 records sequencing reads without variant DNA segment (SV) signals, as sequencing reads that support the presence of variant DNA segment (SV) signals in the sequence after the polishing.

[0047] In a preferred embodiment of the present invention, in S51 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of typing the sequencing data according to the identification results includes:

[0048] S511 If two variant DNA fragments (SVs) come from the same cluster and the supported sequencing reads do not overlap (result of step S3167), then the two variant DNA fragments (SVs) (result of step S41) are set as heterozygous variant DNA fragments (SVs) of different types.

[0049] S512, if the situation in S511 does not occur, detect the ratio of the number of sequencing reads supporting the variant DNA fragment (SV) (the result of step S3167) to the total number of sequencing reads that cross the signal of the variant DNA fragment (SV). If it is greater than 80%, it is set as homozygous; otherwise, it is heterozygous.

[0050] The beneficial effects of the embodiments of the present invention are:

[0051] This invention utilizes sequence alignment to identify variant DNA fragments (SVs). It can quickly identify a portion of simple variant DNA fragments (SVs). For the remaining complex variant DNA fragments (SVs), a sequence assembly method is used to make the final results more accurate without being too slow. Using sequence assembly technology (step S3163), sequences longer than the sequencing reads can be output, which can be used to more accurately identify the type and length of the SV. For sequences containing variant DNA fragments (SVs) (step 3164), a polishing correction is performed to improve the accuracy of the sequences. By using alignment to filter sequencing reads belonging to different variant DNA fragments (SVs) in the same cluster, different types of variant DNA fragments (SVs) can be more accurately distinguished. Attached Figure Description

[0052] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0053] Figure 1 This is a flowchart of the whole-genome structural variation identification method based on third-generation sequencing of the present invention;

[0054] Figure 2 This is a schematic diagram of the complete process of the whole genome structural variation identification method based on third-generation sequencing of the present invention;

[0055] Figure 3 This is a schematic diagram illustrating the SV detection results using ONT sequencing data of different depths based on third-generation sequencing, as described in this invention.

[0056] Figure 4 This is a HiFi data illustration of the whole-genome structural variation identification method based on third-generation sequencing of the present invention;

[0057] Figure 5 This is the comparison result of the whole genome structural variation identification method based on third-generation sequencing of this invention. Detailed Implementation

[0058] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0059] Explanation of terms used in this invention:

[0060] Structural variation (SV) refers to variations in DNA fragments longer than 50 bp. Variations can be categorized by type, such as deletion, duplication, insertion, inversion, and translocation.

[0061] Next-generation sequencing (NGS): High-throughput sequencing, mainly using instruments on the Illumina platform, can sequence hundreds of thousands to millions of DNA molecules at once. The resulting fragments have shorter read lengths, usually only 100-150 bp, and higher accuracy, usually around 99.9%.

[0062] Third-generation sequencing: Single-molecule sequencing technology, mainly referring to sequencing technologies on PacBio or Oxford Nanopore (ONT) platforms, yields longer fragments. Among them, PacBio platform HiFi data has high accuracy, but the length is usually only 10-20kb. ONT data has a higher error rate, but the length is longer (>30kb).

[0063] Polish: Uses third-generation sequencing data to correct single-base or short insertion / deletion errors in assembled genomes.

[0064] Reads (DNA sequencing fragments): A read obtained from sequencing, i.e., a raw sequencing data.

[0065] Assembly: Whether it is short reads in second-generation sequencing or long reads in third-generation sequencing, the resulting sequencing data is still extremely small compared to the entire genome. Sequence assembly is to connect these small fragments through various algorithms and methods to obtain the initial contig.

[0066] Contig: Based on the overlap information obtained from sequencing reads, an existing overlapping fragment is combined to form a contig, which is a sequence that did not have a gap in the initial version.

[0067] Genome: The contig sequence obtained through assembly; currently, the genomes of most species are publicly available.

[0068] Alignment: Using sequence similarity, sequencing reads are located at specific locations in the genome.

[0069] Typing: The normal human genome contains two sets of chromosomes (one from the father and one from the mother). Therefore, the same location in the genome may contain different bases or sequences. Typing in this article determines whether a given SV location contains only one type of SV or two types of SVs.

[0070] Homozygous: The same position contains only one type of SV (same parents).

[0071] Hybrid: The same position contains 2 SVs (one from the parent and one from the mother).

[0072] Cigar: Alignment result information, describing the alignment status of all bases in a read (such as mismatches, insertions, etc.).

[0073] Supplement: A single read cannot be completely aligned to the genome; different parts are aligned to different locations (block alignment).

[0074] Clip: The two ends of a read cannot be aligned to the genome, so they are cut off during alignment. This is usually caused by sequencing errors or SVs.

[0075] Please refer to Figures 1 to 4 The first embodiment of the present invention provides a method for identifying whole-genome structural variations based on third-generation sequencing, comprising:

[0076] Methods for identifying whole-genome structural variations based on third-generation sequencing include:

[0077] S11: Align the sequencing data with the preset reference genome, sort the preset alignment files, and build an index.

[0078] The index was built using the open-source software minimap2 and samtools.

[0079] S21 analyzes the base alignment (cigar string) of each sequencing read to identify SVs and determine whether they contain a signal of a variant DNA fragment (SV).

[0080] S31 corrects erroneous and overlapping variant DNA fragment signals.

[0081] S41 outputs the polished sequence, compares it with the preset reference genome, performs SV identification, and outputs the identification results.

[0082] S51, based on the identification results, the sequencing data is genotyped.

[0083] In a preferred embodiment of the present invention, in S21 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of parsing the base alignment information (cigar string) of each sequencing read to identify whether it contains a variant DNA fragment (SV) includes:

[0084] S211, If ​​the sequencing read contains an insertion or deletion fragment of 30 bp or more, then store the sequencing read and the interval position of the insertion or deletion fragment in a list;

[0085] S212, If the sequencing read contains a clip of 100 bp or more, store the interval positions of the sequencing read and the clip in a list.

[0086] S213, if its cigar string includes more than one different part aligned to different positions (Supplement), then store the sequencing read and multiple alignment positions in a list.

[0087] In a preferred embodiment of the present invention, in S31 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the correction of erroneous and overlapping variant DNA fragment signals includes:

[0088] S311 calculates the alignment depth of the sequencing data in 500bp windows. If the alignment depth of a certain window is much greater than the average depth, the variant DNA fragment (SV) signal caused by the secondary alignment of that window is filtered out to reduce false variant DNA fragment (SV) signals caused by alignment errors.

[0089] S312 clusters overlapping variant DNA fragment (SV) signals.

[0090] S313: For a given cluster, calculate the consistency between signals. If the consistency is high, take the median value and output it, then delete the cluster.

[0091] S314. Further merge clusters. Iterate over any two clusters. If the signals of the two clusters contain overlapping sequencing reads, and there are more than two overlapping sequencing reads, then merge the two clusters into one cluster.

[0092] S315 filters the clusters.

[0093] S316 outputs a variant DNA fragment (SV) of one type for each cluster.

[0094] In a preferred embodiment of the present invention, in S31 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the correction of erroneous and overlapping variant DNA fragment signals further includes:

[0095] S317, Repeat the filtering of clusters until there are no clusters.

[0096] In a preferred embodiment of the present invention, in S312 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of clustering overlapping variant DNA fragment (SV) signals includes:

[0097] S3121 sorts the list of variant DNA fragment (SV) signals in ascending order by chromosome, start coordinate, and end coordinate.

[0098] S3122: The current variant DNA fragment (SV) signal is compared cyclically with the previous variant DNA fragment (SV) signal. If the two signals are on the same chromosome and the distance is within 500 bp, the current signal is added to the cluster of the previous signal; otherwise, it is used as a new clustering starting point.

[0099] In a preferred embodiment of the present invention, in S313 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of calculating the consistency between signals for a cluster, and if the consistency is high, taking the median value and outputting the cluster, includes:

[0100] S3131 checks whether the types of variant DNA fragment (SV) signals are consistent; if they are inconsistent, the cluster is skipped.

[0101] S3132 sequentially detects the variance of the start coordinates, the variance of the end coordinates, and the variance of the length of the variant DNA fragment (SV) signal. If any variance is greater than 50 (parameter), the cluster is skipped.

[0102] S3133: Take the median value of the starting coordinate as the starting position of the mutant DNA fragment (SV), the median value of the ending coordinate as the ending position of the mutant DNA fragment (SV), and the median value of the length as the length of the mutant DNA fragment (SV), and output the mutant DNA fragment (SV).

[0103] S3134, Remove the variant DNA fragment (SV) cluster from the cluster list output above.

[0104] In a preferred embodiment of the present invention, in S315 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the filtering of clusters includes:

[0105] S3151, calculate the minimum start position and maximum end position of the cluster variant DNA fragment (SV) signal, and use them as the signal interval of the cluster.

[0106] S3152, based on the depth information obtained in S311, calculate the average depth of the cluster signal interval.

[0107] S3153, if the number of reads contained in the variant DNA fragment (SV) signal of the cluster is less than 30% of the average depth, or less than 3, it is considered that the variant DNA fragment (SV) signal is caused by an excessively high error rate or an alignment error, and the cluster is deleted.

[0108] In a preferred embodiment of the present invention, in S316 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of outputting a variant DNA fragment (SV) of one type for each cluster includes:

[0109] S3161, calculate the minimum start position and the maximum end position of the cluster variant DNA fragment (SV) signal, and use them as the signal interval of the cluster.

[0110] S3162 checks whether the signal range of the cluster can be crossed by any sequencing read in the cluster.

[0111] S3163: If there is coverage, the sequencing reads are used as the reference sequence for the cluster; if there is no coverage, the sequencing reads of the cluster are paired, the graph is assembled, and the longest contig sequence is output as the reference sequence for the cluster.

[0112] S3164 aligns all sequencing reads to the reference sequence of the cluster and filters out alignments containing variant DNA fragment (SV) signals. The reference sequence is then polished using the filtered consistency alignments.

[0113] S3165 outputs the polished sequence.

[0114] In a preferred embodiment of the present invention, in S316 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of outputting a variant DNA fragment (SV) of one type for each cluster further includes:

[0115] S3166, re-align all sequencing reads in the cluster to the above-output variant DNA fragments, filter out sequencing reads without variant DNA fragment (SV) signals, and delete the signals contained in the filtered sequencing reads from the cluster.

[0116] S3167 records sequencing reads without variant DNA segment (SV) signals, as sequencing reads that support the presence of variant DNA segment (SV) signals in the sequence after the polishing.

[0117] In a preferred embodiment of the present invention, in S51 of the above-described method for identifying whole-genome structural variations based on third-generation sequencing, the step of typing the sequencing data according to the identification results includes:

[0118] S511 If two variant DNA fragments (SVs) come from the same cluster and the supported sequencing reads do not overlap (result of step S3167), then the two variant DNA fragments (SVs) (result of step S41) are set as heterozygous variant DNA fragments (SVs) of different types.

[0119] S512, if the situation in S511 does not occur, detect the ratio of the number of sequencing reads supporting the variant DNA fragment (SV) (the result of step S3167) to the total number of sequencing reads that cross the signal of the variant DNA fragment (SV). If it is greater than 80%, it is set as homozygous; otherwise, it is heterozygous.

[0120] This invention utilizes ONT sequencing data at different depths to perform SV detection results as follows: Figure 3 As shown in the figure, it can be seen that GrandSV's results at different sequencing depths are better than those of similar software.

[0121] High accuracy and sensitivity:

[0122] (1) For sequences containing variant DNA fragments (SV) (step S3164), a polishing correction was performed to improve the accuracy of the sequences (from 90-97% to about 99%). Then, by re-aligning the more accurate sequences, the identified variant DNA fragments (SV) and their boundaries were made more accurate.

[0123] (2) For complex regions, assembly technology (step S3163) can be used to output sequences longer than the sequencing reads. Then, by using longer and more accurate sequence re-alignment, more accurate alignment results can be obtained, which will ultimately lead to more accurate identification results of variant DNA fragments (SVs).

[0124] (3) When identifying different variant DNA fragments (SVs) by typing, the traditional approach is to simply compare the signal coordinate ranges of the variant DNA fragments (SVs). If the coordinate ranges of variant DNA fragments (SVs) of different types are not significantly different, they cannot be distinguished. However, the Grand variant DNA fragment (SV) method uses alignment filtering to identify sequencing reads belonging to different variant DNA fragments (SVs) within the same cluster (step S3166), which more accurately distinguishes different types of variant DNA fragments (SVs).

[0125] (4) By iterative looping, low-frequency variant DNA fragments (SVs) with low support numbers can be identified.

[0126] (5) Figure 4As shown in the HiFi data, due to the limitations of sequencing read length, long variant DNA fragments (SVs) are not covered by HiFi data. Therefore, the traditional approach only reports the breakpoints of the variant DNA fragments (SVs), and cannot report the type and length of the variant DNA fragments (SVs). GrandSV can utilize assembly technology (step S3163) to output the sequence that spans the variant DNA fragment (SV). Using this sequence, the type and length of the variant DNA fragment (SV) can be identified more accurately.

[0127] Please refer to Figures 1 to 5 The second embodiment of the present invention provides a method for identifying whole-genome structural variations based on third-generation sequencing.

[0128] like Figure 5 As shown, the identification steps are as follows, taking ONT third-generation sequencing as an example:

[0129] Figure 5 In the diagram, the numbers on the left represent the Read ID, the blue boxes with numbers indicate the signal for inserting that number length, and the horizontal lines with numbers indicate the signal for deleting that number length.

[0130] (1) Sequence alignment: Aligning sequencing reads to the genome. The alignment results are as follows: Figure 5 As shown.

[0131] (2) Identify the presence of SV signals for each sequencing read. For example... Figure 5 As shown, read 1 contains a 232bp delete signal, read 2 contains a 76bp insert signal, read 3 contains a 302bp delete signal, read 4 does not contain an SV signal, read 5 contains a 305bp delete signal, reads 6-9 do not contain SV signals, read 10 contains a 74bp insert signal, read 11 contains a 70bp insert signal, read 12 contains 223bp and 512bp insert signals, read 13 contains a 66bp insert signal, read 14 does not contain an SV signal, read 15 contains a 75bp insert signal, read 16 contains a 305bp delete signal, read 17 does not contain an SV signal, read 18 contains a 305bp delete signal, read 19 does not contain an SV signal, read 20 contains a 78bp insert signal, read 21 contains 301bp, 69bp, and 367bp insert signals, and reads 22-24 do not contain SV signals. Reads without SV signals will be ignored in subsequent analysis.

[0132] (3) The alignment depth of the genome was calculated according to a window of 500bp. It was found that the alignment depth of this window was not significantly different from the average depth, so no filtering was performed.

[0133] (4) Cluster the SV signals of the same read. For the insert signals contained in Read 12 and 21, since these signals are close together, they can be merged into one signal. The merged Read 12 contains 735bp of insert signals, and Read 21 contains 737bp of insert signals.

[0134] (5) Sort the identified SV signals in ascending order according to chromosome, start coordinate, and end coordinate. Since the SV signals contained in read1,2,3,5,10,11,12,13,15,16,18,20,21 are relatively close, they will be clustered into a large class, named A.

[0135] (6) Calculate the consistency between signals in category A. It is found that the signal types in category A are inconsistent, that is, there are both inserted and deleted signals. Therefore, this category is not output or deleted.

[0136] (7) Filter clusters according to depth information. Since cluster A contains a large number of reads, far exceeding 30% of the average depth, this cluster is not filtered.

[0137] (8) Randomly select a read from A, assuming the selected read is read2. Since the signal range of cluster A can be crossed by read2, and the alignment length of read2 is 91kp, which is much larger than the span of A (about 1.8kb), no assembly is required. Directly use read2 as the reference sequence B1 for this cluster.

[0138] (9) Align the remaining reads 1, 3, 5, 10, 11, 12, 13, 15, 16, 18, 20, 21 to the reference sequence B1. Since reads 1, 3, 5, 12, 16, 18, 21 contain inconsistent SV signals (inconsistent length or category) with the reference sequence B1, their alignment results will show SV signals and are therefore filtered out. Reads 10, 11, 13, 15, 20 show perfect matches, and these perfect match alignment results are used to correct the reference sequence B1 to generate sequence B2.

[0139] (10) Align the remaining reads 1, 3, 5, 10, 11, 12, 13, 15, 16, 18, 20, and 21 back to the reference sequence B2. Since the alignment of reads 10, 11, 13, 15, and 20 does not contain the SV signal, these reads are considered as reads that support the SV contained in the sequence B2. And remove these reads from cluster A.

[0140] (11) Repeat steps 7-10 to obtain the polished sequences B2 (supporting reads 2, 10, 11, 13, 15, 20) and C2 (supporting reads 3, 5, 16, 18). Additionally, read1 forms one cluster, and reads12 and 21 form another cluster. These two clusters contain too few reads and will be filtered out in step 7.

[0141] (12) The polished sequences B2 and C2 are re-aligned to the genome using minimap2 and SVs are identified (the identification process is as in step 2). A 75bp insertion SV1 and a 304bp deletion SV2 are obtained and output.

[0142] (13) Classification. Since SV1 and SV2 come from the same cluster and their supported reads do not overlap, SV1 and SV2 are set as heterogeneous SVs.

[0143] Furthermore, the identification process for pacbio platform samples is the same as that for ONT platform samples, except for the parameter settings used for comparison and polishing in steps 1, 9, 10, and 12.

[0144] The embodiments of this invention aim to protect a method for identifying whole-genome structural variations based on third-generation sequencing, which has the following effects:

[0145] 1. Taking advantage of the speed of sequence alignment method to identify variant DNA fragments (SVs), a portion of simple variant DNA fragments (SVs) can be identified. For the remaining complex variant DNA fragments (SVs), sequence assembly method is used to make the final result more accurate, while not being too slow.

[0146] 2. Using sequence assembly technology (step S3163), sequences longer than the sequencing reads can be output.

[0147] 3. For sequences containing variant DNA fragments (SV) (step S3164), a polishing correction was performed to improve sequence accuracy.

[0148] 4. By comparing and filtering reads belonging to different variant DNA fragments (SVs) in the same cluster (step S3166), different types of variant DNA fragments (SVs) can be distinguished more accurately.

[0149] Finally, it should be noted that the above-described embodiments are merely specific implementations of the present invention, used to illustrate the technical solutions of the present invention, and not to limit it. The scope of protection of the present invention is not limited thereto. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that any person skilled in the art can still modify or easily conceive of changes to the technical solutions described in the foregoing embodiments within the technical scope disclosed in the present invention, or make equivalent substitutions for some of the technical features; and these modifications, changes, or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for identifying whole-genome structural variations based on third-generation sequencing, characterized in that, include: S11, compare the sequencing data with the preset reference genome, sort the preset alignment files, and build an index; S21: For each sequencing read, analyze its base alignment and perform SV identification to determine whether it contains a signal of a mutated DNA fragment. S31, corrects erroneous and overlapping variant DNA fragment signals; In S31, the correction of erroneous and overlapping variant DNA fragment signals includes: S311 calculates the alignment depth of the sequencing data in 500bp windows. If the alignment depth of a certain window is much greater than the average depth, the variant DNA fragment signal caused by the secondary alignment of that window is filtered out to reduce false variant DNA fragment signals caused by alignment errors. S312 clusters overlapping variant DNA fragment signals; S313: For a given cluster, calculate the consistency between signals. If the consistency is high, output the median value and delete the cluster. S314, further merge clusters, iterate over any two clusters, if the sequencing reads contained in the signals of the two clusters overlap and there are more than two overlapping sequencing reads, then merge the two clusters into one cluster. S315, filtering of clusters; S316, for each cluster, outputs a variant DNA fragment of one type; In S31, the correction of erroneous and overlapping variant DNA fragment signals further includes: S317, Repeat the filtering of clusters until there are no clusters; In S312, the clustering of overlapping variant DNA fragment signals includes: S3121, sort the list of mutated DNA fragment signals in ascending order by chromosome, start coordinate, and end coordinate; S3122: The current mutated DNA fragment signal is compared with the previous mutated DNA fragment signal in a loop. If the two signals are on the same chromosome and the distance is within 500 bp, the current signal is added to the cluster of the previous signal; otherwise, it is used as a new clustering starting point. In S313, the step of calculating the consistency between signals for a cluster, and if the consistency is high, taking the median value as the output and deleting the cluster, includes: S3131, detect whether the types of mutated DNA fragment signals are consistent; if they are inconsistent, skip the cluster. S3132, sequentially detect the variance of the start coordinates, the variance of the end coordinates, and the variance of the length of the mutated DNA fragment signal. If any variance is greater than 50, skip the cluster. S3133: Take the median value of the starting coordinate as the starting position of the mutant DNA fragment, the median value of the ending coordinate as the ending position of the mutant DNA fragment, and the median value of the length as the length of the mutant DNA fragment, and output the mutant DNA fragment. S3134, Remove the above-output cluster of variant DNA fragments from the cluster list; In S315, the filtering of clusters includes: S3151, calculate the minimum start position and the maximum end position of the DNA fragment signal for each cluster variation, and use them as the signal interval for that cluster; S3152, Based on the depth information obtained in S311, calculate the average depth of the clustered signal interval; S3153, if the number of reads contained in the variant DNA fragment signal of the cluster is less than 30% of the average depth, or less than 3, it is considered to be a variant DNA fragment signal caused by an excessively high error rate or an alignment error, and the cluster is deleted. In S316, the step of outputting a variant DNA fragment of one type for each cluster includes: S3161, calculate the minimum start position and the maximum end position of the DNA fragment signal for each cluster variation, and use them as the signal interval for that cluster; S3162, detect whether the signal interval of the cluster can be crossed by any sequencing read in the cluster; S3163, if there is coverage, the sequencing read is used as the reference sequence for the cluster; if there is no coverage, the sequencing reads of the cluster are paired, the graph is assembled, and the longest contig sequence is output as the reference sequence for the cluster. S3164 aligns all sequencing reads to the reference sequence of the cluster and filters out alignment results containing variant DNA fragment signals, then uses the filtered consistency alignment to correct the reference sequence. S3165, Output the corrected sequence; In S316, the step of outputting a variant DNA fragment of one type for each cluster also includes: S3166, Realign all sequencing reads in the cluster to the above-output variant DNA fragments, filter out sequencing reads without variant DNA fragment signals, and remove the signals contained in the filtered sequencing reads from the cluster; S3167, records sequencing reads without a signal of a variant DNA fragment, which are considered as sequencing reads that support the variant DNA fragment contained in the sequence after the correction; S41, output the corrected sequence, compare it with the preset reference genome, perform SV identification, and output the identification results; S51, based on the identification results, the sequencing data is genotyped; In S51, the step of genotyping the sequencing data based on the identification results includes: S511 If two variant DNA fragments come from the same cluster and the supported sequencing reads do not overlap, then the two variant DNA fragments are set as heterozygous variant DNA fragments of different types. S512, if the situation in S511 does not occur, detect the ratio of the number of sequencing reads supporting the variant DNA fragment to the total number of sequencing reads that cross the signal of the variant DNA fragment. If it is greater than 80%, it is set as homozygous; otherwise, it is heterozygous.

2. The method for identifying whole-genome structural variations based on third-generation sequencing according to claim 1, characterized in that, In S21, the step of analyzing the base alignment of each sequencing read and identifying whether it contains a mutated DNA fragment includes: S211, If ​​the sequencing read contains an insertion or deletion fragment of 30 bp or more, store the sequencing read and the interval position of the insertion or deletion fragment in a list; S212, If the sequencing read contains an alignment excision fragment of 100 bp or more, store the sequencing read and the interval position of the excision fragment in a list. S213, if the base alignment includes more than one different part aligned to different positions, then the sequencing read and multiple alignment positions are stored in a list.

Images