A multi-mapping correction algal species annotation and abundance estimation method

Abstract

This invention discloses a multi-mapping correction method for algal species annotation and abundance estimation. The invention relates to the fields of bioinformatics processing and environmental monitoring technology, and includes the following steps: acquiring raw metagenomic sequencing data of algae in aquatic environments, performing quality control, and outputting clean reads; constructing a competitive hierarchical reference database containing target, challenge, and bait sequences; aligning clean reads to this database and retaining multiple candidate alignment results; calculating the corrected abundance index for each species using a multi-mapping penalty and genome normalization model; aggregating species-level abundance to multiple taxonomic levels (phylum, class, order, family, genus, species) and independently normalizing them; calculating confidence intervals for relative abundance at each level through bootstrapping resampling; and outputting a final report based on confidence scores. This invention effectively solves the problem of distinguishing between specific and conserved sequences in multi-mapping reads, significantly improving the accuracy of algal species annotation and the reliability of abundance estimation.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics processing and environmental monitoring technology, and in particular to a multi-mapping correction method for algal species annotation and abundance estimation. Background Technology

[0002] The structure of aquatic algal communities is a key indicator for assessing water quality safety and providing early warning of algal blooms. With the development of high-throughput sequencing technology, metagenomics has become a primary method for studying aquatic microorganisms. However, existing techniques for analyzing algae using metagenomic data have the following limitations: First, algal species, especially closely related species, share numerous highly conserved gene regions, such as 18S rRNA, rbcL genes, 16S rRNA, and chloroplast genomes. Sequencing reads often align to multiple species in reference databases simultaneously. Current mainstream software, such as Kraken2 and MetaPhlAn, typically employs strategies like discarding multi-mapped reads or random allocation. The former leads to data waste and underestimation of abundance, while the latter introduces random errors. Second, the genome size of different algae can vary by several orders of magnitude, ranging from 1 Mb to 5 Gb, and the copy number of marker genes, such as rRNA, varies significantly among species. Directly using the number of aligned reads to represent biomass would severely overestimate the abundance of species with large genomes or multiple copy genes. Third, existing methods only provide an abundance value, failing to inform users whether the result originates from specific or conserved sequences, and also lacking a statistical confidence interval for abundance estimation, resulting in a lack of reliable data for environmental decision-making. Furthermore, existing methods mostly focus on species-level reporting, but short readings have limited resolution among closely related species, and species-level results are often unstable due to insufficient database completeness, lacking a mechanism for reporting downgrades to higher taxonomic levels.

[0003] Therefore, it is necessary to propose a multi-mapping correction method for algal species annotation and abundance estimation to solve the above problems. Summary of the Invention

[0004] The purpose of this invention is to provide a multi-mapped correction method for algal species annotation and abundance estimation, in order to solve the problems of multi-mapped classification caused by high homology, bias in genome size and copy number, lack of confidence quantification and statistical uncertainty assessment, and single taxonomic level in existing technologies.

[0005] This invention provides a multi-mapping correction method for algal species annotation and abundance estimation, comprising: Step 1: Obtain raw sequencing data of algal metagenomics in the aquatic environment, perform data preprocessing on the raw sequencing data, and output clean reads; Step 2: Integrate algal genome data, classify reference sequences into three categories according to their roles: target sequences, challenge sequences, and bait sequences, construct a competitive hierarchical reference database, and record metadata such as reference sequence length and marker gene copy number for each species; Step 3: Align the clean read segments to the merged competitive hierarchical reference database, retaining multiple candidate alignment results for each read segment; Step 4: Using the multi-mapping penalty-genome normalization formula, calculate the corrected abundance index of each species based on the candidate alignment results; Step 5: Aggregate the species-level corrected abundance indices upwards to multiple taxonomic levels, including phylum, class, order, family, genus, and species. Each taxonomic level is then independently normalized to obtain the relative abundance. Step 6: Perform sampling with replacement based on the reading segment identifier, and calculate the 95% confidence interval of the relative abundance of each classification level; Step 7: Normalize the corrected abundance index into a percentage, and output the final report by combining the confidence score and 95% confidence interval.

[0006] Furthermore, in step one, the data preprocessing includes removing connectors, low-quality bases, and abnormal reads.

