Method for quality assessment and screening of nucleic acid sequencing data

Abstract

The application provides a nucleic acid sequencing data quality evaluation method and screening method. The difference between a sequencing signal and a corrected signal is used to evaluate base sequencing quality or to screen nucleic acid sequencing data. In particular, the difference between a dephasing corrected signal and an error correction corrected signal is used. This is the first quality evaluation method and screening method specially applicable to error correction code sequencing. The method of the application is also applicable to sequencing data obtained by other sequencing methods, and has a wide range of applications.

Description

TECHNICAL FIELD

[0001] The present application relates to a method for quality evaluation and screening of nucleic acid sequencing data, and belongs to the field of gene sequencing. BACKGROUND

[0002] High-throughput sequencing is a technology that simultaneously determines the sequences of thousands of DNA sequences, which is widely used in basic biology and medical research, as well as in vitro diagnosis. Generally, in addition to outputting the measured DNA sequence, a high-throughput sequencer also gives each base a quality value to represent the measurement accuracy of the base. Usually, this quality value is given in the form of Phred, that is, when the measurement accuracy of the base is D, the given quality value is:

[0003] q = -10 log 10 (1-p)

[0004] For example, a quality value of 10 corresponds to an accuracy of 90%, a quality value of 20 corresponds to an accuracy of 99%, and a quality value of 30 corresponds to an accuracy of 99.9%.

[0005] Currently, the quality value evaluation method commonly used in the field of high-throughput sequencing is the Phred algorithm (Ewing B, Green P. Base-calling of automated sequencer traces using Phred. II. Error probabilities. Genome Res. 8: 186-194 (1998)). The Phred algorithm is widely used in the field of DNA sequencing, but different sequencing technologies use different predictors. For example, IonTorrent uses the following predictors: the position of the base in the sequence; the length of the homopolymer where the base is located; the degree to which the base is close to an integer in the sequencing signal; the degree of phase shift of the base in the sequencing signal, etc. Illumina has also disclosed a quality evaluation method using a neural network, that is, a neural network is trained to predict the quality of a base through the values of the predictors. The quality predictor values of Illumina include online overlap, purity, phasing, start5, hexamer score, motif accumulation, endiness, approximate homopolymer, intensity decay, second-to-last purification, signal overlap with background (SOWB), and / or offset purity G adjustment, etc. (Patent CN202080005431.0).

[0006] A good scoring algorithm should have high discrimination, that is, it can filter out high-quality bases and give them high scores, and filter out low-quality bases and give them low scores. However, when using existing algorithms to evaluate the quality of sequencing data of error correction code sequencing, there is a problem of insufficient / low discrimination. Therefore, it is necessary to develop a base quality evaluation method suitable for error correction code sequencing data and having higher discrimination.

[0007] A high-throughput sequencer can generate a large amount of data in one run. However, due to defects and noise interference in optical imaging, signal acquisition, DNA amplification to produce polyclones, chemical reactions, mechanical systems, algorithm processing, etc., not all generated data is high-quality data. If the data provided to the user after downloading contains low-quality and unreliable sequencing data, it may affect the user's analysis of the sample, causing false positives or false negatives, and in severe cases, it may affect the clinical diagnosis result and delay the disease. Therefore, screening, retaining and outputting high-quality sequencing data is a necessary step for high-throughput sequencing.

[0008] Illumina uses "purity" to screen its sequencing data (US Patent US13006206). Specifically, Illumina labels four bases with different fluorescence, and the base with the strongest fluorescence during sequencing is the measured base, and the purity is defined as the strongest fluorescence intensity divided by the sum of the four fluorescence intensities, and sequences with purity higher than a given threshold are considered high-quality sequences. Ion torrent uses two indicators to screen its sequencing data: SSQ (Euclidean distance between corrected signal and its integer value) and PPF (proportion of non-zero signals) (patent EP19181402, WOUS11067959), and sequences with both indicators less than a given threshold are considered high-quality sequences.

[0009] In practical applications, high-throughput sequencers are generally used to screen sequencing data by the quality values of the bases of the measured DNA sequences. For example, the fastp software will perform a moving average on the quality values of each sequence, and will cut off the sequence at a low average. That is, the screening process must be performed after obtaining the quality values corresponding to the bases, which increases the waiting time, and there is a need to develop a sequencing data screening method that does not depend on the quality values of the bases. SUMMARY

[0010] In one aspect, the present application provides a method for quality evaluation of nucleic acid sequencing data, characterized in that it comprises:

[0011] a) sequencing a reference nucleic acid sample to obtain a set of sequencing signals s; sequencing a nucleic acid sample to be measured using the same sequencing method to obtain a set of sequencing signals s';

[0012] b) signal correction is performed on the sequencing signal s to obtain a corrected signal c, which can be directly converted into a base sequence;

[0013] c) the converted base sequence is aligned to a reference sequence to obtain an alignment result, and bases in the high-quality aligned base sequence are marked as sequencing correct or sequencing error according to the alignment result;

[0014] d) the signal difference of the corresponding part of the sequencing correct or sequencing error base in the sequencing signal s and the corrected sequence c is compared, and the relationship between the signal difference and the base sequencing accuracy is established;

[0015] e) signal correction is performed on the sequencing signal s' in the same way as step b) to obtain a corrected signal c', and then c' is converted into a base sequence;

[0016] f) the signal difference of the corresponding part of each base in the converted base sequence in s' and c' is compared, and the sequencing accuracy of the base is predicted by using the relationship between the signal difference and the base accuracy established in step d).

[0017] According to a preferred embodiment, the type of signal generated by sequencing can be optical signal or electrical signal; the optical signal can be monochromatic or polychromatic.

[0018] According to a preferred embodiment, the nucleic acid to be tested includes DNA, or RNA or other types of nucleic acid molecules.

[0019] According to a preferred embodiment, the genomic sequence of the species to which the reference nucleic acid sample belongs is known; when the reference nucleic acid is DNA, the reference sequence is the genomic sequence of the species to which the reference nucleic acid belongs; when the reference nucleic acid is RNA, the reference sequence is the transcriptome sequence of the species to which the reference nucleic acid belongs.

[0020] According to a preferred embodiment, the converted base sequence is aligned to a reference sequence to obtain an alignment result, and high-quality aligned base sequences are further screened from the alignment result, and the bases in the high-quality aligned base sequences are marked as sequencing correct or sequencing error; the quality value range needs to be selected according to the used alignment software or algorithm.

[0021] According to a preferred embodiment, the sequencing signal is the signal directly collected by the sequencer or the normalized signal; correspondingly, the correction process of the sequencing signal can be the correction of the out-of-phase, and the corrected signal is the signal after the correction of the out-of-phase.

[0022] According to a preferred embodiment, in error-correcting code sequencing or sequencing technologies such as SOLiD with precise chemical readout, the sequencing signal can be a phase-corrected signal; correspondingly, the correction process for the sequencing signal can be an error correction process, and the corrected signal is an error-corrected signal.

[0023] According to a preferred embodiment, in single-molecule real-time sequencing using a circular consensus sequencing mode, or in nanopore sequencing technology that utilizes DNA replication and repeated sequencing to improve accuracy, the sequencing signal can be the original sequence; correspondingly, the correction process for the sequencing signal can be a process of finding a consensus sequence.

[0024] According to a preferred embodiment, the method for establishing the relationship between signal differences and base accuracy is to construct a lookup table between signal differences and base accuracy.

[0025] According to a preferred embodiment, the method for establishing the relationship between signal difference and base accuracy is to divide one or more predictors into several intervals, and to statistically analyze the accuracy of bases in each interval and the quality value corresponding to the accuracy. The evaluation method may be to calculate which predictor interval each base in the tested nucleic acid falls into, and then assign the quality value corresponding to that interval to that base.

[0026] According to a preferred embodiment, the method for establishing the relationship between signal differences and base accuracy, and the corresponding evaluation method, is the Phred algorithm.

