Method, system and device for carrying out SNP / InDel identification and statistics on base sequence

Abstract

The invention discloses a method, a system and a device for carrying out SNP / InDel identification and statistics on a base sequence, and belongs to the technical field of bioinformatics. According to the method, quality control analysis, data splicing and multi-sequence comparison are carried out according to original data obtained after PCR product sequencing, and then SNP / InDel recognition and classification statistics are carried out on the compared data by using a doSnpAndIndel.py code file. According to the method, all conditions of each basic group in the reference sequence after next-generation sequencing can be accurately analyzed and counted, variation conditions are classified, and variation types are judged; the specific variation condition of each sequencing sequence, the initial sequences of the same variation condition and the total number of the initial sequences can be analyzed, the quality of next-generation sequencing and the sequencing preference of the reference sequence can be judged, and the method has good practicability.

Description

Technical Field

[0001] The present invention belongs to the technical field of bioinformatics, and more specifically, relates to a method, system and device for performing SNP / InDel identification and statistics on base sequences. Background Art

[0002] Polymerase Chain Reaction (PCR) is a nucleic acid synthesis technology that uses the principle of DNA double-strand replication to specifically amplify and replicate DNA fragments in vitro. Its principle is to use the parent DNA as a template and specific primers as the extension starting point under the catalysis of DNA polymerase, through repeated "denaturation (Denaturation) - annealing (Annealing) - extension (Extension)" steps, to quickly replicate the daughter DNA complementary to the parent DNA template in vitro, so as to meet the in vitro research and application of nucleic acids. PCR technology is mainly used in gene isolation and cloning, sequence analysis, gene expression regulation, gene polymorphism research and many other aspects. Because the overall or local DNA concentration is too low, direct sampling will cause problems in sequencing, so this experiment will use PCR technology to amplify the weak DNA to increase the concentration and ensure the normal operation of subsequent experiments.

[0003] Next-generation sequencing (NGS), also known as next-generation sequencing, high-throughput sequencing (HTS) or massive parallel sequencing (MPS), has a wide range of applications in biology, medicine and scientific research. Compared with traditional sequencing technology, second-generation sequencing can simultaneously sequence millions of DNA molecules in one run, and its accuracy rate of up to 99% significantly improves the accuracy of PCR product sequencing. The present invention performs NGS on PCR products and performs quality control and sequence alignment on the result data, which can efficiently identify and analyze SNP (single nucleotide polymorphism) and InDel (insertion and deletion) variations, and is helpful for the subsequent screening of phenotypic genes and trait genes of plants and animals.

[0004] Single Nucleotide Polymorphism (SNP) refers to the diversity caused by the variation of a single base in a gene sequence. SNP is the most common type of genetic information variation, accounting for more than 90% of all known genetic variations. Insertion-deletion (InDel) refers to the insertion or deletion of a small sequence fragment at a certain site in the gene sequence. Usually, the length of these variant fragments is less than 50 base pairs (bp), which is also a common type in DNA diversity. At the same time, SNP and InDel markers have high resolution and are more abundant, covering a large range of the genome and are relatively stable genetically. In the human genome, it is estimated that there is an average of 1 SNP per 1,000 base pairs, and the total number is estimated to be more than 3 million. In animal and plant research, it is precisely because of the occurrence of polymorphism and its above-mentioned characteristics that species show different phenotypes, traits, and genes, which helps in breeding, disease resistance and disaster prevention decisions. Therefore, the study of SNP and InDel is a vital part of genetic research, phenotypic research, human gene and animal and plant genetic variation research, and disease research. Generally speaking, the research on InDels is less than that on SNPs, but this experiment completed the simultaneous identification and analysis of SNPs and InDels in the same process.

[0005] Sequence alignment is a core issue in bioinformatics. It is used to compare the similarities and differences of multiple biological sequences (such as DNA, RNA or protein). It can identify the conserved regions, functional regions and mutations in the sequences, and then reveal the genetic relationship and functional differences of organisms. MAFFT (Multiple Alignment Program for Protein and Nucleic Acid Sequences) is a widely used, efficient and accurate multiple sequence alignment tool. It combines a variety of advanced algorithms and technologies and is widely used in genome annotation, protein structure prediction, evolutionary analysis and phylogenetic research. It is an important research tool in the fields of bioinformatics, genomics and proteomics.

