Automatic Sanger sequencing result analysis tool based on peak-valley signal relationship

Abstract

The invention discloses an automatic Sanger sequencing result analysis tool based on a peak-valley signal relationship, which comprises three-stage screening of accurate definition and multi-criterion combination of peak-valley signals, coverage of preliminary screening, a quality score window and iterative denoising, improves the recognition accuracy of Soulder sleeve peak and G peak interference, avoids mutation detection errors caused by missed judgment or misjudgment under special conditions, and improves the accuracy of Sanger sequencing result analysis. In addition, the proportion of effective data in the credible section is remarkably increased, interference of invalid data on follow-up analysis is greatly reduced, automatic high-throughput analysis can be achieved, the method does not depend on personnel experience, misjudgment caused by personnel experience difference, fatigue or negligence can be eliminated, and labor cost is reduced.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology, specifically relating to an automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships. Background Technology

[0002] Sanger sequencing technology, considered the "gold standard" for gene sequence analysis, is widely used in fields such as single-gene disease diagnosis, tumor target mutation verification, and gene editing efficiency testing due to its high accuracy. Its core output is a .ab1 format file containing the base sequence, signal intensity, and sequencing quality score. Efficient and accurate analysis of .ab1 files is crucial for obtaining reliable genetic information. SangerReport is a native tool specifically designed for automated analysis of Sanger sequencing results. It aims to provide a concise, comprehensive, automated, and high-throughput analysis of sequencing results, utilizing peak-valley signal relationships to efficiently and accurately identify insertions / deletions, fragment duplications, independent peak overlap, shoulder overlap, true mutations, G-peak interference, and other issues.

[0003] Existing automated analysis tools for Sanger sequencing results have significant technical deficiencies: 1. Lack of ability to handle special abnormal signals. Taking SangerReport as an example, it can only identify basic mutations and regular peak overlaps, relying on complex CNNs, SVMs, or bimodal probability models, but it cannot handle three key special cases—Shoulder peak overlaps (the base type of the overlap at this site is the same as the main peak at the adjacent site, but it is missed due to signal tailing characteristics), G peak interference (G bases on both sides of this site cause G-type background signals to be misjudged as true overlaps), and long fragment duplications (repetitive fragments >10bp unrelated to the target sequence are missed), resulting in a mutation detection accuracy of less than 70%; 2. Crude judgment of reliable regions. Most tools rely only on a single quality factor. The score threshold screening of sequences does not take into account the clustering effect of suspicious signal sites, resulting in more than 15% invalid data in the reliable segment and a high loss rate of valid sequences; 3. High dependence on manual review, more than 30% of suspected results need to be manually reviewed during high-throughput analysis. Due to differences in personnel experience, the misjudgment rate reaches 15%, and the analysis efficiency is low (≤800 files / hour), which cannot meet the high-throughput needs of clinical and scientific research.

[0004] Therefore, there is a need in this field to develop an automated tool for analyzing Sanger sequencing results, which can analyze Sanger sequencing results in a concise, comprehensive, automated, and high-throughput manner, thereby improving analysis efficiency. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention provides an automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships, comprising the following steps:

[0006] a. Read the .ab1 file generated by Sanger sequencing and extract the base sequence B = {b1, b2, ..., b}. n}, Signal sequence S = {s1, s2, ..., s m} and the quality score sequence Q = {q1, q2, ..., q n};

[0007] b. Determine the reliable section:

[0008] b1. Retain 30-700 bp segments as potential reliable regions;

[0009] b2. Calculate the average quality score of a continuous 6-bp window downstream from 30bp. When the average quality score is greater than 30 for the first time, take the starting position of the window + 2 as the starting point. Calculate the average quality score of a continuous 6-bp window upstream from 700bp. When the average quality score is greater than 30 for the first time, take the starting position of the window - 2 as the ending point.

[0010] b3. Check the number of suspicious signal sites at the start end (Start~Start+40) and the end end (End~End+40). If there are more than 4, adjust the start and end sites. After 15 iterations, determine the final reliable segment [Start_final,End_final].

