Method for analyzing sequencing data of single-cell microbial genome

Abstract

The invention provides a method for analyzing single-cell microorganism genome sequencing data. According to the method, species-level bins with polycellular pollution and background pollution removed and corresponding species assembly results are obtained from a single-cell microbial genome by utilizing an optimized method for assembling a species genome sequence, then strain-level bins are distinguished from the species-level bins according to species SNP characteristics, and corresponding strain assembly results are obtained after assembly. Compared with the prior art, by automatically calculating the number of assembly rounds and simplifying the process, human intervention is reduced, and the operation convenience and efficiency are improved; the pollution degree of the species bin can be effectively controlled, and the accuracy of an assembly result is ensured; the similarity matrix is calculated in groups, so that the memory consumption under the condition of high SAG number is obviously reduced, and meanwhile, the operation speed is improved; the similarity threshold value can be dynamically selected or adjusted, so that the similarity threshold value can adapt to different sample types or library quality.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics, and specifically relates to the field of single-cell microbial genome sequencing, and more specifically 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 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. Science376, 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 efficiently distinguishes similar species based on the reads of microorganisms in each droplet in the single-cell microbial genome sequencing results, distinguishes strains through mutations in individual microorganisms, detects interactions between microorganisms and viruses, and discovers 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] In the above article, a method for assembling species-level genome sequences based on single-cell microbial genome sequencing data is described. Specifically, the collection of SAG (single-amplified genome, referring to the droplets containing microorganisms and their amplified gene sequences) is called a bin, and each SAG forms a separate bin in the initial state, and then the following process is executed in a loop: use SPAde to assemble all reads of the SAG in the bin, and then use sourmash to calculate the similarity between the assembly results of the bin, and merge the bins with similarity greater than the similarity threshold of 5%; until a certain number of cycles are repeated, the species-level genome sequence without filtering contamination is obtained.

