A method of microbial single-cell library chimeric filtration and SAG quality screening

Abstract

The application provides a method for chimeric filtering and SAG quality screening of a microbial single cell library. The method automatically completes chimeric filtering and quality screening of SAGs before assembly based on read-by-read classification and annotation of a metagenome classifier without the need for assembly and known species reference, so that assembly resources are concentrated on high-quality SAGs, and the assembly quality is significantly improved.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics and relates to the analysis of microbial single-cell sequencing data, specifically to a method for microbial single-cell library chimeric filtering and SAG quality screening. Background Technology

[0002] Microbial communities are widely distributed in environments such as the human body, soil, and fermentation systems. Their species composition and abundance structure are closely related to host health and industrial production. Microbe-seq, a single-cell microbial genome technology, utilizes droplet microfluidics combined with bioinformatics analysis to acquire high-throughput genome information of individual microorganisms within a microbial community at single-cell resolution (Wenshan Zheng et al., High-throughput, single-microbegenomics with strain resolution, applied to a human gut microbiome, Science 376, eabm1483, 2022). This technology can obtain tens of thousands of amplified single-cell genomes from complex microbial communities without culturing, providing a data foundation for community diversity analysis and strain identification.

[0003] In single-cell microbial genome sequencing, the sorted individual microbial cells typically undergo whole-genome amplification (such as multiple substitution amplification) to obtain sufficient DNA for library construction and sequencing. The amplified product is the single-cell amplified genome (SAG). However, during isothermal amplification using DNA polymerases such as phi29 with strand substitution and random primer annealing, the extending new strand may dissociate from the current template and continue extending with another template, or undergo rearrangement and recombination in highly branched amplification structures, ultimately forming a chimera molecule composed of two or more discontinuous template fragments. Chimerism is an inherent sequencing effect in various amplification-based library construction processes. In single-cell SAGs, interference from two-cell or multi-cell contamination, reagent contamination, and cross-contamination by closely related species can further exacerbate the problem, resulting in the inclusion of signals from non-target sources in the SAG sequencing reads. The presence of chimeric signals can contaminate downstream genome assembly, species identification, and abundance statistics, and is a significant source of noise in single-celled microbial genome analysis. Therefore, chimeric filtering and quality screening of SAGs in the library are necessary before proceeding to downstream analyses such as assembly.

[0004] However, most existing chimerism detection and quality screening methods are designed for metagenomic library construction and sequencing or amplicon sequencing, and are difficult to directly transfer to single-cell data. Specifically, existing technologies mainly have the following problems:

[0005] (1) Reference-dependent chimerism detection methods, represented by UCHIME, rely on comparing sequences with a known set of parental reference sequences to determine chimerism, which is based on the premise that the species composition of the sample is known; however, the community composition of microbial single-cell libraries is usually unknown, containing a large number of uncultured and uncharacterized species, lacking reliable parental references, so such methods cannot be effective.

[0006] (2) Reference-free chimerism detection methods, such as vsearch and uchime_denovo, rely on pairwise comparisons between a large number of reads to identify chimeric structures spliced ​​from multiple parents, which requires high and uniform sequencing depth. However, after whole-genome amplification of single-cell SAG, the sequencing depth is shallow and the coverage is severely uneven. Furthermore, the chimeric reads introduced by amplification are highly similar to the real reads in sequence. Such methods cannot reliably distinguish between chimeric reads and real reads at shallow depths.

[0007] (3) Before assembly, the existing process only has rough screening dimensions such as read count or sequencing volume for each SAG, and cannot determine whether a SAG has been mixed with chimeric or contaminated signals; the reliable judgment of SAG quality can only rely on indicators such as alignment rate and purity after assembly, while assembly and alignment require a lot of computing resources, and quality control is seriously lagging behind the investment of computing power. Once an unqualified SAG is found, the previously invested computing resources are wasted.

[0008] (4) Sequence-level chimerism detection tools can only output whether a SAG is chimeric at the read level. They cannot quantify the degree of chimerism of a SAG at different classification levels and are difficult to distinguish between species-level close relatives and severe cross-phylum contamination. At the same time, chimerism estimation based on classification ratio will be distorted when the proportion of unclassified reads in the SAG is high. Directly discarding unclassified reads will systematically underestimate the degree of chimerism.

