Real-time SNV and Indel detection method for next-generation sequencing data

Abstract

The invention discloses a real-time SNV and Indel detection method for second-generation sequencing data, and relates to the real-time SNV and Indel detection method for the second-generation sequencing data. The invention aims to solve the problem that a large amount of initial data generated by a modern high-throughput sequencing platform in a sequencing fragment extension and stack information updating process is in an idle state; the problems of serious waste of computing resources and low efficiency of an analysis process are caused, so that the total period from sample sequencing to obtaining of a variation detection result is remarkably prolonged. The method comprises the following steps: acquiring and processing generated sequencing fragment data in real time while a sequencer performs sequencing operation, immediately mapping the fragments onto a reference genome by utilizing an improved sequence alignment algorithm, determining the accurate positioning of the fragments on the genome, and then continuously performing the sequencing process to obtain the accurate positioning of the fragments on the reference genome. And the obtained comparison result is continuously extended and updated, so that the parallel execution of sequencing and variation detection is finally realized, and the variation detection efficiency is remarkably improved.

Description

Technical Field

[0001] The present invention relates to a real-time SNV and Indel detection method for second-generation sequencing data. Background Art

[0002] With the rapid development and wide application of the new generation of high-throughput sequencing technology, variant detection has become an essential basic link and key step in the bioinformatics analysis process. In the fields of genomics research, clinical diagnosis, and precision medicine, various downstream analysis processes, such as genotyping, variant pathogenicity assessment, population genetics analysis, and phenotype association studies, are all premised on and supported by highly accurate variant detection results.

[0003] In genomic variant research, a variant essentially refers to the difference between the genomic sequence of an individual sample and the standard reference genome. According to the widely accepted classification criteria in the international bioinformatics field, small variants on the genome include single nucleotide polymorphism (SNP) and small insertion-deletion variants (Indels) with a length of no more than 50 base pairs (bp). These small variants account for more than 90% of all known genomic variants in humans and have irreplaceable core value and decisive significance in aspects such as the analysis of the pathogenic mechanisms of human genetic diseases, drug target identification, and the formulation of personalized medical treatment plans.

[0004] In response to the variant detection requirements for second-generation sequencing (NGS) data, research institutions and enterprises at home and abroad have developed variant detection tools with various technical routes: traditional variant detection tools represented by GATK (Genome Analysis Toolkit), which analyze and evaluate sequence alignment information through probability statistics methods; emerging deep learning-based variant detection tools represented by Clair3 and DeepVariant, which automatically extract and classify the characteristics of sequencing data through artificial intelligence algorithms such as convolutional neural networks (CNNs) and recurrent neural networks (RNNs).

[0005] For variant detection tools in the existing technical routes, whether they are traditional algorithms based on statistical models and Bayesian inference (such as GATK, SAMtools) or new algorithms based on artificial intelligence technologies such as deep neural networks (such as Clair3, DeepVariant), there is a common defect in the timing architecture: these tools must wait for the sequencer to complete the data generation process of the entire sample and output all sequencing reads before they can start the sequence alignment and subsequent variant detection processes.

[0006] Since modern high-throughput sequencing platforms (such as Illumina NovaSeq, PacBio Sequel, Oxford Nanopore PromethION, etc.) usually take several hours to several days to complete a full sequencing experiment, a large amount of initial data generated during the sequencing fragment extension and stack information update process is left idle and cannot be utilized immediately. This will cause serious waste of computing resources and low efficiency of the overall analysis process, thus significantly extending the overall cycle from sample sequencing to obtaining variant detection results, posing challenges to the timeliness of clinical diagnosis, disease risk assessment, and scientific research projects. Summary of the Invention

[0007] The object of the present invention is to solve the problem that modern high-throughput sequencing platforms usually take several hours to several days to complete a full sequencing experiment, resulting in a large amount of initial data generated during the sequencing fragment extension and stack information update process being left idle and unable to be utilized immediately; causing serious waste of computing resources and low efficiency of the overall analysis process, thus significantly extending the overall cycle from sample sequencing to obtaining variant detection results, and to propose a real-time SNV and Indel detection method for second-generation sequencing data.

[0008] The specific process of a real-time SNV and Indel detection method for second-generation sequencing data is as follows:

[0009] Step 1: Obtain the alignment results of the sequencing reads on the human reference genome;

[0010] Based on the alignment results, obtain the final candidate variant sites;

[0011] Based on the alignment results, extract the unique identifier read_id in the original information of each sequencing read as the index key and store it in the read_info dictionary; write the read_info dictionary into a pickle file through serialization technology; the specific process is as follows:

[0012] Step 11: Obtain the alignment results of the sequencing reads on the human reference genome, divide each chromosome in the aligned human reference genome into non-overlapping blocks, and extract the stack information at each position on the genome within the blocks;

[0013] Step 1-2: Use the mpileup function of samtools to parse the pileup information at each position on the genome within the block extracted in Step 1-1, and obtain the parsing result, which is used as the initial candidate variant sites.

[0014] Step 1-3: Screen the initial candidate variant sites obtained in Step 1-2 to obtain the final candidate variant sites.

[0015] Step 1-4: Call the samtools tool to traverse all the sequencing reads on the chromosome aligned in Step 1-1, and obtain the original information of each sequencing read.

[0016] Extract the unique identifier read_id in the original information of each sequencing read as the index key and store it in the read_info dictionary.

[0017] In the read_info dictionary data structure, store the feature information for each sequencing read.

[0018] Write the read_info dictionary into a pickle file through serialization technology.

[0019] Step 2: After the sequencer outputs a new batch of sequencing reads, perform selective filtering on the new batch of sequencing data to obtain the filtered reads; construct event tuples.

[0020] Traverse each event tuple (pos, event_type, idx), and execute different processing logics according to the event type to obtain the information of each read in the sorted list sorted_events under different processing logics.

[0021] Analyze the differences between each read and the reference sequence site by site on the genome, and update the pileup information at the corresponding sites.

[0022] The specific process is as follows:

[0023] Step 2-1: Load from disk the final candidate variant site information stored in the candidate_positions dictionary saved in Step 1-3 and the sequencing read information stored in the read_info dictionary recorded in Step 1-4.

[0024] After the sequencer outputs a new batch of sequencing reads, perform selective filtering on the new batch of sequencing data, and only retain the reads that have a corresponding relationship with the loaded read_info dictionary.

[0025] The retained new batch of sequencing reads is stored in FASTQ format.

[0026] Step 22: Construct event tuples based on Step 21;

[0027] Step 23: Construct a sorted_events sorted list based on the event tuples constructed in Step 22. Traverse each event tuple (pos, event_type, idx) in the sorted list sorted_events, and execute different processing logics according to the event type to obtain the information of each read segment in the sorted list sorted_events under different processing logics;

[0028] Step 24: Analyze the differences between each read segment in the reads_to_realign set in Step 23 and the reference sequence site by site on the genome, and update the stacking information of the corresponding sites;

[0029] Step 3: Perform variant detection based on Step 2 to obtain a variant detection result, and pass the variant detection result into a VCF file.

[0030] The beneficial effects of the present invention are:

