A method for analyzing single-cell microbial genomic sequencing data

Abstract

The present invention provides a method for analyzing single-cell microbial genome sequencing data. This method uses an optimized method for assembling species genome sequences to obtain species-level bins free of multi-cell contamination and background contamination from single-cell microbial genomes, as well as corresponding species assembly results. Then, according to the species SNP characteristics, strain-level bins are distinguished from the species-level bins, and corresponding strain assembly results are obtained after assembly. Compared with the prior art, the present invention reduces manual intervention, improves the convenience and efficiency of operation, effectively controls the contamination degree of species bins, ensures the accuracy of assembly results, significantly reduces memory consumption in the case of a high number of SAGs, and improves the operation speed by calculating the similarity matrix in groups. It can dynamically select or adjust the similarity threshold to adapt to different sample types or library qualities.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, specifically relates to the field of single-cell microbial genome sequencing, and more specifically relates to a method for analyzing single-cell microbial genome sequencing data. 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 important roles in food production (such as fermentation), medicine (such as the production of antibiotics), environmental protection (such as sewage treatment), etc. Different microorganisms perform different biological functions through their different gene sequences. Through sequencing means, it is currently possible to detect microbial genes on a large scale, so as to deeply understand their action mechanisms and significance.

[0003] Single-cell microbial genome technology (microbe-seq) is a technical means that uses single-cell sequencing to detect the entire gene sequence of a single microorganism with high throughput. 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 combines a variety of droplet microfluidic operation technologies with bioinformatics analysis means, and can obtain the genomic information of thousands of single-cell microorganisms from complex microbial communities without culturing. It 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 of single microorganisms, detect the interaction 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, and can not only quickly identify different strains, but also provide more refined biological characteristic analysis based on the characteristic genes of each strain.

[0004] In the above article, a method for assembling species-level genome sequences based on single-cell microbial genome sequencing data was described. Specifically, a set of SAGs (single-amplified genome, referring to the droplets containing microorganisms and their amplified gene sequences) is called a bin, and initially each SAG forms a separate bin. Then, the following process is repeatedly executed: using SPAde to assemble all reads of the SAGs within the bin, and then using sourmash to calculate the similarity between the assembly results of the bins. Bins with a similarity greater than the similarity threshold of 5% are merged; until a certain number of rounds of cycling, a species-level genome sequence without filtered contamination is obtained.

[0005] However, the above method for assembling species-level genome sequences still has the following problems:

[0006] (1) The analysis process is cumbersome, requiring manual specification of the number of assembly rounds, and each round of analysis requires manual submission of tasks;

[0007] (2) It is impossible to control the contamination level within the bin. If the contamination level is too high, it is difficult to remove the contamination from the assembly results;

[0008] (3) If the number of SAGs contained in the sequencing data is high, the memory consumption is too large, slowing down the operation speed;

[0009] (4) When calculating the bin similarity, too much contamination may be included; and it is impossible to dynamically select the similarity threshold, and different sample types or library qualities require manual adjustment and modification of parameters. Summary of the Invention

[0010] To solve the above problems, the present invention provides an optimized method for analyzing single-cell microbial genome sequencing data. Using the optimized method for assembling species genome sequences, species-level bins free of multi-cell contamination and background contamination and corresponding species assembly results are obtained from single-cell microbial genomes. Then, strain-level bins are distinguished from the species-level bins according to the species SNP characteristics, and the corresponding strain assembly results are obtained after assembly.

[0011] The technical solution adopted by the present invention is: a method for analyzing single-cell microbial genome sequencing data, including:

[0012] S1. Obtain a single-cell microbial genome sequencing library, generate a bin for each SAG, and initialize the bin list;

[0013] S2. Assemble the bins according to the k-mer set for the current assembly round, and evaluate the integrity and contamination level of the assembly results of the bins to obtain assembly information;

[0014] S3. If the assembly information meets the assembly termination condition, execute step S6; otherwise, execute step S4.

[0015] S4. Group all bins, perform hierarchical clustering on the bins within each group according to the similarity between pairs of bins and the species annotation information; then split the hierarchical clustering results into one or more clusters according to a dynamically adjusted similarity threshold until the number of clusters within the group or the similarity threshold meets the preset splitting termination condition; merge the bins in a single cluster within the group into a new bin, and then summarize the new bins of all groups to obtain the bin list for the next assembly round.

[0016] S5. If the bin list for the next assembly round has the same number of bins as the bin list for the current assembly round, execute step S6; otherwise, update the k-mer and then repeat steps S2 to S4.

[0017] S6. Remove polyploid contamination and background contamination within the bins in sequence to obtain species bins, and after assembly, obtain the corresponding species assembly results.

[0018] S7. Cluster the SAGs within the species bins according to the species SNP characteristics to obtain the corresponding strain bins; remove the background contamination within the strain bins, and after assembly, obtain the corresponding strain assembly results.

[0019] Preferably, in step S2, the length of the k-mer is K, and the value range of K is 21 - 61 and K is an odd number.

[0020] Preferably, in step S3, the assembly termination condition is any one of the following (1) - (2):

[0021] (1) The proportion of the number of SAGs corresponding to bins with a completeness greater than 0.5 in the total number of SAGs is greater than the preset completeness threshold.

[0022] (2) The number of assembly rounds is greater than the preset maximum number of assembly rounds.

[0023] Preferably, step S4 includes:

[0024] S4-1. Filter the assembly results of each bin to obtain the filtered assembly results of each bin.

[0025] S4-2. Divide all bins into multiple groups, and calculate the similarity between pairs of bins within each group according to the filtered assembly results of each bin and the k-mer set in the current assembly round.

[0026] S4-3. Modify the similarity according to the integrity and contamination degree of the assembly results of each bin, and perform hierarchical clustering on the bins within the group according to the species annotation information and the modified similarity to obtain the hierarchical clustering results of the group;

[0027] S4-4. Split the hierarchical clustering results of the group into one or more clusters according to the similarity threshold, and calculate the ratio R of the number of all clusters within the group to the number of all bins;

[0028] S4-5. If R is less than the ratio threshold, execute step S4-6; otherwise, lower the similarity threshold. If the lowered similarity threshold is less than the minimum similarity, execute step S4-6; otherwise, repeat step S4-4;

