Error Assembly Detection and Repair Method and System Based on Paired-End Sequencing Data

Abstract

The present invention belongs to the technical field of sequence analysis, and particularly relates to a method and system for detecting and repairing incorrect assemblies based on paired-end sequencing data. By using the alignment information of the paired-end sequencing data with the assembly results, various abnormal signals are extracted and subjected to clustering analysis, which can effectively identify true incorrect assemblies. According to the types of incorrect assemblies, combined with the sequencing data information, the assembly results are corrected, thereby optimizing the assembly results and improving the overall accuracy and integrity of genome assembly.

Description

Technical Field

[0001] The present invention belongs to the technical field of sequence analysis, and particularly relates to a method and system for detecting and repairing incorrect assemblies based on paired-end sequencing data. Background Art

[0002] With the rapid development of genome sequencing technology, various assembly algorithms have emerged one after another, greatly improving the ability to construct genomic sequences from sequencing data. Complete and accurate genome assembly is crucial for analyzing the genetic information of organisms, studying evolutionary mechanisms, and identifying gene variations related to diseases (including single nucleotide polymorphisms, insertions / deletions, and structural variations). This information has important application values in fields such as disease mechanism research, precision medicine, and personalized treatment.

[0003] However, achieving high-quality assembly is challenging. Existing genome assembly algorithms still have certain limitations, which may lead to errors during the assembly process. Incorrect assemblies usually manifest as structural errors, such as insertion / deletion (indel) errors and misjoin errors. These errors not only affect the accuracy of the genome but also have a negative impact on comparative genomics, genetic variation analysis, and downstream bioinformatics research. To solve the problem of incorrect assemblies, various error detection methods have been proposed. These methods are mainly divided into two categories: reference genome-based methods and read alignment-based methods.

[0004] Reference genome-based methods detect errors by aligning the assembly results with the reference genome. For example, QUAST can identify errors such as translocations and inversions, but it is difficult to distinguish between assembly errors and true variations generated by individuals. In addition, due to the lack of accurate reference genomes for many organisms in the real environment, the applicability of this type of method is limited. Read alignment-based methods, on the other hand, use the sequencing data itself to evaluate the assembly quality, and identify repeat regions, incorrect assemblies, and insertion / deletion errors by analyzing features such as changes in coverage depth and abnormal insert fragment lengths. However, existing read alignment-based methods have a high false positive rate, and most tools cannot mark the specific positions where errors occur, let alone have the function of correcting assembly errors, which affects their application in actual scenarios. Summary of the Invention

[0005] Current sequencing technologies and genome assembly tools still face many challenges in the process of complete chromosome reconstruction. These challenges mainly include sequencing errors, limitations in sequencing coverage, and a large number of repetitive sequences present in the genome. Although second-generation sequencing (NGS) technologies provide short-read data with high coverage and high accuracy, errors still commonly exist in genome assembly due to factors such as polyploid structure, genetic diversity, heterozygosity, and repetitive sequences. Therefore, the present invention provides a method and system for detecting and repairing incorrect assemblies based on paired-end sequencing data.

[0006] The technical solution of the present invention is as follows:

[0007] The present invention provides a method for detecting and repairing incorrect assemblies based on paired-end sequencing data, comprising the following steps:

[0008] S1: Based on the assembly result, with a preset suspicious site as the core, extend a first preset length to both sides of the genomic sequence, and use it as a candidate error region;

[0009] Compare the paired-end sequencing data with the assembly result, and based on the comparison information, extract various characteristic signals of the candidate error region. The various characteristic signals are a coverage anomaly signal, an insertion or deletion anomaly signal, a paired direction anomaly signal, an insertion size anomaly signal, a read pair anomaly signal, and a splicing anomaly signal;

[0010] S2: For each candidate error region, based on the various characteristic signals of the extracted candidate error region, perform clustering to obtain the clustering result of each characteristic signal of the corresponding candidate error region; if the clustering result of one characteristic signal is abnormal, then this candidate error region is an incorrect assembly region;

[0011] S3: Based on the incorrect assembly region, determine the type of assembly error. If the type is an insertion or deletion, remove the inserted sequence or the inserted or deleted fragment;