[0031] The present invention aims to propose an innovative real-time variant detection technology architecture. By obtaining and processing the generated sequencing fragment data in real time while the sequencer is performing sequencing operations, using an improved sequence alignment algorithm to immediately map these fragments to the reference genome to determine their precise positions on the genome, and then continuously extending and updating the obtained alignment results as the sequencing process continues, ultimately realizing the parallel execution of sequencing and variant detection, and significantly improving the efficiency of variant detection. Brief Description of the Drawings

[0032] Figure 1 It is a flowchart of the present invention;

[0033] Figure 2 It is a schematic diagram of stack information processing based on block parallelism. Detailed Embodiments

[0034] Detailed Embodiment 1: The specific process of a real-time SNV and Indel detection method for second-generation sequencing data in this embodiment is as follows:

[0035] Step 1: Obtain the alignment result of the sequencing reads on the human reference genome;

[0036] Based on the alignment result, obtain the final candidate variant sites;

[0037] Based on the alignment results, the unique identifier read_id in the original information of each sequencing read is extracted as the index key and stored in the read_info dictionary; the read_info dictionary is written into a pickle file through serialization technology to achieve efficient data storage and fast loading; the specific process is as follows:

[0038] Step 1-1: Obtain the alignment results of the sequencing reads on the human reference genome. Each chromosome in the aligned human reference genome is divided into non-overlapping blocks, and the pileup information at each position on the genome within the block is extracted;

[0039] Step 1-2: Use the mpileup function of samtools to parse the pileup information at each position on the genome within the block extracted in Step 1-1 to obtain the parsing results (perform in-depth analysis on each genomic position within the interval block to identify potential candidate variants (candidate variants) for each sample). The parsing results serve as the initial candidate variant sites;

[0040] Step 1-3: Screen the initial candidate variant sites obtained in Step 1-2 to obtain the final candidate variant sites;

[0041] Step 1-4: Invoke the samtools tool to traverse all the sequencing reads on the chromosome aligned in Step 1-1 to obtain the original information of each sequencing read (recorded in the BAM file);

[0042] The unique identifier read_id in the original information of each sequencing read is extracted as the index key and stored in the read_info dictionary;

[0043] In the read_info dictionary data structure, multi-dimensional feature information is stored for each sequencing read;

[0044] The read_info dictionary is written into a pickle file through serialization technology to achieve efficient data storage and fast loading;

[0045] Step 2: After the sequencer outputs a new batch of sequencing reads (multiple sequencing reads) (50bp or other values set by the user), perform selective filtering on the new batch of sequencing data to obtain the filtered reads; construct event tuples;

[0046] Traverse each event tuple (pos, event_type, idx), and execute different processing logics according to the event type to obtain the information of each read in the sorted list sorted_events under different processing logics;

[0047] Analyze the differences between each read and the reference sequence site by site on the genome, and update the stacking information at the corresponding sites;

[0048] The specific process is as follows:

[0049] Step 2-1: Load from disk the final candidate variant site information stored in the candidate_positions dictionary in Step 1-3 and the sequencing read information stored in the read_info dictionary recorded in Step 1-4;

[0050] After the sequencer outputs a new batch of sequencing reads (multiple sequencing reads) (50bp or other values set by the user), perform selective filtering on the new batch of sequencing data, and only retain the reads that have a corresponding relationship with the loaded read_info dictionary;

[0051] The retained new batch of sequencing reads is stored in FASTQ format;

[0052] The subsequent steps to construct events are based on this loaded result, and the subsequent results need to be updated in candidate_positions and read_info;

[0053] Step 2-2: Construct event tuples based on Step 2-1;

[0054] Step 2-3: Construct a sorted_events sorted list based on the event tuples constructed in Step 2-2. Traverse each event tuple (pos, event_type, idx) in the sorted_events sorted list, and execute different processing logics according to the event type to obtain the information of each read in the sorted_events sorted list under different processing logics;

[0055] Step 2-4: Analyze the differences between each read in the reads_to_realign set in Step 2-3 and the reference sequence site by site on the genome, and update the stacking information at the corresponding sites;

[0056] Step 3: Perform variant detection based on Step 2 to obtain the variant detection result (variant detection quality score), and pass the variant detection result into the VCF file.

[0057] Specific Embodiment 2: The difference between this embodiment and Specific Embodiment 1 is that in Step 1-1, the alignment result of the sequencing reads on the human reference genome is obtained, each chromosome in the aligned human reference genome is divided into non-overlapping blocks, and the stacking information at each position on the genome within the block is extracted;

[0058] The specific process is as follows:

[0059] 1) Use a sequencer to obtain human sequencing reads in real time, take a 150 bp sequencing read as the initial data block, and store it in FASTQ format;

[0060] 2) Use the sequence alignment tool minimap2 to map the sequencing reads on the initial data block to the human reference genome to obtain the original alignment results in SAM format (including aligned and unaligned);

[0061] 3) Call the bioinformatics tool samtools to convert the original alignment results in SAM format into the original alignment results in BAM format (a more efficient binary format);

[0062] 4) Sort the aligned sequencing reads in the original BAM format according to their positions on the human reference genome from front to back and create an index; establish an efficient retrieval structure to lay the foundation for subsequent data processing;

[0063] 5) Based on the position and index of the sequencing reads obtained in 4) on the human reference genome, each chromosome in the aligned human reference genome is divided into M non-overlapping blocks according to a pre-set interval length (for example, there are 22 chromosomes, and the length of each chromosome or the number of bases contained in the chromosome varies, and the number of blocks varies. For example, if the length of chromosome 1 is 248956422 and the set block size is 20000000, a total of 13 blocks are obtained (248956422 / 20000000~=13, rounded up)). Each block is uniquely identified by a triplet of "chromosome:start position:end position"; the value of M is a positive integer;

[0064] 6) Call GNUparallel to parallelize the calculation of different blocks of the same chromosome;

[0065] The mpileup function of samtools is used to extract the stacking information of each position on the genome within the block. The stacking information contains rich sequence features, mainly including the reference sequence name of the current position on the human reference genome, the genomic coordinates of the current position on the human reference genome, the reference base of the current position on the human reference genome, the number of sequencing fragments covering the current position (such as read 1 on the alignment, read 2 on the alignment, read 3 on the alignment, read 4 on the alignment), the bases of the sequencing fragments covering the current position, and the base quality ASCII code sequence of the sequencing fragments covering the current position. The stacking information is all information on the aligned human reference genome;

[0066] SAMTools's MPileup performs quality control filtering based on the preset alignment quality score threshold to eliminate low-quality sequencing fragments and ensure data reliability for subsequent analysis.

[0067] The other steps and parameters are the same as those in the first specific implementation manner.

[0068] Specific implementation manner three: The difference between this implementation manner and the first or second specific implementation manner is that in steps one and two, the mpileup function of samtools is used to analyze the stacking information at each position on the genome within the block extracted in step one, and an analysis result is obtained (performing in-depth analysis on each genomic position within the interval block to identify potential candidate variants for each sample), and the analysis result is used as the initial candidate variant sites;

[0069] The specific process is as follows:

[0070] The mpileup function of samtools is used to analyze the stacking information at each position on the genome within the block extracted in step one, and an analysis result is obtained (performing in-depth analysis on each genomic position within the interval block to identify potential candidate variants for each sample), and the analysis result is complex information such as base match, mismatch, insertion, deletion, alignment strand orientation, and mapping quality;

[0071] The analysis result is used as the initial candidate variant sites;

[0072] The analysis rules for using the mpileup function of samtools to analyze the stacking information at each position extracted in step one in 6) are as follows:

[0073] 1. The '.' symbol represents a base that exactly matches the positive strand of the reference sequence;

[0074] 2. The ',' symbol represents a base that exactly matches the negative strand of the reference sequence;

[0075] 3. The capital letters 'ATCGN' represent the detected base mismatch cases on the positive strand;

[0076] 4. The lowercase letters 'atcgn' represent the detected base mismatch cases on the negative strand;

[0077] 5. The '*' symbol represents a fuzzy base, usually indicating a potential deletion site;

[0078] 6. The '^' symbol identifies the starting position of a read segment. The alignment quality value of the read segment is obtained by subtracting 33 from the ASCII code immediately following the '^'. The '^' symbol modifies the subsequent base, and the base (ATCGatcgNn) immediately following it represents the first base of the read segment.

[0079] 7. The '$' symbol identifies the ending position of a read segment. The '$' symbol modifies the base preceding it.

[0080] 8. The sequence pattern that conforms to the regular expression '\+[0 - 9]+[ACGTNacgtn]+' indicates that an insertion variant is detected at the current site. For example, +3agg indicates that a 3-base insertion event is detected here.

[0081] 9. The sequence pattern that conforms to the regular expression '-[0 - 9]+[ACGTNacgtn]+' indicates that a deletion variant is detected at the current site. For example, -4CTGA indicates that a 4-base deletion event is detected here.

[0082] This precise parsing provides reliable basic data for subsequent variant site screening and quality assessment.

[0083] Other steps and parameters are the same as those in the first or second specific implementation manner.

[0084] Specific implementation manner four: The difference between this implementation manner and one of the first to third specific implementation manners is that in step 13, the initial candidate variant sites obtained in steps 1 and 2 are screened to obtain the final candidate variant sites. The specific process is as follows:

[0085] 1). Calculate the sequencing coverage depth DP of the initial candidate variant sites obtained in steps 1 and 2 and the allele frequency AF of the corresponding variants (assuming that the reference base at a position is A and the obtained stacked information is AAAAACCGGC, then the possible variants are C or G, and the corresponding allele frequencies are 0.3 and 0.2 respectively. These are two candidate variants).

[0086] 2). For SNP variants (the variant type can be obtained from the stacked information), when DP is greater than or equal to min_depth and AF is greater than or equal to min_alt_freq, the corresponding site is determined as a candidate variant site and is recorded in the candidate site index dictionary candidate_positions.

[0087] For inserted deletion mutations (if the reference base is A and the stacking information is AAAAAAA+AT+AT+AT, then there is a candidate mutation of AT insertion with an allele frequency of 0.3; the mutation type can be obtained from the stacking information), when DP is greater than or equal to min_depth and AF is greater than or equal to min_indel_freq, the corresponding site is determined as a candidate mutation site and recorded in the candidate site index dictionary candidate_positions;

[0088] min_depth represents the minimum depth threshold;

[0089] min_alt_freq represents the minimum mutation frequency threshold;

[0090] min_indel_freq represents the minimum inserted deletion mutation threshold;

[0091] 3) In the candidate site index dictionary candidate_positions, using the position coordinate pos on the chromosome as the unique key index, and storing the characteristic information of each candidate mutation site at the same time;

[0092] The characteristic information of each candidate site includes: contig, pos, ref_base, read_name_list, base_list, read_name_dict, has_snp, and has_indel;

[0093] contig, pos, ref_base, read_name_list, base_list, read_name_dict, has_snp, and has_indel; are included in the parsing results of steps one and two;

[0094] contig represents the chromosome identifier it belongs to;

[0095] pos represents the exact position coordinate of the candidate site on the chromosome;

[0096] ref_base represents the base type of the reference genome at the candidate site;

[0097] read_name_list represents the list of unique identifiers of all sequencing reads covering the candidate site;

[0098] base_list represents the list of base compositions of the candidate site, including the base types of each read at the candidate site and potential inserted deletion information;

[0099] The read_name_dict represents a dictionary mapping the names of sequencing reads to indices, which is used for quick retrieval;

[0100] has_snp indicates whether there is a single nucleotide polymorphism at the candidate site;

[0101] has_indel indicates whether there is an insertion or deletion variation at the candidate site;

[0102] After completing the identification and screening of all candidate variant sites, the candidate_positions dictionary is written into a pickle file in binary format through serialization technology; these serialized data will be reloaded and used in subsequent processes.

[0103] Other steps and parameters are the same as those in any one of the specific embodiments one to three.

[0104] Specific embodiment five: The difference between this embodiment and any one of the specific embodiments one to four is that in step one four, the samtools tool is called to traverse all the sequencing reads on the chromosome aligned in step one one to obtain the original information of each sequencing read (recorded in the BAM file);

[0105] The unique identifier read_id in the original information of each sequencing read is extracted and stored as an index key in the read_info dictionary;

[0106] In the read_info dictionary data structure, multi-dimensional feature information is stored for each sequencing read;

[0107] The read_info dictionary is written into a pickle file through serialization technology to achieve efficient data storage and fast loading;

[0108] The specific process is as follows:

[0109] Multi-dimensional feature information is stored for each sequencing read, including: ref_start, ref_end, read_start, read_end, read_length, has_start_soft_clip, has_end_soft_clip, strand_info, start_soft_seq, end_soft_seq, update_round, and alignment_round;

[0110] ref_start represents the starting position of the sequencing read on the reference genome, using a 1-based coordinate system;

[0111] ref_end represents the ending position of the sequencing read on the reference genome;

[0112] read_start represents the starting position of the alignment sequence within the sequencing read, using a 0-based coordinate system;

[0113] read_end represents the ending position of the alignment sequence within the sequencing read;

[0114] read_length represents the total base length of the sequencing read;

[0115] has_start_soft_clip indicates whether there is a soft-clipped region at the start of the sequencing read;

[0116] has_end_soft_clip indicates whether there is a soft-clipped region at the end of the sequencing read;

[0117] strand_info represents the strand direction information, where "+" represents the forward strand and "-" represents the reverse strand;

[0118] start_soft_seq represents the soft-clipped sequence at the start;

[0119] end_soft_seq represents the soft-clipped sequence at the end;

[0120] update_round represents the update round counter, with an initial value of 1, indicating the number of rounds of processing new fastq during the variant analysis process;

[0121] alignment_round represents the alignment round counter, with an initial value of 1, used to record the number of rounds that the sequencing read participates in the alignment;

[0122] Soft-clipped sequences are sequence fragments that are excluded from the aligned region due to local mismatches during the alignment of the sequencing read but are still retained in the original alignment result; these sequences are very important for discovering variants at the read edges and reconstructing complex structural variants; the system ensures the accuracy of extended alignment and variant detection by identifying and storing these soft-clipped sequences;

[0123] When all read information extraction and structuring are completed, the read_info dictionary is written into a pickle file through serialization technology to achieve efficient data storage and fast loading.

[0124] Other steps and parameters are the same as those in any one of the specific embodiments one to four.

[0125] Specific Embodiment 6: The difference between this embodiment and any one of Specific Embodiments 1 to 5 is that in Step 21, the final candidate variant site information stored in the candidate_positions dictionary saved in Step 13 from the disk and the sequencing read information stored in the read_info dictionary recorded in Step 14 are loaded;