[0006] In recent years, with the development and cross-application of biotechnology, intelligence and information technologies and the improvement of relevant theoretical levels, there is an urgent need for a technology that can accurately analyze and classify mutations and ultimately determine the quality of second-generation sequencing and the sequencing preference of reference sequences. Summary of the invention

[0007] In view of the shortcomings of the prior art, the purpose of the present invention is to provide a method for SNP / InDel identification and statistics of base sequences, which is used to identify the specific variation of sequencing sequences. Another purpose of the present invention is to provide a SNP / InDel classification and statistics method for PCR product second-generation sequencing data, which is used to classify the specific variation of sequencing sequences and identify the quality of second-generation sequencing data and the preference of reference sequence sequencing.

[0008] In order to solve the above technical problems, the technical solution adopted by the present invention is as follows:

[0009] A method for SNP / InDel identification and statistics on base sequences, comprising: performing double-end sequence quality control processing on sequence data after PCR product sequencing to obtain data to be analyzed; splicing the data to be analyzed into sequence data with a length of nearly 300 bp; performing BLAST alignment and MAFFT multiple sequence alignment on the spliced ​​sequence data; and using the doSnpAndIndel.py tool to perform SNP / InDel identification and statistics on the sequence data after the alignment.

[0010] Before performing MAFFT multiple sequence alignment, use the INSERT.py tool to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the file to be aligned by MAFFT.

[0011] The specific operation steps of the INSERT.py tool are as follows: converting the result file result.txt obtained by BLAST alignment into the result.xlsx file; storing the content of the result.xlsx file in the cols array by column, taking the data of the 8th and 9th columns of the cols array and storing them as two arrays respectively, taking the respective modes, and obtaining the reference sequence length and starting position; finally, in the reference sequence file ref1.fa, according to the sequence starting position and sequence length, intercepting the fragment of the corresponding sequencing sequence, and inserting the fragment into the first line of the mafftin.fa file.

[0012] The operation process of the doSnpAndIndel.py tool for SNP / InDel identification is as follows: referencing the AlignIO function in the Bio library to read the content of the mafft_align.fa file and save it as a file array, obtaining the sequencing depth, and extracting the first line of the reference sequence and saving it as a string.

[0013] The statistics are as follows: referencing the AlignIO function in the Bio library to read the contents of the mafft_align.fa file and save it as a file array, obtaining the sequencing depth, and extracting the first line of the reference sequence and saving it as a string; creating a site array and a content quantity array for SNP and Deletion, and creating a site array, a quantity array, and a content array for Insertion for vertical base analysis; creating a row array for storing each row of sequence, a site array that is initially "0" to determine whether it is a unique sequence, an array of the number of identical reads, a sequence array of unique sequences, and a sorting array for horizontal sequence analysis.

[0014] The vertical base analysis is specifically as follows: looping the file array columns, if the reference base position is "-", it indicates that an Insertion occurs at the position, the position is stored in the Insertion position array, the column is traversed, the content and quantity of the Insertions in the column are counted, and the positions are stored in the Insertion content array and the quantity array respectively; if the reference base position is not "-", "-" appears when traversing the column, it indicates that a Deletion occurs at the position, the position is stored in the Deletion position array, the traversal continues, the number of "-"s is counted, and the position is stored in the Deletion content quantity array; if the reference base position is not "-", a base different from the reference base appears when traversing the column, it indicates that a SNP occurs at the position, the position is stored in the SNP position array, the traversal continues, the number of each different base is counted, and the positions are stored in the SNP content quantity array respectively; if threshold screening is required, the SNP, Deletion, and Insertion quantity arrays are traversed, and the number of mutations below the threshold is reset to zero.

[0015] The horizontal sequence analysis is specifically as follows: first store each row of sequence into a row array, then traverse the row; if the site array is "0", the row sequence is a unique sequence, and subsequent row traversal is performed; if it is the same sequence row, the site array content is changed to "-1", the quantity array is increased, and the sequence array is updated; if the site array is "-1", the row is a non-unique sequence, and traversal continues to the next row; the quantity array of the same reads is copied and sorted to obtain an array with the largest quantity to the smallest quantity.

[0016] Any of the above methods for SNP / InDel identification and statistics on base sequences comprises the following specific steps:

[0017] 1) Perform NGS sequencing on the PCR product to obtain the original sequencing data;

[0018] 2) Perform double-end sequence quality control on the original sequencing data to obtain the data to be analyzed;

[0019] 3) Splice the data to be analyzed into reads with a length of nearly 300 bp, obtain the spliced ​​amplified fragment sequence file Flash.extendedFrags.fastq, and then convert its format to flash.fa;

[0020] 4) Use the reference sequence ref1.fa to build a database;

[0021] 5) Perform BLAST comparison with the database and the flash.fa file to obtain the BLAST comparison result file result.txt;

[0022] 6) Use the INSERT.py code file to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the mafftin.fa file to obtain a new mafftin.fa file;

[0023] 7) Perform MAFFT multiple sequence alignment on the mafftin.fa file to obtain the mafft_align.fa file;

[0024] 8) Use the doSnpAndIndel.py tool to perform SNP / InDel identification and statistics on the mafft_align.fa file.

[0025] The system for implementing a method for SNP / InDel identification and statistics on base sequences includes the following modules:

[0026] Raw data processing module: performs double-end sequence quality control on raw sequencing data to obtain data to be analyzed;

[0027] Splicing module: Splice the data to be analyzed into reads with a length of nearly 300 bp, obtain the spliced ​​amplified fragment sequence file Flash.extendedFrags.fastq and convert its format to flash.fa;

[0028] Database construction module: construct the database of the reference sequence ref1.fa;

[0029] BLAST comparison module: Use the database and flash.fa file to perform BLAST comparison and obtain the BLAST comparison result file result.txt;

[0030] Insertion module: Use the INSERT.py code file to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the mafftin.fa file to obtain a new mafftin.fa file;

[0031] MAFFT multiple sequence alignment module: Perform MAFFT multiple sequence alignment on the mafftin.fa file to obtain the mafft_align.fa file;

[0032] SNP / InDel identification module, SNP / InDel statistics module: Use the doSnpAndIndel.py tool to perform SNP / InDel identification and statistics on the mafft_align.fa file.

[0033] A device for performing SNP / InDel identification and statistics on a base sequence, comprising at least one processor; and at least one memory for storing at least one program. When the at least one program is executed by the at least one processor, the at least one processor implements any one of the methods for performing SNP / InDel identification and statistics on a base sequence.

[0034] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0035] 1) The present invention performs double-end sequence quality control on the data after PCR product sequencing to obtain the data to be analyzed; splices the data to be analyzed into data with a length of nearly 300 bp; performs BLAST alignment and MAFFT multiple sequence alignment on the spliced ​​data; and uses the doSnpAndIndel.py code file to perform SNP / InDel identification and classification statistics on the data after the alignment.

[0036] 2) The INSERT.py tool of the present invention inserts the reference sequence fragment corresponding to the sequencing sequence into the first line of the file to be compared by MAFFT, and can accurately compare the sequencing sequence with the reference sequence fragment.

[0037] 3) The doSnpAndIndel.py tool for SNP / InDel identification and classification statistics of the present invention can accurately analyze and count all situations that occur to each base in the reference sequence after the second-generation sequencing, and classify the variation to determine whether the variation belongs to SNP or InDel (Insertion or Deletion); it can also analyze the specific variation of each sequencing sequence (including variation site and variation type), as well as the starting sequence and total number of the same variation, to determine the quality of the second-generation sequencing and the sequencing preference of the reference sequence. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 A flow chart of the present invention as a whole;

[0039] Figure 2 It is the system module composition diagram;

[0040] Figure 3 Flowchart for inserting reference sequence fragments into the mafftin.fa file;

[0041] Figure 4 This is the sequence alignment presentation diagram of the MAFFT alignment result file in the MEGA software;

[0042] Figure 5 Flow chart for SNP / InDel analysis;

[0043] Figure 6 It is the vertical base analysis result graph in the SNP / InDel analysis result file;

[0044] Figure 7 This is a horizontal sequence analysis result graph in the SNP / InDel analysis result file. DETAILED DESCRIPTION

[0045] In order to make the purpose, technical scheme and advantages of the present invention clearer, the present invention is further described below in conjunction with specific embodiments. Unless otherwise specified in the following embodiments, the technical means used are conventional means well known to those skilled in the art.

