A drug-resistant group analysis method based on 16S rRNA amplicon sequencing

Abstract

The present application belongs to the field of bioinformatics analysis, and relates to an analysis method for predicting the drug-resistant gene composition and abundance and drug resistance rate of a bacterial community in a sample based on 16S rRNA amplicon sequencing results. The method of the present application greatly expands the application range of 16S rRNA amplicon sequencing, and provides a simple and easy-to-use and low-cost data re-mining method for the study of the spread of drug-resistant genes in the environment.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of bioinformatics analysis, and particularly relates to an analysis method for predicting functional genes based on 16S rRNA amplicon, and inferring the composition and abundance of drug-resistant genes in a sample. BACKGROUND

[0002] Protein function prediction methods rely on databases annotated by experts; one widely used database is the KEGG (Kyoto Encyclopedia of Genes and Genomes) database; this database consists of a comprehensive database of molecular pathways, networks and genes involved in various cellular processes such as metabolism, signal transduction and disease. The KO (KEGG Orthology) database is a database that represents molecular functions in terms of functional orthologs; functional orthologs are defined manually in KEGG molecular networks; KO numbers are obtained from experimental results of genes and proteins in specific organisms, and these numbers are then used to assign homologous genes in other organisms.

[0003] Antibiotic resistance genes or drug-resistant genes (ARGs) are widely concerned as new pollutants. Plasmids, integrins and transposons carrying ARGs in bacteria can be horizontally transferred between strains of the same species and between strains of different species, even if the resistant strains die, the naked DNA carrying ARGs will exist for a long time under the protection of deoxynucleotidase. Currently, the analysis of ARGs mainly relies on high-throughput quantitative polymerase chain reaction (qPCR) and metagenomics methods. Although these methods can produce accurate information about the abundance and distribution of ARGs, they are complex and time-consuming to operate. 16s rRNA-based amplicon sequencing is an important method for studying microbial community structure, and is simple to operate and low in cost, and is widely used in the field of microbial ecology. 16s rRNA-based amplicon sequencing is mainly used to study the structure and composition of microbial communities, while PICRUSt2-based functional prediction realizes the analysis of microbial functions on the basis of amplicon sequencing. If the analysis of drug-resistant genes can be realized on the basis of functional prediction, on the one hand, it can expand the application range of 16s rRNA amplicon sequencing, and on the other hand, it can greatly reduce the analysis cost of ARGs in the environment. SUMMARY

[0004] The technical problem to be solved by the present application is: how to realize the analysis of the composition and abundance of ARGs in the sample based on 16s rRNA amplicon sequencing, thereby expanding the application range of the method and reducing the cost of ARGs analysis in the sample. The analysis method for predicting the composition and abundance of drug-resistant genes based on 16s rRNA established by the present application can realize the abundance and composition of ARGs in the sample based on function prediction, and calculate the drug resistance rate of the sample.

[0005] The purpose of the present application is to provide an analysis method for predicting the composition and abundance of drug-resistant genes in a sample based on 16S rRNA. Specifically, first, obtain the sequence through the amplicon determination of 16S rRNA specific primers, and obtain the ASV sequence information through QIIME2 analysis; then, predict the KEGG functional genes present in the sample through PICRUSt2, establish the number relationship of KEGG in ASV, screen the KO belonging to ARG through the relationship of KO corresponding to ARG in the KEGG database, establish the quantity relationship list of ARG contained in ASV, obtain the abundance of ARG in each sample, and further infer the abundance and composition information of drug-resistant bacteria and drug-resistant genes in the sample.

[0006] The present application adopts the following technical solutions:

[0007] 1. An analysis method for predicting the composition and abundance of drug-resistant genes of bacterial community in a sample based on 16S rRNA amplicon sequencing sequence, characterized in that it comprises the following steps:

[0008] 1) Perform 16S rRNA amplicon sequencing on all bacteria in the sample to obtain 16S rRNA sequences of all bacteria in the sample;

[0009] 2) Obtain the species and abundance of bacteria contained in the sample, which includes analyzing the 16S rRNA sequence using PICRUSt2 software to output the species and abundance of ASV numbered bacteria contained in the sample;