[0029] S4-6. Merge the bins in a single cluster within the group into a new bin, and then summarize the new bins of all groups to update the bin list.

[0030] Preferably, in step S4-1, the method for filtering the assembly results of each bin includes: removing the contigs with lengths less than the length threshold from the assembly results of the bin; the length threshold is one-L of the Contig N50 length in the assembly results.

[0031] Preferably, the value range of L is 5 to 20.

[0032] Preferably, in step S4-3, the method for modifying the similarity includes: if the integrity of a certain bin within the group is not less than 0.9 or the contamination degree is not less than 0.5, modify the similarity between this bin and any bin within the group to 0.

[0033] Preferably, step S4-6 further includes: if the number of bins in a certain cluster within the group is greater than the merging threshold, split this cluster into multiple small clusters, and then merge the bins in a single small cluster into a new bin; wherein, the number of bins in the small cluster is not greater than the merging threshold.

[0034] Preferably, in step S5, the method for updating the k-mer includes: updating the k-mer according to the current assembly round and the k-mer decay step size.

[0035] Preferably, in step S6, the method for removing polyploid contamination within the bin includes:

[0036] S6-1. Align all SAGs within the bin to the conitg of the bin assembly results in step S2 to obtain the SAG-conitg matrix M;

[0037] S6-2. Iterate the following process until a preset termination condition is met: for matrix M or the matrix M obtained in the previous iteration G1 , M G2 perform hierarchical clustering, and then perform binary division on conitg and matrices M, M G1 or M G2 respectively to obtain ContigGroup1, ContigGroup2, M1, and M2; calculate the classification probabilities 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 the SAG according to P1, P2, and a preset classification threshold, and regard the SAG whose P1 or P2 meets the classification threshold as a retained SAG and assign it to SAGGroup1 or SAGGroup2 respectively, and regard the SAG whose P1 and P2 do not meet the classification threshold as a multi-cell SAG and discard it; obtain the SAG-conitg matrices M G1 , M G2 ;

[0038] S6-3. After the termination condition is met, invalidate the classification results of this iteration, merge and assemble the retained SAGs obtained before this iteration to obtain bins without multi-cell contamination and an assembly result without multi-cell contamination.

[0039] Preferably, in step S6, the method for removing background contamination in the bin includes:

[0040] S6-4. Evaluate the assembly result without multi-cell contamination obtained in step S6-3 to obtain a characteristic gene set of the bin, and determine the classification information and 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;

[0041] S6-5. Set the length search boundary, length search step size, coverage search boundary, and coverage search step size of the contigs in the assembly result without multi-cell contamination, and then gradually search for the assumed length threshold and assumed coverage threshold of the contigs;

[0042] S6-6. Combine different assumed length thresholds and assumed coverage thresholds to form several cleaning thresholds. Respectively filter the contigs with characteristic genes in the assembly result after removing polyploid contamination according to different cleaning thresholds to obtain the retained contigs and the characteristic sequences distributed in the retained contigs. Evaluate the integrity and contamination degree of the assembly result after removing polyploid contamination under each cleaning threshold according to the reference characteristic sequences of the bins and the characteristic sequences of the retained contigs. Perform quality rating on the assembly results after removing polyploid contamination under each cleaning threshold with integrity and contamination degree as indicators, and select the best cleaning threshold according to the quality rating.

[0043] S6-7. Re-filter the contigs in the assembly result after removing polyploid contamination according to the best cleaning threshold to obtain the species bins and the corresponding species assembly results.

[0044] Preferably, in step S6-6, the method for evaluating the integrity and contamination degree of the assembly result after removing polyploid contamination under this cleaning threshold according to the reference characteristic sequences of the bins and the characteristic sequences of the retained contigs includes:

[0045] Determine the number of types of characteristic sequences in the retained contigs and the number of types of characteristic sequences that repeatedly appear in the retained contigs;

[0046] Determine the integrity of the assembly result after removing polyploid contamination under this cleaning threshold according to the ratio of the number of types of characteristic sequences in the retained contigs to the number of types of reference characteristic sequences of the bins;

[0047] Determine the contamination degree of the assembly result after removing polyploid contamination under this cleaning threshold according to the ratio of the number of types of characteristic sequences that repeatedly appear in the retained contigs to the number of types of reference characteristic sequences of the bins.

[0048] Preferably, in step S7, the method for clustering the SAGs in the species bins according to the species SNP characteristics to obtain the corresponding strain bins includes:

[0049] S7-1. Align each SAG to the species genome to obtain the SNPs of each SAG relative to the species genome;

[0050] S7-2. Combine the SNPs of all SAGs to obtain a preliminary SNP matrix, filter the SNPs in the preliminary SNP matrix according to preset filtering conditions to obtain a characteristic matrix containing characteristic SNPs;

[0051] S7-3. Hierarchically cluster the SAGs based on the characteristic SNPs, divide the hierarchical clustering results according to different numbers of strains, evaluate each division result using the WSS score as an evaluation index, and determine the optimal number of strains according to the evaluation results;

[0052] S7-4. Divide the hierarchical clustering results according to the optimal number of strains to obtain the SAG set of each strain;

[0053] S7-5. Calculate the coverage of each SAG in the SAG set of each strain when aligned to the species genome;

[0054] S7-6. According to the alignment rate, coverage, relative coverage ratio, and effective coverage ratio of each SAG, use the greedy algorithm to screen out the SAGs for strain assembly from the SAG set to obtain the SAG assembly set, that is, the corresponding strain bin.

[0055] Preferably, step S7-6 includes:

[0056] S7-6-1. Set the initial state of the SAG assembly set. In the first round of the loop, select the p-th SAG with the largest coverage in the current SAG set and add it to the SAG assembly set, merge the coverage regions, and obtain the SAG assembly set in the first round;

[0057] Or in the current non-first round of the loop, select the p-th SAG with the largest effective coverage ratio in the current SAG set and add it to the SAG assembly set obtained in S7-6-3 in the previous round of the loop, merge the coverage regions, and obtain the SAG assembly set in the current round;