[0005] However, the above-mentioned species-level genome sequence assembly methods still have the following problems: (1) The analysis process is cumbersome and requires manual specification of the number of assembly rounds, and each round of analysis requires manual submission of tasks; (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 result. (3) If the sequencing data contains a high number of SAGs, the memory consumption will be too large, slowing down the calculation speed; (4) When calculating bin similarity, too much contamination may be included; and the similarity threshold cannot be selected dynamically. Different sample types or library qualities require manual adjustment of parameters. Summary of the invention

[0006] In order to solve the above problems, the present invention provides an optimized method for analyzing single-cell microbial genome sequencing data, which uses the optimized method for assembling species genome sequences to obtain species-level bins and corresponding species assembly results from single-cell microbial genomes that have removed polycellular contamination and background contamination, and then distinguish strain-level bins from species-level bins based on species SNP characteristics, and obtain corresponding strain assembly results after assembly.

[0007] The technical solution adopted by the present invention is: a method for analyzing single-cell microbial genome sequencing data, comprising: 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 bins and species annotation information; split the hierarchical clustering results into one or more clusters according to the dynamically adjusted similarity threshold until the number of clusters or the similarity threshold in the group meets the preset split termination condition; merge the bins in a single cluster in 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; 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 the background contamination within the strain bin and obtain the corresponding strain assembly result after assembly.

[0008] Preferably, 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.

[0009] Preferably, 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.

[0010] Preferably, step S4 comprises: 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.

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

[0012] Preferably, the value range of L is 5~20.

[0013] Preferably, 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.

[0014] Preferably, 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.

[0015] Preferably, 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 attenuation step length.

[0016] Preferably, in step S6, the method for removing polycellular contamination in a bin includes: S6-1. Align all SAGs in the bin to the assembly result conitg of the bin in step S2 to obtain the SAG-conitg matrix M; S6-2. 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 ; S6-3. 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 bins without polytopic contamination and the assembly results without polytopic contamination.

[0017] Preferably, in step S6, the method for removing background contamination in the bin includes: S6-4. Evaluate the assembly result obtained in step S6-3 to remove polysome contamination, 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 having characteristic sequences; S6-5. Set the contig length search boundary, length search step, coverage search boundary, and coverage search step in the assembly result after removing polysome contamination, and then gradually search for the hypothetical length threshold and hypothetical coverage threshold of the contig; 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 results after polycystic contamination removal 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 polycystic contamination removal under each cleaning threshold according to the reference characteristic sequence of the bin and the characteristic sequence of the retained contig; perform quality rating of the assembly results after polycystic contamination removal under each cleaning threshold using completeness and contamination as indicators, and select the best cleaning threshold according to the quality rating; S6-7. Re-filter the contigs in the assembly results after removing polysome contamination according to the optimal cleaning threshold to obtain species bins and corresponding species assembly results.

[0018] Preferably, in step S6-6, 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.

[0019] Preferably, 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.

[0020] Preferably, step S7-6 includes: S7-6-1. Set the initialization state of the SAG assembly set. In the first cycle, select the SAG with the largest coverage in the current SAG set, add it to the SAG assembly set, merge the coverage areas, and obtain the first round of SAG assembly set; Or in the current non-first round, select the SAG with the largest effective coverage ratio p in the current SAG set and add it to the SAG assembly set obtained in S7-6-3 in the previous round, merge the coverage area, and obtain the SAG assembly set of the current round; S7-6-2. Calculate the first overlap length of the coverage area of ​​each SAG number q that is not currently included in the SAG assembly set and the coverage area of ​​SAG number p, and then calculate the relative coverage ratio of each SAG number q that is not currently included in the SAG assembly set relative to SAG number p based on the first overlap length, the set length of the coverage area of ​​SAG number p, and the read length of SAG number p that is not aligned to the species gene sequence; S7-6-3. Select the q with the largest relative coverage ratio max SAG No. is added to the SAG assembly set obtained in S7-6-1, the coverage area is merged, and the SAG assembly set of the current round is obtained iteratively; S7-6-4. Calculate the second overlap length of the coverage area of ​​each SAG q' that is not currently included in the SAG assembly set and the coverage area of ​​the SAG assembly set obtained in S7-6-3, and then calculate the effective coverage ratio of each SAG q' that is not currently included in the SAG assembly set relative to SAG p based on the second overlap length, the set length of the coverage area of ​​the SAG assembly set obtained in S7-6-3, and the read length of SAG p that is not aligned to the species gene sequence; S7-6-5. Loop through steps S7-6-1 to S7-6-4 until the preset assembly conditions are met to obtain a SAG assembly set, i.e., the corresponding strain bin; wherein the preset assembly conditions are set based on 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.

[0021] Beneficial effects of the present invention: Compared with the method of assembling species-level genome sequences provided in the prior art, the present invention 1) reduces human intervention by automatically calculating the number of assembly rounds and simplifying the process, thereby improving the convenience and efficiency of operation; 2) can effectively control the contamination of species bins to ensure the accuracy of assembly results; 3) by grouping and calculating the similarity matrix, the memory consumption in the case of high SAG numbers is significantly reduced, while the calculation speed is improved, solving the problems existing in large-scale data processing; 4) can dynamically select or adjust the similarity threshold to adapt it to different sample types or library qualities, thereby enhancing the flexibility and applicability of the method. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 A schematic diagram of a process for analyzing single-cell microbial genome sequencing data provided by an embodiment of the present invention.

[0023] Figure 2 The present invention is a flowchart of using steps S1 to S5 of the method of an embodiment of the present invention to co-assemble species genome sequences from a single-cell microbial genome sequencing library. DETAILED DESCRIPTION

[0024] 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.