[0012] If the type is an incorrect connection, generate a consensus sequence through multiple sequence alignment, compare the consensus sequence with the assembly result, and perform repair.

[0013] In the said S1, the specific method for extracting the coverage anomaly signal of the candidate error region is as follows:

[0014] With the candidate error region as the benchmark, extend a second preset length to both sides of it, and use it as an extended region;

[0015] Calculate the average coverage of all sites in the extended region, traverse and detect the coverage of each site in the candidate error region. If the coverage of a certain site is not less than the first preset multiple of the average coverage or not greater than the second preset multiple of the average coverage, then this site is recorded as a coverage anomaly site;

[0016] Statistically count the total number of sites with abnormal coverage within the candidate error region, calculate the ratio of the total number to the total length of the candidate error region, and use this ratio as the coverage anomaly signal.

[0017] For the above S1, extract the insertion or deletion anomaly signal of the candidate error region, specifically:

[0018] Statistically count the number of reads with insertions or deletions at each site in the candidate error region. If the coverage of a certain site is not less than the preset value, calculate the ratio of the number of reads to the coverage of this site, and use this ratio as the insertion or deletion signal of this site;

[0019] Compare the insertion or deletion signals of all sites in the candidate error region, and use the maximum ratio as the insertion or deletion anomaly signal of this candidate error region.

[0020] For the above S1, extract the paired orientation anomaly signal of the candidate error region, specifically:

[0021] Detect the reads with mismatched orientations at each site in the candidate error region. If reads with mismatched orientations are found at a certain site, mark the local region of this site as the mispaired orientation region, statistically count the number of reads with mismatched orientations in the mispaired orientation region and the total number of reads in the mispaired orientation region, and use the ratio of the number of reads to the total number as the paired orientation anomaly signal of this mispaired orientation region;

[0022] Compare the paired orientation anomaly signal magnitudes of all mispaired orientation regions within the candidate error region, and use the paired orientation anomaly signal of the mispaired orientation region with the maximum value as the paired orientation anomaly signal of this candidate error region.

[0023] For the above S1, extract the insertion size anomaly signal of the candidate error region, specifically:

[0024] Detect the reads with abnormal insertion sizes at each site in the candidate error region. If reads with abnormal insertion sizes are found at a certain site, mark the local region of this site as the insertion size abnormal region, statistically count the number of reads with abnormal insertion sizes in the insertion size abnormal region and the total number of reads in the insertion size abnormal region, and use the ratio of the number of reads to the total number as the insertion size anomaly signal of this insertion size abnormal region;

[0025] Compare the insertion size anomaly signal magnitudes of all insertion size abnormal regions within the candidate error region, and use the insertion size anomaly signal of the insertion size abnormal region with the maximum value as the insertion size anomaly signal of this candidate error region.

[0026] For the above S1, extract the read pair anomaly signal of the candidate error region, specifically:

[0027] Isolate reads in the candidate error region are detected site by site. If a paired read that cannot be aligned to the same assembly result is found at a certain site, the local region of this site is marked as an isolated read region. The number of abnormally paired reads in the isolated read region and the total number of reads in the isolated read region are counted, and the ratio of the number of reads to the total number is used as the read pairing abnormal signal in the isolated read region.

[0028] Compare the magnitudes of the read pairing abnormal signals in all isolated read regions within the candidate error region, and use the read pairing abnormal signal with the maximum value in the isolated read region as the read pairing abnormal signal for this candidate error region.

[0029] In step S1, extract the splicing abnormal signal of the candidate error region, specifically:

[0030] Count the number of reads with splicing in the candidate error region site by site. If the coverage of a certain site is not less than the preset value, calculate the ratio of the number of reads to the coverage of this site, and use this ratio as the splicing abnormal signal for this site.

[0031] Compare the splicing abnormal signals of all sites in the candidate error region, and use the maximum ratio as the splicing abnormal signal for this candidate error region.

[0032] In step S2, before clustering, it further includes:

[0033] Randomly select several regions from the assembly results longer than the preset sequence length, and jointly perform abnormal signal clustering analysis with the candidate error region.

[0034] In step S3, if the type is misconnection, generate a consensus sequence through multiple sequence alignment, and compare the consensus sequence with the assembly result for repair, specifically:

[0035] Based on the assembly result of the misconnection type, analyze the alignment characteristic signals and mark the breakpoint region;