[0027] According to a preferred embodiment, the method for establishing the relationship between signal differences and base accuracy is machine learning.

[0028] According to the preferred embodiment, various differences between sequencing signal s and correction signal c can be compared to jointly evaluate base quality using multiple predictors; let s = (s1, s2, ..., sc). n c = (c1, c2, ..., c) n The various differences between the sequencing signal s and the correction signal c include:

[0029] 1)

[0030] 2)

[0031] 3)

[0032] 4) max|s i -c i |

[0033] 5) min|s i -c i|

[0034] 6)

[0035] 7)

[0036] According to preferred embodiments, the difference between the sequencing signal s and the correction signal c can be the difference of a portion of the two signals, including:

[0037] 1) the difference of all sub-signals;

[0038] 2) the difference of the first several sub-signals;

[0039] 3) the difference of the last several sub-signals;

[0040] 4) the difference of the middle several sub-signals;

[0041] 5) the difference of the odd-numbered sub-signals;

[0042] 6) the difference of the even-numbered sub-signals;

[0043] 7) the difference of the sub-signals greater than a given threshold in s or c;

[0044] 8) the difference of the sub-signals less than a given threshold in s or c;

[0045] 9) a combination of the above options, including but not limited to the difference of the first several odd-numbered sub-signals.

[0046] According to preferred embodiments, the local difference between the sequencing signal s and the correction signal c can be compared; the local difference refers to the difference between the sub-signals (s i and the correction signal c i before and after the sub-signals (s i-m , s i-m+1 , s i-m+2 , …, s i+m-1 , s i+m ) and (c i-m , c i-m+1 , c i-m+2 , …, c i+m-1 , c i+m ), i.e. a local difference can be calculated for each sub-signal to obtain a set of local differences; m is an integer less than i.

[0047] According to preferred embodiments, on the basis of the difference between the sequencing signal s and the correction signal c, other predictors can be used to jointly evaluate the base quality, including but not limited to:

[0048] 1) the position of the base in the sequence;

[0049] 2) the length of the homopolymer in which the base is located;

[0050] 3) the position of the base in the homopolymer in which it is located;

[0051] 4) the degree to which the base is located at an integer number of homopolymers;

[0052] 5) the degree of phase shift of the base in the sequencing signal;

[0053] 6) the degree of attenuation of the base in the sequencing signal;

[0054] 7) the parameters estimated in the correction process for the base in the sequencing signal, including unit signal, background signal, attenuation coefficient, advance coefficient, lag coefficient.

[0055] According to a preferred embodiment, in error-correcting code sequencing, the length of the homopolymer in which the base is located is used as a predictor in conjunction with the difference between the sequencing signal s and the corrected signal c to jointly evaluate the base quality.

[0056] According to a preferred embodiment, in error-correcting code sequencing, the number of the more abundant base in the homopolymer in which the base is located is used as a predictor in conjunction with the difference between the sequencing signal s and the corrected signal c to jointly evaluate the base quality.

[0057] The present application provides a nucleic acid sequencing data quality evaluation system, characterized in that it comprises:

[0058] a processor, a memory, and a program for quality evaluation of nucleic acid sequencing data, the program comprising the following instructions:

[0059] a) signal correction of the sequencing signal s obtained by sequencing of a reference nucleic acid sample to obtain a corrected signal c, which can be directly converted into a base sequence;

[0060] b) alignment of the converted base sequence to a reference sequence to obtain an alignment result, and then labeling the bases as correct or incorrect according to the alignment result;

[0061] c) comparing the correct or incorrect bases in the corresponding parts of the sequencing signal s and the corrected signal c, and establishing a relationship between the signal difference and the base accuracy;

[0062] d) signal correction of a new set of sequencing signals s' obtained by sequencing of a nucleic acid sample to be tested using the same correction method as step a) to obtain a corrected signal c', and then converting c' into a base sequence;

[0063] e) for each base in the converted base sequence, comparing the signal difference of the corresponding part in s' and c', using the relationship between the signal difference and the base accuracy established in step c), predicting the sequencing accuracy of the base.

[0064] In another aspect, the present application provides a screening method of nucleic acid sequencing data, characterized in that, comprising:

[0065] a) sequencing the nucleic acid sample to be tested, obtaining a set of sequencing signals;

[0066] b) correcting the sequencing signals to obtain corrected signals, which can be directly converted into base sequences;

[0067] c) comparing the signal difference of the corresponding part in the sequencing signals and the corrected signals for the converted base sequence;

[0068] d) if the signal difference between the sequencing signals and the corrected signals is greater than a given threshold, discarding the sequencing signals, otherwise retaining them.

[0069] According to a preferred embodiment, the type of signal generated by sequencing can be optical signal or electrical signal; the optical signal can be monochromatic or polychromatic.

[0070] According to a preferred embodiment, the nucleic acid to be tested includes DNA, or RNA or other types of nucleic acid molecules.

[0071] According to a preferred embodiment, the sequencing signals are the signals directly collected by the sequencer, or the normalized signals; correspondingly, the correction process of the sequencing signals can be the process of dephasing correction, and the corrected signals are the signals after dephasing correction.

[0072] According to a preferred embodiment, in the sequencing technology of error correction code or SOLiD using precise chemical reading, the sequencing signals can be the signals after dephasing correction; correspondingly, the correction process of the sequencing signals can be the process of error correction, and the corrected signals are the signals after error correction.

[0073] According to a preferred embodiment, in the single molecule real-time sequencing using circular consistency sequencing mode, or the nanopore sequencing technology using DNA replication and repeated sequencing to improve accuracy, the sequencing signals can be the original sequences; correspondingly, the correction process of the sequencing signals can be the process of seeking consistent sequences.

[0074] According to a preferred embodiment, various differences between the sequencing signals and the corrected signals can be compared; let the sequencing signals = (s1, s2,..., sn), and the corrected signals = (c1, c2,..., cn). n ), the corrected signals = (c1, c2,..., cn).n ), then the plurality of different differences between the sequencing signal and the correction signal includes:

[0075] 1)

[0076] 2)

[0077] 3)

[0078] 4) max |s i -c i |

[0079] 5) min |s i -c i |

[0080] 6)

[0081] 7)

[0082] According to a preferred embodiment, the difference between the sequencing signal and the correction signal can be the difference of a portion of the two signals, including:

[0083] 1) the difference of all sub-signals;

[0084] 2) the difference of the first several sub-signals;

[0085] 3) the difference of the last several sub-signals;

[0086] 4) the difference of the middle several sub-signals;

[0087] 5) the difference of the odd-numbered sub-signals;

[0088] 6) the difference of the even-numbered sub-signals;

[0089] 7) the difference of the sub-signals greater than a given threshold in s or c;

[0090] 8) the difference of the sub-signals less than a given threshold in s or c;

[0091] 9) a combination of the above options, including but not limited to the difference of the first several odd-numbered sub-signals.

[0092] According to a preferred embodiment, the local difference between the sequencing signal and the correction signal can be compared; the local difference refers to the sequencing signal s i and the correction signal c i , and the sub-signals before and after the sub-signals (s i-m , s i-m+1 , s i-m+2 , …, s i+m-1, s i+m ) and (c i-m , c i-m+1 , c i-m+2 ,..., c i+m-1 , c i+m ), that is, a local difference can be calculated for each sub-signal to obtain a set of local differences; m is an integer less than i; if a local difference in the set of local differences is greater than a given threshold value, the part with the excessively large local difference is cut off in the sequence, and the truncated sequence is output; if the length of the truncated sequence is less than a preset value, the whole sequence can be discarded.

[0093] According to the preferred embodiment, various differences between the sequencing signal and the correction signal can be compared, and a comprehensive judgment is made to determine whether to discard or retain the sequence to be screened.

[0094] The application also provides a screening system for nucleic acid sequencing data, characterized in that it comprises:

[0095] a processor, a memory, and a program for screening nucleic acid sequencing data, the program comprising the following instructions:

[0096] a) performing signal correction on the sequencing signal obtained by sequencing the nucleic acid sample to be tested to obtain a correction signal, which can be directly converted into a base sequence;

[0097] b) comparing the signal difference between the converted base sequence and the corresponding part of the sequencing signal and the correction signal;

[0098] c) if the signal difference between the sequencing signal and the correction signal is greater than a given threshold value, discarding the sequencing signal, otherwise retaining it.

[0099] Advantages of the application

[0100] The quality evaluation method and screening method for sequencing data disclosed by the application take the signal difference between the sequencing signal and the correction signal as an important judgment standard, and have the following advantages:

[0101] 1. The application selects the difference between the sequencing signal and the correction signal as a predictor, which contains the largest amount of information representing base accuracy and is suitable for error-correcting code sequencing.

[0102] 2. Compared with the existing method, the method of the application has a high degree of discrimination, and the highest quality value evaluated can reach 45, which exceeds 41 of Illumina and 40 of Roche 454.

[0103] 3. The sequencing data screening method of the present application does not need to obtain the quality value of each base in advance, and only needs to use the difference between the sequencing signal and the correction signal to screen the sequencing data, thereby saving the sequencing time.

[0104] 4. The sequencing data screening method of the present application can use the difference between the out-of-phase corrected signal and the error correction corrected signal to evaluate or screen the sequencing data, which is the first sequencing data evaluation and screening method specially suitable for error correction code sequencing. The present method has a wide range of applications and is suitable for sequencing data obtained by various sequencing methods including error correction code sequencing. BRIEF DESCRIPTION OF DRAWINGS

[0105] In order to more clearly illustrate the embodiments of the present application, the drawings required to be used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments consistent with the present application, and those skilled in the art can obtain the drawings corresponding to other embodiments without creative labor on the basis of these drawings.

[0106] Figure 1 . The difference between the out-of-phase correction signal and the error correction correction signal and the corresponding relationship diagram of the base quality value.

[0107] Figure 2 . The evaluation result distinction degree diagram of the quality evaluation table for the error correction code sequencing data of the nucleic acid to be tested.

[0108] Figure 3 . The evaluation diagram of the accuracy of the scoring result of the evaluation method of the present application.

[0109] Figure 4 . The signal difference density distribution diagram of the sequencing signal and the correction signal. DETAILED DESCRIPTION

[0110] In high-throughput sequencing technology, in order to amplify the sequencing signal, the DNA molecules to be detected are generally amplified before sequencing to form DNA clusters. Ideally, all molecules in the clusters formed by amplification of the same DNA molecule should be synthesized into new strands of the same length in each round of sequencing, so that the sequencing signal corresponds to the synthesized new strand one by one, and sequence information can be obtained. However, in actual sequencing, an inevitable limiting factor is dephasing, that is, the new strands of DNA molecules of the same sequence lose synchronism in the synthesis process. This is mainly caused by the incorporation of unexpected nucleotides by polymerase or incomplete extension. Dephasing includes two types: due to the presence of mismatched or impurity nucleotides, the extension length of a part of molecules in the DNA cluster is higher than the expected value, causing "lead"; a part of molecules are not extended by polymerase due to incomplete reaction, and the extension length is lower than the expected value, causing "lag". Lead and lag make the molecules in the DNA cluster have different extension lengths, and this difference will increase with the progress of the sequencing reaction, making the relationship between the sequencing signal and the actual sequence weaker and weaker, causing sequencing errors and severely limiting the increase of read length. In high-throughput sequencing technology based on amplification, signal attenuation and dephasing are the main reasons for limiting the sequencing read length. In addition, defects in optical imaging, signal acquisition, DNA amplification, mechanical systems, algorithm processing and the interference of noise can all lead to final sequencing errors, or in other words, low-quality sequencing data. Before using the sequencing data for subsequent analysis, quality assessment and screening of the sequencing data need to be completed first, and sequences with high sequencing accuracy are screened out for subsequent analysis and calculation.

[0111] In the method for evaluating the quality of the sequencing data, the difference between the sequencing signal and the correction signal is used as a predictor to evaluate the quality of the base sequencing, and the method has high distinguishability.

[0112] Unless otherwise defined, all technical terms used herein have the same meaning as commonly understood by one of ordinary skill in the art. In order to better disclose the method and content of the present application, the key terms in the present application are explained in detail.

[0113] 1. Terminology

[0114] Error Correcting Code (ECC) sequencing and error correction

[0115] The sequencing method is that the sequencing reaction of a nucleic acid sequence to be detected is composed of more than two rounds (for example, three rounds) of double-base injection flow, the three rounds of flow are orthogonal to each other, each round of flow provides half of the information amount of the DNA to be detected, and the signals of the three rounds of flow are integrated, the information redundancy is used, the sequence of the DNA to be detected can be inferred, and the sequencing error that may exist can be detected and corrected, so that the gene sequencing accuracy is greatly improved. The correction process is the error correction.

[0116] Phase correction

[0117] In high-throughput sequencing, thousands to tens of thousands of DNA molecules are contained in each sequencing unit, which are copied from the same DNA template and have the same sequence. The sequencing signal is amplified to a level that can be accurately detected. However, some DNA molecules cannot ensure full reaction in each round during the sequencing process, resulting in lag phenomenon, and some will react in advance due to substrate doping or mismatch, resulting in lead phenomenon. The extension of these DNA molecules will gradually become unsynchronized as the sequencing proceeds, and this phenomenon is called "out-of-phase" in sequencing. The process of fitting and correcting the out-of-phase sequencing signal according to the model of the sequencing reaction is called "out-of-phase correction" in algorithm.

[0118] Precise chemical readout

[0119] The exact call chemistry of the SOLiD next-generation sequencing platform uses connection sequencing to obtain a SOLiD color-coded sequence based on the "double-base coding principle". Subsequent data analysis compares the original color sequence with the reference sequence converted into a color-coded sequence, locates the SOLiD color sequence on the reference, and corrects the sequencing error.

[0120] Cyclic uniformity sequencing

[0121] The circular consensus sequencing (CCS) mode of PacBio company is to repeatedly sequence the same molecule, calculate the consensus sequence, and then improve the accuracy. The circular consensus sequencing (CCS) read is generated by aligning subreads from a single ZMW. Note that the generated CCS read does not include or need to be aligned with a reference sequence. The generated CCS reads use the CCS algorithm to require at least two rounds of reading subreads from the insert.

[0122] Nanopore sequencing technology to improve accuracy

[0123] Technologies represented by intramolecular-ligated nanopore consensus sequencing (INC-Seq) first require circularizing template DNA molecules, and then use rolling circle amplification to amplify the molecules to generate circular molecules composed of multiple repeat units, which are then broken into linear DNA chains before sequencing. Sequencing DNA reads from the same original sequence multiple times can obtain the accurate read of the original DNA sequence by combining these information through calculation. In addition, nanopore sequencing technologies that use DNA replication and repeated sequencing to improve accuracy include Circ-seq, R2C2 and other technologies.

[0124] Alignment

[0125] Align or alignment: "align" is a common concept in bioinformatics, and in bioinformatics, alignment is often used to compare the similarity between different nucleic acids or between different proteins. The alignment in the present application refers to comparing the sequencing base sequence with the reference sequence to determine whether the sequencing base sequence is correct. Commonly used sequence alignment algorithms and software include, but are not limited to, for example, Smith-Waterman algorithm, Bowtie, BWA, SOAP, Needleman-Wunch algorithm, Bowtie2, BLAST, ELAND, TMAP, MAQ, minimap2, SHRiMP, etc.

[0126] Alignment quality

[0127] The sequencing base sequence is aligned to the reference sequence using the foregoing commonly used alignment software or algorithm. For example, bowtie2 uses a scoring mechanism, rather than a complete mismatch to measure the alignment quality, and the higher the score, the greater the probability of correct alignment. The high-quality alignment described in the present application needs to be specifically selected according to the quality value range of the alignment software or algorithm used; for example, when using BWA for sequence alignment, the high-quality aligned base sequence refers to the base sequence with an alignment quality greater than 0, or greater than or equal to 10, or greater than or equal to 20, or greater than or equal to 30, or greater than or equal to 40, or greater than or equal to 50, or greater than or equal to 60.