[0058] S7-6-2. Calculate the first overlapping length of the coverage region of each current q-th SAG not included in the SAG assembly set and the coverage region of the p-th SAG respectively, and then calculate the relative coverage ratio of each current q-th SAG not included in the SAG assembly set relative to the p-th SAG according to the first overlapping length, the set length of the coverage region of the p-th SAG, and the read length of the p-th SAG not aligned to the species gene sequence;

[0059] S7-6-3. Select the q max -th SAG with the largest relative coverage ratio and add it to the SAG assembly set obtained in S7-6-1, merge the coverage regions, and iteratively obtain the SAG assembly set in the current round;

[0060] S7-6-4. Calculate the second overlapping length between the coverage region of each q'-numbered SAG that has not been incorporated into the SAG assembly set and the coverage region of the SAG assembly set obtained in S7-6-3 respectively. Then, calculate the effective coverage ratio of each q'-numbered SAG that has not been incorporated into the SAG assembly set relative to the p-numbered SAG respectively according to the second overlapping length, the set length of the coverage region of the SAG assembly set obtained in S7-6-3, and the read length of the p-numbered SAG that has not been aligned to the species gene sequence.

[0061] S7-6-5. Loop through steps S7-6-1 to S7-6-4 until the preset assembly conditions are met to obtain the SAG assembly set, that is, the corresponding strain bin. Among them, the preset assembly conditions are set according to the total coverage length of the SAG assembly set, the number of SAGs in the SAG assembly set, and the alignment rate of the SAGs incorporated into the SAG assembly set.

[0062] Advantages of the present invention: Compared with the method for assembling species-level genomic sequences provided in the prior art, the present invention 1) reduces human intervention and improves the convenience and efficiency of operation by automatically calculating the number of assembly rounds and simplifying the process; 2) can effectively control the contamination degree of species bins and ensure the accuracy of the assembly results; 3) significantly reduces the memory consumption in the case of a high number of SAGs and improves the operation speed by calculating the similarity matrix in groups, solving the problems existing in large-scale data processing; 4) can dynamically select or adjust the similarity threshold to adapt to different sample types or library qualities, enhancing the flexibility and applicability of the method. Description of the Drawings

[0063] Figure 1 It is a flowchart of a method for analyzing single-cell microbial genomic sequencing data provided by an embodiment of the present invention.

[0064] Figure 2 It is a flowchart of co-assembling the species genomic sequence of a single-cell microbial genomic sequencing library using method steps S1 to S5 of the embodiment of the present invention. Detailed Embodiments

[0065] The following specific embodiments illustrate the implementation manners 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 implementation manners. 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, without conflict, the following embodiments and the features in the embodiments can be combined with each other.

[0066] An embodiment of the present invention provides a method for analyzing single-cell microbial genome sequencing data. After obtaining a single-cell microbial genome sequencing library, this method initializes the bin list. Then, according to the set k-mer, each bin is assembled and the integrity and contamination degree of the assembly result are evaluated. The assembly result is filtered to remove shorter contigs to control the contamination degree. Then, the bins are grouped and the similarity between each pair of bins within the group is calculated. Based on the species annotation information and similarity, the similar bins are hierarchically clustered, and the hierarchical clustering result is split into one or more clusters according to the similarity threshold. During this process, the similarity threshold is dynamically adjusted to adapt to different sample types or library qualities until the number of clusters within the group or the similarity threshold meets the preset splitting termination condition. Next, the similar bins are merged into new bins according to the clusters, and the k-mer value is updated for the next round of assembly. The above assembly process is repeated until the bin list is no longer updated or the assembly meets the preset assembly termination condition, and the bins in the list are used for further assembly of the subsequent species genome. The above steps reduce memory consumption and improve the efficiency and accuracy in the process of assembling the single-cell microbial species genome sequence through an automated process, controlling the contamination degree, and dynamically adjusting the similarity threshold. Subsequently, the obtained bins are successively removed of polyploid contamination and background contamination to obtain species-level bins, and the corresponding species assembly results are obtained after assembly. Finally, according to the species SNP characteristics, the SAGs within the species-level bins are clustered to obtain the corresponding strain-level bins; the background contamination within the strain bins is removed, and the corresponding strain assembly results are obtained after assembly. The specific steps of this method are as follows:

[0067] S1. Obtain a single-cell microbial genome sequencing library, and set the current assembly round counter m to 0, generate a bin for each SAG, and initialize the bin list:

[0068] ;

[0069] wherein,

[0070] SAG i is the i th SAG;

[0071] bin m,i is the m th bin in the i th assembly round;

[0072] binList m is the bin list in the m th assembly round;

[0073] SAGCount is the total number of SAGs in the sequencing library.

[0074] S2. Obtain bin m,i the corresponding fastq file fasta m,i , and then use SPAde. According to the kmer m set in the current m-th assembly round, bin m,i perform assembly, and use checkm to evaluate bin m,i the completeness of the assembly result of Completeness m,i 、 contamination Contamination m,i , and obtain assembly information. Among them, the length of k-mer is K, and the value range of K is 21 - 61 and K is an odd number.

[0075] S3. If the assembly information meets the assembly termination condition, execute step S6-1; otherwise, execute steps S4-1 to S4-6; the assembly termination condition is any one of the following (1) - (2):

[0076] (1) The proportion of the number of SAGs corresponding to bins with a completeness greater than 0.5 in the total number of SAGs is greater than the preset completeness threshold C;

[0077] (2) The assembly round number m is greater than the preset maximum assembly round number MaxM;

[0078] S4-1. Filter the assembly result of bin m,i , and remove the contigs in the assembly result of bin m,i that are less than one L-th of the ContigN50 length, and obtain bin m,i the filtered assembly result of FilteredFasta m,iAmong them, the value range of L is 5 to 20. In a preferred embodiment, L takes the value of 10. Contig N50 is an important statistical indicator in the field of genome sequencing and assembly, used to evaluate the quality of genome assembly. Specifically, Contig N50 refers to arranging all the contigs of the genome in descending order of length, and then accumulating the lengths of these contigs until the accumulated total length reaches half of the total length of all contigs. The length of the last contig accumulated is Contig N50. For example, if the size of a genome is 1M, after the reads obtained by sequencing are assembled into contigs, these contigs are arranged from long to short and then accumulated in turn until the accumulated length reaches 500k (i.e., 50% of 1M), then the length of the last contig accumulated is Contig N50. The larger the value of Contig N50, the better the continuity of the assembly, the more long fragments are assembled, and thus it is closer to the true genome structure. Therefore, Contig N50 is generally considered an important criterion for measuring the quality of genome assembly results. Therefore, in the present invention, according to the lengths of the contigs in the bin and the length of Contig N50, the assembly results of the bin are filtered, and the shorter contigs with poor assembly continuity or integrity are discarded, so as to control the contamination degree of the bin at a low level for subsequent analysis.