[0010] 3) Obtain the drug-resistant genes and their abundance contained in various bacteria, which includes analyzing the 16S rRNA sequence using PICRUSt2 software to output the functional genes and their number contained in various ASV numbered bacteria, and screening and retaining the drug-resistant genes and their number contained in various ASV numbered bacteria through a drug-resistant gene list;

[0011] 4) Based on the species and abundance of ASV numbered bacteria and the drug-resistant genes and their number contained in various ASV numbered bacteria obtained in steps 2) and 3), obtain the species and abundance of drug-resistant genes in the sample by calculation;

[0012] The calculation in step 4) is as follows: the abundance of the bacteria of the ASV number obtained in step 2) is multiplied by the number of drug-resistant genes contained in the bacteria of the ASV number obtained in step 3), to obtain the abundance of the drug-resistant genes contained in the bacteria of the ASV number in the sample; and the abundances of the drug-resistant genes contained in all bacteria of ASV numbers in the sample are added up to obtain the abundance of the drug-resistant genes in the sample, and the abundances of all drug-resistant genes in the sample are further obtained in the same way.

[0013] 2. The method of item 1, further comprising the step of obtaining ASVs of bacteria using QIIME2 software analysis, and obtaining annotations of bacteria by Greengene2 database alignment.

[0014] 3. The method of item 1 or 2, wherein the 16S rRNA sequence is a full-length or V3-V4 region or V4 region of 16S rRNA sequence.

[0015] 4. The method of item 1, wherein steps 2) and 3) further comprise assigning ASV numbers of bacterial species to ASV order.

[0016] 5. The method of item 1, wherein the list of drug-resistant genes in step 3) is obtained from the KEGG database.

[0017] 6. A method for predicting the drug resistance rate of a bacterial community in a sample, comprising the steps of any one of items 1-5, obtaining the species and abundance of bacteria of ASV numbers containing drug-resistant genes in the sample, and adding up the abundances of all bacteria of ASV numbers containing drug-resistant genes in the sample and multiplying by 100% to obtain the drug resistance rate of the sample.

[0018] Definitions

[0019] KO (KEGG Orthology): KO is part of KEGG, which provides a functional classification of genes and proteins; KO helps understand the function and regulation of biological systems by linking genes with biological processes, metabolic pathways, etc. In the classification system of KEGG proteins or enzymes, proteins with similar sequences and similar functions in the same pathway are grouped into a group and defined as a KO, written as "K + number", such as K03585, K18138, K12340, etc.

[0020] OTU: operational taxonomic unit.

[0021] ASV: Amplicon sequence variant, refers to the ASV number of each microorganism obtained by determining and analyzing the DNA sequence in the microbial community through high-throughput sequencing technology; compared with the traditional OTU, ASV can more accurately identify the microbial species, distinguish the microbial species at the subspecies level, and increase the accuracy and sensitivity of the microbial diversity analysis.

[0022] ARB: Antibiotic-resistant bacteria.

[0023] VGT: Vetical gene transfer.

[0024] HGT: Horizontal gene transfer.

[0025] The present application has the following advantages due to the above technical scheme:

[0026] 1. Compared with qPCR and metagenomic methods, the present application provides prediction for the composition of ARG in the sample based on 16S rRNA, which can greatly simplify and reduce cost.

[0027] 2. Based on the 16S rRNA amplicon sequencing sequence, the drug-resistant bacteria and drug-resistant groups can be output, and the resistance genes are calculated and predicted by bacterial community reverse calculation, which provides a new analysis method for the occurrence of ARG in the environment.

[0028] 3. The present application can accurately identify the relationship between ARG and bacterial community, and clearly calculate the drug resistance rate of bacterial community in the sample, realize the comprehensive analysis of drug resistance, and provide rich data for the research of drug resistance. BRIEF DESCRIPTION OF DRAWINGS

[0029] Figure 1 is the VGT and HGT network diagram of each sample, the upper row of diagrams is the VGT network diagram, wherein the small dots are ARG, the large dots are ARB, and the lines represent the relationship; the lower row of diagrams is the HGT network diagram, wherein the dots are ARB, and the lines represent the relationship. DETAILED DESCRIPTION

[0030] The present application will be further described below in combination with the drawings and specific embodiments, so that those skilled in the art can

[0031] better understand the present application and can implement it, but the embodiments are not limiting to the present application.