[0046] The method for SNP / InDel identification and statistics of base sequences of the present invention is as follows: Figure 1 As shown. First, primer design, PCR amplification and second-generation sequencing were performed on the gene sequence to be analyzed to obtain the original sequence data (rawreads); then, double-end sequence quality control (FASTP quality control) and sequence splicing (FLASH splicing) were performed on the original sequence data to obtain sequence data close to 300bp (long clean reads); then, Seqtk conversion was performed (.fastq files were converted to .fasta files), BLAST alignment was performed to obtain the starting site and sequencing length of the sequencing fragment in the reference sequence, INSERT was performed to insert the reference sequence (the sequencing fragment was inserted into the mafftin.fa file), and MAFFT multiple sequence alignment was performed to obtain the mafft_align.fa alignment result file; finally, SNP and InDel analysis was performed on the mafft_align.fa file to obtain the final .xlsx result file.

[0047] The module composition of the system used to implement the above-mentioned method for SNP / InDel identification and statistics on base sequences is shown in the figure below. Figure 2 As shown, it includes the following modules:

[0048] Raw data processing module: performs double-end sequence quality control on raw sequencing data to obtain data to be analyzed;

[0049] Splicing module: Splice the data to be analyzed into reads with a length of nearly 300 bp, obtain the spliced ​​amplified fragment sequence file Flash.extendedFrags.fastq and convert its format to flash.fa;

[0050] Database construction module: construct the database of the reference sequence ref1.fa;

[0051] BLAST comparison module: Use the database and flash.fa file to perform BLAST comparison and obtain the BLAST comparison result file result.txt;

[0052] Insertion module: Use the INSERT.py code file to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the mafftin.fa file to obtain a new mafftin.fa file;

[0053] MAFFT multiple sequence alignment module: Perform MAFFT multiple sequence alignment on the mafftin.fa file to obtain the mafft_align.fa file;

[0054] SNP / InDel identification module, SNP / InDel statistics module: Use the doSnpAndIndel.py tool to perform SNP / InDel identification and statistics on the mafft_align.fa file.

[0055] A device equipped with the above system includes at least one processor; at least one memory for storing at least one program (system); when the at least one program is executed by the at least one processor, the at least one processor implements the above method for SNP / InDel identification and statistics of base sequences.

[0056] The following embodiments are all carried out using the above method, system or device.

[0057] Example 1

[0058] 1. Design primers, perform PCR amplification and next-generation sequencing on the gene sequence to obtain raw reads. According to the wild-type gene sequence of "Nanlin 895 Poplar" (NL895) (SEQ ID NO.1), design primers at both ends of the editing site (SEQ ID NO.2), with the forward and reverse primers in the range of 20-25bp in length and the amplified product in the range of 150-300bp in length; paste the gene sequence to be cloned into the Primer 5 software, set the primer search parameters, select primers with high scores and no hairpin structure, and synthesize PCR primers with adapter primers. The sequences of primers without adapters and with adapters are shown in Table 1, where the underlined sequence represents the adapter sequence of the primer.

[0059] Table 1 Primer sequences without adapter and with adapter

[0060]

[0061] According to the specific conditions of the sequence, appropriate parameters such as annealing temperature and primer concentration (as shown in Table 2) are set, and PCR amplification is performed through the three steps of "thermal denaturation-annealing-primer extension" under the mediation of the above-mentioned PCR primers to achieve amplification of the target fragment. The PCR reaction system is shown in Table 3, and the amplification program is shown in Table 4.

[0062] Table 2 Specific parameters of primers used for PCR amplification

[0063]

[0064]

[0065] Table 3 PCR reaction system

[0066] Reagents Dosage template 2μL Forward and reverse primers 0.5 μL each dNTPs(2mM) 2μL r×Taq(1.0μ / μL) 0.5μL 10×buffer 2μL <![CDATA[ddH 2 O(double distilled water)]]> up to 25 μL

[0067] Table 4 PCR amplification program

[0068]

[0069] After PCR amplification, 2 μL of the product was taken for agarose electrophoresis detection to obtain a single target band; then, the Illumina platform was used for sequencing with a sequencing depth of 8000× to obtain the raw sequencing data (rawreads) shown in Table 5.

[0070] Table 5 Raw sequencing data obtained by second generation sequencing

