A gene sequence comparison method

Abstract

The present invention discloses a gene sequence comparison method, which can quickly and accurately compare degenerate base sequences obtained by degenerate sequencing. The sequence to be tested and its corresponding reference sequence are represented as degenerate polymer sequences by the same coding method, and then a dictionary is constructed, candidate sequences are screened, and penalty points are calculated to obtain the comparison result of the sequence to be tested. The comparison method of the present invention overcomes many problems that the existing methods are not suitable for degenerate sequencing, and greatly improves the determination efficiency of the comparison error type and the comparison accuracy.

Description

Technical Field

[0001] The invention relates to a gene sequence comparison method, belonging to the field of biological information technology. Background Art

[0002] Degenerate sequencing is a new form of sequencing-by-synthesis technology. In each odd-numbered round of sequencing, a degenerate substrate consisting of two base substrates is sampled simultaneously, and in the even-numbered rounds, a degenerate substrate consisting of another two base substrates is sampled. The two degenerate substrates are sampled alternately. Except for the first round of sample injection, each round can have 1 to more base reactions, and there are no empty rounds in sequencing (referred to as 2+2 degenerate sequencing). Unlike other sequencing methods, a 2+2 degenerate sequencing does not directly obtain an accurate sequence consisting of 4 bases, but a degenerate sequence consisting of 2 degenerate bases. The information contained in the degenerate sequence is slightly less than that of the 4-base sequence, but the sequencing speed is greatly improved.

[0003] Existing alignment methods are not suitable for 2+2 degenerate sequencing and have the following defects:

[0004] First, the existing sequence alignment methods represented by BWA are developed for 4-base sequences and cannot be directly used on degenerate sequences. If the existing alignment methods are to be used, the genome needs to be degenerate.

[0005] Secondly, the existing alignment method is optimized to adapt to Illumina's base-by-base sequencing method, and has many inadaptability to the data sequenced by degenerate polymers in 2+2 degenerate sequencing. The penalty algorithm for evaluating alignment accuracy scores bases, but fails to consider the correlation between degenerate polymer length and error rate, which makes it difficult to ensure alignment accuracy.

[0006] Finally, when using existing alignment methods to align degenerate sequences, many local error types will be misjudged. For example, the existing alignment methods always tend to judge alignment errors as base substitution errors, but in 2+2 degenerate sequencing, base substitution errors rarely occur, and most of them are insertion or deletion errors. These errors will cause many misjudgments using existing alignment methods, seriously affecting the effect of sequence alignment.

[0007] In view of the many problems existing in the existing technology, it is necessary to develop a sequence alignment method suitable for 2+2 degenerate sequencing. Summary of the invention

[0008] The present invention discloses a gene sequence comparison method, which is characterized by comprising:

[0009] A. obtaining the sequencing signal intensity result corresponding to the base of the sequence to be tested by degenerate sequencing, and expressing the intensity result as a degenerate polymer sequence as the sequence to be compared; the degenerate polymer sequence is a sequence composed of the length arrangement of degenerate polymers;

[0010] B. Determine a reference sequence, and represent the reference sequence corresponding to the sequence to be tested as an ideal degenerate polymer sequence corresponding to the sequencing method of the sequence to be tested; and then extract a subsequence of length n from the ideal degenerate polymer sequence bit by bit to construct a dictionary;

[0011] C. For each sequence to be compared, extract the first n bits as the subsequence to be compared, compare the subsequence to be compared with the dictionary traversal, and find the candidate dictionary item;

[0012] D. Find the ideal degenerate polymer sequence corresponding to the candidate dictionary item, and extract a sequence having the same length as the sequence to be compared as the candidate sequence;

[0013] E. Construct a penalty function, calculate the alignment penalty of each sequence to be compared with all its corresponding alternative sequences, and compare the calculated minimum penalty value with a preset penalty threshold. If the minimum penalty value is less than the penalty threshold, the sequence to be compared is successfully aligned at the alternative sequence corresponding to the minimum penalty value; if the minimum penalty value is greater than the penalty threshold, the alignment of the sequence to be compared fails.

