A method for detecting microbial structural variations, mutations and indels based on next generation sequencing

Abstract

The present application relates to the technical field of bioinformatics, in particular to a method for detecting microbial structural variation, mutation and insertion deletion based on second-generation sequencing. The present application is based on the short read Illumina second-generation sequencing strategy, through testing, using the current mainstream biological tools and combining the self-developed program, the automatic analysis of microbial metagenome structural variation (SV), single nucleotide variation (SNV / mutation) and small fragment insertion deletion (Indel) is realized, and the problems of low sensitivity, high false positive rate and low efficiency in the prior art for detecting variation under the metagenome background are significantly solved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of bioinformatics, and in particular to a method for detecting microbial structural variation, mutation and insertion deletion based on second-generation sequencing. BACKGROUND

[0002] Microbial structural variation (SV), mutation (SNV) and small fragment insertion deletion (Indel) not only affect microbial genome molecules and cellular processes, but also affect their regulatory functions and gene expression, which may further cause cells from the same species to exhibit complex phenotypes and different growth rates.

[0003] The identification of microbial SV usually requires the use of second-generation sequencing and third-generation sequencing for assembly and then macro-genome analysis, and the existing technology regards all reference genomes as independent entities, which cannot solve the problem of similar genome interference, resulting in the problems of low sensitivity, high false positive rate and low efficiency of the existing technology in the detection of variations in the macro-genome background. SUMMARY

[0004] In order to solve the above problems, the present application provides a method for detecting microbial structural variation, mutation and insertion deletion based on second-generation sequencing. The present application innovatively groups reference genes, implements a hierarchical candidate screening mechanism, and optimizes the data processing process, which significantly solves the problems of low sensitivity, high false positive rate and low efficiency of the existing technology in the detection of variations in the macro-genome background.

[0005] In order to achieve the above purpose, the present application provides the following technical scheme:

[0006] The present application provides a method for detecting microbial structural variation, mutation and insertion deletion based on second-generation sequencing, comprising the following steps:

[0007] Calculate the similarity of the reference genome and group it: according to biological taxonomy, calculate the similarity of each intra-species reference sequence, take 0.95 as the threshold, and group the reference chromosomes with a similarity greater than the threshold as a reference chromosome set to obtain a chromosome grouping annotation file;

[0008] Filter and quality control of raw sequencing sequences: filter and quality control the raw sequencing data of the microorganism to be detected to obtain quality control data;

[0009] Eliminate host genome sequences: align the quality control data with the host reference sequence library to remove host reads in the sequencing data to obtain target microbial sequencing data;

[0010] Alignment: align the target microorganism sequencing data with the microorganism reference sequence library to obtain a preliminary alignment result; align the target microorganism sequencing data with the microorganism reference sequence library using BWA-MEM and adding -a -P parameters to obtain a final alignment result;

[0011] Structural variation candidate read identification: screen the preliminary alignment result to retain potential structural variation related reads to obtain a structural variation related read set; the potential structural variation related reads include: soft clipping reads, split alignment reads, abnormal insert fragment size, hard clipping reads, insert reads, deletion reads, and reads whose mapping positions of paired reads are not on the same chromosome;

[0012] Chromosome grouping and candidate de-duplication: input the final alignment result and the chromosome grouping annotation file, remove high homologous redundant sequences, retain representative candidate chromosomes, and output a TSV file including: remove chromosomes with total read counts below a threshold; select the chromosome with the highest read support; traverse the remaining chromosomes in descending order of read count, and check whether the reads have been covered by all current candidate chromosomes;

[0013] Structural variation chromosome pair summary and confidence screening: input the structural variation related read set and the chromosome grouping annotation file, aggregate potential structural variation related reads in the structural variation related read set into high-confidence chromosome-level or chromosome-group-level structural variation related read events, and calculate the corresponding read support numbers in the final alignment result, output the SV chromosome pair file;

[0014] Accurate parallel detection of structural variations and insertion sequence tracing: the input includes: the SV chromosome pair file, the original sequencing data, the microorganism reference sequence library, and the final alignment result; the first processing outputs a VCF file, which includes: VCF and enhanced VCF for each chromosome pair or detection group; the first processing process includes:

[0015] Align only using reads related to the target chromosome; perform SV detection on the alignment result; align the insertion sequence to the same sub-reference again to locate its source interval on the original chromosome;

[0016] Structural variation VCF file merging and de-duplication: the input is the VCF file, and the output is the main merged file after global de-duplication and the high-quality subset containing only FILTER=PASS;

[0017] Point mutation and small fragment insertion and deletion candidate site summary: input the final alignment result and the chromosome grouping annotation file, output the chromosome-level candidate list and group-level statistical file;

[0018] High-precision parallel detection of point mutations and indels: its input is the chromosome-level candidate list, the original sequencing data and the microbial reference sequence library, and the output after the second processing is the GATK VCF of each chromosome independently; the flow of the second processing includes:

[0019] Only use the reads related to the target chromosome for alignment; perform point mutation and indel detection on the alignment results;

[0020] Merge and deduplicate point mutation and indel VCF files: input is the GATK VCF of each chromosome independently, and output is the main merged file after global deduplication and the high-quality subset containing only FILTER=PASS.

[0021] Preferably, in the chromosome grouping and candidate de-duplication, the coverage rules include:

[0022] a. The same read has the same CIGAR in the candidate chromosome;

[0023] b. Or the candidate chromosome CIGAR is 150M, and the corresponding read of any CIGAR can be covered;

[0024] If all reads of a chromosome are covered by existing candidates, it is marked as redundant and removed; otherwise, it is upgraded to a new candidate.

[0025] Preferably, the process of structure variant chromosome pair summary and confidence screening includes:

[0026] (1) Input analysis

[0027] Read the structure variant related read set, and aggregate all aligned chromosomes according to the read name;

[0028] If `--groups` is provided, load the grouping file and establish `chrom → group` and `group → chroms` bidirectional mapping; if not provided, the chromosome itself is used as the group by default;

[0029] (2) Structure variant type division

[0030] For each read:

[0031] Single-chromosome structure variant: only aligned to the same chromosome or the same group → recorded as (chr, chr) self-pairing;

[0032] Chromosome-spanning structure variant: aligned to ≥2 different chromosomes or different groups → recorded as all ordered group pairs;

