Method for detecting microbial structure variation, mutation and insertion and deletion based on next-generation sequencing

Abstract

The invention relates to the technical field of bioinformatics, in particular to a method for detecting microbial structure variation, mutation and insertion and deletion based on next-generation sequencing. On the basis of a short-read-length Illumina next-generation sequencing strategy, a current mainstream biological tool is tested, and automatic analysis of microbial metagenome structure variation (SV), single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel) is realized by applying and combining a self-developed program; the problems of low sensitivity, high false positive and low efficiency during variation detection under the metagenome background in the prior art are obviously solved.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, and in particular to a method for detecting structural variations, mutations, and insertions / deletions in microorganisms based on next-generation sequencing. Background Technology

[0002] Structural variations (SVs), mutations (SNVs), and small fragment insertions and deletions (Indels) in microorganisms not only affect microbial genome molecules and cellular processes, but also their regulatory functions and gene expression, which may further lead to cells from the same species exhibiting complex phenotypes and different growth rates.

[0003] Identification of microbial SVs typically requires assembly using second-generation and third-generation sequencing before metagenomic analysis. Furthermore, current technologies treat all reference genomes as independent entities, failing to address the interference problem caused by similar genomes. This results in low sensitivity, high false positives, and low efficiency when performing variant detection in a metagenomic context. Summary of the Invention

[0004] To address the aforementioned problems, this invention provides a method for detecting structural variations, mutations, and insertions / deletions in microorganisms based on next-generation sequencing. This invention significantly solves the problems of low sensitivity, high false positives, and low efficiency faced by existing technologies in variation detection within a metagenomic context by innovatively grouping reference genes, implementing a stratified candidate screening mechanism, and optimizing the data processing workflow.

[0005] To achieve the above objectives, the present invention provides the following technical solution: This invention provides a method for detecting structural variations, mutations, and insertions / deletions in microorganisms based on next-generation sequencing, comprising the following steps: Calculate the similarity of reference genomes and group them: According to the biological classification method, calculate the similarity of reference sequences within each species, and use 0.95 as the threshold. Those with similarity greater than the threshold are grouped into a reference chromosome group, and the chromosome grouping annotation file is obtained. Raw sequencing data filtering and quality control: The raw sequencing data of the microorganisms to be tested are filtered and quality controlled to obtain quality control data; Removing host genome sequences: The quality control data is compared with the host reference sequence library to remove host reads from the sequencing data, thereby obtaining the target microorganism sequencing data; Alignment: The sequencing data of the target microorganism is aligned with the microbial reference sequence library to obtain preliminary alignment results; BWA-MEM is used with the -a and -P parameters appended to align the sequencing data of the target microorganism with the microbial reference sequence library to obtain final alignment results; Structural variation candidate read identification: The preliminary alignment results are filtered to retain reads related to potential structural variations, resulting in a set of reads related to structural variations. The reads related to potential structural variations include: soft-cut reads, split alignment reads, abnormal insertion fragment sizes, hard-cut reads, insertion reads, deleted reads, and reads whose mapping positions of paired reads are not on the same chromosome. Chromosome grouping and candidate redundancy removal: Input the final alignment result and the chromosome grouping annotation file, remove highly homologous redundant sequences, retain representative candidate chromosomes, and output a TSV file, including: removing chromosomes with a total number of reads below a threshold; selecting the chromosome with the highest read support number; traversing the remaining chromosomes in descending order of read count and checking whether their reads have been covered by all current candidate chromosomes; Structural variation chromosome pair aggregation and confidence filtering: Input the set of reads related to structural variations and the chromosome grouping annotation file, aggregate the potentially structural variation-related reads in the set of reads related to structural variations into high-confidence chromosome-level or chromosome-group-level structural variation-related read events, calculate the number of reads supported in the final alignment results, and output the SV chromosome pair file; Precise parallel detection of structural variations and insertion sequence tracing: The input includes: the SV chromosome pair file, the raw sequencing data, the microbial reference sequence library, and the final alignment results; after the first processing, a VCF file is output, which includes: the VCF and enhanced VCF for each chromosome pair or detection group; the first processing flow includes: Only reads related to the target chromosome are used for alignment; the alignment results are subjected to SV detection; the inserted sequence is aligned again to the same sub-reference to locate its source region on the original chromosome; Structural variation VCF file merging and deduplication: The input is the VCF file, and the output is the main merged file after global deduplication and a high-quality subset containing only FILTER=PASS; Summary of point mutations and small fragment insertion / deletion candidate sites: The input is the final alignment result and the chromosome grouping annotation file, and the output is a chromosome-level candidate list and a group-level statistical file; High-precision parallel detection of point mutations and insertions / deletions: Its inputs are the chromosome-level candidate list, the raw sequencing data, and the microbial reference sequence library; the output after the second processing is a chromosome-independent GATK VCF; the second processing workflow includes: Only reads relevant to the target chromosome are used for alignment; point mutations and insertions / deletions are detected in the alignment results. Merging and deduplication of point mutation and insertion / deletion VCF files: The input is the GATK VCF file for each chromosome, and the output is the main merged file after global deduplication and a high-quality subset containing only FILTER=PASS.

[0006] Preferably, in chromosome grouping and candidate redundancy removal, the covering rules include: a. The same CIGAR exists in the same read on the candidate chromosome; b. Or the candidate chromosome CIGAR is 150M, which can cover the corresponding reading of any CIGAR; If all reads of a chromosome are covered by existing candidates, they are marked as redundant and removed; otherwise, they are promoted to new candidates.

[0007] Preferably, the process of summarizing and screening the structurally varied chromosome pairs includes: (1) Input parsing Read the set of reads related to the structural variation, and aggregate all aligned chromosomes by read name; If `--groups` is provided, the grouping file will be loaded, and a bidirectional mapping between `chrom → group` and `group → chroms` will be established; if not provided, the chromosome itself will be used as the group by default. (2) Classification of structural variation types For each reading segment: Structural variations within a single chromosome: only aligned to the same chromosome or the same group → denoted as (chr, chr) self-pairing; Cross-chromosomal structural variations: those that align to ≥2 different chromosomes or different groups are recorded as all ordered pairings; (3) Support counting and threshold filtering Count the number of read segments supported for each pair; Filter low-confidence events using `--threshold`, retaining only pairs with a support count greater than or equal to the threshold; Simultaneously output two group-level statistical files: "Before Filtering" and "After Filtering"; (4) Expand into specific chromosome pairs For each pair of groups that are retained: If the number of chromosomes in a group is 1, then output (chr, chr) directly. If the number of chromosomes in a group is ≥2, generate all specific chromosome combinations, along with the source group information; Additional filtering is applied to single-chromosome self-pairing: if the chromosome has already appeared in a cross-chromosome pair, its self-pairing record is deleted.