[0071]

[0072]

[0073] 2. Perform quality control on raw reads to obtain clean reads, and splice the clean reads into one long read.

[0074] In the process of second-generation sequencing, affected by a series of objective factors, there may be low-quality or too many N sequences in the raw reads. In order to ensure that the input of downstream analysis is clean, reliable, high-quality, and adapter-free data, the raw reads need to be quality controlled and preprocessed before starting the formal analysis. The present invention uses FASTP software (version-0.23.4, Shenzhen Hipros Co., Ltd.) to perform double-end data quality control on raw reads.

[0075] The process is: the number of bases with quality values ​​lower than Q20 in a sequence accounts for at most 40% of the total number of bases, and a sequence can have at most 10% N. Use -i, -I, -o, -O to perform double-end sequencing on rawR1.fq and rawR2.fq files to obtain the result fpR1.fq and fpR2.fq files (the specific program language used to implement this function can be: fastp-q 20-u40-n 10-i rawR1.fq-I rawR2.fq-o fpR1.fq-OfpR2.fq).

[0076] The output results of the FASTP software are as follows: (1) the total number of reads, total number of bases, number of Q20 and Q30 bases and their proportions for each .fastq file before and after filtering; (2) the filtering results of the overall quality control, including the number of reads that passed the filter, the number of low-quality reads, the number of reads with too many Ns, the number of reads that are too short, the number of reads with the adapter removed, and the number of bases removed due to the adapter; (3) the duplication rate. The results after the FASTP software quality control are shown in Tables 6 and 7.

[0077] Table 6 Results after quality control by FASTP software

[0078]

[0079] Table 7 Screening results after FASTP software quality control

[0080] Number of sequences that passed the filter 12968 Number of low-quality sequences 142 Too many N sequences 0 Too short sequence number 0 Number of sequences to remove adapters 200 Number of bases removed due to adapters 16068 Repetition rate 74.9962％

[0081] Since the second generation sequencing can only sequence reads of 150bp at most, a double-end sequencing method is used for reads close to 300bp that need to be sequenced. Therefore, in order to conduct subsequent experiments, the clean reads after quality control need to be spliced ​​to form a read of close to 300bp in length. In this embodiment, FLASH software (version-1.2.11) is used to splice reads.

[0082] The process is: set the number of threads to 6, the output file prefix to Flash, and input fpR1.fq and fpR2.fq files for splicing (the program language to implement this function can be: flash--min-overlap 10--max-mismatch-density0.25-t 6fpR1.fq fpR2.fq--output-prefix=Flash).

[0083] The output results of the software operation are: a set of statistical results of read splicing, including the total number of read pairs, the number of spliced ​​read pairs, the number of unspliced ​​read pairs, and the splicing rate, as shown in Table 8. The software operation will obtain a Flash.extendedFrags.fastq file (amplified fragment sequence file after splicing), which is the file used in the subsequent experiments.

[0084] Table 8 Statistical results after splicing by FLASH software

[0085]

[0086] Note: r1 represents sequence 1 involved in splicing; r2 represents sequence 2 involved in splicing; r represents the spliced ​​sequence.

[0087] 3. BLAST alignment and MAFFT multiple sequence alignment

[0088] In order to understand the position, length and other information of the sequencing sequence in the reference sequence, as well as the sequence alignment results between the sequencing sequence and the reference sequence, this example performed BLAST alignment and MAFFT multiple sequence alignment on the sequencing sequence and the reference sequence respectively.

[0089] Since all subsequent operations require the operation of .fasta files, this embodiment uses a fast, lightweight processing software that can perfectly process .fasta and .fastq files, seqtk software (version-1.4-r122), to convert .fastq files into .fasta files.

[0090] The process is: use the conventional conversion parameter seq-a to convert the .fastq file into a .fasta file (the program language to implement this function can be: seqtk seq-a Flash.extendedFrags.fastq>flash.fa).

[0091] The software running results: The .fastq file was perfectly converted to the .fasta file, as shown in Table 9.

[0092] Table 9 flash.fa file converted by seqtk software

[0093]

[0094]