[0036] Based on the paired-end sequencing data, extract the paired-end reads that completely cover the breakpoint region. After filtering, obtain the local read sequence;

[0037] The local read sequence generates a consensus sequence through multiple sequence alignment. Compare the consensus sequence with the assembly result. If the successfully aligned length of the repaired sequence accounts for more than 90% of the breakpoint region, replace the misconnected region in the assembly result with the consensus sequence.

[0038] The present invention also provides an error assembly detection and repair system based on paired-end sequencing data, including:

[0039] Feature signal extraction module: Based on the assembly result, with a preset suspicious site as the core, extend a first preset length to both sides of the genomic sequence to obtain a candidate error region;

[0040] The paired-end sequencing data is aligned with the assembly result. Based on the alignment information, various characteristic signals of the candidate error regions are extracted. The various characteristic signals are coverage anomaly signals, insertion or deletion anomaly signals, paired orientation anomaly signals, insertion size anomaly signals, read pair anomaly signals, and splicing anomaly signals.

[0041] Error assembly region determination module: For each candidate error region, after clustering based on the various characteristic signals of the candidate error region extracted, the clustering results of each characteristic signal of the corresponding candidate error region are obtained; if the clustering result of one characteristic signal is abnormal, then the candidate error region is an error assembly region.

[0042] Repair module: Based on the error assembly region, the type of assembly error is determined. If the type is insertion or deletion, the inserted sequence or the inserted or deleted fragment is removed.

[0043] If the type is incorrect connection, a consensus sequence is generated through multiple sequence alignment, and the consensus sequence is aligned with the assembly result for repair.

[0044] Beneficial effects

[0045] The present invention utilizes the alignment information of the paired-end sequencing data with the assembly result, extracts various abnormal signals, and performs clustering analysis, which can effectively identify true error assemblies; according to the type of error assembly, combined with the sequencing data information, the assembly result is corrected, thereby optimizing the assembly result and improving the overall accuracy and integrity of genome assembly. Specific embodiments

[0046] The following embodiments are intended to illustrate the present invention rather than further limit the present invention.

[0047] The present invention provides a method for detecting and repairing error assemblies based on paired-end sequencing data, including the following steps:

[0048] S1: Based on the assembly result, with a preset suspicious site as the core, after extending a first preset length to both sides of the genomic sequence, it is used as a candidate error region.

[0049] The paired-end sequencing data is aligned with the assembly result. Based on the alignment information, various characteristic signals of the candidate error region are extracted. The various characteristic signals are coverage anomaly signals, insertion or deletion anomaly signals, paired orientation anomaly signals, insertion size anomaly signals, read pair anomaly signals, and splicing anomaly signals.

[0050] S2: For each candidate error region, after clustering based on multiple characteristic signals of the extracted candidate error region, the clustering results of each characteristic signal of the corresponding candidate error region are obtained; if the clustering result of one characteristic signal is abnormal, then this candidate error region is an incorrect assembly region.

[0051] S3: Based on the incorrect assembly region, determine the type of assembly error. If the type is insertion or deletion, then remove the inserted sequence or the inserted and deleted fragment.

[0052] If the type is incorrect connection, then through multiple sequence alignment, generate a consensus sequence, compare the consensus sequence with the assembly result, and perform repair.

[0053] The present invention utilizes the alignment information of paired-end sequencing data with the assembly result (contig), extracts multiple abnormal signals, and performs clustering analysis, which can effectively identify true incorrect assemblies; according to the type of incorrect assembly, combined with the sequencing data information, correct the assembly result, thereby optimizing the assembly result and improving the overall accuracy and integrity of genome assembly.

[0054] The following describes each step of the operation in detail.

[0055] S1: Based on the assembly result, with a preset suspicious site as the core, extend a first preset length to both sides of the genomic sequence, and use it as the candidate error region (for example, the total length of the candidate error region is 200bp).

[0056] Align the paired-end sequencing data with the assembly result. Based on the alignment information, extract multiple characteristic signals of the candidate error region. The multiple characteristic signals are coverage abnormal signals, insertion or deletion abnormal signals, paired direction abnormal signals, inserted size abnormal signals, read pair abnormal signals, and shear abnormal signals.