[0008] Preferably, the process of summarizing and screening the structural variation chromosome pairs further includes: (5) Generation of test groups When `--group-size` is specified: high-confidence pairs are sorted in descending order of support count; every N pairs are merged into one test group, and a unique chromosome list is recorded; the output includes group ID, chromosome pair details, number of unique chromosomes, and total support count.

[0009] Preferably, the structural variation VCF file merging and deduplication process includes a three-layer precise deduplication algorithm and PASS filtering: The three-layer precise deduplication algorithm includes: Hash grouping: Calculate the hash using `{chr}:{pos}:{svtype}:{end}:{alt}` as the key, and group the records into buckets with the same value; Exact matches are compared within the hash bucket: chromosome, svtype, and alt are completely identical; coordinate tolerance ≤ 100 bp; Quality priority: When there is a true duplicate, the highest record of `QUAL` is retained; if the keys are the same but not similar, a unique key with the source file suffix is ​​generated.

[0010] Preferably, the processing flow for summarizing the point mutation and small fragment insertion / deletion candidate sites includes: (1) BAM reading and filtering Do not filter MAPQ=0: retain all alignment records to ensure that no low-quality but real variant signals are missed; Optional primary alignment filtering: Use `--primary_only` to retain only the primary alignment, reducing noise from multiple alignments; Aggregate all chromosomes matched by `read_name` to form a mapping dictionary; (2) Chromosome grouping mapping (3) Group-level support statistics Ungrouped mode: Each chromosome forms its own group; Grouping mode: The chromosome matched by the read is replaced with the corresponding group name, and each read counts each group only once to avoid duplicate counting; (4) Threshold filtering and expansion Expand the groups that meet the threshold back to the specific chromosome: Each chromosome inherits the `support_reads` of its group; Output fields: `chromosome`, `group_name`, `support_reads`, `read_names`.

[0011] Preferably, the processing flow for high-precision parallel detection of point mutations and insertions / deletions includes: sequence extraction, indexing and dictionary, BWA index, alignment-sorting pipeline, BAM index, and GATK HaplotypeCaller.

[0012] Preferably, the TSV file includes: the candidate chromosomes retained in each group and their order, a list of redundant chromosomes that have been covered and removed by the candidate chromosomes, and a new grouping table containing only the candidate chromosomes; The SV chromosome pair file includes: a list of chromosome-level structural variation pairs, statistics on chromosome-group-level structural variations, and a file of the group to be tested.

[0013] Preferably, the microbial reference sequence library includes RefSeq whole genome sequences of known microorganisms; the host reference sequence library includes the human reference genome hg19.

[0014] Preferably, the filtering includes dynamically removing bases with a quality value lower than Q20 or containing linkers from the 5' and 3' ends, with a maximum removal length of 20 bp by default, and retaining reads with a length ≥ 50 bp; the tool for calculating the similarity of each intraspecific reference sequence includes Mash; the tool used for the preliminary alignment includes BWA-MEM.

[0015] Beneficial effects: This invention provides an automated analysis method for microbial metagenomic structural variations (SVs), single nucleotide variants (SNVs / mutations), and small insertions / deletions (Indels). The method includes: calculating and grouping similarities between reference genomes in a microbial dataset; reading sequencing data to be analyzed; performing adapter removal and quality control on the sequencing data; removing host reads from the data; aligning the data with all reference genomes in the microbial dataset based on the data after removing host reads to obtain preliminary and final (all) alignment results; extracting SV-related reads based on the preliminary alignment results; obtaining a subset of SV-related data from the analyzed data based on read names; calculating the number of SV read support for each chromosome set based on the final (all) alignment results and filtering according to a threshold; expanding the filtered chromosome sets to obtain all SV candidate chromosomes and their corresponding read support numbers; performing secondary alignment and SV detection on each group of SV candidate chromosomes; calculating the read support number for each chromosome set based on the final (all) alignment results; expanding chromosome sets with values ​​above the threshold to obtain candidate chromosomes; and performing secondary alignment and single nucleotide variant (SNV / mutation) and small insertion / deletion (Indel) detection on each candidate chromosome. This invention is based on the short-read Illumina next-generation sequencing strategy. By testing and applying current mainstream biological tools and combining them with a self-developed program, it realizes the automated analysis of microbial metagenomic structural variations (SVs), single nucleotide variations (SNVs / mutations), and small fragment insertions and deletions (Indels). It significantly solves the problems of low sensitivity, high false positives, and low efficiency faced by existing technologies in variation detection in the context of metagenomics. Attached Figure Description

[0016] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the embodiments will be briefly described below.

[0017] Figure 1 This is a schematic flowchart of the detection method of the present invention. Detailed Implementation