[0128] Normalized signal

[0129] The normalized signal, i.e. the more accurate unit signal after decay correction of the original sequencing signal, the ratio of the decay-corrected sequencing signal of each sequencing site to the unit signal of this site is the normalized signal of this site in each sequencing round.

[0130] Predictor

[0131] The predictors are the features of the base, including, for example, the difference between the sequencing signal s and the correction signal c, the position of the base in the sequence, the length of the homopolymer in which the base is located, the position of the base in the homopolymer, the number of the more abundant base in the degenerate polymer in which the base is located, the degree of proximity to an integer of the sequencing signal in which the base is located, the degree of phase loss of the sequencing signal in which the base is located, the degree of decay of the sequencing signal in which the base is located, and the parameters estimated in the correction process, such as the unit signal, the background signal, the decay coefficient, the leading coefficient, and the lagging coefficient. The above-mentioned features (predictors) associated with the base call are used to assign an error probability to the base call. In this method, a training data set is first used to match the features in the training data set with a known error rate (which is usually determined by alignment with a reference genome). Then, for a target data set, the above-mentioned method can be applied to map the features to an error probability (or a quality value, which is a logarithmic transformation of the error probability). That is, the predictors do not help to improve the quality of the base call, but can help to distinguish between low-quality and high-quality base calls. In the prior art, various methods are used to assign an error probability to the corresponding base call, such as position-specific error probability, Phred, Pyrobayes, logistic regression, SVM, vector quantization, etc.

[0132] Length of degenerate multimer

[0133] Degenerate sequencing is a multi-base sequencing, which is different from single-base sequencing in that each round of reaction only extends one nucleotide molecule, and multi-base sequencing may extend multiple nucleotides in each round of reaction. The intensity of the fluorescent signal released in the sequencing reaction is positively correlated with the number of fluorescent groups released, and under ideal conditions without decay and phase loss, the fluorescent signal of each round of reaction reflects the number of bases extended in this round, which is referred to as the length of the degenerate polymer (DPL).

[0134] Degenerate base

[0135] In the present application, the degenerate bases are represented by the letters in Table 1 according to the IUPAC symbol naming rules (nucleic acid notation), for example, the letter M represents A and / or C.

[0136] Table 1

[0137] Letter Base represented M A / C K G / T R A / G Y C / T W A / T S C / G B C / G / T D A / G / T H A / C / T V A / C / G

[0138] Base quality value

[0139] The value can be an integer mapping of the probability of base recognition error. The commonly used representation is: formula q = -10 log (1-p), wherein: p is the measurement accuracy of the base, and the higher the value indicates that the base recognition is more reliable and the possibility of base error is smaller. For example, the quality value is 20 (usually written as Q20), P = 0.99 is calculated, and the sequencing accuracy is 99%. The quality value is 30 (usually written as Q30), P = 0.999 is calculated, and the sequencing accuracy is 99.9%. At present, the highest quality value that can be obtained by Illumina sequencer is 41, and the highest quality value that can be obtained by Roche454 is 40. It can be understood that the quality value is only a representation of the base sequencing accuracy, and the representation itself is not important, and the accuracy can also be directly represented. 10 (1-p), wherein: p is the measurement accuracy of the base, and the higher the value indicates that the base recognition is more reliable and the possibility of base error is smaller. For example, the quality value is 20 (usually written as Q20), P = 0.99 is calculated, and the sequencing accuracy is 99%. The quality value is 30 (usually written as Q30), P = 0.999 is calculated, and the sequencing accuracy is 99.9%. At present, the highest quality value that can be obtained by Illumina sequencer is 41, and the highest quality value that can be obtained by Roche454 is 40. It can be understood that the quality value is only a representation of the base sequencing accuracy, and the representation itself is not important, and the accuracy can also be directly represented. 2. DETAILED DESCRIPTION

[0141] Specifically, the first aspect of the present application provides a quality evaluation method of nucleic acid sequencing data, characterized by comprising:

[0142] a) sequencing a reference nucleic acid sample to obtain a set of sequencing signals s; and sequencing a to-be-tested nucleic acid sample by using the same sequencing method to obtain a set of sequencing signals s';

[0143] b) performing signal correction on the sequencing signals s to obtain a corrected signal c, wherein the corrected signal c can be directly converted into a base sequence;

[0144] c) comparing the converted base sequence to a reference sequence to obtain a comparison result, and then marking the base as correct or incorrect according to the comparison result;

[0145] d) comparing the correct or incorrect base in the corresponding part of the sequencing signals s and the corrected sequence c, and establishing a relationship between the signal difference and the base sequencing accuracy;

[0146] e) performing signal correction on the sequencing signals s' by using the same correction method as step b) to obtain a corrected signal c', and then converting c' into a base sequence;

[0147] f) comparing the corresponding part of the signal difference in s' and c' for each base in the converted base sequence, and predicting the sequencing accuracy of the base by using the relationship between the signal difference and the base accuracy established in step d).

[0148] According to a preferred embodiment, the sequencing method includes but is not limited to: dideoxynucleotide termination method (Sanger), pyrosequencing method (454), semiconductor sequencing method (Ion torrent), reversible cycle termination method (Illumina), sequencing by synthesis method (Solexa), error-correcting code sequencing method (ECC), fuzzy sequencing method (patent CN201611088606.0), joint probe anchor ligation method (Huada), oligonucleotide ligation assay sequencing method (SOLiD), single molecule fluorescent sequencing method (Helicos), single molecule real-time sequencing method (PacBio), nanopore sequencing method (Nanopore), etc.; also includes upgrading technologies based on various sequencing methods, such as the circular consensus sequencing method (CCS) improved by PacBio on the basis of single molecule real-time sequencing method, that is, by repeatedly sequencing the same molecule, the consistent sequence can be calculated, and the accuracy is improved; intramolecular-ligated nanopore consensus sequencing (INC-Seq) based on nanopore sequencing method, etc.

[0149] According to a preferred embodiment, the nucleic acid to be tested includes DNA, or RNA or other types of nucleic acid molecules, such as Xenonucleic acids (or XNAs), peptide nucleic acid (PNA), locked ribonucleic acid (LNA), etc. The type of nucleic acid to be tested does not change the core of the present application.

[0150] In some embodiments, the reference nucleic acid sample and the nucleic acid sample to be tested are sequenced at the same time, and in optional embodiments, the reference nucleic acid sample and the nucleic acid sample to be tested can be sequenced at different times, for example, the reference nucleic acid sample can be sequenced first, and then the nucleic acid sample to be tested is sequenced using the same sequencing method; the order of sequencing is not important, the condition that must be met is that the above-mentioned nucleic acid samples are sequenced using the same sequencing method, for example, both use error-correcting code sequencing.

[0151] According to a preferred embodiment, the signal generated by the sequencing can be an optical signal or an electrical signal; the optical signal can be monochromatic or polychromatic. For example, for semiconductor sequencing method, the signal generated by sequencing is an electrical signal; for pyrosequencing method, reversible cycle termination method, sequencing by synthesis method, etc., the sequencing signal generated is a fluorescent signal, i.e. an optical signal.

[0152] In a specific embodiment, the base sequence converted from the corrected signal is aligned to a reference sequence to obtain an alignment result, and then the bases are marked as correct or incorrect according to the alignment result; preferably, the high-quality aligned base sequences are further screened from the alignment result, and then the bases in the high-quality aligned base sequences are marked as correct or incorrect, and the bases that cannot be determined (i.e. the bases that cannot be successfully aligned to the reference sequence or the bases with low alignment quality) are ignored. According to the alignment result, the bases with a "match" result are marked as "correct", and the bases with a "mismatch", "insertion" or "deletion" result are marked as "incorrect". The high-quality alignment in the present application needs to be specifically selected according to the quality value range of the used alignment software or algorithm; for example, when BWA is used for sequence alignment, the high-quality aligned base sequence refers to the base sequence with an alignment quality greater than 0, or greater than or equal to 10, or greater than or equal to 20, or greater than or equal to 30, or greater than or equal to 40, or greater than or equal to 50, or greater than or equal to 60.

