A method for analyzing soil biological complex pollution using metagenomics

Abstract

The present invention discloses a method for analyzing soil biological complex pollution using metagenomics, which mainly includes: (1) constructing a gene set and a protein set; (2) extracting the sequence on the protein set after comparing it to the SARG database to construct an ARGs species traceability library for drug resistance gene traceability analysis; (3) constructing a personalized species database for the gene set for metagenomic microbial abundance calculation; (4) collecting pathogenic bacteria species information to construct a pathogenic bacteria species database, and obtaining the abundance of pathogens in the sample (5) calculating the biological complex pollution index. The present invention realizes the traceability of drug resistance genes, can easily track the situation of drug resistance genes carried by pathogens, and also reintegrates the pathogenic bacteria species in the public database and constructs a pathogenic bacteria database for extracting pathogenic bacteria information in the analysis results, proposes a new soil biological complex pollution calculation method, and provides a reference basis for the quantitative calculation of complex pollution.

Description

Technical Field

[0001] The present invention relates to the fields of soil biology and bioinformatics, and in particular to a method for analyzing soil biological complex pollution using metagenomes. Background Art

[0002] Soil biocontamination, consisting of antimicrobial resistance genes (ARGs) and plant and animal pathogens, has garnered widespread attention worldwide. With the advancement of intensive agricultural production and the increasing overuse of pesticides and antibiotics, soil ARGs have accumulated significantly, severely damaging soil ecosystems and leading to the frequent occurrence of soil-borne diseases. The interaction between soil pathogens and ARGs further exacerbates the risk of biocontamination in soil-plant systems, posing a serious threat to food safety and human health. Currently, the analysis of ARGs and pathogens in environmental samples primarily relies on next-generation sequencing (NGS) technology. In recent years, metagenomics (mNGS) technology has been increasingly applied to the detection of ARGs and pathogens in the environment, and mNGS-based databases and analysis tools have gradually matured and diversified. For example, analytical pipelines such as SARG-OAP and PathoLive enable the detection and analysis of antimicrobial resistance genes and pathogens at the mNGS read level.

[0003] At present, research on pathogens and ARGs is often in an independent and parallel state. Researchers often ignore the interactive effects of the two, and there is a lack of exploration of the interaction mechanism, transmission rules and influencing factors of the two pollution factors in the soil-plant system. At present, the existing analysis method is to perform high-throughput sequencing of samples (such as mNGS), obtain the species abundance of the samples and the abundance of ARGs, and then obtain potential high-risk carriers of ARGs through correlation networks. However, the analysis method is complicated, and the correlation network-based association method is far less reliable than reads comparison, and it is impossible to obtain the true abundance of ARGs carried by pathogens. At present, there is no way to evaluate the biological complex pollution caused by the interactive effect of the two. Therefore, it is urgent to establish a method based on pathogen analysis and drug-resistant gene tracing to provide theoretical and technical support for the research and regulation of soil biological complex pollution. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for analyzing soil biological complex pollution using metagenomics. By constructing a soil (pathogen) microbial species database and a non-redundant drug-resistant gene set database, pathogen species analysis and drug-resistant gene traceability annotation are performed based on metagenomic quality control data, thereby realizing rapid analysis of soil biological complex pollution.

[0005] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0006] A method for analyzing soil biological complex pollution using metagenomes, comprising the following steps:

[0007] Step 1: Construct gene sets and protein sets;

[0008] Step 2: After comparing the protein set to the SARG database, the aligned sequences were extracted to construct an ARGs species traceability library for drug resistance gene traceability analysis;

[0009] Step 3: Gene set construction of personalized species database for metagenomic microbial abundance calculation;

[0010] Step 4: Collect pathogenic bacteria species information to build a pathogenic bacteria species database and obtain the abundance of pathogens in the sample;

[0011] Step 5: Calculate the biological composite pollution index (SBI) based on the abundance of pathogenic bacteria and ARGs;