[0033] (3) Support count and threshold filtering

[0034] Count reads support for each group pair

[0035] Filter low confidence events by `--threshold`, only keep group pairs with support ≥ threshold;

[0036] Output both "before filtering" and "after filtering" group-level statistics files;

[0037] (4) Expand to specific chromosome pairs

[0038] For each group pair kept:

[0039] If the number of chromosomes in the group is 1, directly output (chr, chr);

[0040] If the number of chromosomes in the group is ≥ 2, generate all specific chromosome combinations and attach the source group information;

[0041] Additional filtering for single-chromosome self-pairing: if the chromosome has appeared in a cross-chromosome pair, delete its self-pairing record.

[0042] Preferably, the process of structural variation chromosome pair aggregation and confidence filtering further comprises:

[0043] (5) Group generation for detection

[0044] When `--group-size` is specified: arrange high-confidence group pairs in descending order of support; merge every N group pairs into one group for detection, and record the unique chromosome list; output the results including group ID, chromosome pair details, unique chromosome number, and total support.

[0045] Preferably, the process of structural variation VCF file merging and deduplication comprises a three-layer accurate deduplication algorithm and PASS filtering:

[0046] The three-layer accurate deduplication algorithm comprises:

[0047] Hash grouping: calculate hash with `{chr}:{pos}:{svtype}:{end}:{alt}` as the key, and divide records into the same value bucket;

[0048] Accurate matching within the hash bucket: chromosome, svtype, and alt are completely identical; coordinate tolerance ≤ 100 bp;

[0049] Quality priority: keep the record with the highest `QUAL` when there is a real duplicate; generate a unique key with a source file suffix when the keys are the same but not similar.

[0050] Preferably, the process of point mutation and small fragment insertion / deletion candidate site aggregation comprises:

[0051] (1) BAM reading and filtering

[0052] No filtering MAPQ=0: keep all alignment records, make sure no low-quality but real variant signal is missed;

[0053] Optional primary alignment filtering: keep only primary alignments by `--primary_only`, reduce multi-alignment noise;

[0054] Aggregate all chromosomes a read is aligned to by `read_name`, form a mapping dictionary;

[0055] (2) Chromosome grouping mapping

[0056] (3) Group-level support counting

[0057] No grouping mode: each chromosome forms a group by itself;

[0058] Grouping mode: replace the chromosomes a read is aligned to with the corresponding group name, count each read only once for each group, avoid duplicate counting;

[0059] (4) Threshold filtering and unfolding

[0060] Unfold groups that meet the threshold back to specific chromosomes:

[0061] Each chromosome inherits the `support_reads` of the group it belongs to;

[0062] Output fields: `chromosome`, `group_name`, `support_reads`, `read_names`.

[0063] Preferably, the processing procedure for high-precision parallel detection of point mutations and indels includes sequence extraction, indexing and dictionary, BWA indexing, alignment-sorting pipeline, BAM indexing, and GATK HaplotypeCaller.

[0064] Preferably, the TSV file includes candidate chromosomes retained in each group and their order, details of redundant chromosomes excluded by candidate chromosomes, and a new group table containing only candidate chromosomes.

[0065] The SV chromosome pair file includes a list of chromosome-level structural variation pairs, chromosome-group-level structural variation statistics, and a group file to be detected.

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

[0067] Preferably, the filtering comprises dynamically removing bases with quality values below Q20 or adapters from both 5' and 3' ends, with a maximum removal length of 20 bp by default, and retaining read lengths with lengths ≥ 50 bp; the tool used for calculating the similarity of each intra-species reference sequence comprises Mash; and the tool used for the preliminary alignment comprises BWA-MEM.

[0068] Advantages:

[0069] The present application provides a kind of microbial macro genome structural variation (SV), single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel) automatic analysis method, comprising: calculating the similarity of reference genome in microbial data set and grouping;Read the sequencing data to be analyzed;The data to be analyzed is de-ligated and quality controlled;Remove host reads in data;Based on the data of removing host reads, all reference genomes in microbial data set are aligned to obtain preliminary alignment result and final (all) alignment result;According to preliminary alignment result, extract SV related reads;According to reads name, obtain SV related data subset from analysis data;According to final (all) alignment result, calculate the number of reads supporting chromosome SV, and filter according to threshold value;After filtering, chromosome is unfolded to obtain all SV candidate chromosomes and corresponding reads support number;Each group of SV candidate chromosomes is subjected to secondary alignment and SV detection;According to final (all) alignment result, the number of reads supporting each chromosome is calculated;The chromosome group higher than threshold value is unfolded to obtain candidate chromosome;Each candidate chromosome is subjected to secondary alignment and single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel) detection.The present application is based on short read Illumina second generation sequencing strategy, by testing, using current mainstream biological tools and combining self-developed program, realizes microbial macro genome structural variation (SV), single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel) automatic analysis, significantly solves the problem of low sensitivity, high false positive and low efficiency when the prior art detects variation in the background of macro genome. BRIEF DESCRIPTION OF DRAWINGS

[0070] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced as follows.

[0071] Figure 1 The flowchart of the detection method of the present application is shown. DETAILED DESCRIPTION

[0072] The present application provides a method for detecting microbial structural variation, mutation and insertion and deletion based on second generation sequencing, comprising the following steps:

[0073] Similarity of reference genomes is calculated and grouped: download reference genomes of microorganisms, calculate the similarity of each reference sequence within a species according to biological taxonomy, take 0.95 as the threshold, and take the reference chromosome set with similarity greater than the threshold as a reference chromosome set, to obtain a chromosome grouping annotation file;

[0074] Filtering and quality control of raw sequencing sequences: filter and quality control the raw sequencing data of the microorganism to be tested to obtain quality control data;

[0075] Eliminate host genome sequences: align the quality control data with the host reference sequence library to remove host reads in the sequencing data to obtain target microorganism sequencing data;

[0076] Alignment: align the target microorganism sequencing data with the microorganism reference sequence library to obtain a preliminary alignment result; use BWA-MEM and add -a -P parameters to align the target microorganism sequencing data with the microorganism reference sequence library to obtain a final alignment result;

[0077] Structural variation candidate read identification: screen the preliminary alignment result to retain potential structural variation related reads to obtain a structural variation related read set; the potential structural variation related reads include: soft cutting reads, split alignment reads, abnormal insertion fragment size, hard cutting reads, insertion reads, deletion reads, and reads whose mapping positions of paired reads are not on the same chromosome;