[0095] Since the length of the reference sequence is much longer than that of the sequencing sequence, it is necessary to know the length of the sequencing sequence and its position in the reference sequence in order to perform more accurate sequence comparison and analysis. In this embodiment, a database is constructed for the reference sequence ref1.fa (doi.org / 10.1038 / s41467-020-19559-2) and a BLAST comparison is performed between the database and the spliced ​​flash.fa to obtain the start position, end position and matching length of the sequencing sequence in the reference sequence, so that the reference sequence fragment corresponding to the sequencing sequence is cut out from the reference sequence according to the above information.

[0096] The process of building a database is as follows: use makeblastdb to build a nucleic acid library, the database type is nucl, and the database name is ref1_db (the command to build the database can be: makeblastdb-in ref1.fa-dbtype nucl-outref1_db).

[0097] The BLAST sequence alignment process used in this embodiment is: setting the e-value threshold to 1e-5, the output file format to 6, and using the ref1_db database and the flash.fa file to perform nucleotide sequence alignment (the BLAST sequence alignment command may be: blastn-query flash.fa-db ref1_db-out results1.txt-evalue 1e-5-outfmt 6).

[0098] The above operation obtains the starting position and length of the flash.fa file in the reference sequence. The file result.txt is a file that records the BLAST comparison results, in which the 7th and 8th columns respectively represent the starting and ending positions of the comparison region on the query sequence, which can indirectly derive the length of the sequencing sequence; the 9th and 10th columns respectively represent the starting and ending positions of the comparison region on the reference sequence. Columns 7 to 10 of the result.txt file are shown in Table 10.

[0099] Table 10 Columns 7 to 10 in the result.txt file

[0100]

[0101] Before performing MAFFT multiple sequence alignment, the reference sequence fragment corresponding to the sequencing sequence needs to be inserted into the first line of the mafftin.fa file to be used for MAFFT alignment. Only in this way can the sequencing sequence and the reference sequence fragment be accurately aligned to obtain the desired result.

[0102] This embodiment uses the INSERT.py tool, and its workflow diagram is as follows: Figure 3 The main process is as follows: First, convert the result file result.txt obtained by BLAST alignment into the result.xlsx file; then, save the data in the 8th and 9th columns of result.xlsx into two arrays respectively, and take their respective modes, so as to ensure the accuracy of the starting site and matching length; finally, according to the sequence starting site, length and format characteristics of the .fasta file, the real reference fragment is intercepted and stored in the mafftin.fa file. The first 3 lines of the mafftin.fa file are shown in Table 11.

[0103] Table 11 The first 3 lines of the mafftin.fa file

[0104]

[0105]

[0106] Perform MAFFT multiple sequence alignment on the mafftin.fa file to align all sequenced sequences with the inserted reference sequence fragments.

[0107] The process is: use automatic selection algorithm parameters to perform MAFFT alignment on the mafftin.fa file to obtain the mafft_align.fa result file (the program language to implement this function can be: mafft--auto mafftin.fa>mafft_align.fa).

[0108] There are no gaps in the sequences in the mafftin.fa file and the lengths are different. The MAFFT software compares each sequencing sequence with the first line of the reference sequence. For the missing bases that are not aligned, "-" is used to fill them, ultimately ensuring that all sequences are of the same length and the alignment is accurate. The MEGA software (Molecular Evolutionary Genetics Analysis, MegaLimited, version-11.0.13) displays the bases of all sequences one by one by marking different bases with different colors, making it easier for observers to observe base variations. The results of multiple sequence alignment using the MEGA software are shown below. Figure 4 shown.

[0109] 4. Identify, classify and count the SNPs and InDels in the aligned mafft_align.fa file.

[0110] According to the above steps, operations and analysis of the classification of SNPs and InDels, this embodiment makes a classification of SNPs / InDels as shown in Table 12, and uses the doSnpAndIndel.py tool to realize the identification, classification and statistics of SNPs / InDels.

[0111] The doSnpAndIndel.py tool performs SNP and InDel identification, classification, and statistics on the aligned mafft_align.fa file. The workflow is as follows: Figure 5 The specific operation steps are as follows:

[0112] 1) Read the MAFFT comparison results

[0113] By referencing the AlignIO function in the Bio library, the contents of the mafft_align.fa file are read and saved as a file array, the sequencing depth is obtained (if variations caused by objective reasons such as the machine need to be screened out, a variation threshold can be set), and the first line of the reference sequence is extracted and saved as a string.

[0114] 2) Vertical base analysis