[0025] The embodiment of the present invention provides a method for analyzing single-cell microbial genome sequencing data. The method initializes the bin list after obtaining the single-cell microbial genome sequencing library. Then, each bin is assembled according to the set k-mer and the integrity and contamination of the assembly result are evaluated. The assembly result is filtered to remove shorter contigs to control the contamination. Then, the bins are grouped and the similarity between the bins in the group is calculated. The similar bins are hierarchically clustered according to the species annotation information and the similarity, and the hierarchical clustering results are split into one or more clusters according to the similarity threshold. In this process, the similarity threshold is dynamically adjusted to adapt to different sample types or library quality until the number of clusters in the group or the similarity threshold meets the preset split termination condition. Then, similar bins are merged into new bins according to the clustering, and the k-mer value is updated for the next round of assembly. Repeat the above assembly process until the bin list is no longer updated, or the assembly meets the preset assembly termination condition, and the bins in the list will be used for further assembly of subsequent species genomes. The above steps reduce memory consumption and improve the efficiency and accuracy of the genome sequence assembly process of single-cell microbial species by automating the process, controlling contamination, and dynamically adjusting the similarity threshold. Subsequently, the obtained bins are successively removed from polycellular 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: S1. Obtain a single-cell microbial genome sequencing library and set the current assembly round counter m is 0, each SAG generates a bin, and initializes the bin list: ; in, SAG i For the i No. SAG; bin m,i For the m Assembly wheel i No.bin; binList m For the m The bin list of the assembled wheel; SAGCount is the total number of SAGs in the sequencing library.

[0026] S2. Acquisition bin m,i Corresponding fastq file fasta m,i , and then use SPAde, according to the current m-th assembly round set kmer m ,right bin m,i Assemble and evaluate using checkm bin m,i The completeness of the assembly result Completeness m,i 、 Pollution Contamination m,i , and obtain assembly information. The length of k-mer is K, the value range of K is 21-61 and K is an odd number.