[0032] Example 1: Method for obtaining ARG abundance and drug resistance rate of bacterial community in sample

[0033] The method for obtaining ARG composition, abundance and resistance rate based on 16S rRNA amplicon sequencing sequence information comprises the following steps:

[0034] 1. Sequencing step: based on amplicon sequencing, obtain the 16S rRNA hypervariable region sequence of the bacterial genome in the sample. In this embodiment, sludge and water samples collected from a sewage treatment plant in Beijing are taken as examples for illustration.

[0035] The V3-V4 region of the 16S gene sequence is a hypervariable region, which is usually used for bacterial diversity analysis because the sequence difference in this region can well reflect the classification information of microorganisms; this part of the region can provide sufficient variability to distinguish different microbial communities, facilitating comparison and annotation; in this embodiment, the V4 region can be used to distinguish different types of bacteria.

[0036] 2. Analyze the 16S rRNA sequence by PICRUSt2 (http: / / github.com / picrust / picrust2) to output the types of bacteria with ASV numbers contained in the sample, which is used to predict the functional characteristics of the bacterial community contained in the sample, wherein A1, A2, A3, B1 and B2 are different samples for illustration, and the specific steps are as follows (or refer to the instructions for using the PICRUSt2 software):

[0037] According to the data name, make table A_Bac_manifest.tsv, as shown in Table 1.

[0038] Table 1

[0039]

[0040] cd / 01X_amplicon

[0041] conda activate base

[0042] conda activate qiime2_amplicon_2023.9

[0043] Use QIIME2 tools (www.qiime2.org) to transform and import sequence data

[0044] time qiime tools import \

[0045] --type 'SampleData[PairedEndSequencesWithQuality]' \

[0046] --input-path A_Bac_manifest.tsv \

[0047] --output-path A_Bac_manifest.qza \

[0048] --input-format PairedEndFastqManifestPhred33V2

[0049] Show basic information of imported sequences, determine parameters for DADA2 according to Q value of sequences

[0050] time qiime demux summarize \

[0051] --i-data A_Bac_manifest.qza \

[0052] --o-visualization 01A_Bacsample_demux_PEinfor.qzv

[0053] Sequence quality control and denoising: Use DADA2 to perform sequence quality control and denoising

[0054] time qiime dada2 denoise-paired \

[0055] --i-demultiplexed-seqs A_Bac_manifest.qza \

[0056] --p-trim-left-f 19 \

[0057] --p-trim-left-r 20 \

[0058] --p-trunc-len-f 220 \

[0059] --p-trunc-len-r 220 \

[0060] --p-n-threads 50 \

[0061] --o-table A_Bac_table.qza \

[0062] --o-representative-sequences A_Bac_rep-seqs.qza \

[0063] Show the obtained feature table information

[0064] time qiime feature-table summarize \

[0065] --i-table A_Bac_table.qza \

[0066] --o-visualization 03A_Bac_table.qzv \

[0067] Show basic information about the rep sequences obtained

[0068] time qiime feature-table tabulate-seqs \

[0069] --i-data A_Bac_rep-seqs.qza \

[0070] --o-visualization 04A_Bac_rep-seqs.qzv

[0071] Find the number of minimum sample valid reads from 03A_Bac_table.qzv as `MIN_READS` and rarefy:

[0072] time qiime feature-table rarefy \

[0073] --i-table A_Bac_table.qza \

[0074] --p-sampling-depth MIN_READS \

[0075] --o-rarefied-table A_Bac_rarefied-table.qza

[0076] Species annotation: Training of the classifier - Greengene2 database

[0077] Change min-length and max-length from 04A_Bac_rep-seqs.qzv

[0078] time qiime feature-classifier extract-reads \

[0079] --i-sequences 03Greengene2 / Greengene2.qza \

[0080] --p-f-primer GTGYCAGCMGCCGCGGTAA

[0081] --p-r-primer GGACTACNVGGGTWTCTAAT

[0082] --p-min-length 215

[0083] --p-max-length 290

[0084] --p-n-jobs 50

[0085] --o-reads 03Greengene2 / Greengene2-ref-seqs.qza

[0086] time qiime feature-classifier fit-classifier-naive-bayes

[0087] --i-reference-reads 03Greengene2 / Greengene2-ref-seqs.qza

