Method, system and device for removing contamination in single-cell microbial genome species assembly results based on classification feature gene sets

Abstract

The present invention provides a method, system and device for removing contamination in the results of single-cell microbial genome species assembly based on classification feature gene sets. Based on the regularity of the appearance of feature gene sets in species genomes, the present invention gradually searches for the best combination of length threshold and coverage threshold, which is used to filter the results of single-cell microbial genome species assembly, thereby removing contaminated contigs assembled from background reads. The method provided by the present invention not only overcomes the defect of the prior art that the length threshold cannot be automatically determined, but also does not need to assume that the coverage of the contigs of the species is normally distributed or multi-peaked when determining the contig coverage threshold, thereby achieving high-throughput detection and automatic removal of multi-cell contamination and background contamination in the results of single-cell genome assembly, and obtaining the best contamination removal results as much as possible.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and specifically relates to a method, system and device for removing contamination in the assembly results of single-cell microbial genome species based on a classification feature gene set. Background Art

[0002] The relationship between microorganisms and humans is complex and close. Microorganisms can be pathogens that cause diseases, such as bacteria and viruses. However, most microorganisms are beneficial to humans, and they play an important role in food production (such as fermentation), medicine (such as the production of antibiotics), and environmental protection (such as sewage treatment). Different microorganisms play different biological functions through their different gene sequences. Through sequencing, microbial genes can now be detected on a large scale, so as to deeply understand their mechanism and significance.

[0003] Single-cell microbial genomics technology (microbe-seq) is a technical means of high-throughput detection of the entire gene sequence of a single microorganism using a single-cell sequencing method. It was published in Science in 2022 (Wenshan Zheng et al., High-throughput, single-microbe genomics with strain resolution, applied to a human gut microbiome. Science 376, eabm1483 (2022). DOI: 10.1126 / science.abm1483). This technology uses a variety of droplet microfluidic operation technologies combined with bioinformatics analysis methods to obtain the genome information of thousands of single-cell microorganisms from complex microbial communities without cultivation, and can efficiently distinguish similar species based on the reads of microorganisms in each droplet in the single-cell microbial genome sequencing results, distinguish strains through mutations in individual microorganisms, detect interactions between microorganisms and viruses, and discover horizontal gene transfer events in microbial populations. Compared with previous microbial genome detection technologies, single-cell microbial genome technology has significantly improved in throughput and resolution. It can not only quickly identify different strains, but also provide more detailed biological characteristics analysis based on the characteristic genes of each strain.

[0004] Single-cell transcriptome sequencing based on microfluidics can use micro-scale pipelines and control systems to efficiently and accurately operate and analyze single cells. For example, the 10x Genomics Chromium platform, which is based on the principle of an oil-in-water micro-reaction system, is widely used. It performs transcriptome analysis of single cells by capturing, dissolving, reacting, and sequencing single cells. In terms of workflow, the 10x Genomics Chromium platform first uses a microfluidic chip to capture single cells, ensuring that each cell is separated and fixed in a specific position, then lyses each single cell, extracts RNA, and synthesizes cDNA. The synthesized cDNA is amplified and built into a library, and barcode tags and unique molecular identifier (UMI) sequences are added. High-throughput sequencing is then performed to assemble and analyze the reads in SAG (single-amplified genome, referring to a collection of single microbial sequencing sequences in a droplet). The analysis includes gene expression levels, cell type identification, and population structure analysis.

[0005] However, due to its technical characteristics, single-cell sequencing methods based on microfluidics inevitably introduce two types of contamination: (1) background contamination, which is manifested by the reads from a certain species being uniformly distributed in all droplets without specificity due to factors such as flow channel contamination and droplet rupture. This part of the background reads will be assembled into contigs with shorter length and lower coverage; (2) multi-cell contamination, which is manifested by two or more cells being encapsulated in a droplet. Both types of contamination will affect the assembly effect and need to be removed as much as possible during the analysis process.

[0006] In the prior art, Dr. Zheng Wenshan et al. also provided an algorithm for removing multicellular contamination and background contamination based on single-cell microbial sequencing data in the article (Wenshan Zheng et al., High-throughput, single-microbe genomics with strain resolution, applied to a human gutmicrobiome. Science 376, eabm1483 (2022). DOI: 10.1126 / science.abm1483). Among them, because the contigs assembled from background reads usually have lower coverage and length, this method removes contigs below a certain length according to a fixed empirical threshold. The probability density distribution of the coverage of the remaining contigs is then calculated, and the probability density distribution is fitted with a multi-normal distribution to calculate the main peak from the non-contaminated contigs and the contaminated peaks of the contaminated contigs. The coverage threshold is calculated according to the position of the peak, and the contigs with coverage below the threshold are removed.

[0007] However, the above-mentioned method for removing background contamination still has the following problems. First, the method filters the length of contig according to a fixed empirical threshold, so it is impossible to automatically calculate the contig length threshold, and thus it is impossible to perform high-throughput filtering. Secondly, the coverage probability density of contig may not show a normal distribution, so the probability density of coverage may show multiple unclear peaks, so it may not be possible to calculate the coverage threshold according to the position of the peak, making the coverage threshold unreliable. Summary of the invention

[0008] In order to solve the problems in the prior art of removing background contamination in the results of single-cell microbial genome species assembly that (1) the contig length threshold cannot be automatically calculated; and (2) the contig coverage threshold is not reliable, the present invention provides a method, system and device for removing contamination in the results of single-cell microbial genome species assembly based on a classification feature gene set.

[0009] The technical solution adopted by the present invention is: a method for removing contamination in the assembly results of single-cell microbial genome species based on a classification feature gene set, comprising: removing multi-cell contamination and background contamination in a bin in turn;

[0010] Among them, the methods for removing background pollution in the Bin include:

[0011] Assembling a Bin with polysomal contamination removed, evaluating the obtained assembly result with polysomal contamination removed, obtaining a characteristic gene set of the Bin, and determining classification information and a reference characteristic sequence of the Bin based on the characteristic gene set; wherein the characteristic gene set is composed of characteristic genes with characteristic sequences;

[0012] Set the contig length search boundary, length search step, coverage search boundary, and coverage search step in the assembly results after removing polysome contamination, and then gradually search for the contig's hypothetical length threshold and hypothetical coverage threshold;

[0013] Combine different hypothetical length thresholds and hypothetical coverage thresholds to form several cleaning thresholds, and filter the contigs with characteristic sequences in the assembly results after removing polycellular contamination according to different cleaning thresholds to obtain retained contigs; evaluate the completeness and contamination of the assembly results after removing polycellular contamination under each cleaning threshold according to the reference characteristic sequence of Bin and the characteristic sequence distributed in the retained contigs; use the completeness and contamination as indicators to rate the quality of the assembly results after removing polycellular contamination under each cleaning threshold, and select the best cleaning threshold according to the quality rating;

[0014] According to the optimal cleaning threshold, the contigs in the assembly results after removing polycellular contamination are filtered again to obtain the assembly results after removing polycellular contamination and background contamination.

[0015] A marker gene set refers to a set of non-redundant and conservative marker genes in a species genome, and a marker gene has multiple feature sequences. In the complete and uncontaminated single-cell microbial genome assembly result sequence, the feature genes all appear and only appear once. Therefore, the present invention uses the appearance rule of feature genes in the species genome to propose a new method for removing background contamination in the single-cell microbial genome species assembly result. First, the checkm software is used to evaluate the assembly result after removing multicellular contamination, obtain the feature gene set of the Bin, and obtain the species classification information of the Bin and the feature sequence of the corresponding species according to the feature gene set. Different combinations of assumed length thresholds and assumed coverage thresholds are used to filter the contigs in the assembly result after removing multicellular contamination, and then the integrity and contamination of the assembly results obtained under different combinations are evaluated according to the reference feature sequence type of the Bin, the type of feature sequence retained in the contig, and the type of repeats. The assumed length threshold and assumed coverage threshold combination used by the assembly result with the best integrity and contamination are used as the best cleaning threshold. In this process, the acquisition of the optimal cleaning threshold does not require manual judgment based on inherent experience, nor does it require the assumption that the coverage of the species' contigs is normally distributed or multi-peaked, thereby achieving high-throughput detection and automated removal of background contamination in single-cell genome assembly results. Combined with multi-cell contamination removal methods, the optimal contamination removal results may be obtained. In addition, for assembly results without coverage information or inaccurate information, trusted-contig and / or untrusted-contig are generally added to fill gaps and / or determine repetitive sequences during the assembly process, and the method provided by the present invention also supports the removal of contamination from such assembly results.

[0016] In some embodiments, one L of the contig length is used as the contig length search boundary, and L is in the range of 5 to 20; one C of the contig coverage is used as the contig coverage search boundary, and C is in the range of 5 to 20.

[0017] In some embodiments, the method for evaluating the integrity and contamination of the assembly result after removing polysome contamination under the cleaning threshold according to the reference characteristic sequence of Bin and the characteristic sequence of retained contig includes:

[0018] Determine the number of characteristic sequences to be retained in the contig and the number of characteristic sequences that appear repeatedly in the contig;

[0019] According to the ratio of the number of characteristic sequences in the retained contig to the number of reference characteristic sequences in the Bin, the completeness of the assembly result after removing polysome contamination under the cleaning threshold is determined;

[0020] The contamination degree of the assembly result after removing polysome contamination under the cleaning threshold is determined according to the ratio of the number of species of repeated characteristic sequences in the retained contig to the number of species of the reference characteristic sequences in the Bin.

[0021] Since in the complete and pollution-free single microorganism assembly result sequence, the characteristic genes all appear and only appear once. Therefore, the present invention determines the completeness of the assembly result of removing polycellular contamination under the cleaning threshold value according to the ratio of the number of species of characteristic sequences in the retained contig to the number of species of the reference characteristic sequences of the Bin, that is, the larger the ratio, the closer the number of species of characteristic sequences in the retained contig is to the number of species of the reference characteristic sequences of the Bin, which represents the more complete the assembly result. The present invention determines the contamination degree of the assembly result of removing polycellular contamination under the cleaning threshold value according to the ratio of the number of species of characteristic sequences that repeatedly appear in the retained contig to the number of species of the reference characteristic sequences of the Bin, that is, the larger the ratio, the more species of characteristic sequences that repeatedly appear in the retained contig, and the repeated characteristic sequences may come from the contaminated contigs assembled by the background reads, which represents the more serious the contamination of the assembly result.

[0022] In some embodiments, the assembly results of removing polytopic contamination at each cleaning threshold are rated for quality using completeness and contamination as indicators, and the method for selecting the best cleaning threshold according to the quality rating includes:

[0023] According to the preset classification conditions, the quality of the assembly results after removing polysomal contamination at each cleaning threshold is graded with completeness and contamination as indicators, and the quality rating includes high quality, medium quality and low quality;

[0024] Taking completeness and contamination as indicators, the distances between the assembly results after polytopic contamination removal at each cleaning threshold and the preset ideal points, high-quality points, and medium-quality points are calculated.

[0025] The cleaning threshold of the assembly result with the smallest distance to the ideal point and the highest quality rating for removing polytopic contamination is selected as the optimal cleaning threshold.

[0026] In some embodiments, the distances between the assembly results after polytopic contamination removal at each cleaning threshold and the preset ideal points, high-quality points, and medium-quality points are calculated as follows:

[0027] When the cleaning threshold has the assumed length threshold when the search step is i and the assumed coverage threshold when the search step is j, the completeness of the assembly result after removing polytopic contamination under this cleaning threshold is Completeness i,j , Contamination i,j ;