[0009] Therefore, how to automatically and accurately perform chimeric filtering and SAG quality screening on microbial single-cell libraries with unknown community composition before assembly, without the need for assembly or reference to known species, is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0010] To address the aforementioned problems in the existing technology, this invention provides a method for chimeric filtering and SAG quality screening of microbial single-cell libraries. Without requiring assembly or known species references, it automatically completes chimeric filtering and quality screening of each SAG before assembly based on read-by-read classification annotation of a metagenomic classifier.

[0011] The technical solution adopted in this invention is: a method for chimeric filtration and SAG quality screening of microbial single-cell libraries, comprising:

[0012] S1. Input all SAG sequencing reads in the microbial single-cell sequencing library into a metagenomic classifier for classification to obtain k-mer classification annotations for each read; the k-mer classification annotations include the classification identifier taxid of each k-mer hit of the read and the corresponding k-mer count.

[0013] S2. Split the k-mer classification annotations of each read according to the SAG, and count the number of classified reads and the number of unclassified reads for each SAG;

[0014] S3. For each SAG, calculate the chirp ratio at two or more classification levels. For each target classification level: for each classifying read of the SAG, group each k-mer of the read according to the ancestor classifiers of the target classification level based on the hit taxid. Use the ratio of the k-mer count of the classifier with the largest k-mer count after grouping to the total number of k-mers with classification annotations in the target classification level as the purity of the read in the target classification level. Subtract the purity from 1 as the chirp contribution of the read in the target classification level. Use the mean of the chirp contributions of the classifying reads of the SAG as the measured chirp rate. Use the measured chirp rate to estimate the chirp contribution of the unclassified reads of the SAG. Use the ratio of the sum of the chirp contributions of the classifying reads of the SAG to the sum of the estimated chirp contributions of the unclassified reads to the total number of reads in the SAG as the chirp ratio of the SAG in the target classification level.

[0015] S4. Perform two-level screening for each SAG: The first level of screening is to remove SAGs whose unclassified proportion exceeds a first preset threshold, where the unclassified proportion is the ratio of the number of unclassified reads of an SAG to its total number of reads; The second level of screening is to remove SAGs that pass the first level of screening and whose genus-level chirp ratio exceeds a second preset threshold, wherein SAGs with missing genus-level chirp ratio values ​​are considered to have passed the second level of screening.

[0016] S5. Output the filtering result table of SAGs that have passed the two-level filtering.

[0017] The difficulty in chimerism detection stems from two inherent characteristics of single-cell SAG data: unknown community composition and shallow, uneven sequencing depth after amplification. To address this, this invention does not directly identify the splicing structure of chimeric reads using sequence alignment. Instead, it indirectly reflects the degree of chimerism by utilizing k-mer-level classification signals naturally generated by metagenomic classifiers during the classification process. A non-chimeric read from a single real template should have most of its k-mers concentrated in the same taxonomic unit; however, a chimeric read composed of two or more templates will have its k-mers dispersed across different taxonomic units. Based on this observation, this invention groups the k-mers of each read according to the ancestral taxonomic unit at the target taxonomic level based on the matched taxid. The proportion of k-mers in the dominant taxonomic unit is used as the purity of the read at that level, and its complementarity value is the chimerism contribution. This metric does not rely on any known parental references, nor does it require pairwise comparisons between reads, thus circumventing the respective applicability prerequisites of reference-dependent and reference-free methods. Furthermore, it calculates chimerism ratios at multiple taxonomic levels, allowing for the differentiation of species-level close-related contamination and cross-phylum contamination at different levels. The granularity of chimerism characterization is significantly finer than read-level chimerism determination. Secondly, considering that unclassified reads often constitute a significant proportion in single-cell data, directly discarding them would systematically underestimate the degree of chimerism. This invention estimates the chimerism contribution of unclassified reads using the measured chimerism rate of classified reads, ensuring the chimerism ratio remains robust across the entire read range. Building upon this, this invention further employs a two-stage cascaded screening strategy: the first stage characterizes the taxonomic effectiveness of SAGs using the proportion of unclassified reads, eliminating empty cells, background, and low-quality SAGs; the second stage uses the genus-level chimerism ratio as a quality gate, eliminating SAGs with severe chimerism contamination. The reason for choosing the genus level as the gating level is that species-level annotation is easily mislabeled due to the high similarity of sequences between closely related species within the same genus or the absence of such species in the database. In contrast, genus-level annotations are mostly included in taxonomic databases, and the sequence differences between genera are far greater than those between species within the same genus. This results in more robust taxonomic signals, effectively eliminating true cross-generic chimerism and minimizing the risk of incorrectly deleting valid SAGs due to species-level mislabeling. These steps work synergistically, enabling the present invention to perform SAG-by-SAG quality determination at extremely low computational cost before assembly, without requiring any species references, thus concentrating valuable assembly computational power on high-quality SAGs.