[0011] c. Identifying anomalous signals based on peak-valley signal relationships: Calculating the peak s at base sites. p With troughs t After eliminating background noise, independent peak overlap and Shoulder peak overlap are distinguished by multiple criteria, and G peak interference is eliminated.

[0012] d. Detection of fragment duplication and insertion / deletion: Identify fragment duplications >10bp within a reliable region and align them with the target sequence, allowing for insertions / deletions of ≤2 consecutive sites;

[0013] e. Output an analysis report containing the location and type of the anomaly.

[0014] Specifically, in step a, the .ab1 file is parsed using BioPython's SeqIO module, and each b of the base sequence B is analyzed. i The corresponding core signal site s in the signal sequence S j The quality score sequence Q originates from the "QUAL" data field in the .ab1 file.

[0015] Specifically, in step b2, the average quality score of the 6 consecutive bp windows is calculated as follows:

[0016] For the starting end: Starting from the base position corresponding to 30bp, let the starting calculation position be i. Select six consecutive bases downstream, including the current position i and the next five bases (i.e., the bases corresponding to positions i, i+1, i+2, i+3, i+4, i+5, i+6). Extract the quality score value corresponding to each of these six bases. Add the six values ​​together and divide by 6 to obtain the average quality score of the consecutive 6bp window.

[0017] For the termination end: Starting from the base position corresponding to 700bp, let the starting calculation position be i. Select six consecutive bases upstream, including the current position i and the five bases before it (i.e., the bases corresponding to positions i-6, i-5, i-4, i-3, i-2, i-1, and i). Extract the quality score value corresponding to each of these six bases. Add the six values ​​together and divide by 6 to obtain the average quality score of the consecutive 6bp window.

[0018] Specifically, in step b3, the suspected signal site meets any of the following conditions: qi < 20, or the peak-to-valley ratio of the signal sequence corresponding to the site is < 1.05; when adjusting the start and end sites, the start end is the next site of the suspected signal site closest to the middle section, and the end end is the previous site of the suspected signal site closest to the middle section.

[0019] Specifically, in step c, peak s p With troughs t The calculation method is as follows:

[0020] peaks p =max(s j-3 ,s j-2 ,s j-1 ,s j ,s j+1 ,s j+2 ,s j+3 );

[0021] troughs t =max(min(s) j-3 ,s j-2 ,s j-1 ,s j ),min(s j+1 ,s j+2 ,s j+3 ));

[0022] Where j is the position of base site i in the signal sequence.

[0023] Specifically, in step c, the condition for background noise removal is: base secondary peak b. i，top2 <0.15×base main peak b i，top1 .

[0024] Specifically, in step c, the determination of independent peaks satisfies any of the following conditions:

[0025] (1)b i，top2 >0.3×b i，top1 ;

[0026] (2) The peak S of the signal secondary peak p,top2 ≥0.3×b i，top1 And S p,top2 >1.15×S t,top2 .

[0027] Specifically, in step c, the determination of Shoulder peaks satisfies any of the following conditions:

[0028] (1)S p,top2 ≥0.3×b i，top1 Furthermore, the base type of the secondary signal peak is similar to that of b. i-1 or b i+1 The main peak base types are consistent;

[0029] (2)S p,top2 >45、S t,top2 >40、b i，top2 >3×b i，top3 And the minimum strength of the Shoulder segment of the signal is >10.

[0030] Specifically, in step c, the condition for eliminating G-peak interference is: b i-1 Type = G, b i+1 Type = G, Signal secondary peak type = G.

[0031] Specifically, in step d, the length of the repeated fragment L>10bp, and the allowable range of insertions and deletions is ≤2 consecutive sites; if there are ≥3 consecutive insertions and deletions, it is marked as an aberrant insertion or deletion.

[0032] The present invention has the following beneficial effects:

[0033] 1. The automated tool for Sanger sequencing result analysis provided by this invention improves the accuracy of identifying Shoulder peak overlap and G peak interference from 60%–70% in existing tools to over 95% through precise definition of peak and valley signals and combination of multiple criteria, avoiding mutation detection errors caused by missed or misjudged cases. Through three-level screening of initial screening, quality score window, and iterative noise reduction, it significantly increases the proportion of valid data in reliable regions, greatly reduces the interference of invalid data on subsequent analysis, and improves analysis efficiency.