[0018] This invention provides a method for detecting structural variations, mutations, and insertions / deletions in microorganisms based on next-generation sequencing, comprising the following steps: Calculate the similarity of reference genomes and group them: Download the reference genomes of microorganisms, calculate the similarity of reference sequences within each species according to the biological taxonomy, use 0.95 as the threshold, and group the similarity sequences greater than the threshold as a reference chromosome group to obtain the chromosome grouping annotation file; Raw sequencing data filtering and quality control: The raw sequencing data of the microorganisms to be tested are filtered and quality controlled to obtain quality control data; Removing host genome sequences: The quality control data is compared with the host reference sequence library to remove host reads from the sequencing data, thereby obtaining the target microorganism sequencing data; Alignment: The sequencing data of the target microorganism is aligned with the microbial reference sequence library to obtain preliminary alignment results; BWA-MEM is used with the -a and -P parameters appended to align the sequencing data of the target microorganism with the microbial reference sequence library to obtain final alignment results; Structural variation candidate read identification: The preliminary alignment results are filtered to retain reads related to potential structural variations, resulting in a set of reads related to structural variations. The reads related to potential structural variations include: soft-cut reads, split alignment reads, abnormal insertion fragment sizes, hard-cut reads, insertion reads, deleted reads, and reads whose mapping positions of paired reads are not on the same chromosome. Chromosome grouping and candidate redundancy removal: Input the preliminary alignment results and the chromosome grouping annotation file, remove highly homologous redundant sequences, retain representative candidate chromosomes, and output a TSV file; the TSV file includes: the candidate chromosomes retained in each group and their order, details of redundant chromosomes removed by being covered by candidate chromosomes, and a new grouping table containing only candidate chromosomes; the removal of highly homologous redundant sequences and retention of representative candidate chromosomes includes: (1) removing chromosomes with a total number of reads below a threshold; (2) candidate chromosome screening and redundancy removal at each level: 1) First candidate: select the chromosome with the highest read support number; 2) Subsequent candidate loop: traverse the remaining chromosomes in descending order of read number and check whether their reads have been covered by all current candidate chromosomes; Structural Variation Chromosome Pair Summarization and Confidence Filtering: Input the final (all) alignment results, the set of reads related to structural variations, and the chromosome grouping annotation file. Aggregate the reads related to potential structural variations in the set of reads related to structural variations into high-confidence chromosome-level or chromosome-group-level structural variation-related read events, calculate the corresponding read support number in the final alignment results, and output the SV chromosome pair file. The SV chromosome pair file includes: a list of chromosome-level structural variation pairs, chromosome-group-level structural variation statistics, and a file of the groups to be detected. Precise parallel detection of structural variations and insertion sequence tracing: The input includes: the SV chromosome pair file, the raw sequencing data, the microbial reference sequence library, and the final alignment results; after the first processing, a VCF file is output, which includes: the VCF and enhanced VCF for each chromosome pair or detection group; the first processing flow includes: 1) Target sequence extraction: Extract the specified chromosome into a temporary FASTA file; 2) BWA Index: Creates an index for this sub-reference; 3) BWA-MEM alignment: Alignment is performed using only reads relevant to the target chromosome; 4) Sorting and Indexing: Sort and index the generated alignment results; 5) Manta detection: Perform SV detection on the current alignment results and generate detection results in VCF format; 6) Insertion sequence tracing: Parse the records in the VCF; re-align the inserted sequence to the same sub-reference using BWA to locate its source region on the original chromosome; generate a new record and write it to the additional fields of the VCF; Structural variation VCF file merging and deduplication: The input is the VCF file, and the output is the main merged file after global deduplication and a high-quality subset containing only FILTER=PASS; High-precision parallel detection of point mutations and insertions / deletions: Its inputs are the chromosome-level candidate list, the raw sequencing data, and the microbial reference sequence library; the output after the second processing is a chromosome-independent GATK VCF; the second processing workflow includes: 1) Target sequence extraction: Extract the specified chromosome into a temporary FASTA file; 2) BWA Index: Creates an index for this sub-reference; 3) BWA-MEM alignment: Alignment is performed using only reads relevant to the target chromosome; 4) Sorting and Indexing: Sort and index the generated alignment results; 5) GATK detection: Perform point mutation and insertion / deletion detection on the current alignment results and generate detection results in VCF format; Merging and deduplication of point mutation and insertion / deletion VCF files: The input is the GATK VCF file for each chromosome, and the output is the main merged file after global deduplication and a high-quality subset containing only FILTER=PASS.

[0019] The core of this invention lies in proposing a method for detecting structural variations (SVs), mutations (SNVs), and small insertions / deletions (Indels) in microbial metagenomics based on next-generation sequencing, reference genome similarity grouping, and hierarchical screening. This method significantly solves the problems of low sensitivity, high false positives, and low efficiency faced by existing technologies in variation detection in a metagenomic context by innovatively grouping reference genes, implementing a hierarchical candidate screening mechanism, and optimizing the data processing workflow.

[0020] This invention introduces a genome clustering and grouping mechanism based on sequence similarity into the metagenomic detection process, which can detect microbial metagenomic structural variations (SVs), single nucleotide variations (SNVs / mutations), and small fragment insertions and deletions (Indels) using only second-generation sequencing data, without relying on second-generation sequencing and third-generation sequencing assembly.

[0021] Because there are a large number of similar sequences in microbial reference sequences, existing technologies treat all reference genomes as independent entities, which cannot solve the problem of interference from similar genomes. This invention first performs genome clustering and grouping, and then performs SV detection in a hierarchical manner inside and outside the group, which can more accurately detect microbial structural variations, mutations and insertions and deletions, and improve the calculation speed.

[0022] Existing SV detection methods only identify the insertion time for 50-100 base (bp) length insert fragments originating from other microorganisms, without further determining the source of the fragment. This invention can further determine the microorganism from which the fragment originates.

[0023] This invention first uses the Mash tool to calculate the similarity of all reference genomes in the microbial dataset and groups them to obtain a group list.

[0024] Then, the raw second-generation sequencing data, FastQ1 and FastQ2, are read; adapter removal (using FASTP software) and quality control (using FASTQC software) are performed on the sequencing data; host reads are removed from the sequencing data using Bowtie2 software, generating new Clean_FastQ1 and Clean_FastQ2; based on Clean_FastQ1 and Clean_FastQ2, the BWA tool is used to align with all reference genomes in the microbial dataset to obtain preliminary and final alignment results. Generally, SV detection methods directly use this preliminary result for SV detection. However, because the reference genomes of many microorganisms are very similar, the alignment quality in the results is very low, and the alignment position cannot be accurately identified during SV detection. Therefore, direct SV detection is ineffective. The method provided in this invention can accurately identify the alignment position, and its SV detection performance is superior to existing technologies.

[0025] As one implementation method, the present invention extracts SV-related reads based on preliminary alignment results: regardless of alignment quality, the following conditions are considered SV-related reads: soft-clipped reads, split-alignment reads, abnormal insert size, hard-clipped reads, insert reads, delete reads, and reads whose mapped positions of paired reads are not on the same chromosome.

[0026] In one implementation, the present invention obtains SV-related fq data subsets sv_fastq1 and sv_fastq2 from clean_fastq1 and clean_fastq2 based on the SV-related read names; calculates the number of SV read support for each chromosome pair in the group list based on the major and minor alignment positions of the SV-related reads, and filters them according to a threshold (5 reads); the filtered chromosome pairs are expanded to obtain all SV candidate chromosome pairs.