[0018] Furthermore, in step S3, the two or more classification levels include seven classification levels: domain, phylum, class, order, family, genus, and species. These seven levels cover the complete taxonomic lineage from domain to species, allowing the degree of chimerism to be quantified at any granularity, facilitating the selection of appropriate screening levels for different analytical purposes.

[0019] Further, in step S1, the k-mer classification annotation is:

[0020] ;

[0021] Where A(r) is the set of k-mer category annotations for read r; t i For the taxid of the i-th k-mer hit of read r, c i For the corresponding k-mer count, K r The total number of k-mers contained in read r;

[0022] In step S3, the k-mer count G of classification unit a at the target classification level ρ after grouping is calculated. ρ (r,a) represents the total number of k-mer values ​​K that the read can classify at the target classification level. ρ (r) The chimerism contribution χ of this read at the target classification level ρ ρ (r) is calculated using the following formula:

[0023] ;

[0024] Wherein, ρ(t) i ) for taxid t i The taxid located at the target classification level ρ in the ancestor chain;

[0025] When K ρ When (r)=0 or the read has no resolvable class unit at the target classification level, χ ρ (r)=1.0. Treating reads with no unresolved k-mer as fully nested can prevent such reads with missing information from being misclassified as pure reads and thus lowering the nesting ratio.

[0026] Further, in step S3, the measured splicing rate chimera_ratio ρ (s) are calculated using the following formulas:

[0027] ;

[0028] Where, χ ρ (r) represents the chimerism contribution of the SAG classifier read r at the target classification level ρ. These represent the categorized read set and the uncategorized read set of the SAG s, respectively. The first term in the numerator is the measured chirping contribution, the second term is the estimated chirping contribution of the uncategorized reads, and the denominator is the total number of reads in the SAG. The chirping ratio is calculated using the total number of reads to ensure that SAGs with a high proportion of uncategorized reads are not systematically underestimated in terms of chirping due to missing information.

[0029] Specifically, in step S4, for paired-end sequencing data, the unclassified ratio is the ratio of the number of unclassified single-end reads to the 2x read logs; the first preset threshold ranges from 0.3 to 0.7, preferably 0.7. The unclassified ratio reflects the classification effectiveness of SAG, and the threshold can be configured within the range of 0.3 to 0.7 according to the library noise level; libraries with higher noise levels can use a more stringent threshold.

[0030] Specifically, in step S4, the second preset threshold ranges from 0.1 to 0.3, preferably 0.2. Distributions of genus-level chirp ratios in multiple real-world projects show that the genus-level chirp ratios of high-quality libraries are mostly below 0.2, while the chirp ratios of high-noise libraries generally exceed 0.2. A value of 0.2 can effectively distinguish between high-quality and low-quality libraries. Furthermore, both the quality indicators and thresholds for the second-level screening are configurable. When not configured, the genus-level chirp ratio is used as the quality indicator to accommodate different analytical needs.

[0031] Specifically, in step S1, the metagenomic classifier is Kraken2. All SAG sequencing reads in the library are decompressed in a multi-threaded parallel manner and concatenated into a single read stream, which is then input into a single Kraken2 classification process for classification. The classification database for this process is loaded only once. Merging all SAG reads into a single read stream and inputting it into the single classification process reduces the number of times the classification database is loaded from multiple times for SAGs to once. This streamlined parallel decompression and classification pipeline avoids I / O bottlenecks and significantly reduces the time overhead of the classification process. It should be noted that any metagenomic classifier capable of outputting k-mer to taxid annotations for each read can be used in this invention, and is not limited to Kraken2.