[0079] S4-2. Divide all bins into M groups according to N bins in each group; according to bin m,i the filtered assembly results FilteredFasta m,i and the k-mer m set in the current m-th assembly round bin m,i calculate the minhash features of bin m,i using the MinHash function of sourmash, and then calculate the bin m,ii (that is, any bin in the n-th group of the m-th assembly round except bin m,i ) similarity with

[0080]

[0081] ;

[0082] Among them,

[0083] Gm,n is the set of bins for the nth group of the mth assembly round;

[0084] mh m,n,i is the minhash feature of the ith bin in the nth group of the mth assembly round;

[0085] S m,n,i,ii is the similarity between the ith bin and the iith bin in the nth group of the mth assembly round.

[0086] In this step, all bins are divided into M groups with N bins in each group. For example, when N = 5000 and n = 2, all bins are divided into M groups with 5000 bins in each group. Then the 2nd group includes SAGs numbered from 5000 to 9999. Next, calculate the similarity between every two of the 5000 bins within the group, and then perform the subsequent merging steps in units of groups, namely steps S4-3 to S4-6. Compared with the prior art method of calculating the similarity of all bins together and then merging, the present invention adopts a grouping method to calculate the similarity matrix, which can reduce the memory consumption when the SAG number is high and improve the operation speed. Moreover, after multiple rounds of co-assembly, it has no impact on the assembly result.

[0087] S4-3. If the integrity of a certain bin within the group is not less than 0.9 or the contamination degree is not less than 0.5, then modify the similarity between this bin and any bin within the group to 0; according to the modified similarity, use the linkage function of scipy to perform hierarchical clustering on the bins within the group to obtain the hierarchical clustering result of the group;

[0088] ;

[0089] wherein,

[0090] C m,n,i is the assembly result of the ith bin in the nth group of the mth assembly round.

[0091] In this step, the purpose of modifying the similarity is to prevent excessive contamination from being incorporated by excluding bins with relatively high integrity (i.e., integrity ≥ 0.9) or relatively high contamination (i.e., contamination ≥ 0.5) after assembly from subsequent merging. For bins with already high integrity, further merging will significantly increase their contamination levels and the difficulty of subsequent decontamination. Bins with relatively high contamination may have aggregated multi-cellular SAGs or SAGs with high backgrounds, and further merging will also increase the difficulty of subsequent decontamination. Therefore, in this step, the similarity between a bin with relatively high integrity or relatively high contamination after assembly and any bin within the group is modified to 0, so that this bin will not be split into the same cluster as the remaining bins during subsequent hierarchical clustering and splitting, and thus will no longer be included in the subsequent merging process. In addition, the method of performing hierarchical clustering on the bins within the group in units of small groups can effectively control memory consumption and improve the operation speed.

[0092] S4-4. Set the current split round counter mm to 0, and set the similarity threshold for the current mm-th split round to SH mm ; According to the similarity threshold SH mm , use the fcluster function of scipy to split the hierarchical clustering results of the small group into one or more clusters, and calculate the ratio R of the number of all clusters within the small group to all bins. In this step, the similarity threshold can be dynamically adjusted according to the split results. Since different sequencing libraries have different coverage levels, and the coverage directly affects the similarity, it is necessary to dynamically adjust the similarity threshold for bin merging in order to 1) adapt to libraries of different plasmids; 2) be able to merge bins of some species with fewer SAGs. This solves the problem in the prior art that the similarity threshold cannot be dynamically selected, and the similarity threshold needs to be manually adjusted and modified for different sample types or library qualities.

[0093] S4-5. If R is less than the ratio threshold, then execute step S4-6; otherwise, according to the preset similarity step StepSH lower the similarity threshold SH mm , , to obtain the similarity threshold SH mm+1 for the next split round, i.e., the (mm + 1)-th split round. If the lowered similarity threshold SH mm+1 is less than the minimum similarity minSH , then execute step S4-6; otherwise, repeat step S4-4.

[0094] S4-6. Combine the bins in a single cluster within a group into a new bin, and then summarize the new bins of all groups to obtain the bin list for the next assembly round (the (m + 1)-th assembly round). binList m+1 ; where, if the number of bins in a cluster within a group is greater than the merging threshold maxC, then split the cluster into multiple small clusters, and then combine the bins in a single small cluster into a new bin, and the number of bins in the small cluster is not greater than the merging threshold maxC. For example, the first group in the first assembly round has 5000 bins. After hierarchical clustering, and then according to the similarity threshold SH 1 set in the first splitting round, 3000 clusters are obtained. And the ratio R of the number of 3000 clusters to the number of 5000 bins is 0.6, which is less than the ratio threshold. Then combine multiple bins within a single cluster into a new bin. So, a total of 3000 new bins are obtained from the 3000 clusters in the first group. These 3000 new bins will enter the second assembly round together with the new bins obtained from other groups.

[0095] S5. If the bin list for the next assembly round (the (m + 1)-th assembly round) has the same number of bins as the bin list for the current assembly round (the m-th assembly round), that is binList m+1 = binList m , then execute step S6-1. Otherwise, update m according to the current assembly round number StepkMer and the k-mer decay step size kMer m to obtain kMer m+1 for the (m + 1)-th assembly round, and then repeat steps S2 to S4;

[0096] ;

[0097] where

[0098] StepkMerLen is the number of round intervals required for each decay; StepkMer is the step size for each decay, and the decay step size must be an odd number.