[0153] In the present application, when the reference nucleic acid is DNA, the reference sequence is the genomic sequence of the species to which the reference nucleic acid belongs; when the reference nucleic acid is RNA, the reference sequence is the transcriptome sequence of the species to which the reference nucleic acid belongs. It needs to be noted that the reference nucleic acid sample must be of a species whose genomic sequence is known, and a low base mutation rate should be selected to facilitate the determination of whether the difference between the sequencing sequence and the reference sequence is caused by sequencing error or normal base mutation. For example, the standard nucleic acid sample can be selected from the nucleic acids of lambda phage, Escherichia coli, Saccharomyces cerevisiae, etc.

[0154] According to the preferred embodiment, the sequencing signal can be the signal directly collected by the sequencer or the signal after normalization; correspondingly, the correction process of the sequencing signal can be the correction of the out-of-phase, and the corrected signal is the signal after the correction of the out-of-phase; the signal correction of the sequencing signal s is the same as the correction of the sequencing signal s'.

[0155] According to the preferred embodiment, in the single molecule real-time sequencing using circular consensus sequencing, or some nanopore sequencing using DNA replication, repeated sequencing to improve accuracy (such as Circ-seq, INC-seq, R2C2, etc.), the sequencing signal can be the original sequence; correspondingly, the correction process of the sequencing signal can be the process of obtaining the consensus sequence; the signal correction of the sequencing signal s is the same as the correction of the sequencing signal s'.

[0156] According to the preferred embodiment, in the error correction code sequencing, SOLiD sequencing technology using exact call chemistry technology and other sequencing technologies, the sequencing signal can be the signal after phase correction; Correspondingly, the correction process of the sequencing signal can be the error correction process. The signal correction of the sequencing signal s is the same as the correction of the sequencing signal s'. For error correction code sequencing, the specific error correction process can be referred to in patent CN201510944878.5, and the error correction process itself is not the focus of the present application. The signal difference between the corrected signal after error correction and the phase-corrected signal before correction is the focus of the present application. The representation of the signal difference can be divided into two categories, the first category is the general difference, which can be the difference of all or part of the sub-signals, and the second category is the local difference, which will be described below. It should be understood that the representation of the signal difference is only the specific application of the mathematical formula, and does not affect the core content of the present application.

[0157] According to the preferred embodiment, a plurality of different differences between the sequencing signal s and the corrected signal c can be compared as a plurality of predictors to jointly evaluate the base quality; let s = (s1, s2,..., sn), c = (c1, c2,..., cn), then the plurality of different differences between the sequencing signal s and the corrected signal c include: n n

[0158] 1)

[0159] 2)

[0160] 3)

[0161] 4)max|s i -c i |

[0162] 5)min|s i -c i |

[0163] 6)

[0164] 7)

[0165] It can also be that a part of the sub-signals of the sequencing signal s and the corrected signal c are used to calculate the signal difference by using the above formula, and the sub-signals include the following types:

[0166] 1) the difference of all sub-signals;

[0167] 2) the difference of the first several sub-signals; ​​

[0168] 3) the difference between the last several sub-signals;

[0169] 4) the difference between the middle several sub-signals;

[0170] 5) the difference between the odd-numbered sub-signals;

[0171] 6) the difference between the even-numbered sub-signals;

[0172] 7) the difference between the sub-signals in s or c that are greater than a given threshold;

[0173] 8) the difference between the sub-signals in s or c that are less than a given threshold;

[0174] 9) a combination of the above, including but not limited to the difference between the first several odd-numbered sub-signals.

[0175] In preferred embodiments, the local difference between the sequencing signal s and the corrected signal c can be compared; the local difference refers to the difference between the sequencing signal s and the corrected signal c for the sub-signals numbered i and i-m, i.e. for each sub-signal, a local difference can be calculated, resulting in a set of local differences; m is an integer less than i. i and the corrected signal c i the difference between the sub-signals (s i-m , s i-m+1 , s i-m+2 , …, s i+m-1 , s i+m ) and (c i-m , c i-m+1 , c i-m+2 , …, c i+m-1 , c i+m ) before and after, i.e. for each sub-signal, a local difference can be calculated, resulting in a set of local differences; m is an integer less than i.

[0176] In some implementations, a variety of different differences between the sequencing signal s and the corrected signal c can be compared as multiple predictors to collectively assess base quality.

[0177] In some implementations, other reported predictors can be used to collectively assess base quality based on the difference between the sequencing signal s and the corrected signal c, including but not limited to:

[0178] 1) the position of the base in the sequence;

[0179] 2) the length of the homopolymer in which the base is located;

[0180] 3) the position of the base in the homopolymer in which it is located;

[0181] 4) the degree to which the sequencing signal in which the base is located is close to an integer;

[0182] 5) the degree of phase shift of the sequencing signal in which the base is located;

[0183] 6) the degree of attenuation of the sequencing signal at the base;

[0184] 7) the parameters estimated in the correction process for the sequencing signal at the base, such as unit signal, background signal, attenuation coefficient, pre- coefficient, post-coefficient.

[0185] In a preferred embodiment, the method of establishing the relationship between signal difference and base accuracy can be to construct a table of correspondence between signal difference and base accuracy. For example, when the signal difference is 0.05, the corresponding accuracy is 99.9%; when the signal difference is 0.1, the corresponding accuracy is 99%; and so on, to construct a table of correspondence. In the f) step of the basic steps, only the accuracy corresponding to the signal difference in the table of correspondence needs to be queried to predict the accuracy of the base being measured. In some implementations, the signal difference in the table of correspondence can be expressed in the form of an interval, for example, when the signal difference is between 0.05 and 0.08, the corresponding accuracy is 99.9%.

[0186] In the basic steps, the method of establishing the relationship between signal difference and base accuracy can be to divide one or more predictors into several intervals, and to count the accuracy of the base and the quality value corresponding to the accuracy in each interval; the evaluation method can be to calculate which interval of the predictor each base in the nucleic acid being measured falls into, and then to assign the quality value corresponding to the interval to the base. Specifically, the quality evaluation table can be queried row by row, and if the predictor value of the base to be evaluated is less than the predictor threshold value of a record in the table, the quality value of the record is assigned to the base, otherwise the next record is queried.

[0187] In the basic steps, the method of establishing the relationship between signal difference and base accuracy, and the corresponding evaluation method, can be the Phred algorithm (Ewing B, Green P. Base-calling of automated sequencer traces using Phred. II. Error probabilities. Genome Res. 8: 186-194 (1998)). This method is divided into two steps of training and evaluation, and the training step is briefly described as follows:

[0188] 1. Use a certain supervised method to mark each base as correct or incorrect.

[0189] 2. Use the sequencing signal to calculate the value of each base for a number of predictors. Let there be m predictors (p1, p2,..., pm), and the value of each predictor monotonically increases with the increase of error rate. m

[0190] ​3. Divide each predictor into several thresholds. Suppose each predictor is divided into n1, n2,..., n m interval. m

[0191] 4. Calculate the error rate of the base corresponding to each interval. Suppose the thresholds of an interval are (t1, t2,..., tn), then the base corresponding to this interval is the base whose predictor values satisfy p1≤t1, p2≤t2,..., pn≤tn. m m m

[0192] 5. Select the interval with the lowest error rate, and record the predictor thresholds, corresponding error rate and quality value of this interval into the quality evaluation table.

[0193] 6. Delete the selected interval and its corresponding base in the last step.

[0194] 7. If all intervals or all bases have been deleted, end the training and output the quality evaluation table, otherwise go back to step 4.

[0195] The evaluation steps are briefly described as follows:

[0196] 1. Calculate the value of the predictor for the base to be evaluated.

[0197] 2. Query the quality evaluation table row by row. If the predictor value of the base to be evaluated is less than the predictor threshold of a record in the table, assign the quality value of this record to the base, otherwise continue to query the next record.

[0198] 3. If the quality evaluation table is queried and the base to be evaluated is not assigned a quality value, a uniform quality value such as 0, 10 or 20 can be assigned.

[0199] In the basic steps, the method of establishing the relationship between signal difference and base accuracy can be machine learning. Specifically, the value of one or more predictors is taken as input, and whether the base is correct is taken as output, and a classifier between input and output is trained by machine learning method. The evaluation method can be to calculate the predictor value of the base in the DNA to be tested, input into the trained classifier, and obtain the accuracy of the base and its corresponding quality value according to the output of the classifier. Calculating the classification accuracy according to the output of the classifier is a common technique in machine learning. Common methods are, for example, soft-max, that is, if the outputs representing "correct" and "incorrect" in the output of the classifier are a and b respectively, then the accuracy is:

[0200] ​​​​

[0201] In some preferred embodiments, for error-correcting code sequencing, the difference between the sequencing signal s and the corrected signal c can be used in combination with the length of the degenerate polymer that the base is in as a predictor to jointly assess the base quality. For example, in the sequence ACTTGAAATC, the 5th base G is in three degenerate polymers TTG, GAAA, G, and the length of the corresponding degenerate polymers are 3, 4, 1 respectively.

[0202] In some preferred embodiments, for error-correcting code sequencing, the difference between the sequencing signal s and the corrected signal c can be used in combination with the number of the more abundant base in the degenerate polymer that the base is in as a predictor to jointly assess the base quality. For example, in the sequence ACTTGAAATC, the 5th base G is in three degenerate polymers TTG, GAAA, G, and the more abundant base in each of the three degenerate polymers is T, A, G, with the number of 2, 3, 1 respectively.

[0203] The present application also provides a nucleic acid sequencing data quality assessment system, characterized in that it comprises a processor, a memory, and a program for quality assessment of nucleic acid sequencing data, the program comprising the following instructions:

[0204] a) signal correction of the sequencing signal s obtained by sequencing of a reference nucleic acid sample to obtain a corrected signal c, which can be directly converted into a base sequence;

[0205] b) alignment of the converted base sequence to a reference sequence to obtain an alignment result, and then marking the bases as correct or incorrect according to the alignment result;

[0206] c) comparison of the signal difference between the corresponding parts of the sequencing signal s and the corrected signal c for the correct or incorrect bases, and establishment of the relationship between the signal difference and the base accuracy;

[0207] d) signal correction of a new set of sequencing signals s' obtained by sequencing of a nucleic acid sample to be tested using the same correction method as step a) to obtain a corrected signal c', and then converting c' into a base sequence;

[0208] e) comparison of the signal difference between the corresponding parts of s' and c' for each base in the converted base sequence, and prediction of the sequencing accuracy of the base using the relationship between the signal difference and the base accuracy established in step c).

[0209] Each of the features discussed in the first aspect of the detailed part of the present application is equally applicable to the nucleic acid sequencing data quality assessment system of the present application. As indicated above, all other features are not repeated here and should be considered as repeated by reference. A person of ordinary skill in the art would understand how the features identified in these implementations can be readily combined with the basic feature groups identified in other implementations.

[0210] According to the general practice in the field, after the base quality assessment is completed, the measured nucleic acid sequence can be further screened according to the base quality values. For example, if the average quality value of the bases on a sequence is lower than a given threshold, the sequence can be discarded and not output; if the average quality value of the bases at the tail of a sequence is lower than a given threshold, the tail can be cut off and only the high-quality head part can be output. There are many tools that can complete this step, such as FASTP software, FASTQC software, etc. However, this screening method must be performed after the base quality values are obtained, which will inevitably waste a lot of waiting time, and for users who need sequencing, it is most desirable to obtain high-quality sequencing results after screening in the shortest possible time. Therefore, it is necessary to develop a method for screening high-quality sequencing data without scoring base quality.

[0211] Another aspect of the present application provides a nucleic acid sequencing data screening method, characterized in that it comprises:

[0212] a) sequencing a nucleic acid sample to be measured to obtain a set of sequencing signals;

[0213] b) correcting the sequencing signals to obtain corrected signals, the corrected signals can be directly converted into base sequences;

[0214] c) comparing the signal difference of the converted base sequences in the corresponding parts of the sequencing signals and the corrected signals;

[0215] d) if the signal difference between the sequencing signals and the corrected signals is greater than a given threshold, the sequencing signals are discarded, otherwise they are retained.

[0216] According to the preferred embodiment, the method for sequencing the nucleic acid sample to be tested includes but is not limited to: dideoxynucleotide termination method (Sanger), pyrophosphate sequencing method (454), semiconductor sequencing method (Ion torrent), reversible cycle termination method (Illumina), fluorescent generation sequencing method (Sena), error correction code sequencing method (ECC), fuzzy sequencing method (patent CN201611088606.0), joint probe anchor ligation method (Huada), oligonucleotide ligation detection sequencing method (SOLiD), single molecule fluorescent sequencing method (Helicos), single molecule real-time sequencing method (PacBio), nanopore sequencing method (Nanopore) and the like. It also includes upgrading technologies based on various sequencing methods, such as the circular consensus sequencing method (CCS) improved by PacBio based on the single molecule real-time sequencing method, that is, by repeatedly sequencing the same molecule, the consistent sequence can be calculated, and the accuracy is improved; intramolecular-ligated nanopore consensus sequencing (INC-Seq) based on nanopore sequencing method, etc.

[0217] According to the preferred embodiment, the nucleic acid to be tested includes DNA, or RNA or other types of nucleic acid molecules, such as Xenonucleic acids (or XNAs), peptide nucleic acid (PNA), locked ribonucleic acid (LNA) and the like. The type of nucleic acid molecule to be tested does not change the core content of the present application.

[0218] According to the preferred embodiment, the type of signal generated by sequencing can be optical signal or electrical signal; the optical signal can be monochromatic or polychromatic. For example, for semiconductor sequencing method, the signal generated by sequencing is electrical signal; for pyrophosphate sequencing method, reversible cycle termination method, fluorescent generation sequencing method and the like, the sequencing signal generated is fluorescent signal, i.e. optical signal.

[0219] According to the preferred embodiment, the sequencing signal can be the signal directly collected by the sequencer, or the signal after normalization; correspondingly, the correction process of the sequencing signal can be the out-of-phase correction, and the corrected signal is the signal after out-of-phase correction.

[0220] According to the preferred embodiment, in the single molecule real-time sequencing with circular consensus sequencing, or some nanopore sequencing (such as Circ-seq, INC-seq, R2C2, etc.) using DNA replication and repeated sequencing to improve accuracy, the sequencing signal can be the original sequence; correspondingly, the correction process of the sequencing signal can be the process of obtaining the consensus sequence.

[0221] According to the preferred embodiment, in the error correction code sequencing, SOLiD sequencing technology using exact call chemistry technology, the sequencing signal can be the signal after phase error correction; correspondingly, the correction process of the sequencing signal can be the error correction process. For error correction code sequencing, the specific error correction process can be referred to in patent CN201510944878.5, and the error correction process itself is not the focus of the present application. The signal difference between the corrected signal after error correction and the signal before phase error correction is the focus of the present application. The representation of the signal difference can be divided into two categories, the first category is the general difference, which can be the difference of all or part of the sub-signals, and the second category is the local difference, which will be described below. It should be understood that the representation of the signal difference is only the specific use of mathematical formula, and does not affect the core content of the present application.

[0222] According to the preferred embodiment, various differences between the sequencing signal and the corrected signal can be compared, assuming that the sequencing signal = (s1, s2,..., sn), and the corrected signal = (c1, c2,..., cn), then the various differences between the sequencing signal and the corrected signal include: n n