[0057] Considering that: (1) Correctly assembled sequences usually exhibit uniform coverage. Deviations from the expected coverage (such as significantly higher or lower coverage) may respectively indicate folding or expanded repeats. (2) The presence of a large number of insertions and deletions (indels) in certain regions may indicate assembly errors because correctly assembled contigs should show consistency with reads during alignment. (3) Paired reads derived from the same double-stranded DNA molecule should be aligned within a reasonable distance and have the correct direction on the same contig. Incorrect assembly may lead to abnormal read alignments, such as paired reads being aligned to different contigs, unreasonable inserted sizes of paired reads aligned to the same contig, or incorrect paired directions of reads. Therefore, the present invention first extracts multiple characteristic signals of the candidate error region.

[0058] 1. Regarding coverage abnormality

[0059] Generally, reads can be uniquely aligned to contigs, while reads of repetitive sequences may have multiple alignment positions. Coverage is the number of reads aligned at a position. Correct assembly is usually accompanied by uniform coverage with a small fluctuation range. If there are folded or extended repeats in the candidate error region, abnormal coverage higher or lower than the average coverage of adjacent regions may occur.

[0060] Therefore, extract the abnormal coverage signal of the candidate error region, specifically:

[0061] Taking the candidate error region ( ), extend it by a second preset length (such as 400bp) on both sides to form an extended region ( , and the total length of the extended region is 1000bp);

[0062] Calculate the average coverage of all sites in the extended region ( ), traverse and detect the coverage of each site in the candidate error region. If the coverage of a certain site i ( ) is not less than the first preset multiple (such as 2 times) of the average coverage, or not greater than the second preset multiple (such as 0.5 times) of the average coverage, then this site is marked as an abnormal coverage site (that is, the site type is marked as 1);

[0063] Count the total number of abnormal coverage sites in the candidate error region, calculate the ratio of the total number to the total length of the candidate error region, and take this ratio as the abnormal coverage signal of this candidate error region ( ).

[0064] It can be implemented through the following formula:

[0065] ;

[0066] ;

[0067] ;

[0068] ;

[0069] ;

[0070] In the formula, , , , They are respectively the start position of the candidate error region, the end position of the candidate error region, the start position of the extended region, and the end position of the extended region; i refers to a specific base site in the region to be detected. According to the BAM file, the alignment information in the region to be detected is preprocessed and saved in the array Base. Then, the corresponding information is extracted from the Base array base by base in the region to be detected and the corresponding abnormal signals are calculated.

[0071] 2. Regarding insertions or deletions (Indels) anomalies

[0072] There may be a large number of insertions or deletions in the candidate error region. Insertion-deletion events mean that the aligned reads show insertions or deletions relative to the contig. For example, if the contig has the sequence CGCA, but the aligned reads at the corresponding position indicate an insertion of the TG sequence (CTGGCA), then the contig may contain a deletion error.

[0073] Therefore, the insertion or deletion abnormal signals in the candidate error region are extracted, specifically as follows:

[0074] Count the number of reads with insertions or deletions at each site in the candidate error region. If the coverage of a certain site is not less than the preset value, then calculate the ratio of the number of reads ( ) to the coverage of this site ( ), and take this ratio as the insertion or deletion signal at this site ( );

[0075] Compare the insertion or deletion signals of all sites in the candidate error region, and take the maximum ratio as the insertion or deletion abnormal signal of this candidate error region ( ).

[0076] It can be achieved through the following formula:

[0077] ;

[0078] ;

[0079] Among them, for a certain low-coverage site, when its coverage is less than the preset value, such as less than 0.2 times the average coverage, it indicates that the number of effective reads at this position is insufficient, which may lead to insufficient analysis results and further increase the false positive risk. To avoid misjudgment, such sites are not considered in subsequent analyses.

[0080] In actual assembly situations, interference from heterozygous genes in diploid organisms may also generate certain indel signals, leading to false positives. At a certain locus on a contig, when a large number of insertion or deletion events occur in the corresponding reads, these reads can be identified as "error-supporting" reads (Indel-supporting Reads), that is, reads that show significant indel differences from the contig in the alignment. Compared with genetic variations, the expected proportion of "error-supporting" reads in assembly errors is higher ( ratio indel > 0.5). By setting a threshold, filtering operations are performed on the insertion-deletion signals ( ), filtering out regions with low indel intensity, which do not participate in subsequent clustering analysis, reducing the occurrence of false positives.