[0034] 2. The Sanger sequencing result analysis automation tool provided by this invention can support the sharing of more than 1,000 .ab1 files per hour, enabling automated high-throughput analysis. The analysis results of different batches of data are consistent, and it does not rely on human experience. It can also eliminate misjudgments caused by differences in human experience, fatigue or negligence, thus reducing labor costs. Attached Figure Description

[0035] Figure 1 A flowchart illustrating the automated tool for Sanger sequencing result analysis provided by this invention;

[0036] Figure 2 The chromatogram shows independent peaks.

[0037] Figure 3 Chromatogram of shoulder peaks;

[0038] Figure 4 The chromatogram shows interference from the G peak. Detailed Implementation

[0039] The technical solution of the present invention will be clearly and completely described below. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0040] In the following examples, all reagents used are commercially available unless otherwise specified. Where specific techniques or conditions are not described in the examples, they should be understood as being performed in accordance with techniques or conditions described in the prior art, or in accordance with the instructions for use of commercial products.

[0041] Example 1

[0042] This embodiment illustrates the analysis process of an automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships. Please refer to [link / reference]. Figure 1 Specifically, it includes the following steps:

[0043] Trusted segment analysis:

[0044] Step a: Read the sequencing results file in .ab1 format and retain a 30-700bp segment as a potential reliable segment.

[0045] Step b: Based on the quality score provided in the .ab1 file, further remove segments with low quality scores at the upstream 5' end and the downstream 3' end. To define the start site, calculate the average quality score of a 6-bp window from base i to i+6 downstream from 30bp to see if it is greater than 30. If it is, take base i+2 as the start site. To define the end site, calculate the average quality score of a 6-bp window from base i-6 to i upstream from 700bp to see if it is greater than 30. If it is, take base i-2 as the end site.

[0046] Step c: When there are more than 10 suspicious signal sites within a consecutive 100bp range at the start / end, even if the quality score meets the requirements of step b, it still indicates a poor signal-to-noise ratio. Therefore, it is necessary to further narrow down the range of the reliable segment. Specifically, the next site (for the start end) or the previous site (for the end end) of the suspicious signal site closest to the middle segment is taken as the start / end site.

[0047] After the above steps are completed, anomaly type analysis will be performed to propose peak-valley signal relationships based on local regions. This analysis will be used to identify whether the suspicious signal location belongs to the category of independent peak overlap, Shoulder peak overlap, or G-peak interference. The analysis includes the following steps:

[0048] Step 1: The Sanger sequencing chromatogram data contains two types of sequences: one type is the base sequence B = {b1, b2, ..., b}. n}, where b i The base representing the i-th position includes four types: A, T, C, and G. Another type is the actual signal sequence S = {s1, s2, ..., s...} measured by the instrument. m The base sequence is a subset of the signal sequence, and the signal sequence reflects the actual waveform of each base site.

[0049] Step 2: Assume that the i-th site of the base sequence corresponds to the j-th site of the signal sequence. To determine the actual waveform of the base sequence at the i-th site, it is necessary to analyze the signal intensity of the signal sequence within the window [j-3, j+3].

[0050] Define the base sequence b i peaks p :

[0051] s p =max(s j-3,s j-2 ,s j-1 ,s j ,s j+1 ,s j+2 ,s j+3 );

[0052] Define the base sequence b i troughs t :

[0053] s t =max(min(s) j-3 s j-2 s j-1 s j ), min(s j+1 s j+2 s j+3 ));

[0054] Step 3, compare the secondary base peak b. i，top2 With the main base peak b i，top1 If the height of the peak is less than 0.15 of the height of the peak, then the secondary peak and others are considered to be background noise and no further peak matching is performed.

[0055] Step 4: Determine independent peak overlap. Independent peak overlap is defined as the presence of a distinct secondary peak below the main peak at the same base sequence site. Specifically, independent peak overlap is determined to exist if either of the following two criteria is met:

[0056] Criterion 4.1, base subpeak b i，top2 Exceeding the main base peak b i，top1 30%, that is:

[0057] b i，top2 >0.3*b i，top1 ;

[0058] Criterion 4.2: The peak height of the secondary signal peak within the window exceeds that of the main base peak b. i，top1 30% of the peak value, and its peak value is 1.15 times higher than its trough value, that is:

[0059] s p，top2 ≥0.3*b i，top1 and

[0060] s p，top2 >1.15*s t，top2 ;

[0061] Step 5: Determine Shoulder overlap. Shoulder overlap is defined as a peak at this site that belongs to the same base type as the main peak at an adjacent site. In this case, the shape and height of the overlap peak may be ignored because it does not meet the requirements of Step 4, so additional criteria are needed.

[0062] Criterion 5.1: The peak height of the secondary signal peak within the window exceeds that of the main base peak b. i，top1 30% of the peaks are of the same base type as the main peaks at adjacent sites, i.e.:

[0063] s p，top2 ≥0.3*b i，top1 and

[0064] s j，type =b i-1，type or b i+1，type ;

[0065] Criterion 5.2: The peak height of the signal secondary peak within the window is greater than 45, the trough height is greater than 40, and the base secondary peak b... i，top2 Higher than the third base peak b i，top3 Three times that of the signal, and the height of the Shoulder segment is always greater than 10, that is:

[0066] s p，top2 >45, s t，top2 >40, and

[0067] b i，top2 >3*b i，top3 and

[0068] min(s j，top2 ~s k，top2 )>10, for s k =b i-1 or b i+1 ;

[0069] Step 6, determine G-peak interference. G-peak interference is defined as a site flanked by G bases, resulting in a higher background signal of type G at that site. This is not a true overlay and should be excluded from steps 4 and 5.

[0070] b i-1，type =G,b i+1，type =G,s j，top2，type =G;

[0071] For any abnormal signals (independent peak overlap, shoulder peak overlap, and G-peak interference), report them to the user for verification. Please refer to the following: Figure 2-4 The figures below are chromatographic diagrams of independent peaks, shoulder peaks, and G-peak interference, respectively. The horizontal axis represents base position, and the vertical axis represents signal intensity. An independent peak refers to a site where, in addition to the main peak, there is a distinct secondary peak. Figure 2 In the middle, an independent secondary peak C appeared below the main peak G; Figure 3In this context, "Shoulder overlap" refers to a signal whose overlapping peak signal is identical to the adjacent main peak signal, affecting this location like a "shoulder." Figure 4 In this context, G-peak interference refers to a strong G-background signal that may appear at a location sandwiched between two main G peaks. However, this background signal is not a true overlay peak and should therefore be excluded from the error category. It can be seen that the analytical tools provided in this application are capable of correctly distinguishing and judging these interferences.

[0072] Analyze the .ab1 file for duplicate segments that do not match the target sequence, including:

[0073] 1. Check for insertions / deletions in the .ab1 file: The reliable regions extracted from the .ab1 file are aligned to the target sequence. This alignment allows for insertions / deletions at a maximum of two consecutive sites. When an insertion / deletion is found, this method will report the .ab1 data where the insertion / deletion occurred and its location; when there are at least three consecutive insertions / deletions, this method will report that the .ab1 file is erroneous and request manual review.

[0074] 2. Check for erroneous segment duplication in the .ab1 file: To analyze whether there are erroneous segment duplications in the .ab1 file, this method traverses the sequences in the .ab1 file, records all segment duplications longer than 10bp within the confidence region, and compares them with the target sequence. If the target sequence does not have the corresponding segment duplication, it means that the experimental data corresponding to .ab1 has erroneous segment duplication. In this case, this method will report the .ab1 data with segment duplication and the location of the error for manual verification.