[0027] As one implementation method, this invention performs secondary alignment and SV detection on each set of SV candidate chromosome pairs. Specifically, the corresponding candidate chromosomes are extracted from the microbial dataset. An alignment index is established using bwa mem2, and alignment is performed using sv_fastq1 and sv_fastq2 to obtain the alignment results. Finally, SV detection is performed using Manta software and a designed Python program, and the SV detection results are merged. The number of reads supported by each chromosome set is calculated based on the final alignment results. Chromosome sets with reads exceeding a threshold are expanded to obtain candidate chromosomes.

[0028] In one implementation method, this invention performs secondary alignment and single nucleotide variant (SNV / mutation) and small fragment insertion / deletion (Indel) detection on each candidate chromosome. Specifically, the corresponding candidate chromosome is extracted from the microbial dataset. An alignment index is built using bwa mem2, and alignments are performed using clean_fastq1 and clean_fastq2 to obtain the alignment results. Then, GATK software is used to detect single nucleotide variant (SNV / mutation) and small fragment insertion / deletion (Indel) detection.

[0029] To further illustrate the present invention, the following detailed description, in conjunction with embodiments and accompanying drawings, describes a method for detecting microbial structural variations, mutations, and insertions / deletions based on next-generation sequencing, but these descriptions should not be construed as limiting the scope of protection of the present invention.

[0030] Example A method for detecting structural variations, mutations, and insertions / deletions in microorganisms based on next-generation sequencing is illustrated in the flowchart below. Figure 1 This includes the following steps: 1. Calculate and group reference genomes: Download reference genomes of microorganisms (bacteria, viruses, fungi, protozoa, and archaea) from the NCBI RefSeq database using the "ncbi-genome-download" tool. Based on biological taxonomy, calculate the similarity of reference sequences within each species using the Mash tool (default parameters). Using 0.95 as a threshold, sequences with similarity greater than the threshold are grouped into a single reference chromosome group, resulting in the chromosome grouping annotation file (all_consolidated_groups.tsv).

[0031] 2. Raw Sequencing Sequence Filtering and Quality Control: The raw sequencing sequences are FASTQ format data generated from paired-end next-generation sequencing, with a single sequence length typically ranging from 75 to 350 bp. The input data consists of a pair of paired-end files (fastq1 and fastq2). This step first uses Fastp software (default parameters) to filter the raw sequencing sequences for low-quality bases and adapter sequences: bases with quality values ​​below Q20 or containing adapters are dynamically removed from the 5' and 3' ends, with a maximum removal length of 20 bp by default, retaining reads ≥50 bp. The filtered sequence set is what this invention refers to as the "preprocessed paired-end sequencing sequence." This sequence set can be all paired-end sequencing data from the sample or a randomly sampled subset, as long as its coverage (coverage depth 5× or higher) is sufficient to detect variations in the target microorganism. The preprocessed paired-end sequencing sequences have removed adapters and low-quality fragments introduced during sequencing and can be directly used for subsequent alignment.

[0032] Then, FastQC software (with default parameters) was used to perform quality control on the preprocessed paired-end sequencing sequences and generate quality control reports in HTML and JSON formats for subsequent data quality assessment.

[0033] 3. Select reference sequence libraries, including a microbial reference sequence library and a host reference sequence library. The microbial reference sequence library is composed of RefSeq whole genome sequences of all known microorganisms (bacteria, viruses, fungi, protozoa, and archaea), with bacteria having complete genomes. An alignment reference database is created using the `bwa index` command. The host reference sequence library is composed of the human reference genome hg19 (GRCh37) and an alignment reference database is created using the `bowtie2-build` command. Both the microbial and host reference sequence libraries can be freely configured according to actual needs.

[0034] 4. Removal of Host Genome Sequences. To focus on the microbial components of the sample, this step aims to align and remove host (human) genome sequences from the sample. The Bowtie2 alignment tool is used to align the quality-controlled sequences from step 2 with the host reference genome sequence library from step 3. The alignment strategy is optimized by setting `--very-sensitive` to high sensitivity mode and utilizing the `--no-mixed`, `--no-discordant`, and `--no-unal` parameters to accurately identify sequences originating from the host genome. The `--un-conc` parameter is used to output unconcordant pairs of sequences that failed to align to the host genome to a new FASTQ file (unaligned human sequences). These unaligned sequences are considered to be non-host nucleic acid sequences from the sample, i.e., the target microbial sequences for subsequent analysis.

[0035] Since a large proportion of the sequencing data from second-generation sequencing comes from the host genome, comparing it with the host genome first can quickly eliminate human contamination and significantly reduce the amount of downstream computation.

[0036] 5. Alignment of "unaligned human sequences" with the microbial sequence library. Using BWA-MEM, the "unaligned human sequences" obtained in step 4 are aligned to the microbial reference sequence library from step 3. A two-level strategy is employed: 1) The default parameter alignment (-M) generates the default alignment result (04_final.sorted.bam), which only retains unique and completely matching read pairs; 2) Add the parameter -aP to perform a "final (all) alignment" to generate results for all possible alignment positions (04_all_final.sorted.bam). This file contains all alignment information, including supplementary alignments, for more in-depth analysis.

[0037] The BAM index was then built using samtools index (default parameters) for fast random access. This comparison confirmed whether the "unaligned human sequence" truly originated from microorganisms and removed any potentially residual unaligned human reads, thereby improving detection accuracy.

[0038] 6. Identification of candidate reads for structural variations: The input is the BAM file (denoted as `04_final.sorted.bam`, which is a custom marker defined in this invention, and the same applies below) generated after alignment to the microbial reference genome by BWA-MEM in step 5. The output is the TSV file (denoted as `06.01_final_sv.tsv`) containing potential structural variation (SV) signals.

[0039] This step is implemented using a Python script, and the core logic is as follows: (1) Input format: The script accepts standard SAM or BAM format files and automatically selects the reading mode according to the extension (.sam is text mode r, .bam is binary mode rb).