[0014] According to a preferred embodiment, step B is performed first, and then step A is performed.

[0015] According to a preferred embodiment, the degenerate sequencing is a sequencing reaction with an unblocked 3' end, and the degenerate sequencing includes two different sequencing reagents: a first sequencing reagent and a second sequencing reagent; the two sequencing reagents are added cyclically; wherein the first sequencing reagent comprises two different nucleotide monomers with detectable labels; the second sequencing reagent comprises two nucleotide monomers with detectable labels, and the two nucleotide monomers are different from the nucleotide monomers present in the first sequencing reagent, and wherein the second sequencing reagent is provided after the first sequencing reagent is provided, and the signal generated by the detectable label is detected after the nucleotide monomers are incorporated into the nucleic acid sequence to be detected.

[0016] According to a preferred embodiment, in the dictionary, the length of each degenerate polymer is taken as one bit, and the degenerate polymers are divided into three categories according to their lengths: length 1, length 2, and length greater than or equal to 3.

[0017] According to a preferred implementation, n is an integer, and the value range of n is 4-10.

[0018] According to a preferred embodiment, the reference sequence is a genome, or a transcriptome, or a part of a genome, or a part of a transcriptome, and the genome includes a forward sequence and a reverse sequence.

[0019] According to a preferred embodiment, the sequence of candidate dictionary items has the same length as the subsequence to be compared, and differs in at most one digit.

[0020] According to a preferred embodiment, the penalty function is

[0021] Where b is the degenerate polymer length at that position in the reference sequence, c is the degenerate polymer length at the corresponding position in the sequence to be compared, and a has no special value requirement and defaults to the constant e.

[0022] According to a preferred embodiment, the penalty function calculates the score using each degenerate polymer length in the sequence as a unit; the error penalty for the sequence to be aligned is the average of the scores of all degenerate polymer units.

[0023] According to a preferred embodiment, the preset penalty threshold is set to 0.1 by default, and the open threshold can be adjusted.

[0024] A method for performing sequence alignment on multiple input sequences obtained by degenerate sequencing, wherein each input sequence is aligned according to the aforementioned method, and then the sequence alignment result is output.

[0025] According to a preferred embodiment, the sequence alignment result is an abc format file, and the file content includes the name of the sequence to be aligned, the penalty value, the direction of alignment to the reference genome, the chromosome name and position information, the phase information, the degenerate polymer length measurement value, the degenerate polymer length, the alignment error information of each site, etc.

[0026] Advantages of the present invention

[0027] The gene sequence comparison method disclosed in the present invention overcomes the problems of existing comparison methods that are not suitable for 2+2 degenerate sequencing, and combines sequencing principles and performance, and has the following advantages:

[0028] 1. When designing the penalty function, the correlation between the length of degenerate polymers and the error rate is fully considered, thus greatly improving the alignment accuracy;

[0029] 2. It can effectively distinguish insertion or deletion errors, and improve the efficiency of determining the error type;

[0030] 3. It is not necessary to obtain an accurate 4-base sequence. The number of reaction rounds required to obtain a degenerate sequence is half of that of base-by-base sequencing and 1 / 6 of that of Ion Torrent's 1*4 sequencing, saving a lot of sequencing time. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] In order to more clearly illustrate the embodiments of the present application, a brief description of the drawings used in the embodiments will be given below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0032] Figure 1 . Flowchart of the sequence alignment method applicable to degenerate sequencing.

[0033] Figure 2 .Schematic diagram of the output file of gene sequence alignment results. DETAILED DESCRIPTION

[0034] The gene sequence comparison method of the present invention is described in detail below with reference to the accompanying drawings and specific examples. It should be understood that these examples should not be interpreted as limiting the scope of the present teaching in any way.

[0035] The existing gene sequence alignment method is an alignment method that uses a 4-base sequence as input information, and 2+2 degenerate sequencing obtains a degenerate base sequence, not specific base sequence information. Alignment using the existing alignment method will result in a high error rate. In order to solve this problem, the present invention provides an alignment method suitable for gene sequences obtained by 2+2 degenerate sequencing, which improves the efficiency of determining the type of alignment errors and improves the accuracy.