[0088] --i-reference-taxonomy 03Greengene2 / Greengene2_tax.qza

[0089] --o-classifier 03Greengene2 / Greengene2-classifier.qza

[0090] Species annotation: Species annotation according to Greengene2 database

[0091] time qiime feature-classifier classify-sklearn

[0092] --i-classifier 03Greengene2 / Greengene2-classifier.qza

[0093] --i-reads A_Bac_rep-seqs.qza

[0094] --p-n-jobs 50

[0095] --o-classification 03Greengene2 / A_Bac_rep-seqs.qza

[0096] Species annotation: Greengene2 output

[0097] time qiime metadata tabulate \

[0098] --m-input-file 03Greengene2 / A_Bac_rep-seqs.qza \

[0099] --o-visualization 03Greengene2 / 01A_Bac_rep-seqs.qzv

[0100] time qiime taxa barplot \

[0101] --i-table A_Bac_rarefied-table.qza \

[0102] --i-taxonomy 03Greengene2 / A_Bac_rep-seqs.qza \

[0103] --m-metadata-file A_Bac_manifest.tsv \

[0104] --o-visualization 03Greengene2 / 02A_Bac_table.qzv

[0105] Species annotation: Output biom-formatted file and table of ASVs

[0106] time qiime tools export \

[0107] --input-path A_Bac_rarefied-table.qza \

[0108] --output-path 03A_Bac_table

[0109] biom convert -i 03A_Bac_table / feature-table.biom \

[0110] -o 03A_Bac_table / feature-table.txt \

[0111] --to-tsv

[0112] qiime tools export \

[0113] --input-path A_Bac_rep-seqs.qza \

[0114] --output-path 06A_Bac_rep-seqs

[0115] conda deactivate

[0116] conda activate picrust2

[0117] picrust2_pipeline.py -s 06A_Bac_rep-seqs / dna-sequences.fasta -i 03A_Bac_table / feature-table.biom -o 06Picrust2_stra -p 50 --stratified

[0118] 3. Analyze 16S rRNA sequences by PICRUSt2 (http: / / github.com / picrust / picrust2) to output the abundance of ASV-numbered bacteria contained in the samples

[0119] The "feature-table.txt" obtained in Step 2 is the read count of OTUs in all collected samples; add an ASV column, and sort the OTU IDs by ASV number; list the abundance of ASVs contained in each sample, as shown in Table 2.

[0120] Table 2

[0121]

[0122] 4. Obtain the correspondence of KOs and ARGs, for example, refer to KEGG Antimicrobial Resistance Genes (https: / / www.kegg.jp / kegg-bin / show_brite?ko01504.keg).

[0123] Obtain the drug-resistant genes contained in various bacteria and their abundance, including the following steps:

[0124] 5. Obtain the relationship between ASV-numbered bacteria and the drug-resistant genes ARGs they contain

[0125] The file "KO_predicted" in the folder 06Picrust2_stra obtained in step 2 shows the number of KOs contained in each ASV; the ASVs in the sample are located, and according to the correspondence of KOs and ARGs obtained in step 4, the KOs belonging to ARGs are screened or located, and the other KOs are deleted, and the ASVs with a value of 0 are deleted (i.e. the columns with a value of 0 and the rows with a value of 0 are deleted), to obtain the number of ARGs of the ARB (antibiotic-resistant bacteria, i.e. bacteria containing drug-resistant genes ARG) in each sample, as shown in Table 3. For example, K03585, K18138, K12340, K18299, K17836, K18139 and K08641 are exemplary ARGs with KO numbers.

[0126] Table 3

[0127]

[0128] 6. The abundance of ARB (antibiotic-resistant bacteria, i.e. bacteria containing drug-resistant genes ARG) in each sample

[0129] In Table 2 of step 3, the ASV numbers of ARB in step 5 are selected, and tabulated again to obtain the abundance of ASV numbers of ARB in each sample, as shown in Table 4. Specifically, for example, since the bacteria with ASV numbers in Table 3 in step 5 all contain ARGs, the bacteria with ASV numbers listed in Table 3 are all ARB, such as ASV11, ASV9370, ASV10030, etc., and then the rows of these ARB are selected in Table 2, or the rows that are not these ARB are deleted, i.e. Table 4 is obtained.