[0081] 3. Abnormal pairing direction

[0082] When reads are aligned with a contig, the reads should be paired in the correct direction. Mismatches between different strands will result in a large number of reads with incorrect directions in the alignment, and these reads usually cluster in a local area of the candidate error region.

[0083] Therefore, extract the abnormal pairing direction signal of the candidate error region, specifically:

[0084] Detect the reads with mismatched directions at each locus in the candidate error region one by one. If reads with mismatched directions are found at a certain locus, mark the local area of this locus as the error pairing direction region, count the number of reads with mismatched directions in the error pairing direction region and the total number of reads in the error pairing direction region, and use the ratio of the number of reads ( ), to the total number ( ) as the abnormal pairing direction signal of this error pairing direction region ( );

[0085] Compare the sizes of the abnormal pairing direction signals of all error pairing direction regions within the candidate error region. Preferably, select the "error pairing direction region" where the proportion of all reads with abnormal pairing directions is greater than 20%, and use the abnormal pairing direction signal of the error pairing direction region with the maximum value as the abnormal pairing direction signal of this candidate error region ( ).

[0086] It can be achieved through the following formula:

[0087] ;

[0088] .

[0089] 4. Regarding abnormal insert size

[0090] The insert size refers to the distance between the two paired ends of paired reads. Under normal circumstances, this distance should fluctuate within a certain range. If the insert size is too short or too long, the insert size is unreasonable. The distance of the insert size of reads can be considered to follow a normal distribution. If the insert size of reads is abnormal and deviates from the limit of normal reads (i.e., less than the minimum value of the insert size of normal reads or greater than the maximum value of the insert size of normal reads ), assembly errors may occur.

[0091] To this end, extract the insert size abnormal signal of the candidate error region, specifically:

[0092] Detect the reads with abnormal insert size at each site in the candidate error region one by one. If reads with abnormal insert size are found at a certain site, mark the local area of this site as the insert size abnormal region, count the number of reads with abnormal insert size in the insert size abnormal region and the total number of reads in the insert size abnormal region, and take the ratio of the number of reads ( ), to the total number ( ) as the insert size abnormal signal of this insert size abnormal region ( );

[0093] Compare the insert size abnormal signals of all insert size abnormal regions within the candidate error region. It is preferably the "insert size abnormal region" where the proportion of all insert size abnormal reads is greater than 20%. Take the insert size abnormal signal of the insert size abnormal region with the maximum value as the insert size abnormal signal of this candidate error region ( ).

[0094] It can be achieved through the following formula:

[0095] ;

[0096] ;

[0097] ;

[0098] ;

[0099] In the formula, , are the mean and standard deviation of the sequencing data respectively.

[0100] 5. Regarding abnormal read pairing

[0101] Paired reads of the same DNA molecule usually align to the same contig, and incorrect assembly of the contig may result in a large number of orphan reads aligning to different contigs.

[0102] Therefore, extract the read pair abnormal signals in the candidate error regions, specifically:

[0103] Detect the orphan reads in the candidate error regions site by site. If paired reads that cannot be aligned to the same assembly result are found at a certain site, the local region of this site is recorded as the orphan read region. Count the number of abnormal paired reads in the orphan read region and the total number of reads in the orphan read region. Take the ratio of the number of reads ( ) to the total number ( ) as the read pair abnormal signal in the orphan read region ( );

[0104] Compare the read pair abnormal signals of all orphan read regions within the candidate error region. Preferably, select the "orphan read region" where the proportion of all abnormal paired reads is greater than 20%. Take the read pair abnormal signal with the maximum value in the orphan read region as the read pair abnormal signal in this candidate error region ( ).

[0105] It can be achieved through the following formula:

[0106] ;

[0107] .

[0108] 6. Regarding clipping anomalies

[0109] During the alignment process, reads of some sequences that cannot be fully aligned to the contig usually occur in complex regions of the genome, such as insertion, deletion, or duplication regions. Clipping means that when the aligned reads are being aligned, there are unaligned sequences at one or both ends.