[0036] Specifically, the present invention discloses a gene sequence comparison method, characterized by comprising:

[0037] A. obtaining the sequencing signal intensity result corresponding to the base of the sequence to be tested by degenerate sequencing, and expressing the intensity result as a degenerate polymer sequence as the sequence to be compared; the degenerate polymer sequence is a sequence composed of the length arrangement of degenerate polymers;

[0038] B. Determine a reference sequence, and represent the reference sequence corresponding to the sequence to be tested as an ideal degenerate polymer sequence corresponding to the sequencing method of the sequence to be tested; and then extract a subsequence of length n from the ideal degenerate polymer sequence bit by bit to construct a dictionary;

[0039] C. For each sequence to be compared, extract the first n bits as the subsequence to be compared, compare the subsequence to be compared with the dictionary traversal, and find the candidate dictionary item;

[0040] D. Find the ideal degenerate polymer sequence corresponding to the candidate dictionary item, and extract a sequence having the same length as the sequence to be compared as the candidate sequence;

[0041] E. Construct a penalty function, calculate the alignment penalty of each sequence to be compared with all its corresponding alternative sequences, and compare the calculated minimum penalty value with a preset penalty threshold. If the minimum penalty value is less than the penalty threshold, the sequence to be compared is successfully aligned at the alternative sequence corresponding to the minimum penalty value; if the minimum penalty value is greater than the penalty threshold, the alignment of the sequence to be compared fails.

[0042] In summary, the method of the present invention uses an alignment algorithm based on a hash table data structure, divides a reference sequence into several reference subsequences of the same length, and stores them in a hash table. When performing sequence alignment, the hash table of the reference subsequence data set is searched, and then a subsequence that can completely match or substantially match the subsequence to be aligned is found, and then the entire sequence to be aligned is aligned to each candidate site and the alignment penalty is calculated, and the penalty is compared with a threshold to obtain the alignment result of the entire sequence to be aligned.

[0043] For step A of the present method, in the process of sequencing, a common multi-base sequencing method is first used to sequence a gene sequence, and a sequencing signal corresponding to the gene sequence is obtained by reading data by a sequencer. This sequencing signal is generally converted from the intensity value or current value of the image. The process of sequencing signal processing does not belong to the specific scope of the present invention. The applicant's previous patents CN201510944878.5, or CN202010061629.2, or CN201610899880.X, etc., describe in detail the method for obtaining a more accurate sequencing signal or the processing method of the sequencing signal. If necessary, the aforementioned patents can be added to this patent in the form of citation. Here, the sequencing signal processing process is simply described: after sequencing to obtain the original sequencing data, first determine whether the data is compliant. For compliant data, extract the sequencing signal intensity result, subtract background noise and normalize it. Generally, the result is an integer or floating-point vector arranged according to the acquisition time. For example, the sequencing signal of a single base is X, then when the sequencing signal obtained is 2X, it can be known that this sequencing may extend 2 bases.

[0044] Degenerate sequencing is a highly efficient sequencing method, and its 3' end is not closed. In 2+2 degenerate sequencing, two different sequencing reagents are included: a first sequencing reagent and a second sequencing reagent; the two sequencing reagents are added cyclically; wherein the first sequencing reagent contains two different nucleotide monomers with detectable labels; the second sequencing reagent contains two nucleotide monomers with detectable labels, and the two nucleotide monomers are different from the nucleotide monomers present in the first sequencing reagent, and wherein the second sequencing reagent is provided after the first sequencing reagent is provided, and the signal generated by the detectable label is detected after the nucleotide monomer is incorporated into the nucleic acid sequence to be tested. For example, the odd-numbered rounds of sequencing reactions are passed through the degenerate substrate A / C, and the even-numbered rounds of sequencing reactions are passed through the degenerate substrate G / T. This method can ensure that base extension will occur in each round of reaction, and there will be no empty reaction, so the reaction speed is faster, and the sequence obtained by the sequencing reaction is not an accurate base sequence, but a degenerate base sequence. For example, MMKMKKKMMMKKKM..., where M represents A and / or C, and K represents G and / or T.