[0078] Chromosome grouping and candidate redundancy removal: input the preliminary alignment result and the chromosome grouping annotation file, remove high homologous redundant sequences, and retain representative candidate chromosomes, and output a TSV file; the TSV file includes: each retained candidate chromosome in the group and its order, details of redundant chromosomes removed by the candidate chromosomes, and a new group table containing only candidate chromosomes; the removal of high homologous redundant sequences and the retention of representative candidate chromosomes include: (1) remove chromosomes with total read counts below a threshold; (2) candidate chromosome step-by-step screening and redundancy removal: 1) first candidate: select the chromosome with the highest read support; 2) subsequent candidate cycle: traverse the remaining chromosomes in descending order of read counts, and check whether their reads have been covered by all candidate chromosomes;

[0079] Structure variant chromosome pair filtering: input the final (all) alignment result, the structure variant related reads set and the chromosome grouping annotation file, aggregate the reads in the structure variant related reads set that are potentially related to a structure variant into high-confidence chromosome-level or chromosome-group-level structure variant related reads events, and calculate the corresponding reads support number in the final alignment result, output the SV chromosome pair file; the SV chromosome pair file includes: a chromosome-level structure variant pair list, a chromosome-group-level structure variant statistics and a to-be-detected group file;

[0080] Structure variant accurate parallel detection and insertion sequence tracing: the input content includes: the SV chromosome pair file, the original sequencing data, the microbial reference sequence library and the final alignment result; the first processing outputs a VCF file, which includes: a VCF and an enhanced VCF of each chromosome pair or detection group; the first processing process includes:

[0081] 1) Target sequence extraction: extract the specified chromosome to a temporary FASTA file;

[0082] 2) BWA indexing: index the sub-reference;

[0083] 3) BWA-MEM alignment: only use the reads related to the target chromosome for alignment;

[0084] 4) Sorting and indexing: sort and index the generated alignment result;

[0085] 5) Manta detection: perform SV detection on the current alignment result to generate a VCF format detection result;

[0086] 6) Insertion sequence tracing: parse the records in the VCF; re-BWA align the insertion sequence to the same sub-reference to locate its source interval on the original chromosome; generate new records and write them to the additional fields of the VCF;

[0087] Structure variant VCF file merging and deduplication: the input is the VCF file, and the output is a globally deduplicated main merged file and a high-quality subset containing only FILTER=PASS;

[0088] Point mutation and insertion-deletion high-precision parallel detection: its input is the chromosome-level candidate list, the original sequencing data and the microbial reference sequence library, and the second processing outputs a GATK VCF independent of each chromosome; the second processing process includes:

[0089] 1) Target sequence extraction: extract the specified chromosome to a temporary FASTA file;

[0090] 2) BWA index: index the sub-reference;

[0091] 3) BWA-MEM alignment: alignment using only reads associated with the target chromosome;

[0092] 4) Sorting and indexing: sorting and indexing the generated alignment results;

[0093] 5) GATK detection: point mutation and insertion-deletion detection on the current alignment results, generating VCF format detection results;

[0094] Point mutation and insertion-deletion VCF file merging and deduplication: input is the GATK VCF of each chromosome independently, and output is the main merged file after global deduplication and high-quality subset containing only FILTER=PASS.

[0095] The core of the present application is to propose a microbial macrogenome-oriented, second-generation sequencing-based, reference genome similarity grouping and hierarchical screening of structural variation (SV), mutation (SNV) and small fragment insertion-deletion (Indel) detection method. The method innovatively groups the reference genes, implements hierarchical candidate screening mechanism, and optimizes the data processing process, which significantly solves the problems of low sensitivity, high false positive rate and low efficiency of the prior art in the context of macrogenome-based variation detection.

[0096] The present application introduces a sequence similarity-based genome clustering grouping mechanism in the macrogenome detection process, which can detect microbial macrogenome structural variation (SV), single nucleotide variation (SNV / mutation) and small fragment insertion-deletion (Indel) using only second-generation sequencing data, without relying on second-generation sequencing and third-generation sequencing assembly.

[0097] Due to the existence of a large number of similar sequences in microbial reference sequences, the prior art regards all reference genomes as independent entities, which cannot solve the problem of similar genome interference. The present application first groups the genomes, and then processes the SV detection within and outside the group in a hierarchical manner, which can more accurately detect microbial structural variation, mutation and insertion-deletion, and improve the calculation speed.

[0098] The existing SV detection method only judges the insertion fragment of 50-100 base pairs (bp) length derived from other microorganisms as an insertion time, without further judging the source of the fragment. The present application can further judge the source of the fragment.

[0099] The present application first calculates the similarity of all reference genomes in the microbial data set using the Mash tool and groups them to obtain a grouping list (group list).

[0100] Then, read the original second-generation sequencing data: fastq1 and fastq2; perform de-adapter (fastp software) and quality control (fastqc software) on the sequencing data; use bowtie2 software to remove host reads in the sequencing data to generate new clean_ fastq1 and clean_ fastq2; based on clean_ fastq1 and clean_ fastq2, use the bwa tool to align all reference genomes in the microbial data set to obtain preliminary alignment results and final alignment results. Generally, the SV detection method is to directly detect SV based on the preliminary results, but because the reference genomes of many microorganisms are very similar, the alignment quality in the alignment results is low, and the alignment position cannot be accurately recognized in the SV detection process, so direct detection of SV is poor, and the method provided by the application can accurately recognize the alignment position, and the effect of detecting SV is better than that of the prior art.

[0101] As an embodiment, the application extracts SV-related reads according to the preliminary alignment results: without considering the alignment quality, if the following conditions exist, it is considered that the reads are SV-related reads: soft-clipped reads, split-alignment reads, abnormal insert size, hard-clipped reads, insert reads, delete reads, and reads whose mapping positions of paired reads are not on the same genome.

[0102] As an embodiment, the application obtains SV-related fq data subsets sv_ fastq1 and sv_ fastq2 from clean_ fastq1 and clean_ fastq2 according to the names of SV-related reads; according to the primary alignment position and the secondary alignment position of the SV-related reads, the SV read support number of each chromosome pair in the group list is calculated, and filtering is performed according to a threshold (5 reads); and the filtered chromosome pairs are unfolded to obtain all SV candidate chromosome pairs.