[0032] Specifically, step S2 further includes: for each SAG, counting the number of reads covered by each taxid, using the taxonomic unit with the most covered reads as the dominant taxonomic identifier (OTU_ID) of the SAG, and using the ratio of the number of reads covered by the dominant taxonomic unit to the total number of classified reads as the purity of the SAG, generating an SAG-level abundance table; in step S3, the parsing results of each taxid parsing along the parent chain to each taxonomic level are cached and reused, and the chirp ratio of each SAG at each taxonomic level is calculated using concurrent tasks, and the calculation results are merged back into the SAG-level abundance table. After the phylogenetic parsing results are cached and reused, the parsing overhead of a single taxid is amortized, and concurrent calculation further shortens the chirp ratio calculation time in multi-SAG scenarios.

[0033] Furthermore, the method also includes step S6: aggregating the SAGs that have passed the two-level screening at the library granularity, calculating the number of retained SAGs for each library, the purity of each SAG, the mean and standard deviation of the chirp ratio at each taxonomic level, the Shannon diversity index at each taxonomic level, and the dominant taxa and their proportions, and outputting a library-level statistical table. The library-level statistical table aggregates the quality control results of each SAG into comparable statistics for each library, which can be directly used for downstream community diversity analysis.

[0034] In a second aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the above-described method for chimeric filtering of microbial single-cell libraries and SAG quality screening.

[0035] The beneficial effects of this invention are as follows: This invention utilizes k-mer-level classification annotation of metagenomic classifiers to reflect the chimeric contribution of reads at two or more classification levels using the complementary value of the purity of the dominant taxonomic unit. It eliminates the need for known species references, assembly, and pairwise comparisons between reads, enabling chimeric filtering and quality screening of single-cell libraries with unknown community composition before assembly. By estimating the chimeric contribution of unclassified reads based on the measured chimeric rate, it avoids the systematic underestimation of chimeric degree in cases with a high proportion of unclassified reads. In the two-stage screening, the first stage removes empty cells, background, and low-quality SAGs, while the second stage uses the chimeric ratio at the genus level with robust classification signals as a quality gate, considering missing index values ​​as passable, thus avoiding erroneous deletions due to data sparsity. Verified in multiple practical projects, this invention can effectively distinguish between high and low-quality libraries. After filtering, the proportion of medium- and high-quality species in the SAGs entering assembly increases from 8.7% to 64%, and the detection of high-quality species is significantly increased, concentrating assembly resources on high-quality SAGs and significantly improving the overall assembly quality. Attached Figure Description

[0036] Figure 1A simplified flowchart of the microbial single-cell library chimeric filtration and SAG quality screening method provided in the embodiments of the present invention;

[0037] Figure 2 A detailed flowchart of the microbial single-cell library chimeric filtration and SAG quality screening method provided in the embodiments of the present invention;

[0038] Figure 3 This is a distribution diagram of the hierarchical splicing ratio of the Project A (high noise) library in Example 1;

[0039] Figure 4 This is a distribution diagram of the hierarchical splicing ratio of the Project C (medium noise) library in Example 1;

[0040] Figure 5 This is a distribution diagram of the hierarchical splicing ratio of the Project D (noise heterogeneous) library in Example 1;

[0041] Figure 6 This is a distribution diagram of the hierarchical chirping ratio of the Project B (low noise) library in Example 1. Detailed Implementation

[0042] The following specific embodiments illustrate the implementation of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, unless otherwise specified, the following embodiments and features described therein can be combined with each other.

[0043] First, the terminology involved in this invention will be explained:

[0044] Single-cell amplified genome (SAG): refers to the amplified products and sequencing data obtained by amplifying the whole genome of a single microbial cell obtained from sorting (such as multiple displacement amplification (MDA)).

[0045] A read refers to a single sequencing sequence output by the sequencer; in paired-end sequencing, a pair of reads is denoted as R1 and R2.

[0046] k-mer: refers to a nucleotide substring of length k. Metagenomic classifiers segment reads into k-mers and compare them with classification databases.

[0047] Tax ID (taxid): A unique number assigned to each taxon in a taxonomic database; taxon levels are arranged from lowest to highest along the taxonomic hierarchy, such as species, genus, family, order, class, phylum, and domain.

[0048] Chimera: refers to an amplified product molecule spliced ​​from two or more discontinuous template fragments; in this invention, chimera ratio refers to the quantitative index of the degree of chimerism of SAG at a certain taxonomic level, and purity refers to the degree to which the k-mer / read of read or SAG is concentrated in the dominant taxonomic unit.