[0045] 2+2 degenerate sequencing is multi-base sequencing. Different from single-base sequencing, each round of reaction only extends one nucleotide molecule. In multi-base sequencing, each round of reaction may extend multiple nucleotides. The intensity of the fluorescence signal released by the sequencing reaction is positively correlated with the number of released fluorescent groups. Under ideal conditions without attenuation and dephasing, the fluorescence signal of each round of reaction reflects the number of bases extended in that round, which is called degenerate polymer length (DPL). Taking the MK process of the sequence AAGCTGTCCAGG as an example, the bases extended in each round of reaction are (AA, G, C, TGT, CCA, GG), so the DPL is (2, 1, 1, 3, 3, 2), which is a degenerate polymer sequence. A degenerate polymer sequence refers to a sequence composed of degenerate polymer length arrangements. The actual measured degenerate polymer length is often not an integer, and this value can also be called the degenerate polymer length measurement value. In the method of the present invention, the sequence to be compared obtained by 2+2 degenerate sequencing is a degenerate polymer sequence, which refers to a sequence composed of degenerate polymer length arrangements. The degenerate polymer sequence also needs to be marked in the form of a phase, such as a sequence with odd bits K and even bits M, and its phase is recorded as KM. The phase refers to the base combination order corresponding to each digit marked when the sequence is recorded as a degenerate polymer sequence. The corresponding relationship between the base sequence, the degenerate base sequence and the degenerate polymer sequence can be seen in the example shown in Table 1.

[0046] Table 1. Representation of degenerate polymer sequences

[0047] Base sequence GTGACAGCAGAGCTGCGTAATCT Degenerate base sequence KKKMMMKMMKMKMKKMKKMMKMK Degenerate polymer sequence (Phase KM)3--3--12-11112-12-2-111

[0048] Existing sequence alignment methods are mainly alignment methods with 4-base sequences as input sequences, which are not applicable to 2+2 degenerate sequencing methods. When using existing alignment methods to align degenerate sequences, there will be many misjudgments of local error types. For example, in the example shown in Table 2, the existing alignment method always tends to base substitution errors in the judgment of alignment errors, and in fact, base substitution errors rarely occur in 2+2 degenerate sequencing, and the vast majority are insertion or deletion errors, therefore, the accuracy of alignment can be seriously reduced. In order to adapt to the principle of 2+2 degenerate sequencing methods, the alignment method of the present invention does not use base sequence as the input sequence of alignment, but uses the degenerate polymer sequence rewritten by the degenerate base sequence as the sequence to be aligned. It is understandable that the sequence to be aligned described in the present invention refers to the part of the joint sequence that has been removed that is unrelated to the target sequence.

[0049] Table 2. Error judgments of degenerate polymer sequence alignments and base sequence differences

[0050]

[0051] In the process of rewriting the degenerate base sequence obtained by 2+2 degenerate sequencing into a degenerate polymer sequence, it is also necessary to pay attention to the concept of phase. Degenerate sequences are written as M / K, R / Y, W / S according to the standard degenerate nucleotide identifiers, where M is A / C, K is G / T; R is A / G, Y is C / T, W is A / T, and S is G / C; the above MK or KM, RY or YR, WS or SW, etc. are the phases of the degenerate polymer sequence. The phase is the base combination order corresponding to each digit marked when the sequence is recorded as a degenerate polymer sequence, which is determined by the type and order of the sequencing substrates added during sequencing.