[0028] The distance between the assembly result after removing polycellular contamination under the cleaning threshold and the preset ideal point is DistanceIdeal i,j for:

[0029] ;

[0030] The distance between the assembly result after removing polycellular contamination under this cleaning threshold and the preset high-quality point is DistanceHigh i,j for:

[0031] ;

[0032] The distance between the assembly result after removing polycellular contamination under this cleaning threshold and the preset medium mass point DistanceMedium i,j for:

[0033] .

[0034] In some embodiments, the method for removing polycellular contamination in a Bin comprises:

[0035] Align all SAGs in Bin to the original assembly result conitg of Bin to obtain the SAG-conitg matrix M;

[0036] Iterate the following process until the preset termination condition is met: for the matrix M or the matrix M obtained in the previous iteration G1 、M G2 Perform hierarchical clustering and perform clustering on conitg and matrices M, M in turn according to the hierarchical clustering results. G1 or M G2Perform binary division to obtain ContigGroup1, ContigGroup2, M1, and M2 respectively; calculate the classification probability P1 and P2 of SAG belonging to ContigGroup1 and ContigGroup2 respectively according to M1, M2, and the total number of reads in each SAG; classify SAG according to P1, P2, and the preset classification threshold, and take SAGs whose P1 or P2 meets the classification threshold as retained SAGs and assign them to SAGGroup1 or SAGGroup2 respectively, and take SAGs whose P1 and P2 do not meet the classification threshold as multi-cell SAGs and discard them; obtain the SAG-conitg matrix M of the retained SAGs belonging to SAGGroup1 or SAGGroup2 from the matrix M G1 、M G2 ;

[0037] After the termination condition is met, the classification results of this round of iteration are invalidated, and the retained SAGs obtained before this round of iteration are merged and assembled to obtain the Bin without polytopic contamination and the assembly result without polytopic contamination.

[0038] In some embodiments, the termination condition is any one of the following (1) to (3):

[0039] (1) The ratio of the number of retained SAGs in SAGGroup1 and SAGGroup2 to the number of all SAGs in the Bin is less than the preset retention threshold;

[0040] (2) The number of reserved SAGs in SAGGroup1 or SAGGroup2 is less than or equal to the preset number threshold;

[0041] (3) The difference between the number of retained SAGs in SAGGroup1 or SAGGroup2 and the number of all SAGs in the Bin is less than or equal to the preset difference threshold.

[0042] The method to remove polysome contamination is to align the reads in the SAG of the Bin to the conitg in the original assembly result to obtain the SAG-conitg matrix. The matrix value is how many reads in the SAG are aligned to the conitg. Then perform hierarchical clustering on the matrix and divide the contig into two groups. If the Bin is clustered due to polysome contamination, there will be a clearer binary result. Most of the non-polysome SAGs in the Bin will be aligned to one of the groups, while the polysome SAGs in the Bin will be aligned to two groups, and then the polysome SAGs that do not meet the classification threshold will be removed from the two groups. This method can be used to remove polysome SAGs. Reassembling the SAG after removing the polysomes and regrouping can remove the contamination caused by polysomes in the results.

[0043] In some embodiments, the method for removing contamination in the assembly results of single-cell microbial genome species based on the classification feature gene set also includes:

[0044] Use checkm to re-evaluate the assembly results obtained after removing polysomal contamination and background contamination. The evaluation content includes completeness and contamination.

[0045] The present invention also provides a system for removing contamination in the results of assembling single-cell microbial genome species based on a classification feature gene set, comprising:

[0046] The module for removing polytopic contamination in Bin is used to align all SAGs in Bin to the original assembly result conitg of Bin to obtain the SAG-conitg matrix M; iterate the following process until the preset termination condition is met: for the matrix M or the matrix M obtained in the previous iteration G1 、M G2 Perform hierarchical clustering and perform clustering on conitg and matrices M, M in turn according to the hierarchical clustering results. G1 or M G2 Perform binary division to obtain ContigGroup1, ContigGroup2, M1, and M2 respectively; calculate the classification probability P1 and P2 of SAG belonging to ContigGroup1 and ContigGroup2 respectively according to M1, M2, and the total number of reads in each SAG; classify SAG according to P1, P2, and the preset classification threshold, and take SAGs whose P1 or P2 meets the classification threshold as retained SAGs and assign them to SAGGroup1 or SAGGroup2 respectively, and take SAGs whose P1 and P2 do not meet the classification threshold as multi-cell SAGs and discard them; obtain the SAG-conitg matrix M of the retained SAGs belonging to SAGGroup1 or SAGGroup2 from the matrix M G1 、M G2 ; After the termination condition is met, the classification results of this round of iteration are invalidated, and the retained SAGs obtained before this round of iteration are merged and assembled to obtain the Bin without polytopic contamination and the assembly result without polytopic contamination;

[0047] A module for removing background pollution in Bins is used to evaluate the assembly results of Bins with polycellular pollution removed, obtain feature information of Bins, and determine classification information and reference feature sequences of Bins based on the feature information; wherein the feature information includes feature sequences, feature genes, and feature gene sets; set the length search boundary, length search step, coverage search boundary, and coverage search step of contigs in the assembly results with polycellular pollution removed, and then gradually search for the assumed length threshold and assumed coverage threshold of contigs; combine different assumed length thresholds and assumed coverage thresholds to form several cleaning thresholds, and filter contigs in the assembly results with polycellular pollution removed according to different cleaning thresholds to obtain retained contigs; evaluate the integrity and contamination of the assembly results with polycellular pollution removed under each cleaning threshold according to the reference feature sequence of Bins and the feature sequence of retained contigs; perform quality rating on the assembly results with polycellular pollution removed under each cleaning threshold using integrity and contamination as indicators, and select the best cleaning threshold according to the quality rating; re-filter contigs in the assembly results with polycellular pollution removed according to the best cleaning threshold to obtain the assembly results with polycellular pollution and background pollution removed.