[0099] In a specific embodiment, the above steps S1 to S5 are used to perform co-assembly of the species genome on a single-cell microbial genome sequencing library with 21914 SAGs. Part of the assembly results are shown in Table 1. Among them, the specific parameter settings are as follows:

[0100] The length of the initial k-mer is set to 51;

[0101] The integrity threshold C (the proportion of the number of SAGs corresponding to bins with a completeness greater than 0.5 in the total number of SAGs) is set to 0.75;

[0102] The maximum number of assembly rounds MaxM is set to 20;

[0103] The length threshold is set to one-tenth of the Contig N50 length in the assembly result;

[0104] For each round of bins, they are divided into groups of 8000 bins per group;

[0105] The similarity threshold SHmm is set to 0.25;

[0106] The similarity step size StepSH is set to 0.01;

[0107] The proportion threshold (the ratio of the number of all clusters within a group to the number of all bins) is set to 0.75;

[0108] The merging threshold maxC (the maximum number of bins in a certain cluster within a group) is set to 15;

[0109] The k-mer decay step size StepkMer is set to 5;

[0110] The number of rounds interval StepkMerLen required for each decay is set to 5.

[0111] Table 1. Some bins that meet the high assembly quality standards in the co-assembly results of the genomes of single-cell microbial genomic sequencing libraries

[0112] bin number Completeness Contamination Genome size N50 8_3079 100 29.36782 4354587 32055 8_3100 100 11.45801 2866009 14853 8_3271 100 5.46595 2687405 25206 8_3499 100 37.22571 8003849 12305 8_3500 100 30.0627 6864230 20754 8_3145 99.91948 3.587963 3043472 28843 8_3148 99.90338 2.586188 3025084 36201 8_3098 99.9002 5.588822 2560968 32014 8_3033 99.60145 21.82574 4332449 15661

[0113] Among them, N50 is the Contig N50 length of the bin.

[0114] In this co-assembly, 21914 SAGs were assembled into 3538 species-level bins. Among them, there are 288 species-level bins with a completeness greater than 50. Their average contamination level is 9.33, and the highest contamination level is 72. And there are 87 species-level bins that meet the high assembly quality standards. Table 1 shows the assembly information of some species-level bins, such as bin 8_3271, 8_3145, 8_3148, 8_3098, which have better species genome assembly effects.

[0115] S6. Sequentially remove the multi-cell contamination and background contamination within the bins obtained in step S5 to obtain species bins, and after assembly, obtain the corresponding species assembly results.

[0116] S7. Cluster the SAGs within the species bin according to the species SNP characteristics to obtain the corresponding strain bins; remove the background contamination within the strain bins and obtain the corresponding strain assembly results after assembly.

[0117] In a specific implementation, step S6 specifically includes the following steps:

[0118] S6-1. For each bin in the bin list obtained in step S5, align all the SAGs within the bin to the assembly result conitg of the bin in step S2 to obtain the SAG-conitg matrix M.

[0119] S6-2. Iterate the following process until the preset termination condition is met: perform hierarchical clustering on matrix M or the matrix M obtained in the previous iteration G1 、M G2 ; perform binary division on conitg and matrices M, M G1 or M G2 respectively to obtain ContigGroup1, ContigGroup2, M1, and M2; calculate the classification probabilities P1 and P2 of the SAGs belonging to ContigGroup1 and ContigGroup2 respectively according to M1, M2, and the total number of reads within each SAG; classify the SAGs according to P1, P2, and the preset classification threshold, and use the SAGs whose P1 or P2 meets the classification threshold as the retained SAGs and assign them to SAGGroup1 or SAGGroup2 respectively, and discard the SAGs whose P1 and P2 do not meet the classification threshold as multi-cellular SAGs; respectively obtain the SAG-conitg matrices M G1 、M G2 of the retained SAGs belonging to SAGGroup1 or SAGGroup2 from matrix M. In some embodiments, the termination condition is any one of the following (1) to (3):

[0120] (1) The proportion of the number of retained SAGs within SAGGroup1 and SAGGroup2 to the total number of all SAGs within the bin is less than the preset retention threshold;

[0121] (2) The number of retained SAGs within SAGGroup1 or SAGGroup2 is less than or equal to the preset quantity threshold;

[0122] (3) The difference between the number of retained SAGs within SAGGroup1 or SAGGroup2 and the total number of all SAGs within the bin is less than or equal to the preset difference threshold.

[0123] The method for removing polyploid contamination is to align the reads in the SAG of the Bin to the contig in the original assembly result to obtain the SAG-contig matrix, where the matrix value indicates how many reads in the SAG are aligned to the contig. Then, hierarchical clustering is performed on this matrix and the contig is bisected. If the Bin is clustered due to polyploid contamination, a relatively clear bisection result will be obtained. Most of the non-polyploid SAGs in the Bin will be aligned to one of the groups, while the polyploid SAGs in the Bin will be aligned to both groups. Then, the polyploid SAGs that do not meet the classification threshold are removed from the two groups. In this way, polyploid SAGs can be removed. Reassembling the SAGs after removing polyploids and regrouping them can remove the contamination caused by polyploids in the result.

[0124] S6-3. After meeting the termination condition, invalidate the classification result of this round of iteration, merge and assemble the retained SAGs obtained before this round of iteration to obtain the bin with polyploid contamination removed and the assembly result with polyploid contamination removed.

[0125] S6-4. Evaluate the assembly result with polyploid contamination removed obtained in step S6-3 to obtain the characteristic gene set of the bin, and determine the classification information and 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.

[0126] S6-5. Set the length search boundary, length search step size, coverage search boundary, and coverage search step size of the contig in the assembly result with polyploid contamination removed, and then gradually search for the assumed length threshold and assumed coverage threshold of the contig. In some embodiments, one L-th of the contig length is used as the length search boundary of the contig, where L ranges from 5 to 20; one C-th of the contig coverage is used as the coverage search boundary of the contig, where C ranges from 5 to 20.