[0049] Coverage: refers to the extent to which the sequencing sequence of SAG covers the reference genome.

[0050] Unclassified ratio: refers to the proportion of SAG reads that have not received classification annotations.

[0051] The method for microbial single-cell library chimeric filtration and SAG quality screening provided in this embodiment has the following overall process: Figure 1 As shown, the detailed process is as follows: Figure 2 As shown, the organization is based on merging and classifying, splitting by SAG, statistical analysis of each SAG, multi-level nested calculation, two-level filtering, and library-level aggregation. The specific details are as follows.

[0052] S1. Merged Streaming Classification: All R1 / R2 sequencing files of SAGs in the library are collected, decompressed in multi-threaded parallel processing, and concatenated into a continuous single-read stream. This stream is then input into a single metagenomic classification process for classification. In this embodiment, the metagenomic classifier used is Kraken2, whose classification database is loaded only once throughout the entire process, reducing the loading times from multiple times for SAGs to once. Parallel decompression and classification processes are seamlessly integrated, avoiding I / O becoming a bottleneck. The classification process outputs merged, read-by-read classification results. Each line includes a classification flag (C / U, representing classified and unclassified, respectively), read name, matched taxid, read length, and a k-mer-by-k taxid:count annotation, i.e., the k-mer classification annotation for each read.

[0053] S2. SAG Splitting and SAG-by-SAG Statistics: The classification results are parsed and merged line by line. SAG identifiers are extracted from read names, and each line is written to the corresponding SAG's classification file. Simultaneously, the number of classified and unclassified reads for each SAG, as well as the read frequency for each taxid, are accumulated. Further, a classification report is generated for each SAG, counting the number of reads covered by each taxid. The taxonomic unit with the most covered reads is used as the dominant taxonomic identifier (OTU_ID) for that SAG. The purity of that SAG is the ratio of the number of reads covered by the dominant taxonomic unit to the total number of classified reads. The unclassified ratio is also calculated. For paired-end sequencing data, the unclassified ratio (unclassified_ratio(s)) is the ratio of the number of unclassified single-end reads to the 2x read logarithm, i.e.:

[0054] ;

[0055] Where, N s U is the number of read logs of SAG s. s This represents the number of unclassified single-ended reads. The above statistics are written into a SAG-level abundance table, with each row corresponding to one SAG.

[0056] S3. Multi-level chimerism ratio calculation: Load the taxid node relationships, names, and hierarchical mappings from the taxonomic database. Parse each appearing taxid along its parent chain into a lineage of seven taxonomic levels: kingdom, phylum, class, order, family, genus, and species. The parsing results are cached and reused. For each SAG classification file, read and parse its k-mer classification annotations line by line. Where A(r) is the set of k-mer category annotations for read r; t i For the taxid of the i-th k-mer hit of read r, c i For the corresponding k-mer count, K r The total number of k-mers contained in read r; At each target classification level ρ, the following is performed: Each k-mer of the read is grouped according to its ancestor classification unit at that level based on the matched taxid. After grouping, the k-mer count of classification unit a at level ρ is calculated. (For all ρ(t) i (sum of i values ​​of a), the total number of k-mers that this read can classify at the target classification level. (Sum over all i), where ρ(t) i ) for taxid t iThe taxid of the ancestral chain at level ρ; the purity of the read at that level is the proportion of the k-mer count of the classifier with the largest k-mer count after grouping to the total number of k-mers in that read, and the chimerism contribution of the read at that level is... When K ρ When (r)=0 or the read has no resolvable classification unit at that level, χ ρ (r) = 1.0, indicating complete chimerism. Subsequently, the mean chimerism contribution of all classifiers of this SAG is used as the measured chimerism rate. (When the number of classified reads is 0, the value is 0). The chirping contribution of unclassified reads is estimated based on the measured chirping rate. Finally, the chirping ratio of the SAG at this level is calculated based on the total number of reads.

[0057]

[0058] in, These are the categorized read set and the uncategorized read set for the SAG s, respectively. The chirp ratio of all SAGs is calculated using concurrent tasks, and the calculated 7-level lineage names, identifiers, and 7-level chirp ratios are merged back into the SAG-level abundance table.