[0048] The present invention also provides a device for removing contamination in the assembly results of single-cell microbial genome species based on classification feature gene sets, comprising a processor and a memory; the processor and the memory are connected via a communication bus; wherein the processor is used to call and execute a program stored in the memory; the memory is used to store a program, and the program is at least used to execute the method for removing contamination in the assembly results of single-cell microbial genome species based on classification feature gene sets.

[0049] Beneficial effects of the present invention: Based on the regularity of the occurrence of characteristic gene sets in the genome of species, the present invention gradually searches for the optimal combination of length threshold and coverage threshold, which is used to filter the assembly results of single-cell microbial genome species, thereby removing the contaminated contigs assembled from background reads. The method provided by the present invention not only overcomes the defect that the length threshold cannot be automatically determined in the prior art, but also does not need to assume that the coverage of the contigs of the species is normally distributed or multi-peaked when determining the contig coverage threshold, thereby achieving high-throughput detection and automatic removal of multi-cell contamination and background contamination in the single-cell genome assembly results, and obtaining the best contamination removal results as much as possible. In addition, for assembly results without coverage information or inaccurate information, trusted-contig and / or untrusted-contig are generally added to fill gaps and / or determine repetitive sequences during the assembly process, and the method provided by the present invention also supports the removal of contamination from such assembly results. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 A schematic flow chart of a method for removing contamination in single-cell microbial genome species assembly results based on a classification feature gene set provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0051] The following describes the embodiments of the present invention through specific embodiments, and those skilled in the art can easily understand other advantages and effects of the present invention from the contents disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and the details in this specification can also be modified or changed in various ways based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that the following embodiments and features in the embodiments can be combined with each other without conflict.

[0052] The embodiment of the present invention provides a method for removing contamination in the results of assembling species of a single-cell microbial genome based on a classification feature gene set, comprising: S1. removing multicellular contamination within a Bin; S2. removing background contamination within a Bin.

[0053] S1. Methods for removing polycellular contamination in bins include:

[0054] S1-1. Use bowtie2 software to align all SAGs in Bin or the retained SAGs in the previous iteration that were assigned to SAGGroup1 or SAGGroup2 to the contigs in the original assembly results of Bin. Use samtools idxstats software to obtain the SAG-contig matrix M of SAG n aligned to contig m. n,m .

[0055] S1-2. Matrix M or the matrix M obtained in the previous iteration G1 、M G2 All values ​​are log10 transformed, such as: Use scipy's linkage function to perform hierarchical clustering on the log10 transformed matrix. Divide conitg into two parts according to the hierarchical clustering result to obtain ContigGroup1 and ContigGroup2. n,m Split and get the matrix M i,j Two sub-matrices MLog1 n,m 、MLog2 n,m .

[0056] S1-3. Calculate the classification probability P1 of SAG n belonging to ContigGroup1 and ContigGroup2 respectively based on the total number of reads in MLog1, MLog2, and SAG n (AllReadCount). n 、P2 n :

[0057] .

[0058] S1-4. According to P1 n 、P2 n , the preset classification threshold (PurityHold), classify the n-number SAG. If P1 n ≥PurityHold, then SAG belongs to SAGGroup1; if P2 n ≥ PurityHold, then SAG n is assigned to SAGGroup2; n and P2 n SAGs that do not meet the classification threshold are regarded as multi-cell SAGs and discarded, and are not included in the subsequent process. n With P2 n From the definition of P1, n With P2 n The sum is 1.

[0059] S1-5. Obtain the SAG-conitg matrix M of the reserved SAG belonging to SAGGroup1 or SAGGroup2 from the matrix M G1 、M G2 ;

[0060] S1-6. Iterate steps S1-2 to S1-5 until the preset termination condition is met; after the termination condition is met, the classification result of this round of iteration is invalidated, and the retained SAGs obtained before this round of iteration are merged and assembled to obtain the Bin with polytopic contamination removed and the assembly result with polytopic contamination removed.

[0061] The termination condition is any one of the following (1) to (3):

[0062] (1) The ratio of the sum of the number of retained SAGs in SAGGroup1 (|P1|) and the number of retained SAGs in SAGGroup2 (|P2|) to the total number of SAGs in the Bin (AllSAGCount) is less than the preset retention threshold (KeepRatio):

[0063] ;

[0064] (2) The number of reserved SAGs in SAGGroup1 is less than or equal to the preset number threshold (NHold1), or the number of reserved SAGs in SAGGroup2 is less than or equal to the preset number threshold:

[0065] ;

[0066] (3) The difference between the number of reserved SAGs in SAGGroup1 and the number of all SAGs in Bin is less than or equal to the preset difference threshold (NHold2), or the difference between the number of reserved SAGs in SAGGroup2 and the number of all SAGs in Bin is less than or equal to the preset difference threshold:

[0067] .

[0068] S2. Methods for removing background contamination in the bin include:

[0069] S2-1. Use SPADE software to assemble Bins that have been free of polycystic contamination, use checkm software to evaluate the assembly results obtained to remove polycystic contamination, obtain the characteristic gene set of Bins, determine the species classification information of Bins and the characteristic sequence of the corresponding species based on the characteristic gene set, and use the characteristic sequence of the corresponding species as a reference characteristic sequence; wherein the characteristic gene set is composed of characteristic genes with characteristic sequences. The characteristic gene set should have the distribution information and attribution information of the characteristic sequences, that is, on which contig each characteristic sequence appears and its specific position, and which characteristic sequences belong to the same characteristic gene.

[0070] S2-2. Take one L of the contig length as the contig length search boundary (MaxLen), then at each search step (Step), the contig length search step (StepLen) is:

[0071] .

[0072] Taking one C of the contig length as the coverage search boundary (MaxCov) of the contig, the coverage search step length (StepCov) of the contig at each search step (Step) is:

[0073] .

[0074] Wherein, in some embodiments, L has a value range of 5 to 20, and C has a value range of 5 to 20. The value of L or C is generally set according to the quality of the original sequencing data of the input sample. If the species reads and background reads in the Bin are not clearly distinguished, the quality of the original sequencing data of the sample is not high, and L or C can select a lower value, such as 5 to 10; if the species reads and background reads in the Bin are clearly distinguished, the quality of the original sequencing data of the sample is high, and L or C can select a higher value, such as 10 to 20. For the original sequencing data obtained by single-cell sequencing of the sample, L or C can generally take a value of 10. Those skilled in the art can select an appropriate value according to actual conditions.

[0075] Within the length search boundary, according to the length search step, the hypothetical length threshold of the contig is searched step by step from 0 (TryLen i ):

[0076] , where i is the number of length search steps.

[0077] Within the coverage search boundary, the hypothetical coverage threshold of the contig is searched step by step starting from 0 according to the coverage search step size (TryCov j ):

[0078] , where j is the number of coverage search steps.

[0079] S2-3. Different hypothesis length thresholds (TryLen i ), hypothesis coverage threshold (TryCov j ) are combined to form several cleaning thresholds, and the contigs with characteristic sequences in the assembly results after removing polycystic contamination are filtered according to different cleaning thresholds to obtain the retained contigs and the characteristic sequences distributed in the retained contigs:

[0080] ;

[0081] in,

[0082] MarkerFeature i,j TryLen is the assumed length threshold i , assuming the coverage threshold TryCov j The set of characteristic sequences distributed in the retained contigs under the cleaning threshold,

[0083] AllMarkerFeature is the collection of all feature sequences in the assembly result after removing polysome contamination.

[0084] ContigCov is the coverage of contigs with characteristic sequences.

[0085] ContigLen is the length of the contig with the characteristic sequence.

[0086] S2-4. Determine the length threshold TryLen i , assuming the coverage threshold TryCov j The number of characteristic sequences retained in the contig under the cleaning threshold (PresentMarkerFeature i,j ) and the number of species of repeated feature sequences in the contig (DuplicateMarkerFeature i,j ):

[0087] ,

[0088] ;

[0089] in,

[0090] The CountUnique function calculates the MarkerFeature i,j The number of types of characteristic sequences that appear after removing duplicates,

[0091] MarkerFeature i, j, n TryLen is the assumed length threshold i , assuming the coverage threshold TryCov j The nth feature sequence distributed in the set of feature sequences in the retained contig under the cleaning threshold,

[0092] MarkerSet n is the total number of characteristic sequences in the characteristic gene to which characteristic sequence No. n belongs,

[0093] The Count function calculates the MarkerFeature i,j , the number of occurrences of feature sequence number n.

[0094] For example, assuming the length threshold TryLen i , assuming the coverage threshold TryCov j The set of feature sequences distributed in the retained contig under the cleaning threshold MarkerFeature i,jThere are three types of feature sequences in total, namely feature sequence A, feature sequence B, and feature sequence C. Feature sequence A appears 3 times, feature sequence B appears 2 times, and feature sequence C appears 1 time. After calculation by the CountUnique function, the number of feature sequences in the retained contig (PresentMarkerFeature i,j ) is 3, namely feature sequence A, feature sequence B, and feature sequence C. After calculation by the Count function, the number of repeated feature sequences in the retained contig (DuplicateMarkerFeature i,j ) is 2, that is, feature sequence A and feature sequence B.

[0095] S2-5. According to the number of characteristic sequences in the retained contig (PresentMarkerFeature i,j ) to the number of species of the reference feature sequence of the Bin (MarkerSize), and determine the completeness of the assembly result after removing polysomal contamination under this cleaning threshold (Completeness i,j ):

[0096] .

[0097] According to the ratio of the number of repeated feature sequences in the retained contig (DuplicateMarkerFeaturei,j) to the number of reference feature sequences in the Bin, the contamination degree of the assembly result after removing polysome contamination under the cleaning threshold is determined. i,j ):

[0098] .

[0099] S2-6-1. According to the preset classification conditions, the completeness i,j ), Contamination i,j ) is an indicator for the length threshold (TryLen i ), assuming a coverage threshold (TryCov j ) cleanup threshold to remove polysomal contamination and perform quality rating on the assembly results. i,j ), the quality ratings include high quality (High), medium quality (Medium), and low quality (Low):

[0100] .

[0101] According to the completeness and contamination, the quality of the assembly results after removing polytopic contamination under different cleaning thresholds is graded. The high quality (High) rating requires the completeness to be greater than or equal to 0.9 and the contamination to be less than or equal to 0.05; the medium quality (Medium) rating requires the completeness to be greater than or equal to 0.5 and the contamination to be less than or equal to 0.1; those that do not meet the above conditions are graded as low quality (Low).

[0102] S2-6-2. Completeness i,j ), Contamination i,j ) is the indicator, calculated under the assumption that the length threshold TryLen i , assuming the coverage threshold TryCov j The distances between the assembly results after removing polytopic contamination and the preset ideal points, high-quality points, and medium-quality points respectively.

[0103] DistanceIdeal from the preset ideal point i,j for:

[0104] .

[0105] The distance from the preset high-quality point DistanceHigh i,j for:

[0106] .

[0107] DistanceMedium from the preset medium mass point i,j for:

[0108] .

[0109] S2-6-3. Select the cleanup threshold of the assembly result with the smallest distance to the ideal point and the highest quality rating to remove polytopic contamination as the optimal cleanup threshold. Specifically, the MinDistance function can be defined to return the coordinates of the minimum value of the set under the specified quality rating condition; if the set is empty, (-1, -1) is returned. If the length search step i and the coverage search step j are not equal to -1, it is considered that the length threshold TryLen is assumed. i TryCov with hypothesis coverage threshold j is the optimal cleaning threshold. Calculate i and j in the order of high-quality (High) assembly results, medium-quality (Medium) assembly results, and low-quality (Low) assembly results in step S2-6-1:

[0110] .

[0111] Table 1. Assembly results of removing polysome contamination at different cleaning thresholds

[0112] Assembly results No. 1 No. 2 No.3 No. 4 No. 5 (i,j) (2,3) (3,2) (4,6) (5,7) (7,8) <![CDATA[Completeness i,j ]]> 0.92 0.92 0.88 0.84 0.42 <![CDATA[Contamination i,j ]]> 0.02 0.03 0.04 0.07 0.12 <![CDATA[Quality i,j ]]> High High Medium Medium Low <![CDATA[DistanceIdeal i,j ]]> 0.0068 0.0073 0.016 0.0305 0.3508 <![CDATA[DistanceHigh i,j ]]> 0.0013 0.0008 0.0005 0.004 0.2353 <![CDATA[DistanceMedium i,j ]]> 0.1828 0.1813 0.148 0.1165 0.0068

[0113] For example, Table 1 illustrates the quality rating process of assembly results for removing polytopic contamination under five different cleaning thresholds in an embodiment. Traverse each assembly result in order and set assembly result No. 1 to the initial value. First, select the assembly result with the highest quality rating, that is, assembly result No. 1 and assembly result No. 2, both rated as High. Then compare the distances between the assembly results with the highest quality rating and the preset ideal points. The distances between assembly result No. 1 and assembly result No. 2 and the preset ideal points are 0.0068 and 0.0073, respectively, so the cleaning threshold of assembly result No. 1 is retained as the optimal cleaning threshold. Since the quality ratings of assembly results 3 to 5 are not as good as assembly result No. 1, they are all directly discarded. Finally, the cleaning threshold of assembly result No. 1 (assuming the length threshold TryLen2, assuming the coverage threshold TryCov3) is taken as the optimal cleaning threshold.

[0114] S2-7. Re-filter the contigs in the assembly results after polysome contamination removal according to the optimal cleaning threshold to obtain the assembly results after polysome contamination and background contamination removal.

[0115] S2-8. Use checkm to re-evaluate the obtained assembly results after removing polysomal contamination and background contamination. The evaluation contents include completeness and contamination.

[0116] In a specific embodiment, the method in the prior art (Wenshan Zheng et al., High-throughput, single-microbe genomics with strain resolution, applied to a human gutmicrobiome. Science 376, eabm1483 (2022). DOI: 10.1126 / science.abm1483) and the method provided in the above embodiment of the present invention (this method) are respectively used to remove contamination from the assembly results obtained from different sample types, and the results shown in Tables 2 and 3 are obtained. Among them, the high-quality Bin must meet the requirements of completeness greater than 90% and contamination less than 5%; the medium-quality Bin must meet the requirements of completeness greater than 50% and contamination less than 10%, and it is not repeated with the high quality. Table 3 also specifically compares the integrity and contamination obtained by evaluating the same Bin after contamination removal using different methods. The results show that in all sample types, after being processed by this method, more high-quality Bins are obtained than those in the article method, while the number of medium-quality Bins is not much different. Among all the Bins, after the contamination is removed by this method, the quality rating of the Bin is effectively improved, or the integrity is greatly improved at the cost of a smaller contamination degree, thereby obtaining a better assembly effect.

[0117] Table 2. Different methods were used to process different sample types.

[0118] Sample Type Decontamination methods High-quality Bin quantity Number of medium quality bins Lake Article Method 2 104 Lake This method 5 100 Intestinal samples Article Method 27 67 Intestinal samples This method 51 75 Intestinal & Vaginal Article Method 33 84 Intestinal & Vaginal This method 79 80 stool Article Method 44 153 stool This method 79 208

[0119] Table 3. Different methods for processing the same Bin

[0120] bin_ID Quality Rating_Article Method Quality Rating_This Method Completeness_Article Method Completeness_This method Pollution Degree_Article Method Pollution degree_This method 8_2813 High high 92.14573 96.09098 2.430556 2.776314 S_2913 High high 92.06731 94.47115 2.884615 3.125 8_2682 medium high 89.89247 95.9319 0.537634 0.537634 8_2615 medium high 88.8515 91.91239 1.965812 3.365385 8_2772 medium high 85.54505 93.24537 3.020134 2.684564 8_2815 medium medium 85.39485 85.69812 1.257432 1.382696 8_2861 medium high 84.76368 94.52736 2.873134 2.723881 8_2840 medium medium 84.66039 96.75423 1.842051 7.077626 8_2544 medium high 84.39873 91.0654 2.320675 4.276947 8_2606 medium high 84.12993 95.39436 1.308901 2.748691 8_2534 medium high 83.60595 95.65366 3.934325 2.67658

[0121] The embodiment of the present invention also provides a system for removing contamination in the results of assembling single-cell microbial genome species based on a classification feature gene set, comprising:

[0122] The module for removing polytopic contamination in Bin is used to align all SAGs in Bin to the original assembly result conitg of Bin to obtain the SAG-conitg matrix M; iterate the following process until the preset termination condition is met: for the matrix M or the matrix M obtained in the previous iteration G1 、M G2 Perform hierarchical clustering and perform clustering on conitg and matrices M, M in turn according to the hierarchical clustering results. G1 or M G2Perform binary division to obtain ContigGroup1, ContigGroup2, M1, and M2 respectively; calculate the classification probability P1 and P2 of SAG belonging to ContigGroup1 and ContigGroup2 respectively according to M1, M2, and the total number of reads in each SAG; classify SAG according to P1, P2, and the preset classification threshold, and take SAGs whose P1 or P2 meets the classification threshold as retained SAGs and assign them to SAGGroup1 or SAGGroup2 respectively, and take SAGs whose P1 and P2 do not meet the classification threshold as multi-cell SAGs and discard them; obtain the SAG-conitg matrix M of the retained SAGs belonging to SAGGroup1 or SAGGroup2 from the matrix M G1 、M G2 ; After the termination condition is met, the classification results of this round of iteration are invalidated, and the retained SAGs obtained before this round of iteration are merged and assembled to obtain the Bin without polytopic contamination and the assembly result without polytopic contamination;