[0007] Furthermore, in step two, the reference database adopts a competitive design and includes three types of sequences: target taxonomic units, closely related confused species, and background interference library, as well as an independent detection library for protected non-fish mammals.

[0008] Furthermore, the reference database adopts a standardized sequence header format: reference sequence identifier, role, taxonomy, species name, and marker gene.

[0009] Furthermore, in step three, sequence alignment is performed using Bowtie2 software, set to a very sensitive mode, and a maximum of 20 candidate alignment results are retained for each read segment.

[0010] Furthermore, step four includes: For each comparison to a species i Reading paragraph j Calculate its effective contribution weight:

[0011] in, To contribute weight effectively; For reading paragraph j With species i The similarity score between the reference sequence and the reference sequence is a value between 0 and 100. For reading paragraph j The total number of species matched in the database. ≥1; α is the penalty intensity coefficient, ranging from 0.5 to 1.0; α is the similarity index. For the i For each algal species, the weights of all reads aligned to that species are summed, and genome correction is performed to obtain the corrected abundance index:

[0012] in, To correct the abundance index; For species i The reference genome length; For species i The copy number of the marker gene; For all comparisons to species i A collection of reading segments; When copy number data is unavailable, it degenerates into a length-normalized form only: = / ,in For species i The total weighted reading segment.

[0013] Furthermore, in step five, the abundance index and relative abundance of taxa within a taxonomic hierarchy are defined as follows:

[0014] in, Abundance index; To correct the abundance index; Relative abundance; r Indicates the classification hierarchy, u This represents a classification unit within this level. T ( u,r ) represents all target reference sequences belonging to this classification unit; each level is normally normalized independently, and the sum of the relative abundance of all classification units within the same level is 100%.

[0015] Furthermore, in step six, sampling with replacement is performed using the segment identifier as the unit, and after repeating B times, the 2.5% and 97.5% quantiles are taken as the 95% confidence interval, with B ranging from 100 to 1000.

[0016] Furthermore, in step seven, the unique support ratio is calculated. ; in, The only support ratio; This is the only hit weighted reading segment; For total weighted readings; It is divided into three confidence levels: high, medium, and low. High: Confidence percentage ≥ 75%; Medium: Confidence percentage between 40% and 75%; Low: Confidence percentage < 40%.

[0017] Furthermore, when the confidence level of the species classification level is low, it automatically regresses to the genus classification level; when the confidence level of the genus classification level is low, it regresses to a higher classification level.

[0018] The beneficial effects of this invention are as follows: By constructing a competitive hierarchical reference database and combining it with a multi-mapped penalty-genome normalization model, this invention achieves refined decoupling processing of multi-mapped reads. The linear penalty mechanism dynamically reduces the weight of a read based on the number of species matched, effectively suppressing interference from conserved sequences on abundance signals and highlighting the true contribution of specific sequences, significantly improving the ability to distinguish closely related species and the accuracy of species identification. Furthermore, the invention introduces dual correction based on reference sequence length and marker gene copy number, eliminating systematic biases caused by differences in genome size and gene copy number variations, making the calculated corrected abundance index closer to the true cell density level. Simultaneously, the joint output mechanism of confidence percentage and Bootstrap resampling confidence interval provides users with dual risk assessment criteria from classification reliability to sampling stability. Combined with an automatic classification hierarchical degradation reporting strategy, it can degrade to the genus or higher level when species-level results are unreliable, effectively avoiding the risk of over-interpretation due to incomplete databases or insufficient sequence resolution. Therefore, the overall accuracy and reliability of metagenomic quantitative analysis of aquatic algal communities are significantly improved. Attached Figure Description

[0019] To more clearly illustrate the technical solution of the present invention, the drawings used in the embodiments will be briefly introduced below. Obviously, those skilled in the art can obtain other drawings based on these drawings without creative effort.

[0020] Figure 1 This is a flowchart of the multi-mapping correction algal species annotation and abundance estimation method of the present invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention. The technical solutions provided by various embodiments of this invention will be described in detail below with reference to the accompanying drawings.

[0022] Please see Figure 1 The multi-mapping correction algal species annotation and abundance estimation method of the present invention includes: Step 1: Obtain raw sequencing data of algal metagenomics in the aquatic environment, perform data preprocessing on the raw sequencing data, and output clean reads.