[0130] Table 4

[0131]

[0132] 7. Calculation of the drug resistance rate of the bacterial community

[0133] Table 4 is the abundance of ARB in each sample, and the sum of the abundance of all ARB in the sample in Table 4 can be calculated to obtain the drug resistance rate of the bacterial community in the sample (the percentage of the total ARB abundance in the sample equals the drug resistance rate). Specifically, in the A1 sample, the A1 column (the abundance of ARB of each ASV number in A1) of Table 4 is added, for example, 0.5484+2.1521+0.0567+1.3836…=36.5872, indicating that the drug resistance rate of the bacterial community in the A1 sample is 36.59%.

[0134] 8. Calculation of the abundance of drug-resistant genes ARG in ARB (antibiotic-resistant bacteria, i.e. bacteria containing drug-resistant genes ARG)

[0135] For sample A1, the ARB abundance of ASV number in column A1 of step 6 table 4 is multiplied by the corresponding ARG number in step 5 table 3 respectively, to obtain the abundance of ARG in different ARB in sample A1 (as shown in table 5). Specifically, for example, the ARB abundance of ASV number in column A1 of step 6 is 0.5484 (ASV11), 0, 0, 2.1521 (ASV7259), 0.0567 (ASV1704)……, and the ARG (K03585, K18138, K12340, K18299, K17836, K18139, K08641) in ASV11 in step 5 is 3, 3, 1, 1, 0, 2, 3……, respectively, 0.5484 is multiplied by 3, 3, 1, 1, 0, 2, 3……, respectively, to obtain the abundance of ARB in ASV11, which is 1.6452, 1.6452, 0.5484, 0.5484, 0, 1.0968……, respectively; similarly, the abundance of ARB in ASV7259 is 8.6084, 6.4563, 2.1521, 2.1521, 0, 2.1521……, respectively; the abundance of each ARG in other ARB is calculated in turn, which is table 5.

[0136] For other samples, the same method can be used to analyze and calculate the abundance of ARG in each sample.

[0137] Table 5

[0138]

[0139] The commonly used method for detecting the drug resistance rate of bacteria in the prior art is the culture medium sensitivity method, which inoculates the sample to be tested on the culture medium, so that it grows on the plate containing different antibiotics to evaluate the sensitivity to antibiotics. There is no significant difference between the drug resistance rate detected by this method and the drug resistance rate predicted by the method described in the present application. Therefore, the method described in the present application can be used to predict the drug resistance rate of bacteria in the sample.

[0140] 9. Calculation of the abundance of ARG in each sample

[0141] The abundance of ARGs in the sample is calculated by summation, and thus the abundance of ARGs in the sample can be calculated. Specifically, for example, according to Table 5, the abundance of ARGs in the ARB of the A1 sample in the ASV number is known, and the values in the K03585 column are added up to obtain the total of K03585 (an ARG) in the A1 sample, which is equal to 17.4551 in this embodiment, i.e., the abundance of K03585 (an ARG) in the A1 sample is 17.4551; the abundance of other ARGs can be calculated in the same way. The abundance of ARGs in other samples can be calculated in the same way to obtain the abundance of each ARG in each sample, as shown in Table 6.

[0142] Table 6

[0143]

[0144] The composition and abundance of drug resistance genes in the sample commonly used in the prior art are detected by metagenomic analysis and high-throughput quantitative polymerase chain reaction (qPCR), and there is no significant difference between the composition and abundance of drug resistance genes detected by this method and the composition and abundance of drug resistance genes predicted by the method described in the present application. Therefore, the method described in the present application can be used to predict the composition and abundance of drug resistance genes in the sample.

[0145] Example 2 Correlation analysis of antibiotic-resistant bacteria (ARB) and antibiotic resistance genes (ARG)

[0146] 1. The abundance information of ARGs in different samples is obtained by the method described in Example 1.

[0147] 2. The edge file (as shown in Table 7) and the point file (as shown in Table 8) are made using the data in Table 5, and the VGT network diagram of the sample is drawn by Gephi, as shown in Figure 1 . The VGT network of each sample (in order to distinguish the contribution of VGT and HGT to a specific ARG, we assume that if the microorganism carrying the most abundant ARG is also the host with the most abundant ARG, then we consider that the host is the main host of the ARG. The contribution of VGT is the proportion of the main host among all hosts, and the rest is attributed to the contribution of HGT).