[0115] Create the site array and content number array of SNP and Deletion, and create the site array, number array and content array of Insertion. Loop through the file array columns. If the reference base position is "-", it means that the position has been inserted. Store the site in the site array of Insertion, traverse the column, count the content and number of Insertions in the column, and store them in the content array and number array of Insertion respectively. If the reference base position is not "-", "-" appears when traversing the column, it means that Deletion has occurred at the position. Store the site in the site array of Deletion, continue traversing, count the number of "-", and store it in the content number array of Deletion; if the reference base position is not "-", a base different from the reference base appears when traversing the column, it means that SNP has occurred at the position. Store the site in the site array of SNP, continue traversing, count the number of each different base, and store them in the content number array of SNP respectively. If threshold screening is required, traverse the number arrays of SNP, Deletion, and Insertion, and reset the number of variants below the threshold to zero.

[0116] For vertical base analysis, the site information, variation, number and proportion of specific variation, number of each variation base and total statistics are recorded (if the number array values ​​of SNP, Deletion, and Insertion of the site are all "0", it means that the column has not mutated, and the information of the column is deleted), and different types of variation are distinguished by different font colors, which has better visibility and interactivity. The results of vertical base analysis are as follows: Figure 6 shown.

[0117] 3) Horizontal sequence analysis

[0118] Create a row array to store each row of sequence, a site array that is initially "0" to determine whether it is a unique sequence, an array of the number of identical reads, a sequence array of unique sequences, and a sorting array. First store each row of sequence in the row array, then traverse the row. If the site array is "0", the row sequence is a unique sequence, and traverse the subsequent rows. If it is the same sequence row, change the site array content to "-1", increase the number array, and update the sequence array; if the site array is "-1", the row is not a unique sequence, and continue to traverse the next row. Copy and sort the array of the number of identical reads to obtain an array with the largest number in descending order.

[0119] For horizontal sequence analysis, combined with the array of vertical base analysis, each unique sequence is statistically analyzed for quantity, proportion, mutation type, mutation site, specific mutation situation, and comprehensive statistics of specific sequence content. For non-unique sequences, the row of information is deleted. The results of horizontal sequence analysis are as follows: Figure 7 shown.

[0120] 4) Create an Excel table

[0121] The variation was recorded according to the column headings of the table (Table 12).

[0122] Table 12 Classification of SNP / InDel analyzed by the present invention

[0123]

Claims

1. A method for SNP / InDel identification and statistics on base sequences, characterized in that: The sequence data after PCR product sequencing was subjected to double-end sequence quality control processing to obtain the data to be analyzed; the data to be analyzed was spliced ​​into sequence data with a length of nearly 300bp; the spliced ​​sequence data was subjected to BLAST alignment and MAFFT multiple sequence alignment; the doSnpAndIndel.py tool was used to perform SNP / InDel identification and statistics on the sequence data after completion of the alignment.

2. The method for SNP / InDel identification and statistics of base sequences according to claim 1, characterized in that: Before performing MAFFT multiple sequence alignment, use the INSERT.py tool to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the file to be aligned by MAFFT.

3. The method for SNP / InDel identification and statistics of base sequences according to claim 2, characterized in that: The specific operation steps of the INSERT.py tool are as follows: converting the result file result.txt obtained by BLAST alignment into the result.xlsx file; storing the content of the result.xlsx file in the cols array by column, taking the data of the 8th and 9th columns of the cols array and storing them as two arrays respectively, taking the respective modes, and obtaining the reference sequence length and starting position; finally, in the reference sequence file ref1.fa, according to the sequence starting position and sequence length, intercepting the fragment of the corresponding sequencing sequence, and inserting the fragment into the first line of the mafftin.fa file.

4. The method for SNP / InDel identification and statistics of base sequences according to claim 1, characterized in that: The operation process of the doSnpAndIndel.py tool for SNP / InDel identification is as follows: referencing the AlignIO function in the Bio library to read the content of the mafft_align.fa file and save it as a file array, obtaining the sequencing depth, and extracting the first line of the reference sequence and saving it as a string.