[0040] (2) SV determination rules: - Rule 1 (Compact Idiosyncratic Gapped Alignment Report (CIGAR) characteristics): If the CIGAR string of a read contains any of the operators I (insertion), D (deletion), S (soft cut), or H (hard cut), it is considered a potential SV signal; Soft shearing: Aligns with the reference sequence, but of lower quality and not sufficient to be considered a reliable match; Hard shearing: This part is completely removed, usually because these parts are not aligned with the reference sequence or are of too low quality; The difference is that during the alignment process, the soft-cut portion participates in the alignment process and is used to calculate the final alignment score of the read, while the hard-cut portion does not participate in the alignment and is not used to calculate the alignment score of the read. - Rule 2 (SA tag): If the read has a supplementary alignment (SA tag), it indicates that it is a split-read, with different segments aligned to different locations in the genome, suggesting possible translocation, inversion or large fragment deletion; - If any rule is met, it is determined to be an "SV candidate read".

[0041] (3) Filtering mechanism: - First, exclude all unaligned reads with "is_unmapped=True"; - Only retain reads from the primary or supplementary alignment ("is_secondary=False" or "is_supplementary=True"); - Record the read name, aligned chromosome, 1-based start position, CIGAR string, alignment quality MAPQ, whether it is a supplementary alignment, SA tag content and XA tag content (alternative alignment positions).

[0042] Step 6 efficiently filters out conventional linear alignment reads from massive alignment results, retaining only candidate reads carrying SV signals. This provides a highly condensed dataset containing potential structural variation signals for accurate SV typing downstream, thereby improving computational speed.

[0043] 7. Chromosome Grouping and Candidate Redundancy Removal: The input is the BAM file generated in step 5 (`04_all_final.sorted.bam`, generated by "bwa mem -a -P") and the chromosome grouping annotation file (`all_consolidated_groups.tsv`). The output is multiple TSV files, including: - `final_candidate_chromosomes.tsv`: The final candidate chromosomes retained within each group and their order; - `covered_chromosomes_analysis.tsv`: Details of redundant chromosomes removed by covering candidate chromosomes; - `new_chromosome_groups.tsv`: A new grouping table containing only candidate chromosomes; This step is accomplished using a written Python script, and the core process is as follows: (1) Establishment of chromosome grouping mapping Read `all_consolidated_groups.tsv` (no header, the first column is the group name, and subsequent columns are the chromosomes contained in that group), and construct a bidirectional dictionary: - `chr_to_group`: Chromosome → Group Name; - `group_to_chrs`: Group name → Chromosome list.

[0044] (2) BAM parsing and read deduplication Iterate through "04_all_final.sorted.bam": - Distinguish between reads1 and reads2 based on the SAM flag or read-name suffix; - Using chromosomes as units and read_name as the key, only the first occurrence of the record is retained (duplicates are removed), forming a nested dictionary: (3) Low-abundance chromosome filtering Based on the user-defined threshold `--min_reads` (recommended value is 5), chromosomes with a total number of reads (reads1+reads2) lower than the threshold are removed to reduce noise.

[0045] (4) Stepwise screening and redundancy removal of candidate chromosomes Perform independently for each group: 1) First candidate: Select the chromosome with the highest number of reads supported; 2) Subsequent candidate loop: Traverse the remaining chromosomes in descending order of read count, and check whether their reads have been "covered" by all current candidate chromosomes; - Override rules: a. The same CIGAR exists in the same read on the candidate chromosome; b. Alternatively, the candidate chromosome CIGAR is "150M" (representing a perfect full-length alignment), which can cover the corresponding read of any CIGAR; - If all reads of a chromosome are covered by existing candidates, they are marked as redundant and removed; otherwise, they are promoted to new candidates.

[0046] (5) Result output - `final_candidate_chromosomes.tsv`: Group name, candidate sort, chromosome, number of reads (reads1 / reads2), total number of chromosomes in the group, total number of reads in the group, coverage; - `covered_chromosomes_analysis.tsv`: Covered chromosomes, candidate sources of coverage, number of covered reads, number of uncovered reads, percentage of coverage; - `new_chromosome_groups.tsv`: Contains only a simplified grouping table of the final candidate chromosomes; Step 7 allows for the rapid removal of highly homologous redundant sequences across the entire microbial genome, preserving the most representative candidate chromosomes. This process eliminates interference from similar reference chromosomes, providing a concise and highly reliable target set for subsequent structural variations.

[0047] 8. Structural Variation (SV) Chromosome Pair Summarization and Confidence Screening. The input is the SV candidate TSV file generated in step 6 (`06.01_final_sv.tsv`), and the optional chromosome grouping annotation file (`all_consolidated_groups.tsv`). The output is: - Chromosome-level SV pair list (`06.02_final_sv_summary.tsv`); - Chromosome-group level SV statistics (`06.02_final_sv_group_summary.tsv`); - When `--group-size` is enabled, an additional "groups to be detected" file (`*.detection_groups.tsv`) is generated.

[0048] This step is accomplished using a Python script, with the core logic as follows: (1) Input parsing - Read `06.01_final_sv.tsv` and aggregate all aligned chromosomes by `read_name`; - If `--groups` is provided, the grouping file will be loaded, and a bidirectional mapping between `chrom → group` and `group →chroms` will be established; if not provided, the chromosome itself will be used as the group by default.

[0049] (2) SV type classification For each read: - Single chromosome SV: Only aligned to the same chromosome (or the same group) → denoted as (chr, chr) self-pairing; - Cross-chromosome SV: Aligned to ≥2 different chromosomes (or different groups) → recorded as all ordered pairs (duplicates removed after sorting).

[0050] (3) Support counting and threshold filtering - Count the number of reads supported for each pair; - Filter low-confidence events by using `--threshold` (default 1), keeping only pairs with a support count ≥ the threshold; - Simultaneously output two group-level statistical files, "before filtering" and "after filtering" (`*unfiltered.tsv` and `*group_summary.tsv`), for easy quality control.

[0051] (4) Expand into specific chromosome pairs For each pair of groups that are retained: - If the number of chromosomes in the group is 1, then output (chr, chr) directly. - If the number of chromosomes in a group is ≥2, generate all specific chromosome combinations (sorted and deduplicated), and include source group information; - Additional filtering for single-chromosome self-pairing: If the chromosome has already appeared in a cross-chromosome pair, its self-pairing record is deleted to avoid redundancy.