[0123] A module for removing background pollution in Bins is used to evaluate the assembly results of Bins with polycellular pollution removed, obtain feature information of Bins, and determine classification information and reference feature sequences of Bins based on the feature information; wherein the feature information includes feature sequences, feature genes, and feature gene sets; set the length search boundary, length search step, coverage search boundary, and coverage search step of contigs in the assembly results with polycellular pollution removed, and then gradually search for the assumed length threshold and assumed coverage threshold of contigs; combine different assumed length thresholds and assumed coverage thresholds to form several cleaning thresholds, and filter contigs in the assembly results with polycellular pollution removed according to different cleaning thresholds to obtain retained contigs; evaluate the integrity and contamination of the assembly results with polycellular pollution removed under each cleaning threshold according to the reference feature sequence of Bins and the feature sequence of retained contigs; perform quality rating on the assembly results with polycellular pollution removed under each cleaning threshold using integrity and contamination as indicators, and select the best cleaning threshold according to the quality rating; re-filter contigs in the assembly results with polycellular pollution removed according to the best cleaning threshold to obtain the assembly results with polycellular pollution and background pollution removed.

[0124] An embodiment of the present invention also provides a device for removing contamination in the assembly results of single-cell microbial genome species based on a classification feature gene set, comprising a processor and a memory; the processor and the memory are connected via a communication bus; wherein the processor is used to call and execute a program stored in the memory; the memory is used to store a program, and the program is at least used to execute the method for removing contamination in the assembly results of single-cell microbial genome species based on a classification feature gene set.

[0125] The embodiments described above are merely descriptions of preferred implementations of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope of the present invention.

Claims

1. A method for removing contamination from the results of assembling single-cell microbial genome species based on a classification feature gene set, characterized in that: include: Remove multiple cell contamination and background contamination in the bin in turn; Among them, the methods for removing background pollution in the Bin include: Assembling a Bin with polysomal contamination removed, evaluating the obtained assembly result with polysomal contamination removed, obtaining a characteristic gene set of the Bin, and determining classification information and a reference characteristic sequence of the Bin based on the characteristic gene set; wherein the characteristic gene set is composed of characteristic genes with characteristic sequences; Set the contig length search boundary, length search step, coverage search boundary, and coverage search step in the assembly results after removing polysome contamination, and then gradually search for the contig's hypothetical length threshold and hypothetical coverage threshold; Combine different hypothetical length thresholds and hypothetical coverage thresholds to form several cleaning thresholds, and filter the contigs with characteristic genes in the assembly results after removing polycellular contamination according to different cleaning thresholds to obtain the retained contigs and characteristic sequences distributed in the retained contigs; evaluate the completeness and contamination of the assembly results after removing polycellular contamination under each cleaning threshold according to the reference characteristic sequence of the Bin and the characteristic sequence of the retained contig; use the completeness and contamination as indicators to rate the quality of the assembly results after removing polycellular contamination under each cleaning threshold, and select the best cleaning threshold according to the quality rating; Re-filter the contigs in the assembly results after polycellular contamination removal according to the optimal cleaning threshold to obtain the assembly results after polycellular contamination and background contamination removal; in, The method for evaluating the integrity and contamination of the assembly result after removing polysome contamination under the cleaning threshold according to the reference characteristic sequence of Bin and the characteristic sequence of retained contig includes: Determine the number of characteristic sequences to be retained in the contig and the number of characteristic sequences that appear repeatedly in the contig; According to the ratio of the number of characteristic sequences in the retained contig to the number of reference characteristic sequences in the Bin, the completeness of the assembly result after removing polysome contamination under the cleaning threshold is determined; The contamination degree of the assembly result after removing polysome contamination under the cleaning threshold is determined according to the ratio of the number of species of repeated characteristic sequences in the retained contig to the number of species of the reference characteristic sequences in the Bin.

2. The method according to claim 1, characterized in that The method for removing polycellular contamination in Bin comprises: Align all SAGs in Bin to the original assembly result conitg of Bin to obtain the SAG-conitg matrix M; Iterate the following process until the preset termination condition is met: for the matrix M or the matrix M obtained in the previous iteration G1 、M G2 Perform hierarchical clustering and perform clustering on conitg and matrices M, M in turn according to the hierarchical clustering results. G1 or M G2 Perform binary division to obtain ContigGroup1, ContigGroup2, M1, and M2 respectively; calculate the classification probability P1 and P2 of SAG belonging to ContigGroup1 and ContigGroup2 respectively according to M1, M2, and the total number of reads in each SAG; classify SAG according to P1, P2, and the preset classification threshold, and take SAGs whose P1 or P2 meets the classification threshold as retained SAGs and assign them to SAGGroup1 or SAGGroup2 respectively, and take SAGs whose P1 and P2 do not meet the classification threshold as multi-cell SAGs and discard them; obtain the SAG-conitg matrix M of the retained SAGs belonging to SAGGroup1 or SAGGroup2 from the matrix M G1 、M G2 ; After the termination condition is met, the classification results of this round of iteration are invalidated, and the retained SAGs obtained before this round of iteration are merged and assembled to obtain the Bin without polytopic contamination and the assembly result without polytopic contamination.