[0012] Step 6: Complete the packaging and construction of the above-mentioned ARGs traceability database, personalized species database, and biological composite pollution index calculation and analysis process;

[0013] Step 7: After the environmental samples are collected, DNA is extracted and a library is constructed. Illumina sequencing technology is used to complete metagenomic sequencing of the environmental samples. The sequencing data is quality controlled and human genome contamination is removed to obtain quality control data as the input source of the process;

[0014] Step 8. Input the quality control data into the biological complex pollution analysis process encapsulated and constructed in step 6 to obtain the results of pathogen abundance and ARGs species traceability abundance, and calculate the statistical table of biological complex pollution index; the results can be visualized.

[0015] Furthermore, in step 1, we collected the coding proteins of more than 100,000 reference genomes in the database, used CD-HIT software to remove redundancy and obtained 13,051 species gene sets; used Prodigal to predict the proteins in the gene sets and construct the corresponding species protein sets. The specific steps are as follows:

[0016] Step 1: Collect bacterial reference genome coding protein sequences from public databases;

[0017] Step 2: Merge all protein-coding sequences to species level to eliminate redundancy and construct a gene set;

[0018] Step 3: Prodigal predicts proteins in the gene set and constructs a protein set.

[0019] Furthermore, in step 1, the coding proteins of more than 100,000 reference genomes in the database were collected, and 13,051 species gene sets were obtained by removing redundancy using the cd-hit software; Prodigal was used to predict the proteins in the gene sets and construct the corresponding species protein sets.

[0020] Furthermore, in the step 2, the protein set is blastp-aligned with the SARG database, and the aligned sequences are extracted to construct an ARGs species traceability database, and finally an ARGs species traceability library containing 20,811 drug-resistant genes is established; then the ARGs-OAP software based on Perl language is used to encapsulate the ARGs species traceability database, and the metagenomic quality control reads are input into the encapsulation software to realize the drug-resistant gene species traceability annotation and quantitative analysis.

[0021] Among them, the quantitative analysis method is: the ARGs-OAP software has used USEARCH+BLASTX to compare the ARGs species traceability database to align and quantify the metagenomic sequencing reads. Specifically: the quality control data of all samples are compared with the ARGs species traceability database using BLASTX to obtain the resistance base and calculate the read abundance of each gene in each sample; the final result's abundance quantitative unit is the number of reads matching the gene per million reads.

[0022] Furthermore, in step three, a personalized species database is established based on the previously constructed gene set. First, the database download function of Kraken2 is used to download the taxonomy classification library, and then the gene set sequence is built using kraken-build, and then repackaged using Kraken2; the metagenomic sample quality control data is used as input, and according to the annotation results of Karken2, the number of read sequences annotated to different microbial species in each sample is counted at the kingdom, phylum, class, order, family, genus and species levels to obtain a species annotation abundance table.

[0023] Furthermore, in the step 4, a pathogenic bacteria species database is constructed based on the species level, and based on the pathogenic bacteria species information, the corresponding pathogenic bacteria species can be quickly extracted from the species abundance table and the drug-resistant gene species traceability abundance table.

[0024] Furthermore, in step 5, the sample's biological composite contamination index (SBI) is calculated based on the following formula:

[0025]

[0026]

[0027] SBI j =SBI patho,j+SBI ARG,j (Formula 3)

[0028] Among them, j represents the j-th sample, k represents the total number of samples, 0 < j < k; i represents the i-th pathogenic bacterium or drug resistance gene identified; P i,j and ARG i,j respectively represent the abundance of the pathogenic bacterium i in sample j and the abundance of the drug resistance gene. SBI patho,j represents the index of pathogenic bacteria in the biological complex pollution of sample j. SBI ARG,j represents the index of drug resistance genes in the biological complex pollution. SBI j represents the biological complex pollution index.

[0029] The beneficial effects of the present invention are as follows:

[0030] (1) The technology of annotating ARGs based on metagenomic sequencing data is mature and there are many methods. However, at present, the species tracing of ARGs only stays at judging potential drug resistance gene-carrying species based on the correlation network, and it is impossible to focus on pathogenic bacteria species in the environment. Since the present invention realizes the tracing of drug resistance genes, it can easily track the situation of drug resistance genes carried by pathogenic bacteria, and will provide a new analysis method for the study of biological complex pollution.

[0031] (2) Most of the database species of metagenomic species taxonomic annotation software focus on the strain level. Due to the relatively redundant taxonomic information in the database, false positives may occur in sequence comparison, resulting in a decrease in accuracy and affecting the abundance results. The present invention constructs a personalized database for non-redundant species gene set sequences and encapsulates it based on the Kraken2 software, successfully realizing the abundance annotation at the non-redundant species level. In addition, the current information of pathogenic bacteria databases is scattered and the species are incomplete. The present invention re-integrates the pathogenic bacteria species in the public database and constructs a pathogenic bacteria database for extracting pathogenic bacteria information in the analysis results.

[0032] (3) At the present stage, the research on pathogenic bacteria and drug resistance genes is often in an independent parallel state, ignoring the interaction effects generated by both. The emergence of the present invention will further promote the research on the interaction relationship between pathogenic bacteria and drug resistance genes. In addition, since the soil biological complex pollution composed of drug resistance genes and pathogenic bacteria is a brand-new pollution index and there is no evaluation standard at present, the novel calculation method of soil biological complex pollution proposed by the present invention will also provide a reference basis for the quantitative calculation of complex pollution. Description of the Drawings

[0033] Figure 1 is the flowchart of biological complex pollution analysis;

[0034] Figure 2This is a statistical table of the biological composite pollution index calculated in a specific example; where: SampleID: sample name, pathogene_total_reads: total number of pathogen reads in the sample, ARG_total_ppm: total number of drug-resistant genes in the sample in ppm, SBI: biological composite pollution index;

[0035] Figure 3 is the number of pathogenic bacteria reads in a specific instance (Top 10 species);

[0036] Figure 4 It is the relative abundance of ARGs types in a specific instance (ARGs Type level). DETAILED DESCRIPTION

[0037] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0038] like Figure 1 The present invention provides a method for analyzing soil biological complex pollution using metagenomics, which mainly includes: (1) collecting bacterial coding gene reference genome sequences in the database to construct a gene set, and performing protein prediction to obtain a protein set; (2) after comparing the protein set with the SARG database, extracting the aligned sequences to construct an ARGs species traceability library for drug resistance gene traceability analysis; (3) constructing a personalized species database from the gene set for metagenomic microbial abundance calculation; (4) collecting pathogenic bacteria species information to construct a pathogenic bacteria species database and obtain the pathogenic bacteria abundance in the sample; (5) calculating the biological complex pollution index. Finally, a rapid quantitative calculation of soil biological complex pollution based on metagenomics is achieved, thereby providing technical support for reducing soil biological complex pollution.

[0039] The detailed description of each step is as follows:

[0040] (1) Construction of gene sets and protein sets

[0041] We collected the protein encoding data from over 100,000 reference genomes in the database and used the CD-HIT software to remove redundancy and obtain 13,051 gene sets. We then used Prodigal to predict the protein encoding data from these gene sets and the corresponding protein sets. The specific steps are as follows:

[0042] Step 1: Collect bacterial reference genome coding protein sequences from public databases;

[0043] Step 2: Merge all protein-coding sequences to species level to eliminate redundancy and construct a gene set;

[0044] Step 3: Prodigal predicts proteins in the gene set and constructs a protein set.

[0045] (2) ARGs species traceability library packaging

[0046] The protein set was blastp-matched with the SARG database, and the sequences on the alignment were extracted to construct an ARGs species traceability database, ultimately establishing an ARGs species traceability library containing 20,811 drug-resistant genes. The ARGs-OAP software based on the Perl language was used to encapsulate the ARGs species traceability database, and the metagenomic quality control reads were input into the encapsulated software to realize the annotation and quantitative analysis of drug-resistant gene species traceability. This method not only provides more comprehensive annotation and more accurate quantification of drug-resistant genes in the metagenomic sequencing information of the environment, but also realizes the species traceability of soil drug-resistant gene contamination. The internal quantitative analysis method of the software is described as follows:

[0047] The ARGs-OAP software internally uses USEARCH+BLASTX to align and quantify metagenomic sequencing reads against the ARGs species traceability database. Specifically, quality control data from all samples are aligned against the ARGs species traceability database using BLASTX, obtaining resistance bases and calculating read abundance for each gene in each sample. The final abundance quantification unit is the number of reads per million reads (ppm) matching the gene.

[0048] (3) Species database construction and abundance calculation

[0049] The personalized species database is constructed based on the gene set sequence. Kraken2 is a software that uses accurate k-mer alignment and the most LCA voting results to quickly annotate metagenomic DNA sequences for species. Although the software itself can use multiple databases to identify species abundance of metagenomic quality control data, it is unable to perform taxonomic abundance annotation at the non-redundant species (Species) level. Therefore, it is necessary to build a personalized species database based on the gene set constructed in the early stage. First, use Kraken2's own database download function to download the taxonomy classification library, then use kraken-build to build a library for the gene set sequence and repackage it using Kraken2. The metagenomic sample quality control data is used as input. According to the annotation results of Karken2, the number of read sequences annotated to different microbial species in each sample is counted at the kingdom, phylum, class, order, family, genus and species levels to obtain a species annotation abundance table.

[0050] (4) Construction of pathogenic bacteria database

[0051] Based on the species level, we constructed a pathogenic bacteria species database. Considering that there may be redundancy in pathogenic bacteria strains, we removed the redundancy and divided it into 1,640 species at the species level. The specific composition of pathogenic bacteria at the Phylum level is as follows: Figure 3As shown in the figure, a Perl language script was written independently to quickly extract the corresponding pathogenic bacteria species from the species abundance table and the drug-resistant gene species traceability abundance table based on the information of pathogenic bacteria species.

[0052] (5) Calculation of biological composite pollution factors

[0053] Currently, there is no algorithm and calculation standard for evaluating biological composite pollution. To evaluate the soil biological composite pollution caused by the interaction between pathogenic bacteria and ARGs, an algorithm based on the following formula was used to calculate the biological composite pollution index (SBI, Soil biopullutome index) of the samples:

[0054]

[0055]

[0056] SBI j = SBI patho,j + SBI ARG,j (Formula 3)

[0057] Among them, j represents the jth sample, k represents the total number of samples, and 0 < j < k. And i represents the ith pathogenic bacterium (Formula 1) or drug-resistant gene (Formula 2) noted. P i,j and ARG i,j respectively represent the abundance of pathogenic bacterium i in sample j and the abundance of drug-resistant gene, SBI patho,j represents the biological composite pollution pathogenic bacterium index of sample j, SBI ARG,j represents the biological composite pollution drug-resistant gene index, SBI j represents the biological composite pollution index.

[0058] The following is further illustrated with specific examples:

[0059] This example uses a total of 40 rhizosphere soil samples collected in situ from tomato fields in Qilin Town, Nanjing. The 40 samples are divided into healthy and diseased treatments, and a new method for detecting biological composite pollution is implemented, including the following specific steps:

[0060] 1. Use Perl language to complete the encapsulation and construction of the above ARGs traceability database, personalized species-level database, and biological composite pollution index (SBI) calculation and analysis process.

[0061] 2. After sample collection, extract DNA from it to build a library, use Illumina sequencing technology to complete the metagenomic sequencing of environmental samples, and perform quality control on the sequencing data and remove human genome contamination to obtain quality control data as the input source of the process.

[0062] 3. Input the quality control data into the biological complex pollution analysis process to obtain the results of pathogen abundance and ARGs species traceability abundance, and calculate the statistical table of biological complex pollution index, such as Figure 2 The specific pathogenic bacteria reads absolute abundance and ARGs relative abundance (ppm) are as follows: Figure 3 and Figure 4 .