[0052] (5) Generation of test groups (optional) When `--group-size` is specified: - Sort high-confidence pairs in descending order of support count; - Every N pairs are merged into one "test group" and a unique list of chromosomes is recorded; - Output `*.detection_groups.tsv`, which includes group ID, chromosome pair details, number of unique chromosomes, and total number of supported chromosomes, facilitating subsequent parallel validation or experimental verification.

[0053] Step 8 allows for the rapid aggregation of hundreds of thousands of SV candidate reads generated in step 6 into high-confidence chromosome-level or chromosome-group-level SV events, providing a streamlined and parallelizable list of "to-detect groups" for downstream detection and improving computational efficiency.

[0054] 9. Precise parallel detection of structural variations (SVs) and tracing of inserted sequences, with the following input: - The SV chromosome pair file generated in step 8 (`06.02_final_sv_summary.tsv` or `*.detection_groups.tsv`); - Raw double-ended FASTQ(`fq1 / fq2`); - Microbial reference genome; - Pre-aligned BAM (`04_all_final.sorted.bam`, for efficient read extraction).

[0055] The output is: - VCF for each chromosome pair or test group (located in `06_sv / vcf / `). - Optional enhanced VCF (*.enhanced.vcf, with inserted sequence source annotations).

[0056] This step is accomplished using a written Python script, and the core process is as follows: (1) Task splitting and parallel scheduling - Parse input files: Supports two formats – “SV pairing” format (`chromosome_1`, `chromosome_2`); – “To be tested” format (`unique_chromosomes` column, which can contain multiple chromosomes at once).

[0057] - Generate independent tasks based on the chromosome list; process in parallel, with each parallel process further subdivided into: sequence extraction → BWA indexing → alignment and sorting → manta detection.

[0058] (2) Reads preprocessing - One-time extraction: The main process extracts all read names related to the target chromosome in batches from `04_all_final.sorted.bam` (without filtering special cases such as MAPQ=0).

[0059] - Build the raw sequence database: Store the raw FASTQ using an SQLite database, with the read name as the primary key, and save it to shared memory.

[0060] (3) Single task execution chain Perform the following for each chromosome / set: 1) Target sequence extraction: Extract the specified chromosome into a temporary FASTA file; 2) BWA Index: Creates an index for this sub-reference; 3) BWA-MEM alignment: Alignment is performed using only reads (a subset in memory) relevant to the target chromosome; 4) Sorting and Indexing: Sort and index the generated alignment results; 5) Manta detection: Use the manta tool (default parameters) to perform SV detection on the current alignment results and generate detection results in VCF format; 6) Insertion sequence tracing: – Parse records in VCF; – Align the inserted sequence again to the same sub-reference using BWA to locate its source region on the original chromosome; – Generate a new record and write it to the VCF, with additional fields.

[0061] By employing the parallel and group detection strategies in step 9, high-precision SV detection results can be obtained, and the internal chromosomal origin of the inserted fragment can be automatically labeled, thereby improving computational efficiency.

[0062] 10. Merging and deduplicating Structural Variation (SV) VCF files, the input is: - All chromosome / detection group VCF files generated in step 9; The output is: - The main merged file `06_sv / merged_final.vcf` after global deduplication; - A high-quality subset `06_sv / merged_final_pass.vcf` containing only FILTER=PASS; This step is accomplished using a written Python script, and the core process is as follows: (1) High-performance file scanning and parallel parsing - Zero-copy traversal: using `mmap` memory mapping + batch reading; - Multi-process parsing: Grouped by `--batch_size` (default 100 files) for parallel execution; (2) Three-layer precise deduplication algorithm Hash grouping: Calculate the hash using `{chr}:{pos}:{svtype}:{end}:{alt}` as the key, and group the records into buckets with the same value; Exact matches are compared within the hash bucket: chromosome, svtype, and alt are completely identical; coordinate tolerance ≤ 100 bp; Quality priority: When there is a true duplicate, the highest record with `QUAL` is retained; if the keys are the same but not similar, a unique key with the source file extension is generated. (3) PASS filtering - FILTER field check: Only keep rows where FILTER == "PASS"; Step 10 yields a high-quality, non-redundant SV complete set.

[0063] 11. Summary of candidate sites for point mutations (SNVs) and small fragment insertions / deletions (indels), the input is: - The microbial alignment BAM file generated in step 5 (`04_final.sorted.bam`, which is sorted and indexed); - Chromosome grouping annotation file (`all_consolidated_groups.tsv`); The output is: - Chromosomal-level mutation / indel candidate list (`07.02_final_mutation_indel_summary.tsv`); - Group-level statistics files (*_all_groups.tsv and *_filtered_groups.tsv).

[0064] This step is accomplished using a written Python script, and the core process is as follows: (1) BAM reading and filtering - Iterate through `04_final.sorted.bam`; - Do not filter MAPQ=0: retain all alignment records to ensure that no low-quality but real variant signals are missed; - Optional primary alignment filtering: Only primary alignments are retained using `--primary_only`, reducing noise from multiple alignments; - Aggregate all chromosomes matched by `read_name` to form a mapping dictionary.

[0065] (2) Chromosome grouping mapping (optional) - If `--group_file` is provided, read headerless TSVs: – The first column is the group name; – The chromosomes of the members of this group are listed below.

[0066] - Construct a bidirectional mapping: – `chromosome_to_group`: Chromosome → Group Name; – `group_to_chromosomes`: group name → set of chromosomes.

[0067] (3) Group-level support statistics - No grouping mode: Each chromosome forms its own group; - Grouping mode: Replace the chromosome matched by the read with the corresponding group name, and count each read for each group only once to avoid duplicate counting; - Generate `group_stats`: – `support_reads`: The unique number of reads that cover this group; – `chromosomes`: The actual set of chromosomes involved in the group.

[0068] (4) Threshold filtering and expansion - Filter high-confidence groups using `--min_support` (default 1, recommended 5); - Expand the groups that meet the threshold back to the specific chromosome: – Each chromosome inherits the `support_reads` of its group; – Output fields: `chromosome`, `group_name`, `support_reads`, `read_names`.

[0069] The output results can be directly used as a chromosome whitelist for step 11 (GATK mutation detection); group-level statistics help determine which chromosomes / regions have high-coverage mutation signals and prioritize in-depth analysis. This invention detects point mutations and small fragment insertions / deletions on a chromosome-set basis, which can improve detection accuracy.