3. The method according to claim 2, characterized in that The termination condition is any one of the following (1) to (3): (1) The ratio of the number of retained SAGs in SAGGroup1 and SAGGroup2 to the number of all SAGs in the Bin is less than the preset retention threshold; (2) The number of reserved SAGs in SAGGroup1 or SAGGroup2 is less than or equal to the preset number threshold; (3) The difference between the number of retained SAGs in SAGGroup1 or SAGGroup2 and the number of all SAGs in the Bin is less than or equal to the preset difference threshold.

4. The method according to claim 1, characterized in that The search boundary of the contig length is one-L of the contig length, and the value range of L is 5~20; One-fifth of the contig coverage, C, was used as the coverage search boundary of the contig, and the value range of C was 5-20.

5. The method according to claim 1, characterized in that The quality of the assembly results after removing polytopic contamination at each cleaning threshold is rated using completeness and contamination as indicators. The method for selecting the best cleaning threshold according to the quality rating includes: According to the preset classification conditions, the quality of the assembly results after removing polysomal contamination at each cleaning threshold is graded with completeness and contamination as indicators, and the quality rating includes high quality, medium quality and low quality; Taking completeness and contamination as indicators, the distances between the assembly results after polytopic contamination removal at each cleaning threshold and the preset ideal points, high-quality points, and medium-quality points are calculated. The cleaning threshold of the assembly result with the smallest distance to the ideal point and the highest quality rating for removing polytopic contamination is selected as the optimal cleaning threshold.

6. The method according to claim 5, characterized in that The distances between the assembly results after polytopic contamination removal under the cleaning thresholds and the preset ideal points, high-quality points, and medium-quality points are calculated as follows: When the cleaning threshold has the assumed length threshold when the search step is i and the assumed coverage threshold when the search step is j, the completeness of the assembly result after removing polytopic contamination under this cleaning threshold is Completeness i,j , Contamination i,j ; The distance between the assembly result after removing polycellular contamination under the cleaning threshold and the preset ideal point is DistanceIdeal i,j for: ； The distance between the assembly result after removing polycellular contamination under this cleaning threshold and the preset high-quality point is DistanceHigh i,j for: ； The distance between the assembly result after removing polycellular contamination under this cleaning threshold and the preset medium mass point DistanceMedium i,j for: 。 7. The method according to claim 1, characterized in that Also includes: Use checkm to re-evaluate the assembly results obtained after removing polysomal contamination and background contamination. The evaluation content includes completeness and contamination.

8. A system for removing contamination from the results of assembling single-cell microbial genome species based on a classification feature gene set, characterized in that: include: The module for removing polytopic contamination in Bin is used to iterate the following process for the input Bin until the preset termination condition is met: all SAGs in Bin or the retained SAGs in the previous iteration are assigned to SAGGroup1 or SAGGroup2, respectively, and aligned to conitg in the original assembly result of Bin to obtain the SAG-conitg matrix M; hierarchical clustering is performed on the matrix M, and conitg and the matrix M are divided into two parts in turn according to the hierarchical clustering results to obtain ContigGroup1, ContigGroup2, M1, and M2 respectively; according to the total number of reads in M1, M2, and each SAG, , calculate the classification possibilities P1 and P2 of SAG belonging to ContigGroup1 and ContigGroup2 respectively; classify SAG according to P1, P2 and the preset classification threshold, take SAG whose P1 or P2 meets the classification threshold as the reserved SAG and assign it to SAGGroup1 or SAGGroup2 respectively, and take SAG whose P1 and P2 do not meet the classification threshold as polymorphic SAG and discard it; after the termination condition is met, invalidate the classification result of this round of iteration, merge and assemble the reserved SAG obtained before this round of iteration, and obtain the Bin without polymorphic contamination and the assembly result without polymorphic contamination; A module for removing background pollution in Bins is used to evaluate the assembly results of Bins after polycellular pollution is removed, obtain the characteristic information of Bins, and determine the classification information and reference characteristic sequence of Bins based on the characteristic information; wherein, the characteristic information includes characteristic sequence, characteristic gene, and characteristic gene set; set the length search boundary, length search step, coverage search boundary, and coverage search step of contigs in the assembly results after polycellular pollution is removed, and then gradually search the assumed length threshold and assumed coverage threshold of contigs; combine different assumed length thresholds and assumed coverage thresholds to form several cleaning thresholds, and filter the contigs in the assembly results after polycellular pollution is removed according to different cleaning thresholds to obtain retained contigs; evaluate the integrity and contamination of the assembly results after polycellular pollution is removed under each cleaning threshold according to the reference characteristic sequence of Bin and the characteristic sequence of retained contigs; perform quality rating on the assembly results after polycellular pollution is removed under each cleaning threshold using integrity and contamination as indicators, and select the best cleaning threshold according to the quality rating; re-filter the contigs in the assembly results after polycellular pollution is removed according to the best cleaning threshold to obtain the assembly results after polycellular pollution and background pollution are removed; in, The method for evaluating the integrity and contamination of the assembly result after removing polysome contamination under the cleaning threshold according to the reference characteristic sequence of Bin and the characteristic sequence of retained contig includes: Determine the number of characteristic sequences to be retained in the contig and the number of characteristic sequences that appear repeatedly in the contig; According to the ratio of the number of characteristic sequences in the retained contig to the number of reference characteristic sequences in the Bin, the completeness of the assembly result after removing polysome contamination under the cleaning threshold is determined; The contamination degree of the assembly result after removing polysome contamination under the cleaning threshold is determined according to the ratio of the number of species of repeated characteristic sequences in the retained contig to the number of species of the reference characteristic sequences in the Bin.

9. A device for removing contamination from the results of assembling single-cell microbial genome species based on a classification feature gene set, characterized in that: It comprises a processor and a memory; the processor and the memory are connected via a communication bus; wherein the processor is used to call and execute a program stored in the memory; the memory is used to store a program, and the program is used to execute at least the method for removing contamination in the assembly results of single-cell microbial genome species based on a classification feature gene set as described in any one of claims 1 to 7.

Images