[0063] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the above embodiments do not limit the scope of protection of the present invention in any form, and all technical solutions obtained by equivalent substitution and the like fall within the scope of protection of the present invention.

[0064] The parts not involved in the present invention are the same as the existing technology or can be implemented by using the existing technology.

Claims

1. A method for analyzing soil biological complex pollution using metagenomics, characterized in that The following steps are involved: Step 1: Construct gene sets and protein sets; details are as follows: Step 1: Collect 10w+ bacterial reference genome coding protein sequences from public databases; Step 2: Use cd-hit software to remove redundancy at the species level of the protein coding sequences and construct a gene set; Step 3: Use Prodigal software to predict proteins on the gene set and construct the corresponding species protein set; Step 2: After comparing the protein set to the SARG database, the sequences on the comparison were extracted to construct the ARGs species traceability library for drug resistance gene traceability analysis; the details are as follows: The protein set was blasted against the SARG database, and the sequences on the alignment were extracted to construct the ARGs species traceability database, and finally the ARGs species traceability library of drug-resistant genes was established; the ARGs-OAP software based on Perl language was used to encapsulate the ARGs species traceability database, and the metagenomic quality control reads were input into the encapsulation software to realize the annotation and quantitative analysis of drug-resistant gene species traceability; The quantitative analysis method is: the ARGs-OAP software uses USEARCH+BLASTX to compare the ARGs species traceability database to compare and quantify the metagenomic sequencing reads. Specifically: the quality control data of all samples are compared with the ARGs species traceability database using BLASTX to obtain the resistance base and calculate the reads abundance of each gene in each sample; the final result's abundance quantitative unit is the number of reads matching the gene per million reads; Step 3: Gene set construction of personalized species database for metagenomic microbial abundance calculation; details are as follows: Based on the gene set constructed in the early stage, a personalized species database was established. First, the taxonomy classification library was downloaded using the database download function of Kraken2, and then the gene set sequence was built using kraken-build, and then repackaged using Kraken2. The metagenomic sample quality control data was used as input, and according to the annotation results of Karken2, the number of read sequences annotated to different microbial species in each sample was counted at the kingdom, phylum, class, order, family, genus and species levels to obtain the species annotation abundance table. Step 4: Collect information on pathogenic bacteria species to build a pathogenic bacteria species database and obtain the abundance of pathogenic bacteria in the sample; Construct a pathogenic bacteria species database based on the species level, and based on the pathogenic bacteria species information, quickly extract the corresponding pathogenic bacteria species from the species abundance table and the drug-resistant gene species traceability abundance table; Step 5: Calculate the biological composite pollution index (SBI) based on the abundance of pathogenic bacteria and ARGs; the formula is as follows: SBI j = SBI patho,j + SBI ARG,j (Formula 3) Among them, j represents the j-th sample, k represents the total number of samples, 0 < j < k; i represents the i-th pathogenic bacterium or drug resistance gene identified; P i,j and ARG i,j respectively represent the abundance of pathogenic bacterium i in sample j and the abundance of drug resistance gene, SBI patho,j represents the biological composite pollution pathogenic bacterium index of sample j, SBI ARG,j represents the biological composite pollution drug resistance gene index, SBI j represents the biological composite pollution index; Step 6: Complete the packaging and construction of the above-mentioned ARGs traceability database, personalized species database, and biological complex pollution index calculation and analysis process; Step 7: After the environmental samples are collected, DNA is extracted and a library is constructed. Illumina sequencing technology is used to complete metagenomic sequencing of the environmental samples. The sequencing data is quality controlled and human genome contamination is removed to obtain quality control data as the process input source; Step 8. Input the quality control data into the biological complex pollution analysis process encapsulated and constructed in step 6, obtain the results of pathogen abundance and ARGs species traceability abundance, and calculate the statistical table of biological complex pollution index.

Images