[0070] 12. High-precision parallel detection of point mutations (SNVs) and insertions / deletions (indels), with the following input: - Chromosome-level candidate list generated in step 11 (`07.02_final_mutation_indel_summary.tsv`); - Raw double-ended FASTQ(`fq1 / fq2`); - Microbial Reference Genome FASTA; The output is: a separate GATK VCF for each chromosome (`07_mutation_indel / *.gatk.vcf`); This step is accomplished using a Python script, and the core process is as follows (all software tools used are configured with default parameters): Single-chromosome workflow (each task is executed serially). 1. Sequence Extraction – `samtools faidx ref.fa chromosome>tmp / {chr} / {chr}.fa`; 2. Indexes and Dictionaries – `samtools faidx {chr}.fa` (index); – `gatk CreateSequenceDictionary -R {chr}.fa -O {chr}.dict` (dictionary); 3. BWA Index – `bwa index {chr}.fa`; 4. Comparison-sorting pipeline – `bwa mem | samtools view -Sb | samtools sort -o {chr}.sorted.bam`; 5. BAM Index – `samtools index {chr}.sorted.bam`; 6. GATK HaplotypeCaller – `gatk --java-options "-Xmx10g" HaplotypeCaller -R {chr}.fa -I{chr}.sorted.bam -O {chr}.gatk.vcf`; Step 12 enables precise SNV / indel detection of tens of millions of microbial sequencing reads within hours on a multi-core server, outputting a standard VCF file. Step 12 of this invention, through parallel and group detection strategies, can accurately detect microbial mutations and insertions / deletions while improving computational efficiency.

[0071] 13. Merging and deduplicating point mutation (SNV) and insertion / deletion (indel) VCF files, the input is: - All chromosome / detection group VCF files generated in step 12; The output is: - The main merged file after global deduplication: `07_mutation_indel / merged_mutation_indel.vcf`; - A high-quality subset `07_mutation_indel / merged_mutation_indel_pass.vcf` containing only FILTER=PASS; This step is accomplished using a written Python script, and the core process is as follows: (1) High-performance file scanning and parallel parsing - Multi-process parsing: Grouped by `--batch_size` (default 100 files) for parallel execution; (2) Three-layer precise deduplication algorithm Hash grouping: Calculate the hash using `{chr}:{pos}:{ref}:{alt}` as the key, and divide the records into buckets with the same value; Exact matching is performed within the hash bucket: chromosome, ref, and alt are completely identical; Quality priority: When there is a true duplicate, the highest record with `QUAL` is retained; if the keys are the same but not similar, a unique key with the source file extension is generated. (3) PASS filtering - FILTER field check: Only keep rows where FILTER == "PASS"; Step 13 yields a complete set of point mutations (SNVs) and insertions / deletions (indels) without redundancy and of high quality.

[0072] The alignment tool BWA involved in step 9 can be replaced with other alignment tools, such as Bowtie, minimap2, etc., and the SV detection tool can be replaced with other detection tools, such as Manta, Lumpy, GRIDSS, etc.

[0073] result: This invention simulates five types of structural variation (SV): deletion (DEL), insertion (INS), duplication (DUP), inversion (INV), and translocation (TRA), with four simulated test data for each type. The simulated data and real sequencing data (accessed in the NCBI SRA database as SRR32394700) were then mixed to simulate real-world sequencing data, achieving a final overall accuracy of 95%. The results show that this invention can detect the corresponding structural variation sites (Table 1).

[0074] Table 1 SV Detection Results

[0075] This invention simultaneously simulates single nucleotide variants (SNVs / mutations) and small fragment insertions / deletions (Indels), and mixes simulated data with real sequencing data (SRR32394700) to simulate sequencing data in a real-world scenario, achieving a final overall accuracy of 97.71%. The results demonstrate that this invention can detect both the corresponding single nucleotide variants and small fragment insertion / deletion sites.

[0076] Table 2 SNV and Indel detection results

[0077] In summary, this invention, through innovative grouping of reference genes, implementation of a stratified candidate screening mechanism, and optimization of data processing procedures, significantly solves the problems of low sensitivity, high false positives, and low efficiency faced by existing technologies in variation detection in a metagenomic context.

[0078] Although the above embodiments have provided a detailed description of the present invention, they are only some embodiments of the present invention, and not all embodiments. People can obtain other embodiments based on these embodiments without creative effort, and these embodiments all fall within the protection scope of the present invention.

Claims

1. A method for detecting structural variations, mutations, and insertions / deletions in microorganisms based on next-generation sequencing, characterized in that, Includes the following steps: Calculate the similarity of reference genomes and group them: According to the biological classification method, calculate the similarity of reference sequences within each species, and use 0.95 as the threshold. Those with similarity greater than the threshold are grouped into a reference chromosome group, and the chromosome grouping annotation file is obtained. Raw sequencing data filtering and quality control: The raw sequencing data of the microorganisms to be tested are filtered and quality controlled to obtain quality control data; Removing host genome sequences: The quality control data is compared with the host reference sequence library to remove host reads from the sequencing data, thereby obtaining the target microorganism sequencing data; Alignment: The sequencing data of the target microorganism is aligned with the microbial reference sequence library to obtain preliminary alignment results; BWA-MEM is used with the -a and -P parameters appended to align the sequencing data of the target microorganism with the microbial reference sequence library to obtain final alignment results; Structural variation candidate read segment identification: The preliminary comparison results are filtered to retain read segments related to potential structural variations, thus obtaining a set of read segments related to structural variations; The reads associated with potential structural variations include: soft-splitting reads, split alignment reads, abnormal insertion size, hard-splitting reads, insertion reads, deletion reads, and reads whose mapping positions of paired reads are not on the same chromosome. Chromosome grouping and candidate redundancy removal: Input the final alignment result and the chromosome grouping annotation file, remove highly homologous redundant sequences, retain representative candidate chromosomes, and output a TSV file, including: removing chromosomes with a total number of reads below a threshold; selecting the chromosome with the highest read support number; traversing the remaining chromosomes in descending order of read count and checking whether their reads have been covered by all current candidate chromosomes; Structural variation chromosome pair aggregation and confidence filtering: Input the set of reads related to structural variations and the chromosome grouping annotation file, aggregate the potentially structural variation-related reads in the set of reads related to structural variations into high-confidence chromosome-level or chromosome-group-level structural variation-related read events, calculate the number of reads supported in the final alignment results, and output the SV chromosome pair file; Precise parallel detection of structural variations and insertion sequence tracing: The input includes: the SV chromosome pair file, the raw sequencing data, the microbial reference sequence library, and the final alignment results; after the first processing, a VCF file is output, which includes: the VCF and enhanced VCF for each chromosome pair or detection group; the first processing flow includes: Only reads related to the target chromosome are used for alignment; the alignment results are subjected to SV detection; the inserted sequence is aligned again to the same sub-reference to locate its source region on the original chromosome; Structural variation VCF file merging and deduplication: The input is the VCF file, and the output is the main merged file after global deduplication and a high-quality subset containing only FILTER=PASS; Summary of point mutations and small fragment insertion / deletion candidate sites: The input is the final alignment result and the chromosome grouping annotation file, and the output is a chromosome-level candidate list and a group-level statistical file; High-precision parallel detection of point mutations and insertions / deletions: Its inputs are the chromosome-level candidate list, the raw sequencing data, and the microbial reference sequence library; the output after the second processing is a chromosome-independent GATK VCF; the second processing workflow includes: Only reads relevant to the target chromosome are used for alignment; point mutations and insertions / deletions are detected in the alignment results. Merging and deduplication of point mutation and insertion / deletion VCF files: The input is the GATK VCF file for each chromosome, and the output is the main merged file after global deduplication and a high-quality subset containing only FILTER=PASS.