[0023] Specifically, the raw paired-end sequencing data (FastQ format) underwent adapter removal, low-quality base filtering (Q<20), and length filtering. The fastp tool was used to complete adapter detection, quality filtering, simultaneous output of paired-end reads, and generation of HTML / JSON quality control reports in a single step. The clean reads after quality control were saved to a unified directory, ensuring that different analysis versions were based on the same batch of clean reads, with differences stemming only from the reference library and algorithm.

[0024] Step 2 involves integrating algal genome data, classifying reference sequences into three categories based on their roles: target sequences, challenge sequences, and bait sequences, constructing a competitive hierarchical reference database, and recording metadata such as reference sequence length and marker gene copy number for each species.

[0025] Specifically, download public algal genomes and extract core marker genes, such as 16S, 18S, rbcL, tufA, and chloroplast genomes. Record two key parameters for each species i: reference sequence length Li (in bp) and average copy number Ci of the marker gene. Classify the reference sequences into three categories based on their role: target taxonomic units, closely related confounding species, and background interference sequences.

[0026] Target: Contains target taxonomic units that are to be included in the formal abundance denominator, such as target algae (cyanobacteria, diatoms, green algae, golden algae, etc.).

[0027] Challenge: Contains closely related species of the target species or genus. Its purpose is to test whether a reading can distinguish the target from closely related groups. Challenge records are not included in the abundance denominator, but are used to reduce the confidence level of species-level determinations or to push the reporting level down from species to genus.

[0028] Decoy: A background interference library that includes humans, livestock, waterfowl, benthic animals, crustaceans, shellfish, chloroplasts of higher plants, common bacteria, fungi, and protozoa. The decoy library is not included in the abundance denominator and is primarily used to absorb non-target background reads.

[0029] Protected non-fish library (optional): Used for independent detection of mitochondrial signals in protected non-fish mammals (such as finless porpoises and closely related dolphins). This library is not part of the fish target library and does not participate in the denominator of fish or algal abundance.

[0030] If a read truly supports only a specific target, it should preferentially and uniquely hit that target. If it hits both the target and related challenges, it indicates that the read cannot consistently support species-level conclusions. If it primarily hits decoys, it suggests that it is more likely background interference. This database structure is more robust than simply collecting a large number of target species. The database uses a standardized sequence header (FASTA header) format: reference_id, role, group, scientific_name, and marker. During the sequence alignment mapping file parsing stage, simply splitting the reference sequence name by vertical lines yields the role, group, and species information for the hit.

[0031] Step 3: Align the clean read segments to the merged competitive hierarchical reference database, and retain multiple candidate alignment results for each read segment.

[0032] Specifically, Bowtie2 is used to align clean reads to the merged reference library. A very sensitive mode is set to improve the detection capability of short reads or reads with a small number of mismatches; and a maximum of 20 candidate alignment results are retained for each read. Multi-hit information is retained as the basis for subsequent multi-mapping penalties and downgrade decisions.

[0033] Step four: Using the multi-mapping penalty-genome normalization formula, the corrected abundance index of each species is calculated based on the candidate alignment results.

[0034] The alignment results were analyzed, and the corrected abundance index for each species was calculated using the multi-mapping penalty-genome normalization model proposed in this invention. This model consists of three stages: read weight calculation, species corrected index calculation, and relative abundance normalization.

[0035] For each comparison to a species i Reading paragraph j Calculate its effective contribution weight:

[0036] in, To contribute weight effectively; For reading paragraph j With species i The similarity score between the reference sequence and the reference sequence is a value between 0 and 100. For reading paragraph j The total number of species matched in the database. ≥ 1; The penalty intensity coefficient ranges from 0.5 to 1.0, with 0.8 being preferred; α is the similarity index, which defaults to 1.0; when α > 1, the weight of high-similarity segments is amplified, while the weight of low-similarity segments is compressed; when 0 < α < 1, the similarity difference is smoothed.

[0037] If the read segment is a unique match ( = 1), the penalty term is 0, and the weight depends only on similarity.

[0038] If the read segment matches multiple species ( >1) The penalty takes effect, and the more matching species there are, the lower the weight. The linear penalty design makes the unique hit read contribute more to the confidence score, while the multiple hit read can still retain some information, but it will not be over-amplified.