[0110] Therefore, extract the clipping abnormal signals in the candidate error regions, specifically:

[0111] Count the number of reads with clipping at each site in the candidate error region. If the coverage at a certain site is not less than the preset value, calculate the ratio of the number of reads ( ) to the coverage at this site ( ). Take this ratio as the clipping abnormal signal at this site ( );

[0112] Compare the clipping abnormal signals of all sites in the candidate error region. Take the maximum ratio as the clipping abnormal signal in this candidate error region ).

[0113] It can be achieved through the following formula:

[0114] ;

[0115] ;

[0116] Among them, for a certain low-coverage site, when its coverage is less than the preset value, such as less than 0.2 times the average coverage, it may lead to untrustworthy analysis results, thereby increasing the false positive risk. To avoid misjudgment, such sites are not considered in subsequent analysis.

[0117] S2: For each candidate error region, based on multiple characteristic signals of the extracted candidate error region, after clustering, the clustering results of each characteristic signal of the corresponding candidate error region are obtained; if the clustering result of one characteristic signal is abnormal, then the candidate error region is an incorrect assembly region.

[0118] In the case where all candidate error regions are assembly errors, since clustering such as K-means may misclassify some real incorrect assemblies as normal. To avoid this problem, before clustering, it also includes:

[0119] Randomly select several regions from the assembly results longer than the preset sequence length and perform abnormal signal clustering analysis together with the candidate error regions. A certain number of normal regions are introduced into the candidate dataset to assist clustering, so as to enhance the discrimination ability of incorrect assembly detection.

[0120] Clustering can group objects with similar characteristics together. K-means clustering aims to divide into two categories: normal regions and incorrect regions according to the abnormal signals extracted from candidate error regions. K-means performs well in low dimensions, and the selected abnormal signals are independent of each other. Therefore, each candidate error region is separately clustered according to its multiple signal characteristics (such as abnormal coverage, indels, incorrect pairing direction, etc.) and the clustering results of the corresponding characteristic signals are output.

[0121] Analyze the clustering results of multiple abnormal signal characteristics of the candidate error region respectively, and the abnormal characteristic information extracted from each region, and conduct a comprehensive evaluation of the assembly errors of multiple signals. For each candidate error region of the genome, independent clustering results of six abnormal signals (coverage, indels, pairing direction, insert size, paired reads, and clipping) of this region are obtained { , , , , , }, where 0 represents normal and 1 represents abnormal. If the clustering results of these abnormal signals are all normal, it is considered that there is no abnormality in this region; compared with the signals in the correctly assembled region, the signals in the misassembled region show significant characteristic differences. Clusters containing real assembly errors often show significant performance in a certain signal feature. If at least one abnormal signal is clustered as abnormal, it is determined that there is an assembly error in this region. By defining an overall evaluation function Fregion to summarize the K-means clustering results of multiple signals and give the final abnormality evaluation, the accuracy of assembly error detection is improved, and the missed detection or misjudgment caused by the misjudgment of a single signal is avoided.

[0122] .

[0123] S3: Based on the misassembled region, determine the type of assembly error. If the type is insertion or deletion, remove the inserted sequence or the inserted and deleted fragment;

[0124] If the type is misconnection, through multiple sequence alignment, generate a consensus sequence, compare the consensus sequence with the assembly result, and perform repair.

[0125] Problems such as insertion, deletion, or misjoin may exist in the assembled genome. These assemblies will damage the gene structure and thus affect downstream functional analysis. Therefore, after accurately locating the misassembled region, different repair strategies are adopted for different error types to improve the consistency and accuracy of the assembled sequence.

[0126] Insertion or deletion errors usually show significant differences from the reads alignment signals. In the BAM file, the reads sequences supporting the insertion or deletion are extracted by parsing the CIGAR signal, the proportion of reads supporting this error is calculated, and the abnormal region is marked. The support ratio of the reads supporting the assembly error is statistically analyzed. If the support ratio at a certain position exceeds the preset threshold (such as >0.5), the correction operation is performed, that is, for insertion errors, the redundant sequence in the contig is removed; for deletion errors, the deleted fragment is inserted to restore integrity. If multiple indel error signals are detected and their support frequencies do not exceed the threshold, no indel error repair process is performed.