[0126] After the sequencer outputs a new batch of sequencing reads (multiple sequencing reads) (50bp or other values set by the user), selective filtering is performed on the new batch of sequencing data, and only the reads corresponding to those existing in the loaded read_info dictionary are retained;

[0127] The retained new batch of sequencing reads is stored in FASTQ format;

[0128] Subsequent steps to construct events are based on this loaded result, and subsequent results need to be updated in candidate_positions and read_info;

[0129] The specific process is as follows:

[0130] 1) Parallel processing is achieved at the chromosome level. For each chromosome, the final candidate variant site information stored in the candidate_positions dictionary saved in Step 13 from the disk and the sequencing read information stored in the read_info dictionary recorded in Step 14 are loaded;

[0131] Among them, the candidate_positions dictionary uses chromosome position coordinates as the only key value and stores multi-dimensional feature information of variant sites;

[0132] The read_info dictionary uses the read unique identifier as the key value and records key information such as the alignment position, strand direction, and soft clipping of each read;

[0133] 2) After the sequencer outputs a new batch of sequencing reads, selective filtering is performed on the new batch of sequencing data, and only the reads corresponding to those existing in the read_info dictionary loaded in 1) are retained; the specific process is as follows:

[0134] Sequentially traverse each sequencing read in the new batch of sequencing reads, extract the unique identifier read_id of each sequencing read. If read_id exists in the loaded read_info dictionary, the corresponding sequencing read information is written into a newly generated temporary FASTQ file; if read_id does not exist in the loaded read_info dictionary, the corresponding sequencing read is deleted;

[0135] After screening all the sequencing reads in the new batch of sequencing reads, the retained sequencing reads are stored in FASTQ format;

[0136] Use the bgzip tool to efficiently compress the stored sequencing reads in FASTQ format, and then call the fqidx function of the samtools tool to create an index for the sequencing reads in FASTQ format to support subsequent random access operations;

[0137] Other steps and parameters are the same as those in any one of the specific embodiments one to five.

[0138] Specific embodiment seven: The difference between this embodiment and any one of the specific embodiments one to six is that in step two two, event tuples are constructed based on step two one; the specific process is as follows:

[0139] Define multiple event type constants, including: NEW_READ_START, SNP_START, INSERTION, DELETION_START, DELETION_END, and NEW_READ_END;

[0140] NEW_READ_START represents the start of each sequencing read retained in step two one (in FASTQ);

[0141] SNP_START represents the start of a single nucleotide variant (in the candidate_positions dictionary);

[0142] INSERTION represents an insertion variant (in the candidate_positions dictionary);

[0143] DELETION_START represents the start of a deletion variant (in the candidate_positions dictionary);

[0144] DELETION_END represents the end of a deletion variant (in the candidate_positions dictionary);

[0145] NEW_READ_END represents the end of each sequencing read retained in step two one (in FASTQ);

[0146] According to the strand direction information (strand_info) of each sequencing read retained in step two one, different processing strategies are adopted (which determines the pos in the triple):

[0147] For reads on the forward strand ("+"), extend downstream from the ref_end position on the reference genome;

[0148] For reads on the reverse strand ("-"), extend upstream from the ref_start position on the reference genome;

[0149] This strand-specific processing ensures the correct direction of read extension.

[0150] To achieve efficient read extension and update of variant site stacking information, the "event tuple" mechanism and the scan line algorithm are introduced;

[0151] Each event tuple is a triple data structure (position, event_type, index); multiple events in one block;

[0152] Among them, position represents the exact position coordinate of the event on the genomic reference sequence;

[0153] event_type represents the event type identifier (NEW_READ_START, SNP_START, INSERTION, DELETION_START, DELETION_END, and NEW_READ_END);

[0154] index represents the index value in the event type (the serial number of NEW_READ_START, SNP_START, INSERTION, DELETION_START, DELETION_END, and NEW_READ_END).

[0155] Other steps and parameters are the same as one of the specific embodiments one to six.

[0156] Specific embodiment eight: The difference between this embodiment and one of the specific embodiments one to seven is that in step 23, a sorted_events sorted list is constructed based on the event tuples constructed in step 22. Traverse each event tuple (pos, event_type, idx) in the sorted list sorted_events, and execute different processing logics according to the event type to obtain the information of each read segment in the sorted list sorted_events under different processing logics. The specific process is as follows:

[0157] After constructing all event tuples, sort the event tuples constructed in step 22 from front to back according to the genomic position coordinates, and judge whether there is only 1 event at each position;

[0158] If so, form a sorted_events sorted list by sorting from front to back according to the genomic position coordinates;

[0159] Otherwise, based on the genomic position coordinates sorted from front to back, sort the positions with more than 1 event in the order of event types to finally form the sorted_events list (first sort by the first position in the triple. If there are multiple event tuples at one position, then sort by the order of event types in the event tuple);

[0160] The order of event types is: NEW_READ_START, SNP_START, INSERTION, DELETION_START, DELETION_END, NEW_READ_END;

[0161] During the event processing, maintain two key dynamic data structures: the active read segment set (active_new_read_set) and the candidate variant activity dictionary (active_candidate_dict). Among them, active_new_read_set is an efficient set data structure used to track all read segment indexes in the active state (that is, not yet ended); while active_candidate_dict is a complex nested dictionary structure used to record the read segment coverage related to each candidate variant site;

[0162] By sequentially traversing each event tuple (pos, event_type, idx) in the sorted list sorted_events, execute different processing logics according to the event type:

[0163] 1) When encountering a NEW_READ_START event, add the read segment identifier read_id to the active_new_read_set, indicating that the read segment has entered the active state;

[0164] 2) When encountering a NEW_READ_END event, first remove the read segment from the active_new_read_set, indicating that the read segment activity ends;

[0165] Subsequently, check whether the read segment has been marked as needing to be realigned again (that is, whether it is in the reads_to_realign set); if it has been marked, do nothing (take no action and continue to process the following); if it has not been marked, perform differential updates according to the strand direction information (strand_info) of the read segment: for forward-strand read segments, update the human reference end position (ref_end) of the forward-strand read segment to extend the reference end position downstream by an appropriate length, and the appropriate length is the length of the current read segment; for reverse-strand read segments, update the human reference start position (ref_start) of the reverse-strand read segment to extend the reference start position upstream by an appropriate length, and the appropriate length is the length of the current read segment;

[0166] Meanwhile, update the update_round counter of the read segment to the current processing round to record the processing status;

[0167] 3) When encountering the SNP_START event, check all the read segments in the active_new_read_set;

[0168] If the active_new_read_set is empty, do nothing (take no action and continue with the subsequent processing);

[0169] If the active_new_read_set is not empty, for each read segment, add the read segment to the reads_to_realign set and mark it as needing to be realigned;

[0170] If the update_round of the read segment is inconsistent with the alignment_round, add the read segment to the middle_fastq_needed set, indicating that the unaligned read segments in the middle round need to be extracted and stored as FASTQ-formatted sequencing read segments to ensure sequence continuity; if the update_round of the read segment is consistent with the alignment_round, the read segment is not added to the middle_fastq_needed set;

[0171] 4) When encountering the INSERTION event, adopt a processing logic similar to that of the SNP_START event, mark all the read segments in the current active_new_read_set as needing to be realigned, and if the update_round of the read segment is inconsistent with the alignment_round, add the read segment to the middle_fastq_needed set; if the update_round of the read segment is consistent with the alignment_round, the read segment is not added to the middle_fastq_needed set;