[0148] Table 7

[0149]

[0150] Table 8

[0151]

[0152] 3. Position ARB in step 2 to genus, as shown in Table 9, make edge file (as shown in Table 10) and point file (as shown in Table 11), further draw VGT and HGT network diagram of each sample by Gephi software (http: / / gehi.org), as shown in Figure 1

[0153] In the study of antibiotic resistance, VGT and HGT network diagram plays an important role, which helps to understand the way of transmission of resistant genes and its impact on microbial population. Specifically, VGT refers to the process of gene transmission from parent to offspring, that is, the transmission of gene information to offspring through reproduction. VGT network diagram shows the transmission mode of resistant genes between generations of bacteria, which is crucial to understand the accumulation of resistance within the same species over time, and helps to track the evolution of strains and assess how drug resistance is gradually enhanced in offspring. HGT refers to the direct exchange of genes between different species or individuals, without relying on reproduction. In bacterial populations, HGT is the main way of rapid transmission of resistant genes, especially through plasmids, transposons or bacteriophages and other media. HGT network diagram can help to clarify the rapid transmission path of resistant genes between different bacterial species, and explain how resistance spreads rapidly between different species and even different ecological niches, leading to difficult-to-control drug-resistant infections.

[0154] Table 9

[0155]

[0156] Table 10

[0157]

[0158] Table 11

[0159]

[0160] Although the embodiments of the present application have been shown and described above, it should be understood that the above-mentioned embodiments are exemplary and should not be construed as limiting the present application, and those skilled in the art can make changes, modifications, replacements and variations to the above-mentioned embodiments within the scope of the present application.​

Claims

1. An analysis method for predicting the resistant gene composition and abundance of a bacterial community in a sample based on 16S rRNA amplicon sequencing sequences, characterized by, The method comprises the following steps: 1) 16S rRNA amplicon sequencing of all bacteria in the sample to obtain 16S rRNA sequences of all bacteria in the sample; 2) obtaining the species and abundance of bacteria contained in the sample, which comprises analyzing the 16S rRNA sequences by using the PICRUSt2 software, outputting the species and abundance of ASV numbered bacteria contained in the sample; 3) obtaining the drug-resistant genes and their abundance contained in various bacteria, which comprises analyzing the 16S rRNA sequences by using the PICRUSt2 software, outputting the functional genes and their number contained in various ASV numbered bacteria, and screening and retaining the drug-resistant genes and their number contained in various ASV numbered bacteria by using a drug-resistant gene list; 4) based on the species and abundance of ASV numbered bacteria and the drug-resistant genes and their number contained in various ASV numbered bacteria obtained in steps 2) and 3), the species and abundance of drug-resistant genes in the sample are obtained by calculation; The calculation in step 4) is: multiplying the abundance of ASV numbered bacteria obtained in step 2) and the number of drug-resistant genes contained in the ASV numbered bacteria obtained in step 3) to obtain the abundance of the drug-resistant gene contained in the ASV numbered bacteria in the sample; and adding up the abundance of the drug-resistant gene contained in all ASV numbered bacteria in the sample to obtain the abundance of the drug-resistant gene in the sample, and further obtaining the abundance of all drug-resistant genes in the sample by the same method.

2. The method of claim 1, wherein, It also comprises the step of using QIIME2 software analysis to obtain ASV of bacteria, and obtaining the annotation of bacteria by Greengene2 database comparison.

3. The method of claim 1 or 2, characterized in that The 16S rRNA sequence is a full-length or V3-V4 region or V4 region 16S rRNA sequence.

4. The method of claim 1, wherein, Steps 2) and 3) further comprise assigning the ASV numbered bacterial species as ASV order.

5. The method of claim 1, wherein, The drug-resistant gene list in step 3) is obtained from the KEGG database.

6. A method of predicting the resistance rate of a bacterial community in a sample, characterized in that, Using the steps of the method of any one of claims 1-5, the species and abundance of ASV numbered bacteria containing drug-resistant genes in the sample are obtained, and the abundance of all ASV numbered bacteria containing drug-resistant genes in the sample is added and multiplied by 100% to obtain the drug resistance rate of the sample.

Images