[0059] S4. Two-level SAG screening: Two-level screening applies to the SAG-level abundance table. The first-level screening removes SAGs whose unclassified proportion exceeds a first preset threshold; that is, the SAG set that passes the first-level screening. , where θ u The first preset threshold is set to a value ranging from 0.3 to 0.7. In this embodiment, 0.7 is used by default. That is, SAGs with an unclassified proportion exceeding 0.7 are considered invalid for classification (empty cells, background, or low-quality SAGs) and are removed. The second-level screening involves removing SAGs that have passed the first-level screening from those with a genus-level chirp ratio exceeding the second preset threshold. In other words, the set of SAGs that have passed both levels of screening is now complete. , where θ f The second preset threshold value ranges from 0.1 to 0.3, and is set to 0.2 by default in this embodiment; the quality index function v(·) is the chimera_ratio of the hierarchical level. genus (·), For safe mapping of missing / non-numerical indicators: when v(s) is numerical When v(s) is missing In other words, SAGs with missing chirp ratio values ​​at the genus level are considered to have passed the second-level screening to avoid accidental deletion due to data sparsity. Furthermore, both the quality metric and threshold for the second-level screening can be configured by the user; when not configured, the chirp ratio at the genus level is used as the quality metric.

[0060] S5. Output Screening Results Table: Output the screening results table of SAGs that have passed the two-level screening, that is, the SAG subset that has passed the screening in the SAG abundance table, for downstream assembly and analysis.

[0061] S6. Library-level aggregation: The SAGs that have passed the two-level screening are aggregated at the library level. The mean and standard deviation of the number of retained SAGs, purity, and chimerism ratio at each taxonomic level, as well as the Shannon diversity index at each taxonomic level and the dominant taxa and their proportions, are calculated for each library. The library-level statistical table is output for reuse in downstream community diversity analysis.

[0062] Example 1: Multi-item hierarchical splicing ratio distribution and threshold validity verification

[0063] In four real-world microbial single-cell libraries, the genus-level chimerism ratio of each SAG was calculated using the method described in this invention to verify the effectiveness of the second preset threshold. Project A was a high-noise project, Project B a low-noise project, Project C a medium-noise project, and Project D had varying noise levels across libraries.

[0064] The distribution of hierarchical splicing ratios for each project is as follows: Figures 3 to 6 As shown. Figure 3 The distribution of Project A (high noise) shows a wide distribution of chirp ratios that are significantly higher than the 0.2 threshold, indicating poor overall library quality. Figure 4 The distribution of Project C (medium noise) shows that the splicing ratio is concentrated in the middle section, and some SAG splicing ratios are relatively high. Figure 5 The distribution of Project D shows significant differences in the chirping ratios among different libraries, resulting in inconsistent quality. Figure 6 For Project B (low noise), the splicing ratio is concentrated in the low value region below 0.15, indicating good library quality. As can be seen from the distribution, the genus-level splicing ratio of high-quality libraries is mostly below 0.2, while the splicing ratio of high-noise libraries generally exceeds 0.2, proving that 0.2 can be used as an effective quality threshold for the genus-level splicing ratio, i.e., the preferred value of the second preset threshold.

[0065] Example 2: Improvement of Assembly Yield by Blended Filtration

[0066] On the four libraries of Project C, a two-level screening was performed using the method described in this invention (in this embodiment, the first preset threshold is 0.3, and the second preset threshold is 0.2, i.e., SAGs with a genus-level chirp ratio greater than 0.2 or an unclassified ratio greater than 0.3 were removed). The SAGs that passed the screening were then assembled, and the results were compared with those of the default assembly (all SAGs were directly assembled). The number of SAGs before and after filtering for each library is shown in Table 1.

[0067] Table 1. Number of SAGs before and after filtering for each library in Project C

[0068] library SAG count before filtering SAG count after filtering Lib-C1 7304 1177 Lib-C2 8163 847 Lib-C3 17675 101 Lib-C4 12968 120