[0172] 5) When encountering the DELETION_START event, create a record for the deletion mutation in the active_candidate_dict, and add the indices of all the read segments in the current active_new_read_set to the "start_covered" list of the mutation, indicating that all the current read segments cover the start position of the deletion mutation;

[0173] 6) When encountering the DELETION_END event, add the indices of all current reads to the "end_covered" list corresponding to the deletion variant, indicating that all reads in the current active_new_read_set cover the end position of the deletion variant. Subsequently, traverse all read indices in the "start_covered" and "end_covered" lists of the deletion variant, add all read indices in the "start_covered" and "end_covered" lists to the reads_to_realign set, and add them to the middle_fastq_needed set. This mechanism ensures that reads that fully cover the deletion variant region are properly processed.

[0174] 7) Finally, convert the read indices in the reads_to_realign set into a list of read identifiers as the input for subsequent batch alignment processing.

[0175] Other steps and parameters are the same as those in any one of the specific embodiments one to seven.

[0176] Specific Embodiment Nine: The difference between this embodiment and any one of the specific embodiments one to eight is that in step 24, the differences between each read in the reads_to_realign set in step 23 and the reference sequence are analyzed site by site on the genome, and the stacking information of the corresponding sites is updated. The specific process is as follows:

[0177] For reads that need to be realigned, a batch processing technique (controlled by the batch_size parameter) is used to improve the calculation efficiency.

[0178] 1) Initialize l = 1;

[0179] 2) Pass the information of the l-th read in the reads_to_realign set in step 23 to the align_and_update_candidates_batch function;

[0180] 3) In the align_and_update_candidates_batch function, if the l-th read in the reads_to_realign set is a forward-strand read, extract a reference sequence of a specific length in the downstream direction starting from the termination position of the forward-strand read reference genome (known according to read_info).

[0181] If the update round of the read is not the same as the alignment round, there is an intermediate read sequence;

[0182] If the update round of the read segment is the same as the alignment round, there is no intermediate read segment sequence;

[0183] When there is an intermediate read segment sequence, the reference sequence of a specific length is the length of the soft-clipped sequence plus the length of the intermediate read segment sequence (determined by middle_fastq_needed) plus the length of the current read segment sequence, with an additional 200 bases added; subsequently, the soft-clipped sequence, the intermediate read segment sequence, and the current read segment sequence are concatenated in order to form the sequence to be aligned;

[0184] When there is no intermediate read segment sequence, the reference sequence of a specific length is the length of the soft-clipped sequence plus the length of the current read segment sequence, with an additional 200 bases added; subsequently, the soft-clipped sequence and the current read segment sequence are concatenated in order to form the sequence to be aligned;

[0185] The length of the intermediate read segment sequence (if any) specifically means: if the update round of the read segment is different from the alignment round, it means that it is necessary to extract the length of the sequencing read segment aligned in the intermediate round to ensure sequence continuity;

[0186] 4) In the align_and_update_candidates_batch function, if the l-th read segment in the reads_to_realign set is a negative-strand read segment, starting from the starting position of the negative-strand read segment reference genome (known from read_info), extract a reference sequence of a specific length in the upstream direction; perform reverse complementary conversion on the extracted reference sequence of a specific length to match the biological characteristics of the negative-strand read segment;

[0187] The calculation process of the reference sequence of a specific length is as follows:

[0188] If the update round of the read segment is different from the alignment round, there is an intermediate read segment sequence (determined by middle_fastq_needed);

[0189] If the update round of the read segment is the same as the alignment round, there is no intermediate read segment sequence;

[0190] When there is an intermediate read segment sequence, the reference sequence of a specific length is the length of the soft-clipped sequence (known from read_info) plus the length of the intermediate read segment sequence (determined by middle_fastq_needed) plus the length of the current read segment sequence, with an additional 200 bases added; subsequently, the soft-clipped sequence, the intermediate read segment sequence, and the current read segment sequence are concatenated in order to form the sequence to be aligned;

[0191] When there is no middle read segment sequence, the reference sequence of a specific length is the sum of the soft-clipped sequence length and the current read segment sequence length, plus an additional 200 bases; subsequently, the soft-clipped sequence and the current read segment sequence are concatenated in order to form the sequence to be aligned.

[0192] The length of the middle read segment sequence (if it exists) specifically refers to: if the update round of the read segment is not the same as the alignment round, it means that the length of the sequencing read segment aligned in the middle round needs to be extracted to ensure sequence continuity.

[0193] 5), For the forward-strand read segments, use the extension mode of the KSW2 algorithm to perform a local sequence alignment on the reference sequence of the specific length extracted in 3) and the sequence to be aligned constructed. The specific process is as follows:

[0194] Set the KSW2 algorithm parameters: the scoring matrix (mat) optimized for DNA sequences, the deletion penalty (q), and the extension penalty (e).

[0195] Use the reference sequence of the specific length extracted in 3) and the sequence to be aligned constructed as the input of the KSW2 algorithm, and the KSW2 algorithm obtains the local sequence alignment result.

[0196] The alignment result is output in the form of a CIGAR string, recording the matching, mismatching, insertion, and deletion events during the alignment process.

[0197] 6), For the reverse-strand read segments, use the extension mode of the KSW2 algorithm to perform a local sequence alignment on the reference sequence of the specific length extracted in 4) and the sequence to be aligned constructed. The specific process is as follows:

[0198] Set the KSW2 algorithm parameters: the scoring matrix (mat) optimized for DNA sequences, the deletion penalty (q), and the extension penalty (e).

[0199] Use the reference sequence of the specific length extracted in 4) and the sequence to be aligned constructed as the input of the KSW2 algorithm, and the KSW2 algorithm obtains the local sequence alignment result.

[0200] The alignment result is output in the form of a CIGAR string, recording the matching, mismatching, insertion, and deletion events during the alignment process.

[0201] 7), Based on the alignment result in the form of a CIGAR string, analyze the differences between the l-th read segment and the reference sequence site by site on the genome.

[0202] When no base mismatches, insertion variations, or deletion variations are detected, no processing is performed.

[0203] When a base mismatch, insertion variation, or deletion variation is detected, it is determined whether the corresponding site exists in the set of candidate variation sites; if it exists, the stacking information corresponding to the site (information update of base mismatch, insertion variation, or deletion variation) and the relevant read segment identifier (the last aligned position) are updated; if it does not exist, the stacking information corresponding to the site is not updated;

[0204] 8), Let l = l + 1, and repeat steps 2)-7) until all read segments in the reads_to_realign set in step two and three are completed (the same operation is performed for each read segment), and the updated read_info and candidate_positions dictionaries are serialized and saved to disk for use in subsequent processing steps. Through this iterative read segment extension and variation information update mechanism, the system can continuously improve the variation detection results as the sequencing data continues to be generated.

[0205] Other steps and parameters are the same as those in any one of the specific embodiments one to eight.

[0206] Specific embodiment ten: The difference between this embodiment and any one of the specific embodiments one to nine is that in step three, variation detection is performed based on step two to obtain a variation detection result (variation detection quality score), and the variation detection result is passed into the VCF file; the specific process is as follows:

[0207] After completing all read segment extensions and updates of candidate variation sites, variation detection can be performed;