[0127] Misjoin usually occurs during the assembly process. Common forms include cross-chromosome splicing, reverse splicing, or overly long-distance splicing. These errors may be caused by the assembly algorithm or sequence complexity, and are manifested as abnormal contig breakpoints, alignment directions, or overly long splicing distances in the alignment results.

[0128] To accurately identify the misconnected regions, first analyze and compare the characteristic signals (such as an increase in clipping events, a drastic change in coverage, etc.), and mark the breakpoint regions (1.5 kb upstream and downstream are used as the analysis range);

[0129] Subsequently, based on the paired-end sequencing data, extract the paired-end reads that completely cover the breakpoint regions. After filtering, obtain the high-precision local read sequences;

[0130] Use the SPOA tool to perform multiple sequence alignment on the local read sequences, generate a consensus sequence, and evaluate its integrity and consistency. Compare the consensus sequence with the assembly result. If the successfully aligned length of the repaired sequence accounts for more than 90% of the breakpoint region, the repair is considered effective. Use the generated PAF file from the alignment to replace the misconnected region in the assembly result with the consensus sequence to ensure the continuity and accuracy of the overall sequence.

Claims

1. A method for detecting and repairing incorrect assemblies based on paired-end sequencing data, characterized in that It includes the following steps: S1: Based on the assembly result, with a preset suspicious site as the core, extend a first preset length to both sides of the genomic sequence, and use it as the candidate error region; Compare the paired-end sequencing data with the assembly result. Based on the comparison information, extract multiple characteristic signals of the candidate error region. The multiple characteristic signals are coverage anomaly signal, insertion or deletion anomaly signal, paired direction anomaly signal, insert size anomaly signal, read pair anomaly signal, and splicing anomaly signal; S2: For each candidate error region, based on the multiple characteristic signals of the candidate error region extracted, after clustering, obtain the clustering result of each characteristic signal of the corresponding candidate error region; if the clustering result of one characteristic signal is abnormal, then this candidate error region is an incorrectly assembled region; S3: Based on the incorrectly assembled region, determine the type of assembly error. If the type is insertion or deletion, remove the inserted sequence or the inserted and deleted fragment; If the type is incorrect connection, generate a consensus sequence through multiple sequence alignment, compare the consensus sequence with the assembly result, and perform repair.

2. The method for detecting and repairing incorrect assembly based on paired-end sequencing data according to claim 1, wherein For the above S1, the specific method for extracting the coverage anomaly signal of the candidate error region is: Taking the candidate error region as the benchmark, extend a second preset length to both sides of it, and use it as the extended region; Calculate the average coverage of all sites in the extended region, traverse and detect the coverage of each site in the candidate error region. If the coverage of a certain site is not less than the first preset multiple of the average coverage or not greater than the second preset multiple of the average coverage, then this site is marked as a coverage anomaly site; Count the total number of coverage anomaly sites in the candidate error region, calculate the ratio of the total number to the total length of the candidate error region, and use this ratio as the coverage anomaly signal.

3. The method for detecting and repairing incorrect assembly based on double-ended sequencing data according to claim 1, wherein For the above S1, the specific method for extracting the insertion or deletion anomaly signal of the candidate error region is: Count the number of reads with insertions or deletions at each site in the candidate error region one by one. If the coverage of a certain site is not less than the preset value, then calculate the ratio of the number of reads to the coverage of this site, and use this ratio as the insertion or deletion signal of this site; Compare the insertion or deletion signals of all sites in the candidate error region, and use the maximum ratio as the insertion or deletion anomaly signal of this candidate error region.

4. The method for detecting and repairing incorrect assembly based on paired-end sequencing data according to claim 1, wherein, For the above S1, the specific method for extracting the paired direction anomaly signal of the candidate error region is: Detect the reads with mismatched directions at each site in the candidate error region one by one. If reads with mismatched directions are found at a certain site, then mark the local region of this site as the incorrectly paired direction region, count the number of reads with mismatched directions in the incorrectly paired direction region and the total number of reads in the incorrectly paired direction region, and use the ratio of the number of reads to the total number as the paired direction anomaly signal of this incorrectly paired direction region; Compare the sizes of the paired direction anomaly signals of all incorrectly paired direction regions in the candidate error region, and use the paired direction anomaly signal of the incorrectly paired direction region with the maximum value as the paired direction anomaly signal of this candidate error region.