[0027] S3. If the assembly information satisfies 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): (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 C; (2) The number of assembly rounds m is greater than the preset maximum number of assembly rounds MaxM; S4-1. Yes bin m,i Filter the assembly results and remove bin m,i The contigs whose length is less than one-Lth of the length of ContigN50 in the assembly result are obtained. bin m,i The filtered assembly results FilteredFasta m,i. Wherein, the value range of L is 5~20, and in a preferred embodiment, L is 10. Contig N50 is an important statistical indicator in the field of genome sequencing and assembly, which is used to evaluate the quality of genome assembly. Specifically, Contig N50 refers to arranging all contigs of the genome from large to small according to length, and then accumulating the length of these contigs until the total accumulated length reaches half of the total length of all contigs, and the length of the last accumulated contig is Contig N50. For example, if the size of a genome is 1M, after the reads obtained by sequencing are spliced ​​into contigs, these contigs are arranged from long to short, and then accumulated in sequence until the accumulated length reaches 500k (i.e. 50% of 1M), then the length of the last accumulated contig is ContigN50. The larger the Contig N50 value, the better the continuity of the assembly, the more long fragments assembled, and thus closer to the real genome structure, so Contig N50 is generally considered to be an important criterion for measuring the quality of genome splicing results. Therefore, the present invention filters the bin assembly results according to the length of the contig in the bin and the length of Contig N50, discards shorter contigs with poor assembly continuity or integrity, thereby controlling the contamination of the bin at a lower level for subsequent analysis.

[0028] S4-2. Divide all bins into M groups according to the number of N bins in each group; bin m,i The filtered assembly results FilteredFasta m,i And the current mth assembly round set k-mer m , calculated using sourmash's MinHash function bin m,i The minhash feature is then calculated using the Compare function of sourmash bin m,i and bin m,ii (i.e., except for the nth group in the mth assembly round bin m,i Similarity of any bin outside the group; ; in, G m,n The set of bins of group n in the mth assembly round; mhm,n,i The minhash feature of bin number i of group number n in the mth assembly round; S m,n,i,ii It is the similarity between bin No. i and bin No. ii of group No. n in the mth assembly round.

[0029] In this step, all bins are divided into M groups according to N bins in each group. For example, when N=5000 and n=2, all bins are divided into M groups according to 5000 bins in each group. Then group 2 includes SAGs 5000-9999. The similarities between the 5000 bins in the group are calculated, and then the subsequent merging steps, i.e. steps S4-3 to S4-6, are performed in groups. Compared with the method in the prior art of putting all bins together to calculate the similarity and then merging them, the present invention uses a grouping method to calculate the similarity moment, which can reduce memory consumption when the number of SAGs is high and improve the calculation speed. Moreover, after multiple rounds of co-assembly, it has no effect on the assembly results.

[0030] S4-3. If the integrity of a bin in a group is not less than 0.9 or the contamination is not less than 0.5, the similarity between the bin and any bin in the group is modified to 0; based on the modified similarity, the linkage function of scipy is used to perform hierarchical clustering on the bins in the group to obtain the hierarchical clustering results of the group; ; in, C m,n,i It is the assembly result of bin No. i in group No. n in the mth assembly round.

[0031] In this step, the purpose of modifying the similarity is to exclude bins with relatively high completeness (i.e., completeness greater than or equal to 0.9) or high contamination (i.e., contamination greater than or equal to 0.5) after assembly from subsequent merging to prevent excessive contamination from being included. For bins with high completeness, continued merging will greatly increase their contamination, increasing the difficulty of subsequent decontamination. Bins with high contamination may gather polysomal SAGs or SAGs with high backgrounds, and continued merging will also increase the difficulty of subsequent decontamination. Therefore, this step modifies the similarity between bins with relatively high completeness or high contamination after assembly and any bin in the group to 0, so that the bin will not be split into the same cluster with the remaining bins during the subsequent hierarchical clustering and splitting process, and will no longer be included in the subsequent merging process. In addition, hierarchical clustering of bins within a group in small groups can effectively control memory consumption and increase computing speed.

[0032] S4-4. Set the current split round counter mm to 0, and set the similarity threshold of the current mm split round to SH mm ; Based on the similarity threshold SH mm , use scipy's fcluster function to split the group hierarchical clustering results into one or more clusters (cluster), and calculate the ratio R of the number of all clusters in the group to the number of all bins. In this step, the similarity threshold can be dynamically adjusted according to the split results. Because different sequencing libraries have different coverages, the level of coverage directly affects the similarity, so 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 SAG numbers. This solves the problem in the prior art that the similarity threshold cannot be dynamically selected, and different sample types or library qualities require manual adjustment and modification of the similarity threshold.

[0033] S4-5. If R is less than the ratio threshold, execute step S4-6; otherwise, according to the preset similarity step length StepSH Lower the similarity threshold SH mm , , get the similarity threshold for the next splitting round, i.e. the mm+1 splitting round SH mm+1 , if the similarity threshold is lowered SH mm+1 Less than the minimum similarity SH , then execute step S4-6, otherwise repeat step S4-4.

[0034] S4-6. 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 (the m+1th assembly round) binList m+1 ; If the number of bins in a cluster within a group is greater than the merge threshold maxC, the cluster is split into multiple small clusters, and then the bins in a single small cluster are merged into a new bin, where the number of bins in the small cluster is not greater than the merge threshold maxC. For example, the first group in the first assembly round has 5000 bins. After hierarchical clustering, according to the similarity threshold set in the first splitting round SH 1The split results in 3000 clusters, and the ratio of 3000 clusters to 5000 bins (R=0.6) is less than the ratio threshold. Then multiple bins in a single cluster are merged into a new bin, and the 3000 clusters in group 1 will get a total of 3000 new bins. These 3000 new bins will enter the second assembly round together with the new bins obtained from other groups.

[0035] S5. If the bin list of the next assembly round (the m+1th assembly round) has the same number of bins as the bin list of the current assembly round (the mth assembly round), then binList m+1 = binList m , then execute step S6-1, otherwise according to the current assembly round number m and k-mer decay step length StepkMer renew kM m , get the value for the m+1th assembly round kM m+1 , and then repeat steps S2 to S4; ; in, StepkMerLen The number of rounds required for each decay; StepkMer For each decay step, the decay step must be an odd number.

[0036] In a specific embodiment, the above steps S1 to S5 are used to perform species genome co-assembly on a single-cell microbial genome sequencing library with 21914 SAGs, and some assembly results are shown in Table 1. The specific parameters are set as follows: The length of the initial k-mer was set to 51; The completeness threshold C (the ratio of the number of SAGs corresponding to bins with completeness greater than 0.5 to the total number of SAGs) is set to 0.75; The maximum number of assembly rounds MaxM is set to 20; The length threshold was set to one tenth of the length of Contig N50 in the assembly result; The bins for each round are divided into 8,000 bins per group; The similarity threshold SHmm is set to 0.25; The similarity step size StepSH is set to 0.01; The ratio threshold (the ratio of the number of all clusters to all bins within a group) was set to 0.75; The merge threshold maxC (the maximum number of bins in a cluster within a group) is set to 15; The k-mer decay step length StepkMer is set to 5; The number of rounds required for each decay, StepkMerLen, is set to 5.

[0037] Table 1. Some bins that meet high assembly quality standards in the genome co-assembly results of species from single-cell microbial genome sequencing libraries bin number Completeness Pollution 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 Among them, N50 is the Contig N50 length of the bin.

[0038] This time, 21,914 SAGs were assembled into 3,538 species-level bins, of which 288 species-level bins had a completeness greater than 50, with an average contamination of 9.33 and a maximum contamination of 72. There were 87 species-level bins that met the high assembly quality standard. Table 1 shows the assembly information of some of the species-level bins, such as bin 8_3271, 8_3145, 8_3148, and 8_3098, which have better species genome assembly effects.

[0039] S6. Remove the polysomal contamination and background contamination in the bin obtained in step S5 in turn to obtain the species bin, and obtain the corresponding species assembly result after assembly.

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

[0041] In a specific implementation, step S6 specifically includes the following steps: S6-1. For each bin in the bin list obtained in step S5, align all SAGs in the bin to the assembly result conitg of the bin in step S2 to obtain the SAG-conitg matrix M.

[0042] S6-2. 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 In some embodiments, 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.

[0043] 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.

[0044] S6-3. 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 bins without polytopic contamination and the assembly results without polytopic contamination.

[0045] S6-4. Evaluate the assembly result obtained in step S6-3 after removing polysome contamination, 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.

[0046] S6-5. Set the length search boundary, length search step, coverage search boundary, and coverage search step of the contig in the assembly result after removing polysome contamination, and then gradually search the assumed length threshold and assumed coverage threshold of the contig. In some embodiments, one L of the contig length is used as the length search boundary of the contig, and the value range of L is 5 to 20; one C of the contig coverage is used as the coverage search boundary of the contig, and the value range of C is 5 to 20.

[0047] 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 results after removing polycellular contamination according to different cleaning thresholds to obtain retained contigs and characteristic sequences distributed in the retained contigs; evaluate the completeness and contamination of the assembly results after removing polycellular contamination at each cleaning threshold based on the reference characteristic sequence of the bin and the characteristic sequence of the retained contig; perform quality rating of the assembly results after removing polycellular contamination at each cleaning threshold using completeness and contamination as indicators, and select the optimal cleaning threshold based on the quality rating.

[0048] In some embodiments, the method for evaluating the completeness and contamination of the assembly result of removing polycellular contamination under the cleaning threshold according to the reference feature sequence of Bin and the feature sequence of the retained contig includes: determining the number of species of feature sequences in the retained contig and the number of species of feature sequences that appear repeatedly in the retained contig; determining the completeness of the assembly result of removing polycellular contamination under the cleaning threshold according to the ratio of the number of species of feature sequences in the retained contig to the number of species of the reference feature sequence of the bin; determining the contamination of the assembly result of removing polycellular contamination under the cleaning threshold according to the ratio of the number of species of feature sequences that appear repeatedly in the retained contig to the number of species of the reference feature sequence of the bin. Since in the complete and unpolluted single microbial assembly result sequence, the feature 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 according to the ratio of the number of species of feature sequences in the retained contig to the number of species of the reference feature sequence of the bin, that is, the larger the ratio, the closer the number of species of feature sequences in the retained contig is to the number of species of the reference feature sequence of the bin, and the more complete the assembly result is. The present invention determines the contamination degree of the assembly result after removing polysome contamination under the cleaning threshold according to the ratio of the number of species of characteristic sequences repeatedly appearing 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 repeatedly appearing in the retained contig, and the repeatedly appearing characteristic sequences may come from the contaminated contig assembled by the background reads, which means that the assembly result is more seriously contaminated.

[0049] 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: 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.

[0050] 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: 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: .

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

[0052] In this embodiment, step S6 gradually searches for the best combination of length threshold and coverage threshold based on the occurrence pattern of characteristic gene sets in the species genome, which is used to filter the single-cell microbial genome species assembly results, thereby removing the contaminated contigs assembled by background reads. This method 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 pollution removal results as much as possible. For the specific content of step S6, refer to CN119068975 A.

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

[0054] S7-1. Align each SAG to the species genome and obtain the SNP of each SAG relative to the species genome.

[0055] S7-2. Merge the SNPs of all SAGs to obtain a preliminary SNP matrix, filter the SNPs in the preliminary SNP matrix according to preset filtering conditions, and 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 sites, and / or the number of mutations of SNPs at different sites, and / or the gene to which the SNPs belong. SNP sites with low detection depths in sequencing results are prone to inaccurate analysis results. Therefore, in one or more specific embodiments, the detection depths of SNPs at different sites are included in the preset filtering conditions, and sites with too low detection depths, such as SNP sites with a detection depth of less than 5, are removed to ensure the reliability of the selected characteristic SNPs. The difference in the number of mutations of SNPs at a certain site 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 sites is included in the preset filtering conditions, and SNP sites with too few mutations, such as SNP sites with less than 1% mutations, are removed, thereby effectively identifying mutations with biological significance and avoiding mistaking random mutations for characteristic SNPs.

[0056] 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 the evaluation index, and determine the optimal number of strains based on the evaluation results.

[0057] In some embodiments, step S7-3 includes: S7-3-1. Calculate the mean and standard deviation of each characteristic SNP in the characteristic matrix; S7-3-2. Sort all characteristic SNPs by their standard deviations, select the first x SNP characteristics with standard deviations decreasing from the largest to form an analysis matrix, and perform hierarchical clustering on SAG 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, all characteristic SNPs are included in the analysis matrix; S7-3-3. Determine the maximum number of strains Q traversed by the hierarchical clustering results; 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 results cannot be improved; S7-3-4. Divide the hierarchical clustering results by the number of n strains, where n is any integer between 1 and Q; S7-3-5. The WSS score is used as an evaluation index to evaluate each division result. The calculation formula of the WSS score is: , in, WSS n is the WSS score when divided into n strain numbers; FeatureCount is the number of SNP features in the analysis matrix; Feature i is the value set of the i-th SNP feature in the analysis matrix; ClusterCenter m,i To analyze the characteristic center of the ith SNP feature of the mth strain when the matrix is ​​divided into n strains, the average value of the ith SNP feature of the SAG set of the mth strain is used as ClusterCenter m,i ; P m is the penalty score of the mth strain when the analysis matrix is ​​divided into n strains. m It 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 mth strain in the analysis matrix.

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

[0059] 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: , in, R p is the coverage of SAG number p; MappedCov p is the coverage length of SAG number p on the gene sequence of the species; UnmappedCov p It is the length of reads that SAG number p has not been aligned to the gene sequence of the species.

[0060] S7-6-1. Set the initialization state of the SAG assembly set. In the first cycle, select the SAG with the largest coverage in the current SAG set, add it to the SAG assembly set, merge the coverage areas, and obtain the first round of SAG assembly set; Or in the current non-first round cycle, select SAG number p 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 cycle, merge the coverage area, and obtain the SAG assembly set of the current round.

[0061] S7-6-2. Calculate the first overlap length of the coverage area of ​​each SAG number q that is not currently included in the SAG assembly set and the coverage area of ​​SAG number p, and then calculate the relative coverage ratio of each SAG number q that is not currently included in the SAG assembly set relative to SAG number p based on the first overlap length, the set length of the coverage area of ​​SAG number p, and the read length of SAG number p that is not aligned to the species gene sequence. In some embodiments, the relative coverage ratio TR p, q The calculation formula is: , in, D p,q is the first superposition length of the set of coverage areas of each SAG number q that is not currently included in the SAG assembly set and the set of coverage areas of the SAG number p; C p is the set of coverage areas of SAG number p.

[0062] S7-6-3. Select the SAG with the largest relative coverage ratio qmax and add it to the SAG assembly set obtained in S7-6-1, merge the coverage area, and iterate to obtain the SAG assembly set of the current round.

[0063] S7-6-4. Calculate the second overlap length of the coverage area of ​​each SAG number q' that is not currently included in the SAG assembly set and the coverage area of ​​the SAG assembly set obtained in S7-6-3, and then calculate the effective coverage ratio of each SAG number q' that is not currently included in the SAG assembly set relative to the SAG number p based on the second overlap length, the set length of the coverage area of ​​the SAG assembly set obtained in S7-6-3, and the read length of the SAG number p that is not aligned to the species gene sequence. In some embodiments, the effective coverage ratio R q’, a The calculation formula is: , in, D q’, a is the second superposition length of the set of coverage areas of each SAG with number q' that is not included in the SAG assembly set in the current round a and the set of coverage areas of the SAG assembly set obtained in S6-3; C a It is the set of coverage areas of the SAG assembly set obtained by S7-6-3 in the current round a.

[0064] S7-6-5. Execute steps S7-6-1 to S7-6-4 repeatedly until the preset assembly conditions are met, and obtain the SAG assembly set, i.e., the corresponding strain bin; wherein 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 comparison rate of the SAGs included in the SAG assembly set. In some embodiments, the preset assembly condition is any one of the following (1) to (4): (1) The coverage area of ​​the SAG assembly set accounts for more than 95% of the gene sequence length of the species; (2) The number of SAGs in the SAG assembly set is greater than a set SAG number threshold, and the SAG number threshold is preferably 50; (3) The alignment rate of the SAGs included in the SAG assembly set in round a+1 to the species genome is less than the set alignment rate threshold, and the alignment rate threshold is preferably 50%; (4) All SAGs are included in the SAG assembly set.

[0065] S7-7. Remove background contamination within the strain bin and obtain the corresponding strain assembly results after assembly.

[0066] Step S7 first finely filters the SNP information of SAG based on the preset filtering conditions such as the detection depth of SNPs at different sites, the number of mutations, and the diversity of genotypes in step S7-2 to obtain a feature matrix containing characteristic SNPs. This preprocessing method can more effectively process SAGs with unclear SNP characteristics or low coverage, thereby improving the accuracy of subsequent hierarchical clustering and reducing the risk of erroneous clustering; on the other hand, it can efficiently focus on specific genes, thereby distinguishing bacterial species with specific genes and providing more accurate strain genotype information. Subsequently, when dividing the SAG hierarchical clustering results in step S7-3, an algorithm that can automatically determine the optimal number of strains is introduced. Based on the sequencing results of the single-cell microbial genome, the strains can be accurately distinguished with high throughput, 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 the process of screening SAGs for assembly in step S7-6, the present invention creatively uses the comparison rate, coverage, relative coverage ratio and effective coverage ratio of each SAG as the basis, and adopts a greedy algorithm to efficiently and accurately select the optimal SAG assembly set, which effectively improves the integrity and accuracy of the strain assembly result in step S7-7, reduces contamination, and reduces information loss in the strain assembly process, thereby obtaining a high-quality strain genome sequence. For the specific content of step S7, refer to CN 118737269 A.

[0067] 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 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; The background contamination in the strain bin was removed, and the corresponding strain assembly results were obtained after assembly.

2. The method according to claim 1, characterized in that 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.

3. The method according to claim 1, characterized in that 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.

4. The method according to claim 1, characterized in that 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.

5. The method according to claim 4, characterized in that 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.

6. The method according to claim 5, characterized in that The value range of L is 5~20.

7. The method according to claim 4, 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.

8. The method according to claim 4, 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.

9. 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.

10. 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