[0208] Step three one: Construct the following information for each site in the candidate variation sites in the candidate_positions dictionary in step two: the chromosome of the candidate variation site on the human reference genome, the position of the candidate variation site on the human reference genome, the reference sequence base of the candidate variation site on the human reference genome, the variant base of the candidate variation site, the number of all covered segments (depth DP, for example, if the reads aligned to the current position on the human reference genome are read segment 1, read segment 2, read segment 3, read segment 4, then DP is 4), the number of segments supporting the reference sequence base N ref (Supporting the reference sequence base means the same as the reference sequence), the number of segments supporting the variant base N alt (Supporting the variant base means different from the reference sequence), the number of positive-strand segments and the number of negative-strand segments of the segments supporting the variant base; all are information on the human reference genome;

[0209] Step three two: Use the following probability model for variation detection, and the process is as follows:

[0210] Use the binomial model to compare all possibilities of genotypes, and the calculation formula of the binomial model is as follows:

[0211]

[0212] Among them, Likelihood represents probability;

[0213] (0 / 00 indicates no mutation;

[0214] (0 / 10 indicates the existence of a heterozygous mutation;

[0215] (1 / 10 indicates the existence of a homozygous mutation;

[0216] θ is the probability that the sequencing data is misaligned (default value is 0.03), and it is obtained independently through experience;

[0217] N ref represents the number of fragments supporting the reference sequence base;

[0218] N alt represents the number of fragments supporting the variant base;

[0219] p is the prior probability of the heterozygous genotype, and the prior distribution is configured to be equal for all three genotypes here, p = 1 / 3;

[0220] Generate a Phred-formatted variant detection quality score for each variant;

[0221] Phred = -10 × log(1 - P varient )

[0222] where P varient is the probability that the current variant exists, that is, the maximum value among Likelihood(0 / 0), Likelihood(0 / 1), and Likelihood(1 / 1);

[0223] Phred represents the variant detection quality score;

[0224] The variant detection result (variant detection quality score) is passed into the VCF file;

[0225] The VCF file includes: reference sequence base, variant base, genomic position, variant detection quality score.

[0226] Other steps and parameters are the same as those in any one of the specific embodiments one to nine.

[0227] The present invention can also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and deformations according to the present invention, but these corresponding changes and deformations should all fall within the protection scope of the appended claims of the present invention.

Claims

1. A real-time SNV and Indel detection method for next-generation sequencing data, characterized in that: The specific process of the method is as follows: Step 1: Obtain the alignment results of the sequencing reads on the human reference genome; Based on the alignment results, obtain the final candidate variant sites; Based on the alignment results, extract the unique identifier read_id in the original information of each sequencing read as the index key and store it in the read_info dictionary; write the read_info dictionary into a pickle file through serialization technology; The specific process is as follows: Step 1.1: Obtain the alignment results of the sequencing reads on the human reference genome, divide each chromosome in the aligned human reference genome into non-overlapping blocks, and extract the stacking information at each position on the genome within the blocks; Step 1.2: Use the mpileup function of samtools to parse the stacking information at each position on the genome within the blocks extracted in Step 1.1 to obtain the parsing results, and the parsing results serve as the initial candidate variant sites; Step 1.3: Screen the initial candidate variant sites obtained in Step 1.2 to obtain the final candidate variant sites; Step 1.4: Call the samtools tool to traverse all the sequencing reads on the chromosome aligned in Step 1.1 to obtain the original information of each sequencing read; Extract the unique identifier read_id in the original information of each sequencing read as the index key and store it in the read_info dictionary; In the read_info dictionary data structure, store the feature information for each sequencing read; Write the read_info dictionary into a pickle file through serialization technology; Step 2: After the sequencer outputs a new batch of sequencing reads, perform selective filtering on the new batch of sequencing data to obtain the filtered reads; Construct event tuples; Traverse each event tuple (pos, event_type, idx), execute different processing logics according to the event type, and obtain the information of each read in the sorted list sorted_events under different processing logics; Analyze the differences between each read and the reference sequence site by site on the genome, and update the stacking information at the corresponding sites; The specific process is as follows: Step 2.1: Load from the disk the information of the final candidate variant sites stored in the candidate_positions dictionary saved in Step 1.3 and the information of the sequencing reads stored in the read_info dictionary recorded in Step 1.4; After the sequencer outputs a new batch of sequencing reads, perform selective filtering on the new batch of sequencing data, and only retain the reads that have a corresponding relationship with the loaded read_info dictionary; The retained new batch of sequencing reads is stored in FASTQ format; Step 2.2: Construct event tuples based on Step 2.1; Step 23: Construct a sorted_events sorted list based on the event tuples constructed in Step 22. Traverse each event tuple (pos, event_type, idx) in the sorted list sorted_events, and execute different processing logics according to the event type to obtain the information of each read segment in the sorted list sorted_events under different processing logics; Step 24: Analyze the differences between each read segment in the reads_to_realign set in Step 23 and the reference sequence site by site on the genome, and update the stacking information of the corresponding sites; Step 3: Perform variant detection based on Step 2 to obtain the variant detection results, and pass the variant detection results into the VCF file.

2. The real-time SNV and Indel detection method for next-generation sequencing data according to claim 1, characterized in that: In Step 11, obtain the alignment results of the sequencing reads on the human reference genome, divide each chromosome in the aligned human reference genome into non-overlapping blocks, and extract the stacking information of each position on the genome within the blocks; The specific process is as follows: 1). Use a sequencer to obtain human sequencing reads in real time, take the sequencing reads with a length of 150bp as the initial data block, and store them in FASTQ format; 2). Use the sequence alignment tool minimap2 to map the sequencing reads on the initial data block to the human reference genome to obtain the original alignment results in SAM format; 3). Call the bioinformatics tool samtools to convert the original alignment results in SAM format to the original alignment results in BAM format; 4). Sort the aligned sequencing reads in the original alignment results in BAM format from front to back according to their positions on the human reference genome, and establish an index; 5). Based on the positions and indexes of the sequencing reads on the human reference genome obtained in 4), divide each chromosome in the aligned human reference genome into M non-overlapping blocks according to a preset interval length. Each block is uniquely identified in the form of a triple of "chromosome: start position: end position"; M is a positive integer; 6). Call GNU parallel to perform parallel computing on different blocks of homologous chromosomes; Use the mpileup function of samtools to extract the stacking information of each position on the genome within the block. The stacking information includes the reference sequence name of the current position on the human reference genome, the genomic coordinates of the current position on the human reference genome, the reference base of the current position on the human reference genome, the number of sequencing fragments covering the current position, the bases of the sequencing fragments covering the current position, and the ASCII code of the base quality of the sequencing fragments covering the current position.

3. The real-time SNV and Indel detection method for next-generation sequencing data according to claim 2, wherein: In Step 12, use the mpileup function of samtools to parse the stacking information of each position on the genome within the block extracted in Step 11 to obtain the parsing results, and the parsing results are used as the initial candidate variant sites; The specific process is as follows: Use the mpileup function of samtools to parse the pileup information at each position on the genome within the block extracted in Step 1, obtaining a parsing result, which includes base matching, mismatching, insertion, deletion, alignment strand orientation, and alignment quality; The parsing result serves as the initial candidate variant sites.