[0075] In summary, the above embodiments are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. An automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships, characterized in that, Includes the following steps: a. Read the .ab1 file generated by Sanger sequencing and extract the base sequence B = {b1, b2, ..., b}. n }, Signal sequence S = {s1, s2, ..., s m } and the quality score sequence Q = {q1, q2, ..., q n }; b. Determine the reliable section: b1. Retain 30-700 bp segments as potential reliable regions; b2. Calculate the average quality score of a continuous 6-bp window downstream from 30bp. When the average quality score is greater than 30 for the first time, take the starting position of the window + 2 as the starting point. Calculate the average quality score of a continuous 6-bp window upstream from 700bp. When the average quality score is greater than 30 for the first time, take the starting position of the window - 2 as the ending point. b3. Check the number of suspicious signal sites at the start end (Start~Start+40) and the end end (End~End+40). If there are more than 4, adjust the start and end sites. After 15 iterations, determine the final reliable segment [Start_final,End_final]. c. Identifying anomalous signals based on peak-valley signal relationships: Calculating the peak s at base sites. p With troughs t After eliminating background noise, independent peak overlap and Shoulder peak overlap are distinguished by multiple criteria, and G peak interference is eliminated. d. Detection of fragment duplication and insertion / deletion: Identify fragment duplications >10bp within a reliable region and align them with the target sequence, allowing for insertions / deletions of ≤2 consecutive sites; e. Output an analysis report containing the location and type of the anomaly.

2. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step a, the .ab1 file is parsed using BioPython's SeqIO module, and each b of the base sequence B is analyzed. i The corresponding core signal site s in the signal sequence S j The quality score sequence Q originates from the "QUAL" data field in the .ab1 file.

3. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step b2, the average quality score for a consecutive 6bp window is calculated as follows: For the starting end: Starting from the base position corresponding to 30bp, let the starting calculation position be i. Select six consecutive bases downstream, including the current position i and the next five bases (i.e., the bases corresponding to positions i, i+1, i+2, i+3, i+4, i+5, i+6). Extract the quality score value corresponding to each of these six bases. Add the six values ​​together and divide by 6 to obtain the average quality score of the consecutive 6bp window. For the termination end: Starting from the base position corresponding to 700bp, let the starting calculation position be i. Select six consecutive bases upstream, including the current position i and the five bases before it (i.e., the bases corresponding to positions i-6, i-5, i-4, i-3, i-2, i-1, and i). Extract the quality score value corresponding to each of these six bases. Add the six values ​​together and divide by 6 to obtain the average quality score of the consecutive 6bp window.

4. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step b3, the suspected signal site satisfies any of the following conditions: q i <20, or the peak-to-valley ratio of the signal sequence corresponding to the site is <1.05; when adjusting the start and end sites, the start end is the next site after the signal doubt site closest to the middle section, and the end end is the previous site of the signal doubt site closest to the middle section.

5. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step c, peak s p With troughs t The calculation method is as follows: peaks p =max(s j-3 s j-2 s j-1 s j s j+1 s j+2 s j+3 ); troughs t =max(min(s) j-3 s j-2 s j-1 s j ), min(s j+1 s j+2 s j+3 )); Where j is the position of base site i in the signal sequence.

6. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step c, the condition for background noise removal is: base secondary peak b i，top2 <0.15×base main peak b i，top1 .

7. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step c, the determination of independent peaks satisfies any of the following conditions: (1)b i，top2 >0.3×b i，top1 ; (2) The peak S of the signal secondary peak p,top2 ≥0.3×b i，top1 And S p,top2 >1.15×S t,top2 .

8. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step c, the determination of Shoulder peaks satisfies any of the following conditions: (1)S p,top2 ≥0.3×b i，top1 Furthermore, the base type of the secondary signal peak is similar to that of b. i-1 or b i+1 The main peak base types are consistent; (2)S p,top2 >45、S t,top2 >40、b i，top2 >3×b i，top3 And the minimum strength of the Shoulder segment of the signal is >10.

9. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step c, the condition for eliminating G-peak interference is: b i-1 Type = G, b i+1 Type = G, Signal secondary peak type = G.

10. The automated tool for analyzing Sanger sequencing results based on peak-valley signal relationships according to claim 1, characterized in that, In step d, the length of the repeated fragment L>10bp, and the allowable range for insertion / deletion is ≤2 consecutive sites; if there are ≥3 consecutive insertion / deletion sites, it is marked as an aberrant insertion / deletion.

Images