[0052] According to a preferred embodiment, according to the species and sample type to which the sequencing sequence belongs, a reference sequence is determined, and the reference sequence is a genome, or a transcriptome, or a part of a genome, or a part of a transcriptome. For example, if the sequence measured is a human lesion tissue sample genome, the genome sequence of the human is selected as the reference sequence. The reference sequence can be downloaded from NCBI, Ensembl and UCSC databases. After determining the reference sequence, it is necessary to re-represent its original base sequence, and represent it as a sequence in the same format as the sequence to be compared, such as if the sequence to be compared is the MK phase, then the reference sequence is rewritten as the MK sequence, and then the above sequence is recompiled to become a degenerate polymer sequence, which is an ideal degenerate polymer sequence. The reference genome sequence includes both a forward sequence and a reverse sequence, that is, the forward sequence and the reverse sequence are represented as corresponding degenerate polymer sequences, respectively, for subsequent dictionary construction.

[0053] According to the size and complexity of the reference sequence, determine the value length n of the reference subsequence and construct a dictionary of appropriate size. n is an integer, and the value range of n is 4 to 10. For smaller reference sequences, n can take a slightly larger integer, such as n=7, to speed up the comparison; for larger and more complex reference sequences, such as the human genome, n can take a smaller integer, such as 4, to reduce the size of the dictionary. The reference subsequence is obtained by extracting a subsequence of length n bit by bit from the ideal degenerate polymer sequence. Extract the n-bit index bit by bit to obtain a one-to-one corresponding position and index; merge all the obtained position index tables, sort them from small to large according to the sum of the n-bit indexes, and merge the same indexes to obtain a degenerate dictionary.

[0054] According to a preferred embodiment, step B is performed first, and then step A. That is, the reference sequence is first determined and a dictionary is constructed based on the information of the sample to be tested, and then 2+2 degenerate sequencing and subsequent sequence alignment are performed on the nucleic acid to be tested. In fact, whether step A or step B is performed first is not critical to the method of the present invention, and the implementation order can be flexibly selected according to actual conditions.

[0055] Optionally, if multiple phases of sequencing reactions are to be performed for a certain sample to be tested, then for a certain reference sequence, n can be pre-selected, and dictionaries for all phases of this reference sequence can be constructed, that is, MK (KM), WS (SW), RY (YR) dictionaries are constructed respectively, and then the corresponding dictionary is selected for subsequent alignment according to the phase information of the sequence to be compared obtained by actual sequencing. For example, if the sequence to be compared is in the MK phase, the constructed MK dictionary is selected for subsequent alignment. The advantage of this design is quite obvious, that is, for a specific reference sequence, only one dictionary construction is required, and subsequent sequencing sequences can select the corresponding dictionary to directly perform sequence alignment without rebuilding the dictionary, saving time.

[0056] The dictionary construction method of the present invention also adapts to the characteristics of 2+2 degenerate sequencing. When constructing the dictionary, each degenerate polymer length is one bit, but each bit does not simply record the degenerate polymer length 1 to ∞, but is divided into three categories, namely 1, 2 and greater than or equal to 3; the average probability of the three categories is expected to be 0.5, 0.25, and 0.25. This dictionary construction method first greatly reduces the number of required dictionary items, and secondly greatly avoids the problems of dictionary item density bias caused by unequal probability of occurrence of degenerate polymer lengths, too low utilization of some dictionary items, and too many repetitions of some dictionary items, which reasonably reduces the space occupied by the dictionary and speeds up the overall speed of comparison.

[0057] According to a preferred embodiment, the dictionary is annotated with the position information of the reference subsequence, which includes the forward or reverse direction of the sequence on the chromosome, the sequence number of the chromosome, and the position of the site on the chromosome.

[0058] According to a preferred embodiment, for each sequence to be compared, the first n bits of the sequence are extracted as the subsequence to be compared, and the extracted subsequence to be compared is traversed and compared with the reference subsequence in the dictionary to find all candidate dictionary items; the selection criteria of the candidate dictionary items are that the degenerate polymer sequence is the same length as the subsequence to be compared, and at most one digit is different, that is, there is no vacancy. For example, when n=4, the subsequence to be compared is 1,2,3,1, and the reference subsequence is 1,2,3,2. The two have the same length, which is 4 digits, and only the 4th digit is different, which meets the selection criteria of the candidate dictionary item, and therefore, the reference subsequence is listed as the candidate dictionary item.