[0069] Assembly result comparison: The default assembly yielded 5 high-quality species and 126 medium-quality species, with medium- and high-quality SAGs accounting for 8.7% (3879 / 44536) of the total SAGs. After filtering, the assembly yielded 3 high-quality species and 16 medium-quality species, with the proportion of medium- and high-quality SAGs increasing to 64%. The average purity of the filtered assembly was 0.72, and the average coverage was 0.13. In the default assembly, over 90% of the SAGs were low-quality species and failed to contribute to the effective assembly results, indicating that issues such as chimerism contamination left most SAGs in a low-quality state. The significant increase in this proportion after filtering demonstrates that the two-stage screening of this invention effectively eliminated highly chimeric, low-quality SAGs, concentrating assembly resources on high-quality SAGs and significantly improving the overall assembly quality.

[0070] The same comparative verification was performed on 11 libraries in Project D. The number of SAGs before and after filtering for each library is shown in Table 2.

[0071] Table 2. Number of SAGs before and after filtering for each library in Project D

[0072] library SAG count before filtering SAG count after filtering Lib-D1 5592 3109 Lib-D2 13072 9421 Lib-D3 23076 8840 Lib-D4 17437 3844 Lib-D5 23992 3679 Lib-D6 3544 2033 Lib-D7 11938 2917 Lib-D8 26430 4933 Lib-D9 6324 1508 Lib-D10 22695 5280 Lib-D11 13992 1508

[0073] Assembly results comparison: The proportion of SAGs belonging to medium-to-high quality species increased from 45% before filtration to 58% after filtration; before filtration, 24 high-quality species and 429 medium-quality species were obtained, while after filtration, 59 high-quality species and 233 medium-quality species were obtained, showing a significant increase in the number of high-quality species detected; the average purity of the assembled samples after filtration was 0.75, and the average coverage was 0.1. The assembly quality indicators before and after filtration for both projects are summarized in Table 3.

[0074] Table 3 Comparison of assembly quality indicators before and after filtration

[0075] project index Before filtering (default assembly) Assembly after filtration Project C SAG percentage of medium and high quality species 8.7% 64% Project C Number of high-quality / medium-quality species 5 / 126 3 / 16 Project C Average purity — 0.72 Project C Average coverage — 0.13 Project D SAG percentage of medium and high quality species 45% 58% Project D Number of high-quality / medium-quality species 24 / 429 59 / 233 Project D Average purity — 0.75 Project D Average coverage — 0.1

[0076] Example 3: Validation of Mock Sample Chimeras

[0077] To verify the existence of chimeras and to ensure that the filtering strategy of this invention does not excessively delete valid SAGs, the verification was performed based on the alignment results of paired reads from mock samples composed of known communities. If the alignment of R1 and R2 of a pair of reads meets any of the following conditions, the read is determined to be a possible chimera: (1) R1 and R2 are aligned to different species; (2) R1 and R2 are aligned to different contigs of the same species; (3) R1 and R2 are aligned to the same contig but the distance is greater than 1000 bp.

[0078] The above-mentioned cases were detected in the mock samples. For example, the R1 of Sample-A was aligned to the vicinity of the 1794492 site on NZ_CP015134.1 of Klebsiellapneumoniae, and the R2 was aligned to the vicinity of the 1312350 site on NZ_CP015134.1 of the same species. The distance between the two was greater than 1000 bp, which belongs to the chimerism of case (3); the R1 of Sample-B was aligned to NZ_LRDO01000047.1 of Escherichia coli, and the R2 was aligned to NZ_CP020102.1 of Bacillus subtilis, which belongs to the cross-species chimerism of case (1).

[0079] The chimerism detection results of the two mock samples are shown in Table 4. The average chimerism rate of Mock-A is 0.3%, and the average chimerism rate of Mock-B is 0.22%, indicating that chimerism events do exist but account for a limited proportion in the mock samples. The filtering strategy of this invention will not excessively delete valid SAGs.

[0080] Table 4. Results of chimerism detection in Mock samples

[0081] sample Average splicing rate Examples of detected chimerism Mock-A 0.3% Long distance (>1000bp) within the same contig Mock-B 0.22% R1 / R2 cross-species

[0082] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made by those skilled in the art to the technical solutions of the present invention without departing from the spirit of the present invention should fall within the protection scope of the present invention.

Claims