[0127] S6-6. Combine different assumed length thresholds and assumed coverage thresholds to form several cleaning thresholds, and filter the contigs with characteristic genes in the assembly result with polyploid contamination removed according to different cleaning thresholds to obtain the retained contigs and the characteristic sequences distributed in the retained contigs; evaluate the integrity and contamination degree of the assembly result with polyploid contamination removed under each cleaning threshold according to the reference characteristic sequence of the bin and the characteristic sequences of the retained contigs; perform quality rating on the assembly result with polyploid contamination removed under each cleaning threshold with integrity and contamination degree as indicators, and select the best cleaning threshold according to the quality rating.

[0128] In some embodiments, the method for evaluating the integrity and contamination degree of the assembly result after removing polyploid contamination at the cleaning threshold based on the reference feature sequence of the Bin and the feature sequence of the retained contig includes: determining the number of types of feature sequences in the retained contig and the number of types of feature sequences that repeatedly appear in the retained contig; determining the integrity of the assembly result after removing polyploid contamination at the cleaning threshold according to the ratio of the number of types of feature sequences in the retained contig to the number of types of reference feature sequences of the bin; determining the contamination degree of the assembly result after removing polyploid contamination at the cleaning threshold according to the ratio of the number of types of feature sequences that repeatedly appear in the retained contig to the number of types of reference feature sequences of the bin. Since in the sequence of a single-microorganism assembly result that is complete and pollution-free, the characteristic genes all appear and only appear once. Therefore, the present invention determines the integrity of the assembly result after removing polyploid contamination at the cleaning threshold according to the ratio of the number of types of feature sequences in the retained contig to the number of types of reference feature sequences of the bin, that is, the larger the ratio, the closer the number of types of feature sequences in the retained contig is to the number of types of reference feature sequences of the bin, indicating that the assembly result is more complete. The present invention determines the contamination degree of the assembly result after removing polyploid contamination at the cleaning threshold according to the ratio of the number of types of feature sequences that repeatedly appear in the retained contig to the number of types of reference feature sequences of the Bin, that is, the larger the ratio, the more the number of types of feature sequences that repeatedly appear in the retained contig, and the repeatedly appearing feature sequences may come from the contaminated contigs assembled from the background reads, indicating that the assembly result is more severely contaminated.

[0129] In some embodiments, the method for quality rating the assembly results after removing polyploid contamination at each cleaning threshold with integrity and contamination degree as indicators and selecting the optimal cleaning threshold according to the quality rating includes:

[0130] According to the preset grading conditions, with integrity and contamination degree as indicators, quality rating is performed on the assembly results after removing polyploid contamination at each cleaning threshold, and the quality rating includes high quality, medium quality, and low quality;

[0131] With integrity and contamination degree as indicators, calculate the distances between the assembly results after removing polyploid contamination at each cleaning threshold and the preset ideal point, high-quality point, and medium-quality point respectively;

[0132] Select the cleaning threshold of the assembly result after removing polyploid contamination that has the smallest distance from the ideal point and the highest quality rating as the optimal cleaning threshold.

[0133] In some embodiments, the calculation method for the distances between the assembly results after removing polyploid contamination at each cleaning threshold and the preset ideal point, high-quality point, and medium-quality point respectively is:

[0134] When the cleaning threshold has a hypothesis length threshold at a search step of i and a hypothesis coverage threshold at a search step of j, the integrity of the assembly result after removing polyploid contamination at this cleaning threshold is Completeness i,j and the contamination level is Contamination i,j ;

[0135] Then the distance DistanceIdeal between the assembly result after removing polyploid contamination at this cleaning threshold and the preset ideal point i,j is:

[0136] ;

[0137] The distance DistanceHigh between the assembly result after removing polyploid contamination at this cleaning threshold and the preset high-quality point i,j is:

[0138] ;

[0139] The distance DistanceMedium between the assembly result after removing polyploid contamination at this cleaning threshold and the preset medium-quality point i,j is:

[0140] .

[0141] S6-7. Re-filter the contigs in the assembly result after removing polyploid contamination according to the optimal cleaning threshold to obtain species bins and the corresponding species assembly results.

[0142] In this embodiment, step S6 is based on the occurrence pattern of the characteristic gene set in the species genome, and gradually searches for the optimal combination of the length threshold and the coverage threshold to filter the species assembly result of the single-cell microbial genome, so as to remove the contaminated contigs assembled from the background reads. This method not only overcomes the defect in the prior art that the length threshold cannot be automatically judged, but also does not need to assume that the coverage of the contigs of the species is normally distributed or multi-modal when determining the contig coverage threshold, thus realizing the high-throughput detection and automatic removal of polyploid contamination and background contamination in the single-cell genome assembly result, and obtaining the optimal contamination removal result as much as possible. For the specific content of step S6, refer to CN119068975 A.

[0143] In a specific implementation, step S7 specifically includes the following steps.

[0144] S7-1. Align each SAG to the species genome to obtain the SNPs of each SAG relative to the species genome.

[0145] S7-2. Combine all SNPs of the SAGs to obtain a preliminary SNP matrix, and filter the SNPs in the preliminary SNP matrix according to preset filtering conditions to obtain a characteristic matrix containing characteristic SNPs. In some embodiments, the preset filtering conditions are set according to the detection depth of SNPs at different loci, and / or the number of mutations of SNPs at different loci, and / or the genes to which the SNPs belong. SNPs with low detection depth in the sequencing results are likely to lead to inaccurate analysis results. Therefore, in one or more specific embodiments, the detection depth of SNPs at different loci is incorporated into the preset filtering conditions to remove loci with too low detection depth, such as SNP loci with a detection depth lower than 5, so as to ensure the reliability of the selected characteristic SNPs. Differences in the number of mutations of a certain locus SNP in the sequencing results may lead to misclassification of strains. Therefore, in one or more specific embodiments, the number of mutations of SNPs at different loci is incorporated into the preset filtering conditions to remove SNP loci with too few mutations, such as SNP loci with a mutation number less than 1%, so as to effectively identify mutations with biological significance and avoid misidentifying random mutations as characteristic SNPs.

[0146] S7-3. Perform hierarchical clustering on the SAGs according to the characteristic SNPs, divide the hierarchical clustering results according to different numbers of strains, evaluate each division result using the WSS score as an evaluation index, and determine the optimal number of strains according to the evaluation results.

[0147] In some embodiments, step S7-3 includes:

[0148] S7-3-1. Calculate the mean and standard deviation of each characteristic SNP in the characteristic matrix;

[0149] S7-3-2. Sort all the characteristic SNPs according to their standard deviations, select the top x SNP characteristics with decreasing standard deviations in turn to form an analysis matrix, and perform hierarchical clustering on the SAGs according to the analysis matrix; more specifically, in step S3-2, the parameter x is preferably not less than 500, more preferably not less than 3000, and further preferably 3000. In one or more specific embodiments, if the total number of characteristic SNPs in the characteristic matrix is less than 500, then all the characteristic SNPs are incorporated into the analysis matrix;

[0150] S7-3-3. Determine the maximum number of strains Q traversed by the hierarchical clustering result; in a preferred embodiment, the maximum number of strains Q is not greater than 20. Experiments have found that when the maximum number of strains P exceeds 20, the accuracy of the division result cannot be better improved;

[0151] S7-3-4. Divide the hierarchical clustering result according to n numbers of strains, where n is any integer between 1 and Q;

[0152] S7-3-5. Evaluate each division result using the WSS score. The calculation formula for the WSS score is as follows:

[0153] ,

[0154] where,

[0155] WSS n is the WSS score when divided into n numbers of strains;

[0156] FeatureCount is the number of SNP features in the analysis matrix;

[0157] Feature i is the value set of the i-th SNP feature in the analysis matrix;

[0158] ClusterCenter m,i is the feature center of the i-th SNP feature of the m-th strain when the analysis matrix is divided into n numbers of strains. The average value of the i-th SNP feature of the SAG set of the m-th strain is used as ClusterCenter m,i ;

[0159] P m is the penalty score of the m-th strain when the analysis matrix is divided into n numbers of strains. The P m is set according to the average value of each SNP feature in the analysis matrix, the number of SAGs in the analysis matrix, and the number of SAGs of the m-th strain in the analysis matrix.

[0160] S7-4. Divide the hierarchical clustering result according to the optimal number of strains to obtain the SAG set of each strain.

[0161] S7-5. Calculate the coverage of each SAG in the SAG set of each strain aligned to the species genome. In some embodiments, the coverage of each SAG aligned to the species genome is defined as:

[0162] ,

[0163] where,

[0164] R p is the coverage of the p-th SAG;

[0165] MappedCov p is the covered length of the p-th SAG on the species gene sequence;

[0166] UnmappedCov p is the reads length of the p-th SAG that is not aligned to the species gene sequence.

[0167] S7-6-1. Set the initialization state of the SAG assembly set. In the first round of loop, select the p-th SAG with the largest coverage in the current SAG set and add it to the SAG assembly set, merge the coverage areas, and obtain the SAG assembly set in the first round;

[0168] Or in the current non-first round of loop, select the p-th SAG with the largest effective coverage ratio in the current SAG set and add it to the SAG assembly set obtained in S7-6-3 in the previous round of loop, merge the coverage areas, and obtain the SAG assembly set in the current round.

[0169] S7-6-2. Calculate the first overlapping length between the coverage area of each current q-th SAG not included in the SAG assembly set and the coverage area of the p-th SAG respectively. Then, according to the first overlapping length, the set length of the coverage area of the p-th SAG, and the reads length of the p-th SAG that is not aligned to the species gene sequence, calculate the relative coverage ratio of each current q-th SAG not included in the SAG assembly set with respect to the p-th SAG. In some embodiments, the relative coverage ratio TR p, q is calculated by the formula:

[0170] ,

[0171] where,

[0172] D p,q is the first overlapping length between the set of the coverage area of each current q-th SAG not included in the SAG assembly set and the set of the coverage area of the p-th SAG;

[0173] C p is the set of the coverage area of the p-th SAG.

[0174] S7-6-3. Select the qmax-th SAG with the largest relative coverage ratio and add it to the SAG assembly set obtained in S7-6-1, merge the coverage areas, and iteratively obtain the SAG assembly set in the current round.

[0175] S7-6-4. Calculate the second overlapping length between the coverage area of each current q'-th SAG not included in the SAG assembly set and the coverage area of the SAG assembly set obtained in S7-6-3 respectively. Then, according to the second overlapping length, the set length of the coverage area of the SAG assembly set obtained in S7-6-3, and the reads length of the p-th SAG that is not aligned to the species gene sequence, calculate the effective coverage ratio of each current q'-th SAG not included in the SAG assembly set with respect to the p-th SAG. In some embodiments, the effective coverage ratio R q’, a is calculated by the formula:

[0176] ,

[0177] Among them,

[0178] D q’, a is the second superposition length of the set of coverage regions of each SAG numbered q' that has not been included in the SAG assembly set in the current a-th round and the set of coverage regions of the SAG assembly set obtained by S6-3;

[0179] C a is the set of coverage regions of the SAG assembly set obtained by S7-6-3 in the current a-th round.

[0180] S7-6-5. Repeat steps S7-6-1 to S7-6-4 until the preset assembly conditions are met to obtain the SAG assembly set, that is, the corresponding strain bin; among them, the preset assembly conditions are set according to the total coverage length of the SAG assembly set, the number of SAGs in the SAG assembly set, and the alignment rate of the SAGs included in the SAG assembly set. In some embodiments, the preset assembly conditions are any one of the following (1) to (4):

[0181] (1) The length of the set of coverage regions of the SAG assembly set accounts for more than 95% of the length of the species gene sequence;

[0182] (2) The number of SAGs in the SAG assembly set is greater than the set SAG number threshold, and the SAG number threshold is preferably 50;

[0183] (3) The alignment rate of the SAGs included in the SAG assembly set in the (a + 1)-th round to the species genome is less than the set alignment rate threshold, and the alignment rate threshold is preferably 50%;

[0184] (4) All SAGs are included in the SAG assembly set.

[0185] S7-7. Remove the background contamination in the strain bin and obtain the corresponding strain assembly result after assembly.