[0103] As an implementation, the present application performs secondary alignment and SV detection on each group of SV candidate chromosomes, specifically: extract the corresponding candidate chromosomes from the microbial dataset. Use bwa mem2 to establish alignment index and use sv_fastq1 and sv_fastq2 for alignment to obtain alignment results, and finally use manta software and design python program for SV detection, and combine the SV detection results; calculate the reads support number of each chromosome group according to the final alignment results; and expand the chromosome group higher than the threshold value to obtain the candidate chromosomes.

[0104] As an implementation, the present application performs secondary alignment and single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel) detection on each candidate chromosome. Specifically: extract the corresponding candidate chromosomes from the microbial dataset. Use bwa mem2 to establish alignment index and use clean_fastq1 and clean_fastq2 for alignment to obtain alignment results, and then use GATK software for single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel) detection.

[0105] In order to further illustrate the present application, the method for detecting microbial structural variation, mutation and insertion and deletion based on second-generation sequencing provided by the present application is described in detail below in combination with examples and drawings, but they should not be understood as limiting the scope of protection of the present application.

[0106] Examples

[0107] A method for detecting microbial structural variation, mutation and insertion and deletion based on second-generation sequencing, the flowchart is shown in Figure 1 , comprising the following steps:

[0108] 1. Calculate the similarity of the reference genome and group: use the "ncbi-genome-download" tool to download the reference genome of microorganisms (bacteria, viruses, fungi, protozoa and archaea) from the RefSeq database of NCBI. According to the biological classification, use the Mash tool (default parameters) to calculate the similarity of each intra-species reference sequence, and use 0.95 as the threshold value, and the similarity greater than the threshold value as a reference chromosome group, to obtain the chromosome grouping annotation file (all_consolidated_groups.tsv).

[0109] 2. Raw sequencing sequence filtering and quality control: The raw sequencing sequence is the FASTQ format data generated by the double-end second-generation sequencing, and the length of a single sequence is usually 75-350 bp. The input data is a pair of double-end files (fastql and fastq2). This step first uses the Fastp software (default parameters) to filter the low-quality bases and adapter sequences of the raw sequencing sequence: dynamically cut off the bases with a quality value lower than Q20 or containing adapters from the 5' and 3' ends, and the maximum cutting length is 20 bp by default, and the read length is greater than or equal to 50 bp. The filtered sequence set is the so-called "pre-processed double-end sequencing sequence" in the present application, which can be the entire double-end sequencing data of the sample or a subset after random sampling, as long as its coverage (coverage depth of 5x and above) is sufficient to detect the variation of the target microorganism. The pre-processed double-end sequencing sequence has removed the adapters introduced for sequencing and low-quality fragments, and can be directly used for subsequent alignment.

[0110] Then the pre-processed double-end sequencing sequence is quality controlled by using the FastQC software (default parameters), and HTML and JSON format quality control reports are generated for subsequent evaluation of data quality.

[0111] 3. Selecting the reference sequence library, including the microbial reference sequence library and the host reference sequence library. The microbial reference sequence library is integrated by the RefSeq whole genome sequences of all known microorganisms (bacteria, viruses, fungi, protozoa and archaea), among which the bacteria are complete genomes, and the 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 the alignment reference database is created using the bowtie2-build command, and the microbial reference sequence library and the host reference sequence library can be freely set according to actual needs.

[0112] 4. Eliminating host genome sequences. In order to focus on the microbial components in the sample, this step aims to align and eliminate the host (human) genome sequences in the sample. The Bowtie2 alignment tool is used to align the sequences quality controlled in step 2 with the host reference genome sequence library in step 3. By setting --very-sensitive to run in high sensitivity mode, and using --no-mixed, --no-discordant and --no-unal parameters to optimize the alignment strategy, the sequences derived from the host genome can be accurately identified. Using the --un-conc parameter, the un-concordant pairs that failed to be successfully aligned to the host genome are output to a new FASTQ file (human unaligned sequences). These unaligned sequences are considered to be nucleic acid sequences from non-hosts in the sample, i.e. the target microbial sequences for subsequent analysis.

[0113] Since the amount of sequencing data from the host genome in the second-generation sequencing off-line data accounts for a large proportion, aligning with the host genome first can quickly exclude human contamination and significantly reduce the downstream calculation amount.

[0114] 5, "Human unaligned sequences" are aligned with the microbial sequence library. The "human unaligned sequences" obtained in step 4 are aligned to the microbial reference sequence library in step 3 using BWA-MEM. A two-stage strategy is adopted:

[0115] 1) Default parameter alignment (-M) generates default alignment results (04_final.sorted.bam), and only unique and completely matched read pairs are retained;

[0116] 2) Additional parameters -a-P are used for "final (all) alignment" to generate results of all possible alignment positions (04_all_final.sorted.bam). This file contains all alignment information, including supplementary alignment, for more in-depth analysis.

[0117] Then, the BAM index is established using samtools index (default parameters) for fast random access. Through this alignment, it can be confirmed whether the "human unaligned sequences" are truly derived from microorganisms, and the possible residual unaligned human reads are removed, thereby improving the detection accuracy.