1. A method for chimeric filtration and SAG quality screening of microbial single-cell libraries, characterized in that, include: S1. Input all SAG sequencing reads in the microbial single-cell sequencing library into a metagenomic classifier for classification to obtain k-mer classification annotations for each read; the k-mer classification annotations include the classification identifier taxid of each k-mer hit of the read and the corresponding k-mer count. S2. Split the k-mer classification annotations of each read according to the SAG, and count the number of classified reads and the number of unclassified reads for each SAG; S3. For each SAG, calculate the chirp ratio at two or more classification levels. For each target classification level: for each classifying read of the SAG, group each k-mer of the read according to the ancestor classifiers of the target classification level based on the hit taxid. Use the ratio of the k-mer count of the classifier with the largest k-mer count after grouping to the total number of k-mers with classification annotations in the target classification level as the purity of the read in the target classification level. Subtract the purity from 1 as the chirp contribution of the read in the target classification level. Use the mean of the chirp contributions of the classifying reads of the SAG as the measured chirp rate. Use the measured chirp rate to estimate the chirp contribution of the unclassified reads of the SAG. Use the ratio of the sum of the chirp contributions of the classifying reads of the SAG to the sum of the estimated chirp contributions of the unclassified reads to the total number of reads in the SAG as the chirp ratio of the SAG in the target classification level. S4. Perform two-level screening for each SAG: The first level of screening is to remove SAGs whose unclassified proportion exceeds a first preset threshold, where the unclassified proportion is the ratio of the number of unclassified reads of an SAG to its total number of reads; The second level of screening is to remove SAGs that pass the first level of screening and whose genus-level chirp ratio exceeds a second preset threshold, wherein SAGs with missing genus-level chirp ratio values ​​are considered to have passed the second level of screening. S5. Output the filtering result table of SAGs that have passed the two-level filtering.

2. The method according to claim 1, characterized in that, In step S3, the classification hierarchy includes kingdom, phylum, class, order, family, genus, and species.

3. The method according to claim 1, characterized in that, In step S1, the k-mer classification annotation is: ； Where A(r) is the set of k-mer category annotations for read r; t i For the taxid of the i-th k-mer hit of read r, c i For the corresponding k-mer count, K r The total number of k-mers contained in read r; In step S3, the k-mer count G of classification unit a at the target classification level ρ after grouping is calculated. ρ (r,a), the total number of k-mers that the read can classify at the target classification level, K. ρ (r) The chimerism contribution χ of this read at the target classification level ρ ρ (r) is calculated using the following formula: ； Wherein, ρ(t) i ) for taxid t i The taxid located at the target classification level ρ in the ancestor chain; When K ρ When (r)=0 or the read has no resolvable class unit at the target classification level, χ ρ (r)=1.

0.

4. The method according to claim 1, characterized in that, In step S3, the measured splicing rate chimera_ratio ρ (s) are calculated using the following formulas: ， ; Where, χ ρ (r) represents the chimerism contribution of the SAG classifier read r at the target classification level ρ. These are the categorized read set and the uncategorized read set of the SAG, respectively.

5. The method according to claim 1, characterized in that, In step S4, for paired-end sequencing data, the unclassified ratio is the ratio of the number of unclassified single-end reads to the 2-fold read logs; the first preset threshold ranges from 0.3 to 0.7; the second preset threshold ranges from 0.1 to 0.

3.

6. The method according to claim 1, characterized in that, In step S1, the metagenomic classifier is Kraken2; the sequencing reads of all SAGs in the library are decompressed in parallel by multiple threads and then concatenated into a single read stream, which is then input into the Kraken2 single classification process for classification.

7. The method according to claim 1, characterized in that, Step S2 further includes: for each SAG, counting the number of reads covered by each taxid, using the taxonomic unit with the most covered reads as the dominant taxonomic unit of the SAG, and using the proportion of the number of reads covered by the dominant taxonomic unit to the total number of classified reads as the purity of the SAG, and generating an SAG-level abundance table.

8. The method according to claim 7, characterized in that, In step S3, the parsing results of each taxid parsed along the parent chain to each classification level are cached and reused, and the chirp ratio of each SAG at each classification level is calculated using concurrent tasks. The calculation results are then merged back into the SAG-level abundance table.

9. The method according to any one of claims 7 or 8, characterized in that, It also includes step S6: aggregating the SAGs that have passed the two-level screening at the library level, calculating the number of retained SAGs, the purity of each SAG, the mean and standard deviation of the chirp ratio of each taxonomic level, the Shannon diversity index of each taxonomic level, and the dominant taxonomic units and their proportions, and outputting a library-level statistical table.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1 to 9.

Images