5. The method for SNP / InDel identification and statistics of base sequences according to claim 1, characterized in that: The statistics are as follows: referencing the AlignIO function in the Bio library to read the contents of the mafft_align.fa file and save it as a file array, obtaining the sequencing depth, and extracting the first line of the reference sequence and saving it as a string; creating a site array and a content quantity array for SNP and Deletion, and creating a site array, a quantity array, and a content array for Insertion for vertical base analysis; creating a row array for storing each row of sequence, a site array that is initially "0" to determine whether it is a unique sequence, an array of the number of identical reads, a sequence array of unique sequences, and a sorting array for horizontal sequence analysis.

6. The method for SNP / InDel identification and statistics of base sequences according to claim 5, characterized in that: The vertical base analysis is specifically as follows: looping the file array columns, if the reference base position is "-", it means that the position has an insertion, storing the position in the insertion position array, traversing the column, counting the content and number of insertions in the column, and storing them in the insertion content array and number array respectively; If the reference base position is not "-", and "-" appears when traversing the column, it means that a Deletion has occurred at this position. The site is stored in the Deletion site array, and the traversal continues. The number of "-"s is counted and stored in the Deletion content quantity array; if the reference base position is not "-", and a base different from the reference base appears when traversing the column, it means that a SNP has occurred at this position. The site is stored in the SNP site array, and the traversal continues. The number of each different base is counted and stored in the SNP content quantity array respectively; If threshold screening is required, traverse the SNP, Deletion, and Insertion quantity arrays and reset the number of variants below the threshold to zero.

7. The method for SNP / InDel identification and statistics of base sequences according to claim 5, characterized in that: The horizontal sequence analysis is specifically as follows: first store each row of sequence into a row array, then traverse the rows; if the site array is "0", the row sequence is a unique sequence, and subsequent row traversal is performed; if it is the same sequence row, the site array content is changed to "-1", the quantity array is increased, and the sequence array is updated; if the site array is "-1", the row is a non-unique sequence, and the traversal continues to the next row; the quantity array of the same reads is copied and sorted to obtain an array with the largest quantity to the smallest quantity.

8. The method for SNP / InDel identification and statistics of base sequences according to any one of claims 1 to 7, characterized in that: The specific steps are: 1) Perform NGS sequencing on the PCR product to obtain the original sequencing data; 2) Perform double-end sequence quality control on the original sequencing data to obtain the data to be analyzed; 3) Splice the data to be analyzed into reads with a length of nearly 300 bp, obtain the spliced ​​amplified fragment sequence file Flash.extendedFrags.fastq, and then convert its format to flash.fa; 4) Use the reference sequence ref1.fa to build a database; 5) Perform BLAST comparison with the database and the flash.fa file to obtain the BLAST comparison result file result.txt; 6) Use the INSERT.py code file to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the mafftin.fa file to obtain a new mafftin.fa file; 7) Perform MAFFT multiple sequence alignment on the mafftin.fa file to obtain the mafft_align.fa file; 8) Use the doSnpAndIndel.py tool to perform SNP / InDel identification and statistics on the mafft_align.fa file.

9. A system for implementing the method for SNP / InDel identification and statistics on base sequences according to claim 1, characterized in that: Includes the following modules: Raw data processing module: performs double-end sequence quality control on raw sequencing data to obtain data to be analyzed; Splicing module: Splice the data to be analyzed into reads with a length of nearly 300 bp, obtain the spliced ​​amplified fragment sequence file Flash.extendedFrags.fastq and convert its format to flash.fa; Database construction module: construct the database of the reference sequence ref1.fa; BLAST comparison module: Use the database and flash.fa file to perform BLAST comparison and obtain the BLAST comparison result file result.txt; Insertion module: Use the INSERT.py code file to insert the reference sequence fragment corresponding to the sequencing sequence into the first line of the mafftin.fa file to obtain a new mafftin.fa file; MAFFT multiple sequence alignment module: Perform MAFFT multiple sequence alignment on the mafftin.fa file to obtain the mafft_align.fa file; SNP / InDel identification module, SNP / InDel statistics module: Use the doSnpAndIndel.py tool to perform SNP / InDel identification and statistics on the mafft_align.fa file.

10. A device for SNP / InDel identification and statistics on base sequences, characterized in that: It comprises at least one processor; at least one memory for storing at least one program; when the at least one program is executed by the at least one processor, the at least one processor implements a method for SNP / InDel identification and statistics on base sequences as described in any one of claims 1 to 7.

Images