5. The method for detecting and repairing incorrect assembly based on paired-end sequencing data according to claim 1, wherein For the above S1, the specific method for extracting the insert size anomaly signal of the candidate error region is: Detect the reads with abnormal insert sizes at each site in the candidate error region. If reads with abnormal insert sizes are found at a certain site, mark the local region of this site as the region with abnormal insert sizes, count the number of reads with abnormal insert sizes in the region with abnormal insert sizes and the total number of reads in the region with abnormal insert sizes, and use the ratio of the number of reads to the total number as the insert size abnormal signal of this region with abnormal insert sizes; Compare the insert size abnormal signals of all regions with abnormal insert sizes in the candidate error region, and use the insert size abnormal signal of the region with the maximum value as the insert size abnormal signal of this candidate error region.

6. The method for detecting and repairing incorrect assemblies based on paired-end sequencing data according to claim 1, wherein In step S1, extract the read pair abnormal signal of the candidate error region, specifically: Detect the orphan reads at each site in the candidate error region. If paired reads that cannot be aligned to the same assembly result are found at a certain site, mark the local region of this site as the orphan read region, count the number of reads with paired anomalies in the orphan read region and the total number of reads in the orphan read region, and use the ratio of the number of reads to the total number as the read pair abnormal signal of the orphan read region; Compare the read pair abnormal signals of all orphan read regions in the candidate error region, and use the read pair abnormal signal of the orphan read region with the maximum value as the read pair abnormal signal of this candidate error region.

7. The method for detecting and repairing incorrect assembly based on paired-end sequencing data according to claim 1, wherein In step S1, extract the splicing abnormal signal of the candidate error region, specifically: Count the number of reads with splicing at each site in the candidate error region. If the coverage of a certain site is not less than the preset value, calculate the ratio of the number of reads to the coverage of this site, and use this ratio as the splicing abnormal signal of this site; Compare the splicing abnormal signals of all sites in the candidate error region, and use the maximum ratio as the splicing abnormal signal of this candidate error region.

8. The method for detecting and repairing misassemblies based on paired-end sequencing data according to claim 1, wherein In step S2, before clustering, it further includes: Randomly select several regions from the assembly results longer than the preset sequence length and perform abnormal signal clustering analysis together with the candidate error region.

9. The method for detecting and repairing incorrect assembly based on paired-end sequencing data according to claim 1, wherein In step S3, if the type is misconnection, generate a consensus sequence through multiple sequence alignment, and compare the consensus sequence with the assembly result for repair, specifically: Based on the assembly result of the misconnection type, analyze the alignment characteristic signals and mark the breakpoint region; Based on the paired-end sequencing data, extract the paired-end reads that completely cover the breakpoint region, and after filtering, obtain the local read sequence; The local read sequence generates a consensus sequence through multiple sequence alignment. Compare the consensus sequence with the assembly result. If the successfully aligned length of the repaired sequence accounts for more than 90% of the breakpoint region, replace the misconnected region in the assembly result with the consensus sequence.

10. An error assembly detection and repair system based on double-ended sequencing data, characterized in that It includes: Feature signal extraction module: Based on the assembly result, with a preset suspicious site as the core, extend a first preset length to both sides of the genomic sequence as the candidate error region; Compare the paired-end sequencing data with the assembly result, and based on the comparison information, extract multiple feature signals of the candidate error region. The multiple feature signals are coverage abnormal signal, insertion or deletion abnormal signal, paired direction abnormal signal, insert size abnormal signal, read pair abnormal signal, splicing abnormal signal; Incorrect assembly region determination module: For each candidate incorrect region, after clustering based on multiple characteristic signals of the extracted candidate incorrect region, the clustering results of each characteristic signal of the corresponding candidate incorrect region are obtained; if the clustering result of one characteristic signal is abnormal, then this candidate incorrect region is an incorrect assembly region; Repair module: Based on the incorrect assembly region, determine the type of assembly error. If the type is insertion or deletion, then remove the inserted sequence or the inserted and deleted fragment; If the type is incorrect connection, then through multiple sequence alignment, generate a consensus sequence, compare the consensus sequence with the assembly result, and perform repair.