4. A real-time SNV and Indel detection method for second-generation sequencing data according to claim 3, characterized in that: In Step 13, screen the initial candidate variant sites obtained in Step 12 to obtain the final candidate variant sites; The specific process is as follows: 1). Calculate the sequencing coverage depth DP of the initial candidate variant sites obtained in Step 12 and the allele frequency AF of the corresponding variants; 2). For SNP variants, when DP is greater than or equal to min_depth and AF is greater than or equal to min_alt_freq, the corresponding site is determined as a candidate variant site and recorded in the candidate site index dictionary candidate_positions; For insertion-deletion variants, when DP is greater than or equal to min_depth and AF is greater than or equal to min_indel_freq, the corresponding site is determined as a candidate variant site and recorded in the candidate site index dictionary candidate_positions; min_depth represents the minimum depth threshold; min_alt_freq represents the minimum variant frequency threshold; min_indel_freq represents the minimum insertion-deletion variant threshold; 3). In the candidate site index dictionary candidate_positions, use the position coordinate pos on the chromosome as the unique key index, and simultaneously store the characteristic information of each candidate variant site; The characteristic information of each candidate site includes: contig, pos, ref_base, read_name_list, base_list, read_name_dict, has_snp, and has_indel; contig represents the chromosome identifier to which it belongs; pos represents the position coordinate of the candidate site on the chromosome; ref_base represents the base type of the reference genome at the candidate site; read_name_list represents the list of unique identifiers of all sequencing reads covering the candidate site; base_list represents the list of base compositions at the candidate site, including the base types of each read at the candidate site and potential insertion-deletion information; read_name_dict represents the mapping relationship dictionary between the sequencing read names and indices; has_snp indicates whether there is a single nucleotide polymorphism at the candidate site; has_indel indicates whether there is an insertion or deletion variant at the candidate site; Write the candidate_positions dictionary into a pickle file in binary format through serialization technology.

5. The real-time SNV and Indel detection method for next-generation sequencing data according to claim 4, characterized in that: In Step 14, call the samtools tool to traverse all the sequencing reads on the chromosome aligned in Step 11 to obtain the original information of each sequencing read; Extract the unique identifier read_id in the original information of each sequencing read as the index key and store it in the read_info dictionary; In the read_info dictionary data structure, store the feature information for each sequencing read; Write the read_info dictionary into a pickle file through serialization technology; The specific process is as follows: Each sequencing read stores feature information, including: ref_start, ref_end, read_start, read_end, read_length, has_start_soft_clip, has_end_soft_clip, strand_info, start_soft_seq, end_soft_seq, update_round, and alignment_round; ref_start represents the starting position of the sequencing read on the reference genome; ref_end represents the ending position of the sequencing read on the reference genome; read_start represents the starting position of the aligned sequence within the sequencing read; read_end represents the ending position of the aligned sequence within the sequencing read; read_length represents the total base length of the sequencing read; has_start_soft_clip indicates whether there is a soft clip region at the starting end of the sequencing read; has_end_soft_clip indicates whether there is a soft clip region at the ending end of the sequencing read; strand_info represents the strand direction information, where "+" represents the forward strand and "-" represents the reverse strand; start_soft_seq represents the soft clip sequence at the starting end; end_soft_seq represents the soft clip sequence at the ending end; update_round represents the update round counter, with an initial value of 1; alignment_round represents the alignment round counter, with an initial value of 1 round.

6. The real-time SNV and Indel detection method for next-generation sequencing data according to claim 5, wherein: In step 21, load the final candidate variant site information stored in the candidate_positions dictionary saved in step 13 from the disk and the sequencing read information stored in the read_info dictionary recorded in step 14; After the sequencer outputs a new batch of sequencing reads, perform selective filtering on the new batch of sequencing data, and only retain the reads that have a corresponding relationship with the loaded read_info dictionary; The retained new batch of sequencing reads is stored in FASTQ format; The specific process is as follows: 1) Load the final candidate variant site information stored in the candidate_positions dictionary saved in step 13 from the disk and the sequencing read information stored in the read_info dictionary recorded in step 14; 2) After the sequencer outputs a new batch of sequencing reads, perform selective filtering on the new batch of sequencing data, and only retain the reads that have a corresponding relationship with the read_info dictionary loaded in 1); The specific process is as follows: Sequentially traverse each sequencing read in the new batch of sequencing reads, extract the unique identifier read_id of each sequencing read. If read_id exists in the loaded read_info dictionary, write the corresponding sequencing read information to a newly generated temporary FASTQ file; if read_id does not exist in the loaded read_info dictionary, delete the corresponding sequencing read. After screening all the sequencing reads in the new batch of sequencing reads, the retained sequencing reads are stored in FASTQ format. Use the bgzip tool to compress the stored sequencing reads in FASTQ format, and then call the fqidx function of the samtools tool to create an index for the sequencing reads in FASTQ format.

7. The real-time SNV and Indel detection method for next-generation sequencing data according to claim 6, wherein: In step 22, event tuples are constructed based on step 21. The specific process is as follows: Define event type constants, including: NEW_READ_START, SNP_START, INSERTION, DELETION_START, DELETION_END, and NEW_READ_END. NEW_READ_START represents the start of each sequencing read retained in step 21. SNP_START represents the start of a single nucleotide variant. INSERTION represents an insertion variant. DELETION_START represents the start of a deletion variant. DELETION_END represents the end of a deletion variant. NEW_READ_END represents the end of each sequencing read retained in step 21. According to the strand direction information of each sequencing read retained in step 21, different processing strategies are adopted: For reads on the forward strand, extend downstream from the ref_end position on the reference genome. For reads on the reverse strand, extend upstream from the ref_start position on the reference genome. Each event tuple is a triple data structure (position, event_type, index). Among them, position represents the position coordinate of the event on the genomic reference sequence. event_type represents the event type identifier. index represents the index value in the event type.

8. A real-time SNV and Indel detection method for next-generation sequencing data according to claim 7, characterized in that: In step 23, a sorted_events sorted list is constructed based on the event tuples constructed in step 22. Traverse each event tuple (pos, event_type, idx) in the sorted list sorted_events, and execute different processing logics according to the event type to obtain the information of each read in the sorted list sorted_events under different processing logics. The specific process is as follows: Sort the event tuples constructed in step 22 in ascending order of genomic position coordinates, and judge whether there is only 1 event at each position. If so, sort them in ascending order of genomic position coordinates to form a sorted_events sorted list. Otherwise, on the basis of sorting the genomic position coordinates from front to back, sort the positions with more than 1 event in the order of event types to finally form the sorted_events sorted list; The order of event types is: NEW_READ_START, SNP_START, INSERTION, DELETION_START, DELETION_END, NEW_READ_END; Traverse each event tuple (pos, event_type, idx) in the sorted list sorted_events and execute different processing logics according to the event type: 1) When encountering a NEW_READ_START event, add the read identifier read_id to the active_new_read_set, indicating that the read has entered the active state; 2) When encountering a NEW_READ_END event, first remove the read from the active_new_read_set, indicating the end of the read activity; Subsequently, check whether the read has been marked as needing to be realigned; If it has been marked, do nothing; If it has not been marked, perform differential updates according to the strand direction information of the read: for forward-strand reads, update the human reference end position of the forward-strand read to extend the reference end position downstream by an appropriate length, and the appropriate length is the length of the current read; for reverse-strand reads, update the human reference start position of the reverse-strand read to extend the reference start position upstream by an appropriate length, and the appropriate length is the length of the current read; At the same time, update the update round counter of the read to the current processing round to record the processing status; 3) When encountering a SNP_START event, check all reads in the active_new_read_set; If the active_new_read_set is empty, do nothing; If the active_new_read_set is not empty, for each read, add the read to the reads_to_realign set and mark it as needing to be realigned; If the update round of the read is inconsistent with the alignment round, add the read to the middle_fastq_needed set; if the update round of the read is consistent with the alignment round, the read is not added to the middle_fastq_needed set; 4) When encountering an INSERTION event, mark all reads in the current active_new_read_set as needing to be realigned. If the update round of the read is inconsistent with the alignment round, add the read to the middle_fastq_needed set; if the update round of the read is consistent with the alignment round, the read is not added to the middle_fastq_needed set; 5) When encountering the DELETION_START event, create a record for the deletion variant in the active_candidate_dict, and add the indices of all reads in the current active_new_read_set to the "start_covered" list of the variant; 6) When encountering the DELETION_END event, add the indices of all current reads to the "end_covered" list of the corresponding deletion variant; subsequently, traverse all read indices in the "start_covered" and "end_covered" lists of the deletion variant, add all read indices in the "start_covered" and "end_covered" lists to the reads_to_realign set, and add them to the middle_fastq_needed set; 7) Finally, convert the read indices in the reads_to_realign set to a list of read identifiers.