[0118] 6, Structure variation candidate reads identification: the input is the BAM file generated after aligning to the microbial reference genome in step 5 using BWA-MEM (denoted as `04_final.sorted.bam`, where ` is a custom mark defined by the present application, and the same applies hereinafter), and the output is a TSV file containing potential structure variation (Structure Variation, SV) signals (denoted as `06.01_final_sv.tsv`).

[0119] This step is implemented through a Python script, and the core logic is as follows:

[0120] (1) Input format: The script accepts standard SAM or BAM format files, and automatically selects the reading mode according to the extension (.sam for text mode r,.bam for binary mode rb).

[0121] (2) SV determination rules:

[0122] - Rule 1 (CIGAR feature): if a read contains any of the operators I (insertion), D (deletion), S (soft clip), H (hard clip) in its CIGAR string, it is considered as a potential SV signal;

[0123] Soft clip: aligned to the reference sequence, but with low quality, not enough to be considered as a reliable match;

[0124] Hard clip: the part is completely removed, usually because these parts have no alignment or too low quality with the reference sequence;

[0125] The difference lies in: in the alignment process, the soft clip part participates in the alignment process, and the final alignment score of the read will be calculated. Hard clip does not participate in alignment and does not participate in the calculation of the alignment score of the read.

[0126] - Rule 2 (SA tag): if the read has a supplementary alignment SA tag, it indicates that it is a split-read, different segments are aligned to different positions in the genome, suggesting that there may be translocation, inversion or large fragment deletion;

[0127] - If any of the rules are met, it is determined to be an "SV candidate read".

[0128] (3) Filtering mechanism:

[0129] - First, exclude all unmapped reads with "is_unmapped=True";

[0130] - Only keep reads with primary or supplementary alignment ("is_secondary=False" or "is_supplementary=True");

[0131] - 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 position).

[0132] Through step 6, conventional linear alignment reads can be efficiently filtered out from massive alignment results, and only candidate reads carrying SV signals are retained, providing a highly concentrated dataset containing potential structural variation signals for downstream SV precise typing, thereby improving computing speed.

[0133] 7. Chromosome grouping and candidate de-duplication: input is the BAM file produced in step 5 (`04_all_final.sorted.bam`, generated by "bwa mem -a -P") and the chromosome grouping annotation file (`all_consolidated_groups.tsv`), output is multiple TSV files, including:

[0134] - `final_candidate_chromosomes.tsv`: the final retained candidate chromosomes within each group and their ordering;

[0135] - `covered_chromosomes_analysis.tsv`: details of redundant chromosomes that are eliminated by candidate chromosomes coverage;

[0136] - `new_chromosome_groups.tsv`: new group table containing only candidate chromosomes;

[0137] This step is completed by a Python script written, the core process is as follows:

[0138] (1) Chromosome grouping mapping establishment

[0139] Read `all_consolidated_groups.tsv` (no header, the first column is the group name, and the subsequent columns are the chromosomes contained in the group), and construct a bidirectional dictionary:

[0140] - `chr_to_group`: chromosome -> group name;

[0141] - `group_to_chrs`: group name -> chromosome list.

[0142] (2) BAM parsing and read de-duplication

[0143] Iterate through "04_all_final.sorted.bam":

[0144] - Distinguish reads1 / reads2 according to SAM flag or read-name suffix;

[0145] - Form a nested dictionary in units of chromosomes and with read_name as the key, only retaining the first occurrence (de-duplication):

[0146] (3) Low-abundance chromosome filtering

[0147] According to the user-set threshold `--min_reads` (recommended as 5), chromosomes with total reads (reads1+reads2) below the threshold are removed to reduce noise.

[0148] (4) Candidate chromosome screening and redundancy removal step by step

[0149] For each group, independently perform:

[0150] 1) First candidate: select the chromosome with the highest number of reads support;

[0151] 2) Subsequent candidate loop: traverse the remaining chromosomes in descending order of reads, check if their reads have been "covered" by all current candidate chromosomes;

[0152] - Coverage rules:

[0153] a. The same read has the same CIGAR in the candidate chromosome;

[0154] b. Or the candidate chromosome CIGAR is "150M" (representing full-length perfect alignment), which can cover any CIGAR corresponding read;

[0155] - If all reads of a chromosome are covered by existing candidates, mark it as redundant and remove it; otherwise, promote it as a new candidate.

[0156] (5) Result output

[0157] - `final_candidate_chromosomes.tsv`: group name, candidate ranking, chromosome, reads number (reads1 / reads2), total number of chromosomes in the group, total reads in the group, coverage rate;

[0158] - `covered_chromosomes_analysis.tsv`: covered chromosomes, coverage source candidates, covered reads number, uncovered reads number, coverage percentage;

[0159] - `new_chromosome_groups.tsv`: simplified grouping table containing only the final candidate chromosomes;

[0160] Through step 7, high homologous redundant sequences can be quickly removed in the whole genome of microorganisms, and the most representative candidate chromosomes can be retained, which can exclude the interference of similar reference chromosomes and provide a simplified and highly reliable target set for subsequent structural variation.

[0161] 8. Structural variant (SV) pair summary and confidence filtering. The input is the SV candidate TSV file generated in Step 6 (`06.01_final_sv.tsv`), and optionally a chromosome grouping annotation file (`all_consolidated_groups.tsv`). The output is:

[0162] - a chromosome-level SV pair list (`06.02_final_sv_summary.tsv`);

[0163] - a chromosome-group-level SV summary (`06.02_final_sv_group_summary.tsv`);

[0164] - when `--group-size` is enabled, an additional "detection groups" file is generated (`*.detection_groups.tsv`).

[0165] This step is accomplished by a Python script written for this purpose, with the core logic as follows:

[0166] (1) Input parsing

[0167] - Read `06.01_final_sv.tsv` and aggregate all aligned chromosomes by `read_name`.

[0168] - If `--groups` is provided, load the grouping file and establish `chrom → group` and `group → chroms` bidirectional mappings. If not provided, use the chromosome itself as the group by default.

[0169] (2) SV type classification

[0170] For each read:

[0171] - Intra-chromosome SV: aligned to the same chromosome (or the same group) only → recorded as (chr, chr) self-pairing.

[0172] - Inter-chromosome SV: aligned to ≥2 different chromosomes (or different groups) → recorded as all ordered group pairs (sorted and de-duplicated).

[0173] (3) Support count and threshold filtering

[0174] - Count the number of reads supporting each group pair.

[0175] - Filter low-confidence events by `--threshold` (default 1), keeping only group pairs with support counts ≥ threshold.

[0176] - Output both `*unfiltered.tsv` and `*group_summary.tsv` for each group, for QC purpose.

[0177] (4) Expand to specific chromosome pairs

[0178] For each group pair that is kept:

[0179] - If the number of chromosomes in the group is 1, output (chr, chr) directly.

[0180] - If the number of chromosomes in the group is ≥ 2, generate all specific chromosome combinations (sorted and de-duplicated) and attach the source group information.

[0181] - Additional filtering for single-chromosome self-pairing: if the chromosome has appeared in a cross-chromosome pair, delete its self-pairing record to avoid redundancy.

[0182] (5) Detection group generation (optional)

[0183] When `--group-size` is specified:

[0184] - Sort the high-confidence group pairs in descending order of support number.

[0185] - Merge every N group pairs into a "detection group" and record the unique chromosome list.

[0186] - Output `*.detection_groups.tsv`, containing group ID, chromosome pair details, unique chromosome number, and total support number, for subsequent parallel verification or experimental verification.

[0187] Through step 8, tens of thousands of SV candidate reads generated in step 6 can be quickly aggregated into high-confidence chromosome-level or chromosome-group-level SV events, and a simplified, parallelizable "detection group" list is provided for downstream detection, improving computational efficiency.

[0188] 9. Precise parallel detection of structural variations (SV) and insertion sequence tracing, whose input is:

[0189] - SV chromosome pair file generated in step 8 (`06.02_final_sv_summary.tsv` or `*.detection_groups.tsv`);

[0190] - Raw paired-end FASTQ (`fq1 / fq2`);

[0191] - Microbial reference genome;

[0192] - Pre-aligned BAM (`04_all_final.sorted.bam` for efficient reads extraction).

[0193] Output:

[0194] - VCF per chromosome pair or panel (in `06_sv / vcf / `);

[0195] - Optional enhanced VCF (`*.enhanced.vcf` with insert sequence source annotation).

[0196] This step is done by a Python script written by the author, and the core process is as follows:

[0197] (1) Task splitting and parallel scheduling

[0198] - Parse input files: support two formats

[0199] - "SV pair" format (`chromosome_1`, `chromosome_2`);

[0200] - "Panel to be detected" format (`unique_chromosomes` column, can contain multiple chromosomes at a time).

[0201] - Generate independent tasks according to the chromosome list; parallel processing, in each parallel process, further subdivided into: sequence extraction → BWA index → alignment sorting → manta detection.

[0202] (2) reads preprocessing

[0203] - Batch extraction: the main process extracts all read names related to the target chromosome from `04_all_final.sorted.bam` (except for special cases such as MAPQ=0).

[0204] - Build original sequence database: use SQLite database to store original FASTQ, use read name as primary key, and save to shared memory.

[0205] (3) Single task execution link

[0206] For each chromosome / group:

[0207] 1) Target sequence extraction: extract the specified chromosome to a temporary FASTA file;

[0208] 2) BWA index: create an index for this sub-reference;

[0209] 3) BWA-MEM alignment: Aligning only reads related to target chromosomes (in-memory subset);

[0210] 4) Sorting and indexing: Sorting and indexing the generated alignments;

[0211] 5) manta detection: Using manta tool (default parameters) to detect SVs from current alignments, generating VCF-formatted detection results;

[0212] 6) Insertion sequence tracing:

[0213] - Parsing records in VCF;

[0214] - Re-aligning insertion sequences to the same sub-reference, locating their source intervals on the original chromosomes;

[0215] - Generating new records and writing them to VCF, adding additional fields.

[0216] Through the parallel strategy and grouping detection strategy of step 9, high-precision SV detection results can be obtained, and the chromosomal internal source of insertion fragments is automatically labeled, improving computational efficiency.

[0217] 10, Structural variation (SV) VCF file merging and deduplication, its input is:

[0218] - All chromosome / detection group VCF files generated in step 9;

[0219] Output is:

[0220] - The main merged file after global deduplication `06_sv / merged_final.vcf`;

[0221] - High-quality subset containing only FILTER=PASS `06_sv / merged_final_pass.vcf`;

[0222] This step is completed by a Python script written by the author, and the core process is as follows:

[0223] (1) High-performance file scanning and parallel parsing

[0224] - Zero-copy traversal: Using `mmap` memory mapping + batch reading;

[0225] - Multi-process parsing: Grouping by `--batch_size` (default 100 files) and executing in parallel;

[0226] (2) Three-layer accurate deduplication algorithm

[0227] Hash grouping: compute hash with `{chr}:{pos}:{svtype}:{end}:{alt}` as key, put records into same-value buckets;

[0228] Exact match comparison within hash bucket: chromosome, svtype, alt are exactly the same; coordinate tolerance is <= 100 bp;

[0229] Quality first: keep the record with highest `QUAL` when there are true duplicates; generate unique keys with source file suffix when keys are the same but not identical;

[0230] (3) PASS filtering

[0231] - FILTER field check: keep only rows with FILTER == "PASS";

[0232] By step 10, we get a non-redundant, high-quality SV full set.

[0233] 11. SNV and indel candidate site summary, input is:

[0234] - Microbial alignment BAM file produced in step 5 (`04_final.sorted.bam`, sorted and indexed);

[0235] - Chromosome grouping annotation file (`all_consolidated_groups.tsv`);

[0236] Output is:

[0237] - Chromosome-level mutation / indel candidate list (`07.02_final_mutation_indel_summary.tsv`);

[0238] - Group-level statistics file (`*_all_groups.tsv` and `*_filtered_groups.tsv`).

[0239] This step is completed by a Python script written by the author, and the core process is as follows:

[0240] (1) BAM reading and filtering

[0241] - Traverse `04_final.sorted.bam`;

[0242] - No filtering MAPQ=0: keep all alignment records to ensure that low-quality but real existing variation signals are not missed;

[0243] - Optional primary-only filtering: retain only primary alignments with `--primary_only`, reducing multi-alignment noise;

[0244] - Aggregate all chromosomes a read maps to by `read_name`, forming a mapping dictionary.

[0245] (2) Chromosome grouping mapping (optional)

[0246] - If `--group_file` is provided, read headerless TSV:

[0247] - The first column is the group name;

[0248] - The following columns are the group member chromosomes.

[0249] - Build bidirectional mappings:

[0250] - `chromosome_to_group`: chromosome -> group name;

[0251] - `group_to_chromosomes`: group name -> chromosome set.

[0252] (3) Group-level support statistics

[0253] - No grouping mode: each chromosome forms a group by itself;

[0254] - Grouping mode: replace the chromosome a read maps to with the corresponding group name, count each read only once for each group to avoid repeated statistics;

[0255] - Generate `group_stats`:

[0256] - `support_reads`: the number of unique reads covering the group;

[0257] - `chromosomes`: the actual set of chromosomes involved in the group.

[0258] (4) Threshold filtering and expansion

[0259] - Filter high-confidence groups with `--min_support` (default 1, recommended 5);

[0260] - Expand groups that meet the threshold back to specific chromosomes:

[0261] - Each chromosome inherits the `support_reads` of the group it belongs to;

[0262] - Output fields:

[0263] `chromosome`, `group_name`, `support_reads`, `read_names`.

[0264] Output can be directly used as chromosome whitelist for step 11 (GATK mutation calling); group-level statistics help to determine which chromosomes / regions have high coverage mutation signals, and are prioritized for further analysis. The present application detects point mutations and small indels at the chromosome group level, which can improve detection accuracy.

[0265] 12. High-precision parallel detection of point mutations (SNV) and indels (indel), the input of which is:

[0266] - Chromosome-level candidate list generated in step 11 (`07.02_final_mutation_indel_summary.tsv`);

[0267] - Raw paired-end FASTQ (`fq1 / fq2`);

[0268] - Microbial reference genome FASTA;

[0269] Output: GATK VCF for each chromosome independently (`07_mutation_indel / *.gatk.vcf`);

[0270] This step is completed by a Python script written by the present application, and the core process is as follows (the software tools used are all default parameters):

[0271] Single-chromosome workflow (each task is executed in series)

[0272] 1. Sequence extraction

[0273] - `samtools faidx ref.fa chromosome>tmp / {chr} / {chr}.fa`;

[0274] 2. Indexing and dictionary

[0275] - `samtools faidx {chr}.fa` (indexing);

[0276] - `gatk CreateSequenceDictionary -R {chr}.fa -O {chr}.dict` (dictionary);

[0277] 3. BWA indexing

[0278] - `bwa index {chr}.fa`;

[0279] 4. Alignment-sort pipeline

[0280] - `bwa mem | samtools view -Sb | samtools sort -o {chr}.sorted.bam`;

[0281] 5. BAM indexing

[0282] - `samtools index {chr}.sorted.bam`;

[0283] 6. GATK HaplotypeCaller

[0284] - `gatk --java-options "-Xmx10g" HaplotypeCaller -R {chr}.fa -I{chr}.sorted.bam -O {chr}.gatk.vcf`;

[0285] Step 12 can be completed in several hours on a multi-core server to accurately detect SNV / indel of tens of millions of microbial sequencing reads, and output standard VCF files. The step 12 of the present application can accurately detect microbial mutations and insertions and deletions by parallel strategy and grouping detection strategy, and improve the calculation efficiency.

[0286] 13. Merge and deduplicate point mutation (SNV) and insertion and deletion (indel) VCF files, the input of which is:

[0287] - All chromosome / detection group VCF files generated by step 12;

[0288] The output is:

[0289] - The main merged file after global deduplication `07_mutation_indel / merged_mutation_indel.vcf`;

[0290] - High-quality subset containing only FILTER=PASS `07_mutation_indel / merged_mutation_indel_pass.vcf`;

[0291] This step is completed by a Python script written by the present application, and the core process is as follows:

[0292] (1) High-performance file scanning and parallel parsing

[0293] - Multi-process resolution: grouped by `--batch_size` (100 files by default) and executed in parallel;

[0294] (2) Three-layer accurate deduplication algorithm

[0295] Hash grouping: calculate the hash with `{chr}:{pos}:{ref}:{alt}` as the key, and put the records into the same value bucket;

[0296] Accurate matching comparison within the hash bucket: chromosome, ref, and alt are completely consistent;

[0297] Quality priority: keep the record with the highest `QUAL` when there is a real duplicate; if the keys are the same but not similar, generate a unique key with a source file suffix;

[0298] (3) PASS filtering

[0299] FILTER field check: only keep the lines where FILTER == "PASS";

[0300] Through step 13, a non-redundant and high-quality point mutation (SNV) and insertion-deletion (indel) set can be obtained.

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

[0302] Results:

[0303] The application simulates five cases of structural variation (Structural Variation, SV), which are deletion (Deletion, DEL), insertion (Insertion, INS), duplication (Duplication, DUP), inversion (Inversion, INV), and translocation (Translocation, TRA). Each type contains 4 simulation test data. And mix the simulation data and the real sequencing data (numbered SRR32394700 in the NCBI SRA database) together to simulate the sequencing data in the real scene, and the final comprehensive accuracy is 95%. The results show that the application can detect the corresponding structural variation sites (Table 1).

[0304] Table 1 SV detection results

[0305]

[0306] The application simultaneously simulates single nucleotide variation (SNV / mutation) and small fragment insertion and deletion (Indel), and mixes the simulation data and real sequencing data (SRR32394700) to simulate the sequencing data in a real scene, and the final comprehensive accuracy is 97.71%. The results show that the application can detect corresponding single nucleic acid and small fragment insertion and deletion sites.

[0307] Table 2 SNV and Indel detection results

[0308]

[0309] In summary, the application innovatively groups reference genes, implements a hierarchical candidate screening mechanism, and optimizes the data processing process, significantly solving the problems of low sensitivity, high false positives and low efficiency faced by the prior art in the detection of variations in the metagenomic background.

[0310] Although the above embodiment has made a detailed description of the application, it is only a part of the embodiments of the application, not all the embodiments, and other embodiments can be obtained according to the embodiment without creativity, which belong to the protection scope of the application.

Claims

1. A method for detecting microbial structural variation, mutation and indel based on second-generation sequencing, characterized in that, The method comprises the following steps: Calculating the similarity of the reference genome and grouping: according to the biological taxonomy, the similarity of each intra-species reference sequence is calculated, and a reference genome with a similarity greater than a threshold of 0.95 is taken as a reference chromosome set to obtain a chromosome grouping annotation file; Filtering and quality control of original sequencing sequences: filtering and quality control of original sequencing data of a microorganism to be detected to obtain quality control data; Removing host genome sequences: aligning the quality control data with a host reference sequence library to remove host reads in the sequencing data to obtain target microorganism sequencing data; Alignment: aligning the target microorganism sequencing data with a microorganism reference sequence library to obtain a preliminary alignment result; using BWA-MEM and adding -a -P parameters to align the target microorganism sequencing data with the microorganism reference sequence library to obtain a final alignment result; Identification of candidate reads of structural variations: screening the preliminary alignment result to retain potential structural variation-related reads to obtain a set of structural variation-related reads; The potential structural variation-related reads include: soft clipping reads, split alignment reads, abnormal insert fragment sizes, hard clipping reads, insert reads, deletion reads, and reads whose mapping positions of paired reads are not on the same chromosome; Chromosome grouping and candidate de-redundancy: inputting the final alignment result and the chromosome grouping annotation file, removing highly homologous redundant sequences, retaining representative candidate chromosomes, and outputting a TSV file, which includes: removing chromosomes with a total read number below a threshold; selecting chromosomes with the highest read support number; traversing the remaining chromosomes in descending order of read number and checking whether the reads have been covered by all candidate chromosomes; the covering rules include: a. the same read exists in the same CIGAR in the candidate chromosomes; b. or the CIGAR of the candidate chromosome is 150M, which can cover the corresponding reads of any CIGAR; If all reads of a chromosome are covered by existing candidates, it is marked as redundant and removed; otherwise, it is upgraded as a new candidate; Summary of structural variation chromosome pairs and confidence screening: inputting the set of structural variation-related reads and the chromosome grouping annotation file, aggregating potential structural variation-related reads in the set of structural variation-related reads into high-confidence chromosome-level or chromosome-group-level structural variation-related read events, and calculating the read support number of the corresponding reads in the final alignment result, and outputting an SV chromosome pair file; Precise parallel detection of structural variations and insertion sequence tracing: the input content includes: the SV chromosome pair file, the original sequencing data, the microorganism reference sequence library, and the final alignment result; the first processing outputs a VCF file, which includes: VCF and enhanced VCF of each chromosome pair or detection group; the first processing process includes: Only using reads related to the target chromosome for alignment; performing SV detection on the alignment result; aligning the insertion sequence to the same sub-reference again to locate its source interval on the original chromosome; Structural variant VCF file merging and deduplication: input is the VCF files, output is the main merged file after global deduplication and the high-quality subset containing only FILTER=PASS; Point mutation and small indel candidate site summary: input is the final alignment result and the chromosome grouping annotation file, output is the chromosome-level candidate list and group-level statistics file; High-precision parallel detection of point mutations and indels: its input is the chromosome-level candidate list, the original sequencing data and the microbial reference sequence library, and the output of the second processing is a GATK VCF independent of each chromosome; the second processing includes: Only use the reads related to the target chromosome for alignment; compare the alignment results to detect point mutations and indels; Point mutation and indel VCF file merging and deduplication: input is the GATK VCF independent of each chromosome, output is the main merged file after global deduplication and the high-quality subset containing only FILTER=PASS.

2. The method of claim 1, wherein, The process of the structural variant chromosome pair summary and confidence screening includes: (1) Input analysis Read the structural variant related read set, aggregate all aligned chromosomes by read name; If `--groups` is provided, load the grouping file and establish `chrom → group` and `group → chroms` bidirectional mapping; if not provided, the chromosome itself is used as the group by default; (2) Structural variant type division For each read: Intra-chromosomal structural variation: only aligned to the same chromosome or the same group → recorded as (chr, chr) self-pairing; Chromosomal structural variation: aligned to ≥2 different chromosomes or different groups → recorded as all ordered group pairs; (3) Support count and threshold filtering Count the number of reads supporting each group pair; Filter low-confidence events by `--threshold`, only keep group pairs with support number ≥ threshold; Output two group-level statistics files before and after filtering; (4) Expand to specific chromosome pairs For each retained group pair: If the number of chromosomes in the group is 1, directly output (chr, chr); If the number of chromosomes in the group is ≥2, generate all specific chromosome combinations and attach the source group information; Additional filtering for single-chromosome self-pairing: if the chromosome has appeared in a chromosomal pair, delete its self-pairing record.

3. The method of claim 2, wherein, The process of the structural variant chromosome pair summary and confidence screening also includes: (5) To-be-detected group generation When `--group-size` is specified: arrange the high-confidence group pairs in descending order of support number; merge every N group pairs into a to-be-detected group and record the unique chromosome list; output the results containing group ID, chromosome pair details, unique chromosome number and total support number.

4. The method of claim 1, wherein, The processing flow of the structural variant VCF file merging and deduplication includes three layers of accurate deduplication algorithm and PASS filtering: The three layers of accurate deduplication algorithm include: Hash grouping: compute hash with `{chr}:{pos}:{svtype}:{end}:{alt}` as key, put records into same-value buckets; Exact match comparison within hash bucket: chromosome, svtype, alt exactly match; coordinate tolerance ≤ 100 bp; Quality first: keep record with highest `QUAL` when truly duplicated; generate unique key with source file suffix when keys are identical but not similar.

5. The method of claim 1, wherein, The processing flow of the point mutation and small fragment insertion-deletion candidate site summary includes: (1) BAM reading and filtering No filtering MAPQ=0: keep all alignment records to ensure that no low-quality but real existing variation signal is missed; Optional primary alignment filtering: keep only primary alignment by `--primary_only` to reduce multiple alignment noise; Aggregate all chromosomes that a read aligns to by `read_name` to form a mapping dictionary; (2) Chromosome grouping mapping (3) Group-level support statistics No grouping mode: each chromosome forms a group by itself; Grouping mode: replace the chromosomes that a read aligns to with corresponding group names, count each group only once for each read to avoid repeated statistics; (4) Threshold filtering and expansion Expand groups that meet the threshold back to specific chromosomes: Each chromosome inherits the `support_reads` of the group it belongs to; Output fields: `chromosome`, `group_name`, `support_reads`, `read_names`.

6. The method of claim 1, wherein, The processing flow of the point mutation and insertion-deletion high-precision parallel detection includes sequence extraction, index and dictionary, BWA index, alignment-sorting pipeline, BAM index and GATK HaplotypeCaller.

7. The method of claim 1, wherein, The TSV file includes: candidate chromosomes retained in each group and their order, redundant chromosome details excluded by candidate chromosomes, and a new group table containing only candidate chromosomes; The SV chromosome pair file includes: chromosome-level structural variation pair list, chromosome-group-level structural variation statistics, and group file to be detected.

8. The method of claim 1, wherein, The microbial reference sequence library includes RefSeq whole genome sequences of known microorganisms; the host reference sequence library includes human reference genome hg19.

9. The method of claim 1, wherein, The filtering includes dynamically cutting off bases with quality value lower than Q20 or containing adapters from the 5' end and the 3' end, with the maximum cutting length being 20 bp by default, and keeping read lengths with a length of ≥50 bp; the tool for calculating the similarity of each intra-species reference sequence includes Mash; the tool used for preliminary alignment includes BWA-MEM.

Images