[0039] For the i For each algal species, the weights of all reads aligned to that species are summed, and genome correction is performed to obtain the corrected abundance index:

[0040] in, To correct the abundance index; For species i The reference genome length; For species i The copy number of the marker gene; For all comparisons to species i A collection of reading segments; When copy number data is unavailable, it degenerates into a length-normalized form only: = / ,in For species i The total weighted reading segment. Divided by Eliminate the random matching advantage of large genomes or long sequence species; divide by Eliminate the signal amplification effect of multi-copy genes. Two-dimensional correction makes the abundance index closer to the true cell density.

[0041] Step 5: Aggregate the species-level corrected abundance index upwards to multiple taxonomic levels, including phylum, class, order, family, genus, and species. Each taxonomic level is then independently normalized to obtain the relative abundance.

[0042] In step five, the abundance index and relative abundance of taxa within a taxonomic hierarchy are defined as follows:

[0043] in, Abundance index; To correct the abundance index; Relative abundance; r Indicates the classification hierarchy, u This represents a classification unit within this level. T ( u,r ) represents all target reference sequences belonging to this classification unit. Indicates the same level r All classification units below, Indicates hierarchy r The set of all taxa is considered; each level is normally normalized independently, and the relative abundance of all taxa within the same level is summed to 100%. Therefore, percentages from different levels cannot be directly added together, but the relative contributions of different taxa can be compared within the same level.

[0044] Step 6: Perform sampling with replacement on a per-segment basis and calculate the 95% confidence interval for the relative abundance of each classification level.

[0045] Sampling with replacement was performed using the segment identifier as the unit. After B repetitions, the 2.5% and 97.5% quantiles were used as the 95% confidence intervals. The value of B ranged from 100 to 1000 to evaluate the sampling stability of the relative abundance estimation.

[0046] To assess the uncertainty of relative abundance, resampling with replacement was performed on a per-read identifier basis. For each sample, assuming the set of read identifiers available for abundance calculation is N, N read identifiers were sampled with replacement each time, and the relative abundance for all taxonomic levels was recalculated. After repeating B times, the 2.5% and 97.5% quantiles were used as the 95% confidence intervals.

[0047]

[0048] in, This represents a 95% confidence interval; This represents the 2.5th percentile (i.e., the lower bound, the 2.5th percentile). Represents the 97.5th percentile (i.e., the upper bound, the 97.5th percentile); RA represents relative abundance. This represents the relative abundance value obtained from the 1st to the Bth resampling; B is the number of repetitions, with a recommended value of 100-1000, preferably 200.

[0049] The bootstrap confidence interval reflects the sampling stability within the current read set, not the absolute correctness of the species identification itself. A narrow confidence interval indicates that the relative abundance estimate is stable under resampling, but it should still be interpreted with caution if the taxonomic unit is mainly supported by multiple mapping reads or if there is strong close competition in the challenge library. The final reliability assessment should consider both the confidence level and the bootstrap confidence interval. A high confidence level combined with a narrow confidence interval indicates that the results are stable and reliable, and can be used as the main conclusion. A moderate confidence level combined with a narrow confidence interval indicates that the value is stable but lacks unique support and should be interpreted with caution. A low confidence level (regardless of the width of the confidence interval) indicates that the results are unreliable, and it is recommended to downgrade or exclude them.

[0050] Step 7: Normalize the corrected abundance index into a percentage, and output the final report by combining the confidence score and 95% confidence interval.

[0051] This invention defines the unique support ratio as the proportion of contribution of a uniquely hit read to the total weighted reads. This value is not a statistical confidence interval, but rather a classification reliability indicator, used to measure how much of a classifier's signal comes from unique support, rather than a large number of multi-mapped reads. The unique support ratio is calculated as follows: ; in, The only support ratio; This is the only hit weighted reading segment; For total weighted readings; The higher this indicator, the higher the proportion of unique support for the classification unit. It is categorized into high, medium, and low confidence levels. High: Confidence percentage ≥ 75%, the results are reliable and can be used as the main conclusion; Medium: Confidence percentages between 40% and 75% indicate that the results are somewhat supported, but should be interpreted with caution; Low: Confidence percentage <40%, results are unreliable, it is recommended to downgrade to a higher classification level.