9. A real-time SNV and Indel detection method for second-generation sequencing data according to claim 8, characterized in that: In step 24, analyze the differences between each read in the reads_to_realign set in step 23 and the reference sequence site by site on the genome, and update the stacking information at the corresponding sites; The specific process is as follows: 1) Initialize l = 1; 2) Pass the information of the l-th read in the reads_to_realign set in step 23 to the align_and_update_candidates_batch function; 3) In the align_and_update_candidates_batch function, if the l-th read in the reads_to_realign set is a forward-strand read, extract a reference sequence of a specific length from the termination position of the forward-strand read reference genome in the downstream direction; The calculation process of the reference sequence of a specific length is as follows: If the update round of the read is inconsistent with the alignment round, there is an intermediate read sequence; If the update round of the read is consistent with the alignment round, there is no intermediate read sequence; When there is an intermediate read sequence, the reference sequence of a specific length is the sum of the soft-clipped sequence length, the intermediate read sequence length, the current read sequence length, plus an additional 200 bases; Subsequently, sequentially splice the soft-clipped sequence, the intermediate read sequence, and the current read sequence to form the sequence to be aligned; When there is no intermediate read sequence, the reference sequence of a specific length is the sum of the soft-clipped sequence length, the current read sequence length, plus an additional 200 bases; Subsequently, sequentially splice the soft-clipped sequence and the current read sequence to form the sequence to be aligned; 4) In the align_and_update_candidates_batch function, if the l-th read in the reads_to_realign set is a reverse-strand read, extract a reference sequence of a specific length from the starting position of the reverse-strand read reference genome in the upstream direction; perform reverse complementary conversion on the extracted reference sequence of a specific length; The calculation process of the reference sequence with a specific length is as follows: If the update round of the read segment is inconsistent with the alignment round, there exists an intermediate read segment sequence; If the update round of the read segment is consistent with the alignment round, there is no intermediate read segment sequence; When there exists an intermediate read segment sequence, the reference sequence with a specific length is the sum of the soft-clipped sequence length, the intermediate read segment sequence length, and the current read segment sequence length, plus an additional 200 bases; Subsequently, the soft-clipped sequence, the intermediate read segment sequence, and the current read segment sequence are concatenated in order to form the sequence to be aligned; When there is no intermediate read segment sequence, the reference sequence with a specific length is the sum of the soft-clipped sequence length and the current read segment sequence length, plus an additional 200 bases; Subsequently, the soft-clipped sequence and the current read segment sequence are concatenated in order to form the sequence to be aligned; 5) For the forward-strand read segments, use the extension mode of the KSW2 algorithm to perform local sequence alignment on the reference sequence with the specific length extracted in 3) and the constructed sequence to be aligned; the specific process is as follows: Set the KSW2 algorithm parameters: scoring matrix, gap penalty, and extension penalty; Use the reference sequence with the specific length extracted in 3) and the constructed sequence to be aligned as the input of the KSW2 algorithm, and the KSW2 algorithm obtains the local sequence alignment result The alignment result is output in the form of a CIGAR string, recording the matching, mismatching, insertion, and deletion events during the alignment process; 6) For the reverse-strand read segments, use the extension mode of the KSW2 algorithm to perform local sequence alignment on the reference sequence with the specific length extracted in 4) and the constructed sequence to be aligned; the specific process is as follows: Set the KSW2 algorithm parameters: scoring matrix, gap penalty, and extension penalty; Use the reference sequence with the specific length extracted in 4) and the constructed sequence to be aligned as the input of the KSW2 algorithm, and the KSW2 algorithm obtains the local sequence alignment result; The alignment result is output in the form of a CIGAR string, recording the matching, mismatching, insertion, and deletion events during the alignment process; 7) Based on the alignment result in the form of a CIGAR string, analyze the differences between the l-th read segment and the reference sequence site by site on the genome; When no base mismatches, insertion mutations, or deletion mutations are detected, no processing is performed; When base mismatches, insertion mutations, or deletion mutations are detected, determine whether the corresponding site exists in the candidate variant site set; If it exists, update the stacking information and related read segment identifiers of the corresponding site; If it does not exist, do not update the stacking information of the corresponding site; 8) Let l = l + 1, and repeat steps 2)-7) until all the read segments in the reads_to_realign set in step two and three are completed, and serialize and save the updated read_info and candidate_positions dictionaries to disk.

10. A real-time SNV and Indel detection method for next-generation sequencing data according to claim 9, characterized in that: In step three, mutation detection is performed based on step two to obtain the mutation detection result, and the mutation detection result is passed into the VCF file; The specific process is as follows: Step 3.1: Construct the following information for each variant candidate site in the candidate_positions dictionary obtained in Step 2: the chromosome of the variant candidate site on the human reference genome, the position of the variant candidate site on the human reference genome, the reference sequence base of the variant candidate site on the human reference genome, the variant base of the variant candidate site, the total number of all coverage segments, the number N of segments supporting the reference sequence base ref , the number N of segments supporting the variant base alt , the number of positive-strand segments and the number of negative-strand segments of the segments supporting the variant base; Step three two: Use the following probability model for mutation detection, and the process is as follows: Use the binomial model to compare all possibilities of genotypes. The calculation formula of the binomial model is as follows: Where Likelihood represents probability; (0 / 0) represents no mutation; (0 / 1) indicates the presence of a heterozygous variant; (1 / 1) indicates the presence of a homozygous variant; θ is the probability that the sequencing data is misaligned; N ref represents the number of segments that support the reference sequence bases; N alt Indicates the number of segments supporting the variant base; p is the prior probability of the heterozygous genotype, p = 1 / 3; Generate a Phred-formatted variant detection quality score for each variant; Phred = -10×log(1 - P varient ) where P varient is the probability of the current mutation existing, i.e., the maximum value among Likelihood(0 / 0), Likelihood(0 / 1), and Likelihood(1 / 1); Phred represents the variant detection quality score; The variant detection results are passed into the VCF file; The VCF file includes: reference sequence base, variant base, genomic position, variant detection quality score.

Images