[0223] 1)

[0224] 2)

[0225] 3)

[0226] 4)max|s i -c i |

[0227] 5)min|s i -c i |

[0228] 6)

[0229] 7)

[0230] ​​The above is to take all the sub-signals of the sequencing signal and the correction signal to calculate the difference, and how to take a part of the sub-signals to calculate the difference can be easily given by imitating the above formula.

[0231] According to a preferred embodiment, the difference between the sequencing signal and the correction signal is compared, which can be the difference of a part of the two signals, including but not limited to:

[0232] 1) the difference of all the sub-signals;

[0233] 2) the difference of the first several sub-signals;

[0234] 3) the difference of the last several sub-signals;

[0235] 4) the difference of the middle several sub-signals;

[0236] 5) the difference of the odd-numbered sub-signals;

[0237] 6) the difference of the even-numbered sub-signals;

[0238] 7) the difference of the sub-signals greater than a given threshold in the sequencing signal or the correction signal;

[0239] 8) the difference of the sub-signals less than a given threshold in the sequencing signal or the correction signal;

[0240] 9) a combination of the above options, including but not limited to the difference of the first several odd-numbered sub-signals.

[0241] According to a preferred embodiment, the local difference between the sequencing signal and the correction signal can be compared; the local difference refers to the difference between the sequencing signal s i and the correction signal c i , and the sub-signals before and after them (s i-m , s i-m+1 , s i-m+2 , …, s i+m-1 , s i+m and (c i-m , c i-m+1 , c i-m+2 , …, c i+m-1 , c i+mthe difference between the sequencing signal and the corrected signal, i.e. a local difference can be calculated for each sub-signal, obtaining a set of local differences; m is an integer smaller than i; if a local difference in the set of local differences is greater than a given threshold, the part of the sequence with the local difference that is too large is cut off, and the truncated sequence is output; the threshold is not a fixed value, and for different nucleic acid sequencing reactions, the difference distribution of the sequencing signal and the correction signal can be different, and the threshold needs to be determined according to the difference distribution obtained for each sequencing reaction, and the sequencing signal greater than the threshold is discarded. In some implementations, if the length of the truncated sequence is less than a preset value, the entire sequence can be discarded, for example, the value can be 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, and a suitable value can be selected according to the specific application scenario.

[0242] According to the preferred embodiment, various differences between the sequencing signal and the correction signal can be compared, and a comprehensive judgment is made to decide whether to discard or retain the sequence to be screened. For example, the sum of the absolute values of the difference between each sequencing signal and correction signal, and the sum of the squared values of the difference between each sequencing signal and correction signal can be calculated, and a difference distribution map is drawn according to the two differences, and high-quality signals with both small differences are retained, and low-quality signals with both large differences are discarded.

[0243] The present application also provides a nucleic acid sequencing data screening system, characterized in that it comprises:

[0244] a processor, a memory, and a program for screening nucleic acid sequencing data, the program comprising the following instructions:

[0245] a) correcting the sequencing signal obtained by sequencing the nucleic acid sample to be tested to obtain a corrected signal, which can be directly converted into a base sequence;

[0246] b) comparing the signal difference between the converted base sequence and the corresponding part of the sequencing signal and the corrected signal;

[0247] c) if the signal difference between the sequencing signal and the corrected signal is greater than a given threshold, the sequencing signal is discarded, otherwise it is retained.

[0248] Each of the features discussed in the detailed description of the nucleic acid sequencing data screening method of the present application is equally applicable to the detailed implementation of the nucleic acid sequencing data screening system. As shown above, all other features are not repeated here and should be considered as repeated by reference. A person of ordinary skill in the art will understand how the features identified in these implementations can be easily combined with the basic feature groups identified in other implementations.

[0249] Example 1

[0250] Genomic DNA of lambda phage was purchased from New England Biolabs, and was sequenced after library construction. The error correction code sequencing included MK (two sets of sequencing substrates were M(A, C), K(G, T)), RY (two sets of sequencing substrates were R(A, G), Y(C, T)), WS (two sets of sequencing substrates were W(A, T), S(G, C)) and so on. In the three rounds of sequencing, the signals (i.e. sequencing signals s) of a raw fluorescence signal after signal normalization and phase correction were respectively The signals (i.e. correction signals c) after error correction were respectively A base to be evaluated for quality was respectively located at Cycle i, Cycle j, Cycle k in the three rounds of sequencing, and the difference between the phase correction signal and the error correction signal of the base was defined as:

[0251]

[0252] The obtained DNA sequences were aligned to the reference genome using BWA-MEM software, and unaligned sequences were ignored. According to the alignment results, the bases with a matching alignment result were marked as "correct", and the bases with a mismatch, insertion or deletion alignment result were marked as "error". The difference between the phase correction signal and the error correction signal was divided into several intervals according to the numerical value, the accuracy p of the bases in each interval was calculated, and the quality value was converted according to the formula q = -10 log 10 (1-p). The corresponding relationship between the difference between the phase correction signal and the error correction signal and the quality value is shown in FIG. 1, in which the highest quality value reaches 43 points, corresponding to an accuracy of 99.995%. Figure 1

[0253] Subsequently, error correction code sequencing was performed on another genomic DNA of lambda phage. For each base, the difference between the phase correction signal and the error correction signal was calculated in the same way, and the quality value of each base was evaluated according to the corresponding relationship shown in FIG. 1. Figure 1

[0254] Example 2

[0255] Genomic DNA of lambda phage was purchased from New England Biolabs, and was sequenced after library construction. For each base, the values of the following four predictors were calculated:

[0256] 1. The position of the base in the sequence; ​​

[0257] 2. the length of the homopolymer in which the base is located;

[0258] 3. the difference between the phase correction signal and the error correction signal calculated by the method of Example 1;

[0259] 4. the Euclidean distance between the phase correction signal and the integer signal of the base;

[0260] The obtained DNA sequence was aligned to the reference genome using BWA-MEM software, and unaligned sequences were ignored. According to the alignment results, the bases with alignment results of "match" were marked as "correct", and the bases with alignment results of "mismatch", "insertion" or "deletion" were marked as "error". Using the above four predictors, a quality assessment table was constructed by the Phred algorithm. The table has a total of 447 rows, of which the first 20 rows are as follows:

[0261] Predictor 1 Predictor 2 Predictor 3 Predictor 4 Base quality value 138 4.5 0.6022 0.1764 44 60 1.5 0.7943 0.1305 45 158 1.5 0.7943 0.1513 44 180 5.5 0.6022 0.1305 44 158 4.5 0.5207 0.1596 44 180 1.5 0.7943 0.1305 44 40 5.5 0.6022 0.1678 44 180 4.5 0.6909 0.2101 43 118 1.5 0.9238 0.1678 43 79 3.5 0.5207 0.2386 43 180 1.5 0.7943 0.1513 43 158 4.5 0.7943 0.1305 43 138 1.5 0.9238 0.1513 43 118 5.5 0.4406 0.1764 43 138 3.5 0.7943 0.1966 42 158 1.5 0.7943 0.1596 43 118 3.5 0.6022 0.2386 42 180 1.5 0.7943 0.1764 42 158 3.5 0.7943 0.1678 42 158 1.5 0.9238 0.1513 42

[0262] The quality assessment table was used to evaluate the error correction code sequencing data of the genomic DNA of another λ phage. Figure 2 The proportion of bases with quality values above a certain threshold in the evaluation results was shown, in which the bases with quality values above 40 accounted for 55%, and the highest quality value reached 45. Figure 3 The accuracy of the scoring in the evaluation results was shown, in which the horizontal coordinate represents the predicted quality value according to the quality assessment table, the vertical coordinate represents the actual measured quality value after alignment, the size of the disc represents the number of bases, the dotted line is the straight line represented by the equation y=x, and the two parallel lines are the straight lines represented by the equation y=x±3, and the gray parallelogram surrounded by them is the fluctuation range allowed by the quality assessment. The results show that the obtained quality assessment table can accurately evaluate the sequencing quality of the bases.