[0186] Step S7 first finely filters the SNP information of SAG based on preset filtering conditions such as the detection depth, number of mutations, and genotype diversity of different-site SNPs in step S7-2 to obtain a feature matrix containing characteristic SNPs. This preprocessing method can, on the one hand, more effectively process SAGs with unobvious SNP characteristics or low coverage, thereby improving the accuracy of subsequent hierarchical clustering and reducing the risk of incorrect clustering; on the other hand, it can efficiently focus on specific genes, further distinguish bacterial species with specific genes, and provide more accurate strain genotype information. Subsequently, when dividing the SAG hierarchical clustering results in step S7-3, an algorithm capable of automatically determining the optimal number of strains is introduced, which can accurately distinguish strains in a high-throughput manner based on single-cell microbial genome sequencing results, avoiding the subjective bias caused by manual judgment of the number of strains. Step S7-4 divides the hierarchical clustering results according to the obtained optimal number of strains, and calculates the coverage of each SAG in the SAG set of each strain after division in step S7-5. In step S7-6, during the process of screening SAGs for assembly, the present invention creatively selects the optimal SAG assembly set efficiently and accurately using the greedy algorithm based on the alignment rate, coverage, relative coverage ratio, and effective coverage ratio of each SAG, effectively improving the integrity and accuracy of the strain assembly results in step S7-7, reducing the contamination degree, reducing information loss during the strain assembly process, and thus obtaining high-quality strain genome sequences. For the specific content of step S7, refer to CN 118737269 A.

[0187] The embodiments described above are only used to describe the preferred embodiments of the present invention and do not limit the scope of the present invention. Without departing from the design spirit of the present invention, various deformations and improvements made by those of ordinary skill in the art to the technical solutions of the present invention shall fall within the protection scope of the present invention.

Claims

1. A method for analyzing single-cell microbial genome sequencing data, characterized in that: include: S1. Obtain a single-cell microbial genome sequencing library, generate a bin for each SAG, and initialize the bin list; S2. Assemble the bins according to the k-mer set in the current assembly round, evaluate the integrity and contamination of the bin assembly results, and obtain assembly information; S3. If the assembly information satisfies the assembly termination condition, execute step S6, otherwise execute step S4; S4. Group all bins, hierarchically cluster the bins in the group according to the similarity between the two bins and the species annotation information; then split the hierarchical clustering results into one or more clusters according to the dynamically adjusted similarity threshold, until the number of clusters in the group or the similarity threshold meets the preset split termination condition; Merge the bins in a single cluster within a group into a new bin, and then aggregate the new bins of all groups to obtain the bin list for the next assembly round; S5. If the bin list of the next assembly round has the same number of bins as the bin list of the current assembly round, execute step S6, otherwise update the k-mer and repeat steps S2 to S4; S6. Remove the polycellular contamination and background contamination in the bin in turn to obtain the species bin, and then assemble to obtain the corresponding species assembly result; S7. Cluster the SAGs within the species bin according to the species SNP characteristics to obtain the corresponding strain bin; Remove background contamination in the strain bin and obtain the corresponding strain assembly results after assembly; In step S2, the length of the k-mer is K, the value range of K is 21-61 and K is an odd number; In step S3, the assembly termination condition is any one of the following (1)-(2): (1) The ratio of the number of SAGs corresponding to bins with integrity greater than 0.5 to the total number of SAGs is greater than the preset integrity threshold; (2) The number of assembly rounds is greater than the preset maximum number of assembly rounds; Step S4 includes: S4-1. Filter the assembly results of each bin to obtain the filtered assembly results of each bin; S4-2. Divide all bins into multiple groups, and calculate the similarity between bins in each group based on the filtered assembly results of each bin and the k-mer set in the current assembly round; S4-3. Modify the similarity according to the completeness and contamination of the assembly results of each bin, perform hierarchical clustering on the bins in the group according to the species annotation information and the modified similarity, and obtain the hierarchical clustering results of the group; S4-4. Split the hierarchical clustering results of the group into one or more clusters according to the similarity threshold, and calculate the ratio R of the number of all clusters in the group to the number of all bins; S4-5. If R is less than the ratio threshold, execute step S4-6; otherwise, lower the similarity threshold, if the lowered similarity threshold is less than the minimum similarity, execute step S4-6, otherwise repeat step S4-4; S4-6. Merge the bins in a single cluster within a group into a new bin, then aggregate the new bins of all groups and update the bin list; In step S4-1, the method for filtering the assembly results of each bin includes: removing contigs in the bin assembly results whose length is less than a length threshold; the length threshold is one L of the length of Contig N50 in the assembly results; and the value range of L is 5 to 20.

2. The method according to claim 1, characterized in that In step S4-3, the method for modifying the similarity includes: if the integrity of a bin in the group is not less than 0.9 or the contamination is not less than 0.5, then modifying the similarity between the bin and any bin in the group to 0.

3. The method according to claim 1, characterized in that Step S4-6 also includes: if the number of bins in a cluster within the group is greater than the merging threshold, the cluster is split into multiple small clusters, and then the bins in a single small cluster are merged into a new bin; wherein the number of bins in the small cluster is not greater than the merging threshold.

4. The method according to claim 1, characterized in that In step S5, the method for updating the k-mer includes: updating the k-mer according to the current assembly round number and the k-mer decay step length.

5. The method according to claim 1, characterized in that In step S7, the method of clustering the SAGs in the species bin according to the species SNP characteristics to obtain the corresponding strain bin includes: S7-1. Align each SAG to the species genome to obtain the SNP of each SAG relative to the species genome; S7-2. Merge the SNPs of all SAGs to obtain a preliminary SNP matrix, filter the SNPs in the preliminary SNP matrix according to a preset filtering condition, and obtain a characteristic matrix containing characteristic SNPs; S7-3. Perform hierarchical clustering of SAGs according to characteristic SNPs, divide the hierarchical clustering results according to different numbers of strains, evaluate each division result using the WSS score as an evaluation index, and determine the optimal number of strains based on the evaluation results; S7-4. Divide the hierarchical clustering results according to the optimal number of strains to obtain the SAG set for each strain; S7-5. Calculate the coverage of each SAG in the SAG set of each strain aligned to the species genome; S7-6. According to the alignment rate, coverage, relative coverage ratio, and effective coverage ratio of each SAG, a greedy algorithm is used to screen out the SAG for strain assembly from the SAG set to obtain the SAG assembly set, i.e., the corresponding strain bin.

Images