[0059] After all candidate dictionary items are found, the ideal degenerate polymer sequence corresponding to the candidate dictionary item is found, and a sequence with the same length as the sequence to be compared is intercepted as the candidate sequence. For example, the sequence to be compared is 1,2,3,1,1,4,2,2,1,2,1,1, when n=4, the subsequence to be compared is the first 4 digits: 1,2,3,1, a certain reference subsequence is 1,2,3,2, this reference subsequence is the candidate dictionary item, and the corresponding ideal degenerate polymer sequence is 1,2,3,2,1,1,4,2,2,1,2,1,1,3,3,1,3,5,4,1,1,1,1,1, and the first 12 digits (same length as the sequence to be compared) are selected as the candidate sequence for subsequent penalty calculation.

[0060] Construct a penalty function and calculate the alignment penalty for each sequence to be aligned and all its corresponding alternative sequences. In the alignment algorithm for degenerate polymer sequences, the penalty function is changed to calculate the score based on the length of each degenerate polymer. Referring to the sequencing principle of 2+2 degenerate sequencing, the shorter the polymer, the lower the probability of random errors in 2+2 degenerate sequencing, so the higher the penalty for errors; the longer the polymer, the higher the probability of errors, so the lower the penalty for errors. The error penalty for a sequence is the average of the scores of all degenerate polymer units. The penalty function is:

[0061]

[0062] Where b is the length of the degenerate polymer at that position in the reference sequence, and c is the length of the degenerate polymer at the corresponding position in the sequence to be compared. When c = b, score = 0. When c ≠ b, the farther c deviates from b, the higher the score. There is no special requirement for the value of a, and the default value is the constant e.

[0063] After calculating the penalty for all candidate sequences, it is necessary to determine whether the sequence to be compared is aligned to the reference sequence. This step requires setting a penalty threshold and comparing the calculated minimum penalty value with the penalty threshold. If the minimum penalty value is less than the penalty threshold, the sequence to be compared is successfully aligned at the candidate sequence corresponding to the minimum penalty value; if the minimum penalty value is greater than the penalty threshold, the sequence fails to be aligned to the reference sequence. The default value of the penalty threshold is 0.1, and like other alignment systems, the open threshold can be adjusted.

[0064] According to a preferred embodiment, multiple input sequences obtained by degenerate sequencing are aligned according to the above method, and then the sequence alignment result is output. The sequence alignment result is an abc format file, and the file content includes the name of the sequence to be aligned, the penalty value, the direction of alignment to the reference genome, the chromosome name and position information, the phase information, the degenerate polymer length measurement value, the degenerate polymer length, the alignment error information of each site, etc. Figure 2 This is an example diagram, showing a line of content in the abc format file. For the sake of convenience, it is split into 9 lines for display, where the first line is the name of the sequence to be aligned; the second line is the average penalty value of the sequence; the third line is the direction information of the sequence alignment to the reference genome, 0 means the reverse direction; the fourth line is the position information of the sequence on the chromosome, 21261 means the sequence starts at the 21261st position of the chromosome; the fifth line is the phase information, MK means the first position of the degenerate polymer sequence of the sequence is the M phase; the sixth line is the degenerate polymer length measurement value, It can be seen that most of the values ​​in this part are decimals; the seventh line is the degenerate polymer length of the corresponding position on the reference genome, which is an integer; the eighth line is the alignment error information for each site. If the Query is longer than the Ref, it is represented by uppercase letters: starting from A and the maximum is Z; if the Query is shorter than the Ref, it is represented by lowercase letters: starting from a and the maximum is z; if the two are the same, it is recorded as ".", if the Query has an insertion, it is represented by "+", and if the Query has a deletion, it is represented by "-"; the last line is the name of the chromosome or genome fragment.

[0065] Example 1

[0066] Construct a 4-bit dictionary for the reference genome lambda genome:

[0067] 1. Starting from the forward sequence of the only chromosome / segment 'NC_001416.1' of the lambda genome, the 5' to 3' sequence is:

[0068] 'GGGCGGCGACCTCGCGGGTTTTCGCTATTTATGAAAATTTTCCGGTTTAAGGCGTTTCCGTTCTTCTTCGTCATAACTTAATGTTTTTATTTAAAATACCCTCTGAAAAGAAAGGAAACGACAGGTGC…'

[0069] First, convert the above nucleotide sequence into MK degenerate base sequence:

[0070] 'KKKMKKMKMMMKMKMKKKKKKKMKMKMKKKMKKMMMMKKKKMMKKKKKMMKKMKKKKMMKKKMKKMKKMKKMMKMMMKKMMKKKKKKKMKKKMMMMKMMMMKMKKMMMMKMMMKKMMMMKMMMKKKKM…'

[0071] 2. Convert the above MK degenerate base sequence into a degenerate polymer sequence:

[0072] (Phase KM) '3 1 2 1 1 3 1 1 1 1 7 1 1 1 1 1 3 1 2 4 4 2 5 2 2 1 4 2 3 12 1 2 1 2 2 1 3 2 2 7 1 3 4 1 4 1 1 2 4 1 3 2 4 1 3 4 1…'

[0073] 3. Record the bits greater than or equal to 3 in the previous step as 3:

[0074] (Phase KM)'3 1 2 1 1 3 1 1 1 1 3 1 1 1 1 1 3 1 2 3 3 2 3 2 2 1 3 2 3 12 1 2 1 2 2 1 3 2 2 3 1 3 3 1 3 1 1 2 3 1 3 2 3 1 3 3 1…'

[0075] 4. Extract the 4-bit index bit by bit to get the one-to-one corresponding position and index:

[0076]

[0077] 5. Process the reverse complementary sequence of 'NC_001416.1' in the same manner as steps 1-4 above to obtain indexes for each position.

[0078] 6. Process other chromosomes / segments in the reference genome in the same way as steps 1-5 above to obtain the index of each position (in this example, the lambda genome contains only one segment).

[0079] 7. All the position index tables obtained by merging are sorted from small to large by the sum of the 4-bit indexes, and the same indexes are merged to obtain the MK degenerate dictionary:

[0080]

[0081] 8. Follow steps 1-7 above to do the same for the other two degenerate sequencing methods RY and WS.

[0082] Example 2

[0083] Align a degenerate sequence to the reference genome in Example 1:

[0084] 1. There is a degenerate sequence seq1 to be compared obtained by 2+2 degenerate sequencing:

[0085] RY 1,1,4,1,2,1,1,5,6,1,2,2,1,1,2,1,3,2,1,4,2,1,1

[0086] Take the reference genome refgenome1 containing the ABC dictionary. First select the corresponding dictionary according to the degenerate base types of the sequence to be compared, so here we select the RY degenerate dictionary:

[0087]

[0088]

[0089] 2. Take the first few bits of the seq1 sequence to form a subsequence according to the number of bits in the selected dictionary. Here the dictionary is 4

[0090] So take the first 4 bits:

[0091] RY 1,1,4,1

[0092] The bits greater than 3 are recorded as 3:

[0093] RY 1,1,3,1

[0094] 3. Search for the subsequence in the dictionary index. Set the alignment configuration to allow 1 bit error tolerance. The following matching dictionary entries are obtained:

[0095]

[0096] 4. Find each position in refgenome1 and extract a reference sequence of the same length as the sequence to be compared seq1 from each position

[0097]

[0098]

[0099] 5. According to the penalty function:

[0100]

[0101] Calculate the penalty for each position of seq1 aligned to the reference sequence and calculate the average, taking the sequence at position 'CHROMOSOME2', forward, 638 as an example:

[0102]

[0103]

[0104] 6. Calculate the average penalty for each candidate sequence according to step 5.

[0105]

[0106]