[0052] When the confidence level of the species classification is low, it automatically regresses to the genus classification level; when the confidence level of the genus classification is low, it regresses to a higher classification level. Example: If a sample group shows that cyanobacteria have the highest relative abundance (approximately 47%-54%) at the phylum level, that cyanobacteria are the main contributing class at the class level, and that genera such as *Chromocystis*, *Anabaena*, *Pseudo-Anabaena*, and *Anabaena* contribute at the genus level, then a consistent conclusion of cyanobacterial dominance at the phylum-class-genus levels can be formed, enhancing the reliability of the results.

[0053] The embodiments of the present invention described above do not constitute a limitation on the scope of protection of the present invention.

Claims

1. A method for multi-mapping-corrected algal species annotation and abundance estimation, characterized in that, include: Step 1: Obtain raw sequencing data of algal metagenomics in the aquatic environment, perform data preprocessing on the raw sequencing data, and output clean reads; Step 2: Integrate algal genome data, classify reference sequences into three categories according to their roles: target sequences, challenge sequences, and bait sequences, construct a competitive hierarchical reference database, and record metadata such as reference sequence length and marker gene copy number for each species; Step 3: Align the clean read segments to the merged competitive hierarchical reference database, retaining multiple candidate alignment results for each read segment; Step 4: Using the multi-mapping penalty-genome normalization formula, calculate the corrected abundance index of each species based on the candidate alignment results; Step 5: Aggregate the species-level corrected abundance indices upwards to multiple taxonomic levels, including phylum, class, order, family, genus, and species. Each taxonomic level is then independently normalized to obtain the relative abundance. Step 6: Perform sampling with replacement based on the reading segment identifier, and calculate the 95% confidence interval of the relative abundance of each classification level; Step 7: Normalize the corrected abundance index into a percentage, and output the final report by combining the confidence score and 95% confidence interval.

2. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, In step one, the data preprocessing includes removing connectors, low-quality bases, and abnormal reads.

3. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, In step two, the reference database adopts a competitive design and includes three types of sequences: target taxonomic units, closely related confused species, and background interference library, as well as an independent detection library for protected non-fish mammals.

4. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, The reference database uses a standardized sequence header format: reference sequence identifier, role, taxonomy, species name, and marker gene.

5. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, In step three, sequence alignment is performed using Bowtie2 software, set to a very sensitive mode, and a maximum of 20 candidate alignment results are retained for each read segment.

6. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, Step four includes: For each comparison to a species i Reading paragraph j Calculate its effective contribution weight: in, To contribute weight effectively; For reading paragraph j With species i The similarity score between the reference sequence and the reference sequence is a value between 0 and 100. For reading paragraph j The total number of species matched in the database. ≥1; α is the penalty intensity coefficient, ranging from 0.5 to 1.0; α is the similarity index. For the i For each algal species, the weights of all reads aligned to that species are summed, and genome correction is performed to obtain the corrected abundance index: in, To correct the abundance index; For species i Reference genome length; For species i The copy number of the marker gene; For all comparisons to species i A collection of reading segments; When copy number data is unavailable, it degenerates into a length-normalized form only: = / ,in For species i The total weighted reading segment.

7. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, In step five, the abundance index and relative abundance of taxa within a taxonomic hierarchy are defined as follows: in, Abundance index; To correct the abundance index; Relative abundance; r Indicates the classification hierarchy, u This represents a classification unit within this level. T ( u,r ) represents all target reference sequences belonging to this classification unit; each level is normally normalized independently, and the sum of the relative abundance of all classification units within the same level is 100%.

8. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, In step six, sampling with replacement is performed using the segment identifier as the unit. After repeating B times, the 2.5% and 97.5% quantiles are taken as the 95% confidence interval, with B ranging from 100 to 1000.

9. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 1, characterized in that, In step seven, calculate the unique support ratio. ; in, The only support ratio; This is the only weighted reading segment that was hit; For total weighted readings; It is divided into three confidence levels: high, medium, and low. High: Confidence percentage ≥ 75%; Medium: Confidence percentage between 40% and 75%; Low: Confidence percentage < 40%.

10. The method for multi-mapping correction of algal species annotation and abundance estimation according to claim 9, characterized in that, When the confidence level of the species classification level is low, it is automatically downgraded to the genus classification level; when the confidence level of the genus classification level is low, it is downgraded to a higher classification level.

Images