2. The method according to claim 1, characterized in that, In chromosome grouping and candidate redundancy removal, the covering rules include: a. The same CIGAR exists in the same read on the candidate chromosome; b. Or the candidate chromosome CIGAR is 150M, which can cover the corresponding reading of any CIGAR; If all reads of a chromosome are covered by existing candidates, they are marked as redundant and removed; otherwise, they are promoted to new candidates.

3. The method according to claim 1, characterized in that, The process of summarizing and screening the structural variation chromosome pairs includes: (1) Input parsing Read the set of reads related to the structural variation, and aggregate all aligned chromosomes by read name; If `--groups` is provided, the grouping file will be loaded, and a bidirectional mapping between `chrom → group` and `group → chroms` will be established; if not provided, the chromosome itself will be used as the group by default. (2) Classification of structural variation types For each reading segment: Structural variations within a single chromosome: only aligned to the same chromosome or the same group → denoted as (chr, chr) self-pairing; Cross-chromosomal structural variations: those that align to ≥2 different chromosomes or different groups are recorded as all ordered pairings; (3) Support counting and threshold filtering Count the number of read segments supported for each pair; Filter low-confidence events using `--threshold`, retaining only pairs with a support count greater than or equal to the threshold; Simultaneously output two group-level statistical files: "Before Filtering" and "After Filtering"; (4) Expand into specific chromosome pairs For each pair of groups that are retained: If the number of chromosomes in a group is 1, then output (chr, chr) directly. If the number of chromosomes in a group is ≥2, generate all specific chromosome combinations, along with the source group information; Additional filtering is applied to single-chromosome self-pairing: if the chromosome has already appeared in a cross-chromosome pair, its self-pairing record is deleted.

4. The method according to claim 3, characterized in that, The process of summarizing and screening the structural variation chromosome pairs with confidence levels also includes: (5) Generation of test groups When `--group-size` is specified: high-confidence pairs are sorted in descending order of support count; every N pairs are merged into one test group, and a unique chromosome list is recorded; the output includes group ID, chromosome pair details, number of unique chromosomes, and total support count.

5. The method according to claim 1, characterized in that, The structural variation VCF file merging and deduplication process includes a three-layer accurate deduplication algorithm and PASS filtering: The three-layer precise deduplication algorithm includes: Hash grouping: Calculate the hash using `{chr}:{pos}:{svtype}:{end}:{alt}` as the key, and group the records into buckets with the same value; Exact matches are compared within the hash bucket: chromosome, svtype, and alt are completely identical; coordinate tolerance ≤ 100 bp; Quality priority: When there is a true duplicate, the highest record of `QUAL` is retained; if the keys are the same but not similar, a unique key with the source file suffix is ​​generated.

6. The method according to claim 1, characterized in that, The processing flow for summarizing point mutations and small fragment insertion / deletion candidate sites includes: (1) BAM reading and filtering Do not filter MAPQ=0: retain all alignment records to ensure that no low-quality but real variant signals are missed; Optional primary alignment filtering: Use `--primary_only` to retain only the primary alignment, reducing noise from multiple alignments; Aggregate all chromosomes matched by `read_name` to form a mapping dictionary; (2) Chromosome grouping mapping (3) Group-level support statistics Ungrouped mode: Each chromosome forms its own group; Grouping mode: The chromosome matched by the read is replaced with the corresponding group name, and each read counts each group only once to avoid duplicate counting; (4) Threshold filtering and expansion Expand the groups that meet the threshold back to specific chromosomes: Each chromosome inherits the `support_reads` of its group; Output fields: `chromosome`, `group_name`, `support_reads`, `read_names`.

7. The method according to claim 1, characterized in that, The high-precision parallel detection process for point mutations and insertions / deletions includes: sequence extraction, indexing and dictionary, BWA index, alignment-sorting pipeline, BAM index, and GATKHaplotypeCaller.

8. The method according to claim 1, characterized in that, The TSV file includes: the candidate chromosomes retained in each group and their order, a list of redundant chromosomes that have been removed and covered by the candidate chromosomes, and a new grouping table containing only the candidate chromosomes; The SV chromosome pair file includes: a list of chromosome-level structural variation pairs, statistics on chromosome-group-level structural variations, and a file of the group to be tested.

9. The method according to claim 1, characterized in that, The microbial reference sequence library includes the RefSeq whole genome sequences of known microorganisms; the host reference sequence library includes the human reference genome hg19.

10. The method according to claim 1, characterized in that, The filtering includes dynamically removing bases with a quality value lower than Q20 or containing linkers from the 5' and 3' ends, with a default maximum removal length of 20 bp, and retaining reads with a length ≥50 bp; the tool for calculating the similarity of each intraspecific reference sequence includes Mash; the tool used for the preliminary alignment includes BWA-MEM.

Images