[0107] Screen sites with an average penalty score less than the penalty threshold of 0.1. If the average penalty scores of all sites are greater than or equal to the threshold, the alignment fails. If the average penalty score of a site is less than the threshold, the site with the lowest penalty score is output as the alignment site. In this example, there is a site "'CHROMOSOME1', forward, 4" that meets the average penalty score of less than 0.1, and the alignment is successful. The reference genome position is 'CHROMOSOME1', forward, 4.

[0108] The above is only a preferred specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed by the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A method for comparing gene sequences, It is characterized in that include: A. obtaining the sequencing signal intensity result corresponding to the base of the sequence to be tested by degenerate sequencing, and expressing the intensity result as a degenerate polymer sequence as the sequence to be compared; the degenerate polymer sequence is a sequence composed of the length arrangement of degenerate polymers; B. Determine a reference sequence, and represent the reference sequence corresponding to the sequence to be tested as an ideal degenerate polymer sequence corresponding to the sequencing method of the sequence to be tested; and then extract a subsequence of length n from the ideal degenerate polymer sequence bit by bit to construct a dictionary; C. For each sequence to be compared, extract the first n bits as the subsequence to be compared, compare the subsequence to be compared with the dictionary traversal, and find the candidate dictionary item; D. Find the ideal degenerate polymer sequence corresponding to the candidate dictionary item, and extract a sequence having the same length as the sequence to be compared as the candidate sequence; E. Construct a penalty function, calculate the alignment penalty of each sequence to be compared with all its corresponding alternative sequences, and compare the calculated minimum penalty value with a preset penalty threshold. If the minimum penalty value is less than the penalty threshold, the sequence to be compared is successfully aligned at the alternative sequence corresponding to the minimum penalty value; if the minimum penalty value is greater than the penalty threshold, the alignment of the sequence to be compared fails.

2. The method according to claim 1, It is characterized in that Perform step B first, then step A.

3. The method according to claim 1, It is characterized in that The degenerate sequencing is a sequencing reaction in which the 3' end is not blocked. In the degenerate sequencing, two different sequencing reagents are included: a first sequencing reagent and a second sequencing reagent; the two sequencing reagents are added cyclically; wherein the first sequencing reagent comprises two different nucleotide monomers with detectable labels; the second sequencing reagent comprises two nucleotide monomers with detectable labels, and the two nucleotide monomers are different from the nucleotide monomers present in the first sequencing reagent, and wherein the second sequencing reagent is provided after the first sequencing reagent is provided, and the signal generated by the detectable label is detected after the nucleotide monomers are incorporated into the nucleic acid sequence to be detected.

4. The method according to claim 1, It is characterized in that In the dictionary, the length of each degenerate polymer is taken as one bit, and the degenerate polymer is divided into three categories according to the length of the degenerate polymer: length 1, length 2, and length greater than or equal to 3.

5. The method according to claim 1, It is characterized in that The n is an integer, and the value range of n is 4 to 10.

6. The method according to claim 1, It is characterized in that The reference sequence is a genome, or a transcriptome, or a part of a genome, or a part of a transcriptome, wherein the genome includes a forward sequence and a reverse sequence.

7. The method according to claim 1, It is characterized in that The sequence of the candidate dictionary items has the same length as the subsequence to be compared, and differs in at most one digit.

8. The method according to claim 1, It is characterized in that The penalty function is Where b is the degenerate polymer length at that position in the reference sequence, c is the degenerate polymer length at the corresponding position in the sequence to be compared, and a has no special value requirement and defaults to the constant e.

9. The method according to claim 8, It is characterized in that The penalty function calculates the score using each degenerate polymer length in the sequence as a unit; the error penalty for the sequence to be compared is the average of the scores of all degenerate polymer units.

10. A method for aligning multiple input sequences obtained by degenerate sequencing, wherein each input sequence is aligned according to the method of claims 1-9, and then the sequence alignment result is output.

11. The method according to claim 10, It is characterized in that The sequence alignment result is an abc format file, and the file content includes the name of the sequence to be aligned, the penalty value, the direction of alignment to the reference genome, the chromosome name and position information, the phase information, the degenerate polymer length measurement value, the degenerate polymer length, and the alignment error information of each site.

Images