[0263] The substitution types of the bases with quality values above 40 in the sequencing data were counted, as shown in Table 2:

[0264] Table 2

[0265]

[0266] The error rate of error correction code sequencing is as low as 10 -5 Even 10 -6 orders of magnitude, which makes it have good performance in detecting genetic variations such as single-nucleotide polymorphism (SNP).

[0267] Example 3

[0268] Genomic DNA of lambda phage was purchased from New England Biolabs, and after library construction, error correction code sequencing was performed. For each base, the values of the following 5 predictors were calculated:

[0269] 1. The position of the base in the sequence;

[0270] 2. The difference between the out-of-phase correction signal and the error correction correction signal calculated by the method of Example 1;

[0271] 3. The number of the more common base in the degenerate polymers in which the base is located in the MK round;

[0272] 4. The number of the more common base in the degenerate polymers in which the base is located in the RY round;

[0273] 5. The number of the more common base in the degenerate polymers in which the base is located in the WS round.

[0274] The resulting DNA sequence was aligned to the reference genome using BWA-MEM software, ignoring unaligned sequences. According to the alignment results, any base with a "match" alignment result was marked as "correct", and any base with a "mismatch", "insertion" or "deletion" alignment result was marked as "error". Using the above 5 predictors, a quality assessment table was constructed using the Phred algorithm. The table has 301 rows, of which the first 20 rows are as follows:

[0275]

[0276]

[0277] The performance of this quality assessment table is close to that of the quality assessment table in Example 2.

[0278] Example 4

[0279] A genomic DNA sample of lambda phage was purchased from New England Biolabs, and a sequencing library was constructed and error correction code sequencing was performed, with 50 cycles in each of the MK, RY, and WS rounds, a total of 1797082 bright spots with complete sequencing signals were collected. Among them, let the MK, RY, and WS round out-of-phase corrected signals in each bright spot sequencing signal be The ECC corrected signals are Two kinds of differences are defined:

[0280]

[0281]

[0282] The difference is calculated for each bright spot, and the distribution of the difference is shown in Figure 4 Obviously, the signals can be divided into two clusters with obvious intervals according to the two differences, and the two differences of the cluster in the lower left are smaller, which are high-quality signals and will be reserved. The two differences of the cluster in the upper right are larger, which are low-quality signals and will be deleted. Thus, the screening of high-quality sequencing signals is completed.

[0283] The above is only the preferred embodiment of the present application, and those skilled in the art can make changes and modifications to the above embodiment, therefore, the present application is not limited to the above specific embodiments, and any obvious improvement, replacement or modification made by those skilled in the art on the basis of the present application belongs to the protection scope of the present application.

Claims

1. A method of quality assessment of nucleic acid sequencing data, characterized in that, The method comprises: a) sequencing a reference nucleic acid sample to obtain a set of sequencing signals s; sequencing a test nucleic acid sample with the same sequencing method to obtain a set of sequencing signals s'; b) signal correction of the sequencing signals s to obtain corrected signals c, which can be directly converted into base sequences; c) alignment of the converted base sequences to a reference sequence to obtain alignment results, and marking the bases as correct or incorrect according to the alignment results; d) comparing the correct or incorrect bases in the corresponding parts of the sequencing signals s and the corrected sequences c, and establishing a correlation between the signal difference and the base sequencing accuracy; e) signal correction of the sequencing signals s' in the same way as step b) to obtain corrected signals c', and then converting c' into base sequences; f) comparing the signal difference of each base in the converted base sequences in the corresponding parts of s' and c', and predicting the sequencing accuracy of the base by using the correlation between the signal difference and the base accuracy established in step d).

2. The method of claim 1, wherein, The genomic sequence of the species to which the reference nucleic acid sample belongs is known; when the reference nucleic acid is DNA, the reference sequence is the genomic sequence of the species to which the reference nucleic acid belongs; when the reference nucleic acid is RNA, the reference sequence is the transcriptome sequence of the species to which the reference nucleic acid belongs.

3. The method of claim 1, wherein, The converted base sequences are aligned to the reference sequence to obtain alignment results, from which high-quality aligned base sequences are further screened, and the bases in the high-quality aligned base sequences are marked as correct or incorrect.

4. The method of claim 1, wherein, The sequencing signals are the signals directly collected by the sequencer, or the normalized signals; correspondingly, the correction process of the sequencing signals is the phase correction, and the corrected signals are the phase-corrected signals.

5. The method of claim 1, wherein, In error correction code sequencing or SOLiD sequencing technology using precise chemical reading, the sequencing signals are the phase-corrected signals; correspondingly, the correction process of the sequencing signals is the error correction process, and the corrected signals are the error-corrected signals.

6. The method of claim 1, wherein, In single molecule real-time sequencing using circular consistency sequencing mode, or nanopore sequencing technology using DNA replication and repeated multiple sequencing to improve accuracy, the sequencing signals are the original sequences; correspondingly, the correction process of the sequencing signals is the process of obtaining consistent sequences.

7. The method of claim 1, wherein, The method for establishing the correlation between the signal difference and the base accuracy is to construct a correlation table between the signal difference and the base accuracy.

8. The method of claim 1, wherein, The method for establishing the correlation between the signal difference and the base accuracy is to divide one or more predictors into several intervals, and to count the accuracy of the bases in each interval and the quality value corresponding to the accuracy; the evaluation method is to calculate in which interval of the predictor each base in the measured nucleic acid falls, and to assign the quality value corresponding to the interval to the base.

9. The method of claim 1, wherein, The method for establishing the correlation between the signal difference and the base accuracy and the corresponding evaluation method are the Phred algorithm.

10. The method of claim 1, wherein, The method for establishing the correlation between the signal difference and the base accuracy is machine learning.

11. The method according to any one of claims 7-10, characterized in that, On the basis of using the difference between the sequencing signal s and the corrected signal c, other predictors are used to jointly evaluate the base quality, which include but are not limited to: 1) the position of the base in the sequence; 2) the length of the homopolymer where the base is located; 3) the position of the base in the homopolymer where the base is located; 4) the degree of the sequencing signal where the base is located approaching an integer; 5) the degree of the out-of-phase of the sequencing signal where the base is located; 6) the degree of the attenuation of the sequencing signal where the base is located; 7) the parameters of the sequencing signal where the base is located estimated in the correction process, including the unit signal, the background signal, the attenuation coefficient, the advance coefficient, the lag coefficient.

12. The method of claim 1, wherein, In error-correcting code sequencing, on the basis of using the difference between the sequencing signal s and the corrected signal c, the length of the degenerate polymer where the base is located is combined as a predictor to jointly evaluate the base quality.

13. The method of claim 1, wherein, In error-correcting code sequencing, on the basis of using the difference between the sequencing signal s and the corrected signal c, the number of the more one base in the degenerate polymer where the base is located is combined as a predictor to jointly evaluate the base quality.

14. A system for quality assessment of nucleic acid sequencing data, characterized in that, It comprises: a processor, a memory, and a program for quality evaluation of nucleic acid sequencing data, the program comprising the following instructions: a) signal correction is performed on the sequencing signal s obtained by sequencing the reference nucleic acid sample to obtain the corrected signal c, which can be directly converted into a base sequence; b) the converted base sequence is aligned to a reference sequence to obtain an alignment result, and then the base is marked as correct or incorrect according to the alignment result; c) the correct or incorrect base is compared for signal difference in the corresponding part of the sequencing signal s and the corrected signal c, and the relationship between the signal difference and the base accuracy is established; d) a new set of sequencing signals s' obtained by sequencing the nucleic acid sample to be tested is signal corrected in the same correction manner as step a) to obtain the corrected signal c', and then c' is converted into a base sequence; e) for each base in the converted base sequence, the signal difference in the corresponding part of s' and c' is compared, and the relationship between the signal difference and the base accuracy established in step c) is used to predict the sequencing accuracy of the base.

Images