CN117174165B - 基于宏基因组的环境耐药组分析方法 - Google Patents
基于宏基因组的环境耐药组分析方法 Download PDFInfo
- Publication number
- CN117174165B CN117174165B CN202311389711.8A CN202311389711A CN117174165B CN 117174165 B CN117174165 B CN 117174165B CN 202311389711 A CN202311389711 A CN 202311389711A CN 117174165 B CN117174165 B CN 117174165B
- Authority
- CN
- China
- Prior art keywords
- sequence
- gene
- sequences
- functional
- contig
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000007613 environmental effect Effects 0.000 title claims abstract description 103
- 238000004458 analytical method Methods 0.000 title claims abstract description 70
- 206010059866 Drug resistance Diseases 0.000 title claims abstract description 35
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 333
- 238000000034 method Methods 0.000 claims abstract description 83
- 238000012546 transfer Methods 0.000 claims abstract description 42
- 230000000813 microbial effect Effects 0.000 claims abstract description 38
- 238000012502 risk assessment Methods 0.000 claims abstract description 13
- 239000000203 mixture Substances 0.000 claims description 31
- 108700026244 Open Reading Frames Proteins 0.000 claims description 26
- 230000003115 biocidal effect Effects 0.000 claims description 20
- 108020004465 16S ribosomal RNA Proteins 0.000 claims description 18
- 210000004027 cell Anatomy 0.000 claims description 18
- 102000004169 proteins and genes Human genes 0.000 claims description 16
- 239000013612 plasmid Substances 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 11
- 238000003908 quality control method Methods 0.000 claims description 11
- 150000007523 nucleic acids Chemical group 0.000 claims description 10
- 101000685990 Homo sapiens Specifically androgen-regulated gene protein Proteins 0.000 claims description 9
- 102100023355 Specifically androgen-regulated gene protein Human genes 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 9
- 230000002068 genetic effect Effects 0.000 claims description 9
- 229910052751 metal Inorganic materials 0.000 claims description 9
- 239000002184 metal Substances 0.000 claims description 9
- 238000012165 high-throughput sequencing Methods 0.000 claims description 8
- 230000007923 virulence factor Effects 0.000 claims description 8
- 239000000304 virulence factor Substances 0.000 claims description 8
- 108091028043 Nucleic acid sequence Proteins 0.000 claims description 6
- 230000002759 chromosomal effect Effects 0.000 claims description 6
- 238000011109 contamination Methods 0.000 claims description 6
- 230000002070 germicidal effect Effects 0.000 claims description 6
- NUFKRGBSZPCGQB-FLBSXDLDSA-N (3s)-3-amino-4-oxo-4-[[(2r)-1-oxo-1-[(2,2,4,4-tetramethylthietan-3-yl)amino]propan-2-yl]amino]butanoic acid;pentahydrate Chemical compound O.O.O.O.O.OC(=O)C[C@H](N)C(=O)N[C@H](C)C(=O)NC1C(C)(C)SC1(C)C.OC(=O)C[C@H](N)C(=O)N[C@H](C)C(=O)NC1C(C)(C)SC1(C)C NUFKRGBSZPCGQB-FLBSXDLDSA-N 0.000 claims description 4
- 235000019409 alitame Nutrition 0.000 claims description 4
- 238000003556 assay Methods 0.000 claims description 4
- 230000008030 elimination Effects 0.000 claims description 4
- 238000003379 elimination reaction Methods 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 4
- 210000000349 chromosome Anatomy 0.000 claims description 3
- 238000011160 research Methods 0.000 abstract description 10
- 230000009897 systematic effect Effects 0.000 abstract description 4
- 208000020990 adrenal cortex carcinoma Diseases 0.000 description 40
- 241000894007 species Species 0.000 description 17
- 244000005700 microbiome Species 0.000 description 14
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 14
- 238000011156 evaluation Methods 0.000 description 12
- 238000004590 computer program Methods 0.000 description 11
- 238000007405 data analysis Methods 0.000 description 10
- 238000012216 screening Methods 0.000 description 10
- 238000012163 sequencing technique Methods 0.000 description 9
- 239000010865 sewage Substances 0.000 description 9
- 239000003814 drug Substances 0.000 description 8
- 229940079593 drug Drugs 0.000 description 8
- 108020004414 DNA Proteins 0.000 description 6
- 239000002351 wastewater Substances 0.000 description 6
- 238000001514 detection method Methods 0.000 description 4
- 238000011002 quantification Methods 0.000 description 4
- 238000004065 wastewater treatment Methods 0.000 description 4
- 238000007400 DNA extraction Methods 0.000 description 3
- 229940126575 aminoglycoside Drugs 0.000 description 3
- 230000001580 bacterial effect Effects 0.000 description 3
- 230000000052 comparative effect Effects 0.000 description 3
- 238000012136 culture method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 239000013049 sediment Substances 0.000 description 3
- 239000002689 soil Substances 0.000 description 3
- 229940124530 sulfonamide Drugs 0.000 description 3
- 150000003952 β-lactams Chemical class 0.000 description 3
- 238000012408 PCR amplification Methods 0.000 description 2
- 239000004098 Tetracycline Substances 0.000 description 2
- 108010020764 Transposases Proteins 0.000 description 2
- 229940125364 angiotensin receptor blocker Drugs 0.000 description 2
- 239000003153 chemical reaction reagent Substances 0.000 description 2
- 239000010842 industrial wastewater Substances 0.000 description 2
- 230000000968 intestinal effect Effects 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 239000012528 membrane Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000002243 precursor Substances 0.000 description 2
- 235000019515 salmon Nutrition 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000004062 sedimentation Methods 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 150000003456 sulfonamides Chemical class 0.000 description 2
- 235000019364 tetracycline Nutrition 0.000 description 2
- 150000003522 tetracyclines Chemical class 0.000 description 2
- 241000972773 Aulopiformes Species 0.000 description 1
- 241000894006 Bacteria Species 0.000 description 1
- 238000007399 DNA isolation Methods 0.000 description 1
- 101000694320 Drosophila melanogaster RuvB-like helicase 2 Proteins 0.000 description 1
- 241000192128 Gammaproteobacteria Species 0.000 description 1
- 241000736262 Microbiota Species 0.000 description 1
- 102000008579 Transposases Human genes 0.000 description 1
- 239000011157 advanced composite material Substances 0.000 description 1
- 239000003570 air Substances 0.000 description 1
- 239000002333 angiotensin II receptor antagonist Substances 0.000 description 1
- 239000003242 anti bacterial agent Substances 0.000 description 1
- 230000000844 anti-bacterial effect Effects 0.000 description 1
- 244000052616 bacterial pathogen Species 0.000 description 1
- 239000003899 bactericide agent Substances 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007622 bioinformatic analysis Methods 0.000 description 1
- 238000003766 bioinformatics method Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 229920002678 cellulose Polymers 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 239000000356 contaminant Substances 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000003344 environmental pollutant Substances 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012634 fragment Substances 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 229910001385 heavy metal Inorganic materials 0.000 description 1
- 230000002519 immonomodulatory effect Effects 0.000 description 1
- 230000008102 immune modulation Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000004899 motility Effects 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000008261 resistance mechanism Effects 0.000 description 1
- 239000010802 sludge Substances 0.000 description 1
- 238000010561 standard procedure Methods 0.000 description 1
- 229960005322 streptomycin Drugs 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- FDDDEECHVMSUSB-UHFFFAOYSA-N sulfanilamide Chemical compound NC1=CC=C(S(N)(=O)=O)C=C1 FDDDEECHVMSUSB-UHFFFAOYSA-N 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 229960002180 tetracycline Drugs 0.000 description 1
- 229930101283 tetracycline Natural products 0.000 description 1
- 229940040944 tetracyclines Drugs 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Landscapes
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本申请涉及生物信息学领域，具体提出了一种基于宏基因组的环境耐药组分析方法，包括：基于环境样品，获得所述环境样品的宏基因组数据；和基于所述宏基因组数据，识别功能基因，并在读长序列、重叠群序列和微生物基因组草图三个层次上分别计算所述功能基因在所述环境样品中的相对丰度；同时，本方法还包括了基于功能基因的宿主分析及侧向基因转移潜力分析。本申请提出的环境耐药组分析方法为环境耐药组研究及风险评估提供了***性的技术解决方案。
Description
技术领域
本申请涉及生物信息学领域，具体涉及一种基于宏基因组的环境耐药组分析方法。
背景技术
环境是抗生素耐药性传播扩散的一个关键节点，抗生素抗性基因（antibioticresistance gene，ARGs）已经在水体、沉积物、土壤、污泥和空气等不同生境中被频繁检出，对生态和人类健康都造成了较大的威胁。ARG由此被视为一种新型的环境污染物，其来源、分布和传播规律越来越受到研究者的重视。相关技术中，许多方法被用来检测ARGs，其中基于高通量测序技术的宏基因组学方法克服了传统方法（如细菌培养方法、基于PCR的分子生物学方法等）的局限性，可以用于解析环境样品中所有微生物的功能基因信息，为全面表征环境耐药基因谱提供了技术支撑。
然而，以Illumina为代表的二代测序技术会产生大量短读长（reads）数据，如何有效地利用生物信息学分析手段组织测序数据仍是环境耐药组研究面临的巨大挑战。相关技术中，针对环境耐药组分析往往仅集中于ARG的组成分析，并且涉及的分析工具和数据库繁多，使得数据分析方法单薄，且尚未有整体性的高效数据分析方法。
由此，亟待开发一种合理高效的生物信息数据分析流程以为环境耐药组研究及风险评估提供***性的技术解决方案。
发明内容
本申请至少从以下方面解决了相关技术的问题之一。
为此，本申请的实施例提供了一种基于宏基因组的环境耐药组分析方法，包括：基于环境样品，获得所述环境样品的宏基因组数据；和基于所述宏基因组数据，识别功能基因，并计算所述功能基因在所述环境样品中的相对丰度，包括：a. 基于所述宏基因组数据，对读长序列（reads）进行比对和注释，以识别所述功能基因，并计算所述功能基因在所述环境样品中的第一丰度；b. 对所述读长序列进行组装以获得重叠群序列（contigs），并基于所述重叠群序列进行比对和注释，以识别所述功能基因，并计算所述功能基因在所述环境样品中的第二丰度，可选地，所述重叠群序列的组装长度大于500 bp；和c. 对所述重叠群序列进行分箱以获得微生物基因组草图（MAGs），并基于所述微生物基因组草图进行比对和注释，以识别所述功能基因，并计算所述功能基因在所述环境样品中的第三丰度。
在一些实施例中，所述功能基因包括以下中的一种或多种：抗生素抗性基因、可移动遗传元件、毒力因子基因、和杀菌剂和金属抗性基因，优选为抗生素抗性基因。
在一些实施例中，步骤b中，基于所述重叠群序列进行比对和注释，以识别所述功能基因具体包括：预测所述重叠群序列的开放阅读框，以获得所述重叠群序列对应基因的核酸序列及其翻译后的蛋白序列；基于所述重叠群序列对应基因的核酸序列及其翻译后的蛋白序列的序列相似度对所述重叠群序列对应基因进行聚类，计算所述重叠群序列对应基因的基因丰度并获得非冗余基因集；和基于功能基因数据库，对所述非冗余基因集进行比对和注释，以识别所述功能基因，并获得所述功能基因的组成信息和丰度信息。在一些实施例中，所述比对中，所述非冗余基因集中e-value≤10-5、序列相似性≥80%、覆盖度≥70%的基因序列被识别为所述功能基因。
在一些实施例中，所述功能基因数据库选自由以下中的一个或多个：SARG、MGEdatabase、VFDB和BacMet。
在一些实施例中，步骤c中，所述基于所述微生物基因组草图进行比对和注释，以识别所述功能基因还包括：对所述微生物基因组草图进行去冗余和质量控制，以获得非冗余高质量的微生物基因组草图；预测所述非冗余高质量的微生物基因组草图的开放阅读框，以获得所述非冗余高质量的微生物基因组草图中各对应基因的核酸序列及其翻译后的蛋白序列；和基于功能基因数据库，对所述各对应基因的翻译后的蛋白序列进行比对和注释，以识别所述功能基因，并获得所述功能基因的组成信息和丰度信息。在一些实施例中，基于完整度和污染率对所述微生物基因组草图进行所述质量控制，其中要求所述完整度≥50%并且所述污染率≤10%。
在一些实施例中，所述方法还包括：基于所述功能基因进行宿主分析，包括：i. 识别所述读长序列中的功能基因相似序列，对所述功能基因相似序列进行物种注释，以获得功能基因相似序列在各物种的序列数作为第一宿主丰度，其中所述功能基因相似序列与所述功能基因的序列相似度为≥80%，比对长度≥25 aa，可选地，e值≤10-7；ii. 根据上述任一实施例中所述的方法，识别所述重叠群序列中含有功能基因序列的重叠群序列，对所述含有功能基因序列的重叠群序列进行物种注释，以获得所述含有功能基因序列的重叠群序列在各物种的序列数作为第二宿主丰度；和iii. 根据上述任一实施例中所述的方法，识别所述微生物基因组草图中含有功能基因序列的微生物基因组草图，对所述含有功能基因序列的微生物基因组草图进行物种注释，以获得所述功能基因序列的微生物基因组草图在各物种的序列数作为第三宿主丰度。
在一些实施例中，所述方法还包括：基于所述功能基因进行侧向基因转移潜力分析，包括：x. 根据上述任一实施例中所述的方法，识别所述重叠群序列中的染色体来源序列和质粒来源序列，并获得所述染色体来源序列和质粒来源序列所携带的所述功能基因的组成信息；y. 根据上述任一实施例中所述的方法，识别所述重叠群序列中可能发生侧向基因转移的重叠群序列，并获得所述可能发生侧向基因转移的重叠群序列所携带的所述功能基因的组成信息；和z. 根据上述任一实施例中所述的方法，基于所述微生物基因组草图和所述含有功能基因序列的微生物基因组草图，预测所述含有功能基因序列的微生物基因组草图之间、和所述含有功能基因序列的微生物基因组草图与不含有所述功能基因序列的微生物基因组草图之间的侧向基因转移过程，构建“供体-受体”配对关系的侧向基因转移网络。
在一些实施例中，所述功能基因在所述环境样品中的第一丰度的计算基于16SrRNA基因序列数和细胞数量，计算公式为：
（1）
（2）
（3）
其中Abundance16S表示用16S rRNA基因序列进行标准化的功能基因丰度，Abundancecell表示用细胞数标准化的功能基因丰度，Cnumber表示每个环境样品的细胞数量；Ni(mapped reads)是宏基因组reads中比对到某条特定基因参考序列的序列条数；Lreads是高通量测序获得的序列长度；Li(reference reads)是某条特定基因参考序列的长度；N16S是宏基因组reads中被鉴定为16S rRNA基因的序列条数；L16S是所有16S rRNA基因序列的平均长度；N16S copy number是微生物细胞中16S rRNA基因序列的平均拷贝数；n是所述环境样品中被注释到属于某个基因的参考序列个数。
在一些实施例中，步骤ii还包括：利用软件CoverM计算所述第二宿主丰度，其中取所述含有功能基因序列的重叠群序列在各物种的序列数的平均序列数（mean）为所述第二宿主丰度；并且步骤iii还包括：使用软件metaWRAP (v1.2.1) 的quant_bins模块计算所述第三宿主丰度，其中取每百万序列中基因组拷贝数（Genome copies/ppm reads）表示所述第三宿主丰度。
在一些实施例中，步骤x还包括：对所述重叠群序列进行过滤，以去除序列长度短于1000 bp的所述重叠群序列；和使用Plasflow（v1.1）软件对过滤后的重叠群序列进行分型，以识别出所述过滤后的重叠群序列中的染色体来源序列和质粒来源序列，并且所述步骤y还包括：使用WAAFLE软件识别所述重叠群序列中可能发生侧向基因转移的重叠群序列，其中基于序列长度大于500 bp使用WAAFLE软件筛选出不同宏基因组数据集中发生侧向基因转移的所述重叠群序列。
本申请实施例还提出了一种基于宏基因组的环境耐药组风险评价方法，所述方法包括：根据本申请第一方面任一实施例所述的基于宏基因组的环境耐药组分析方法获得环境样品的功能基因分别在读长序列、重叠群序列和微生物基因组草图水平上的组成信息和丰度信息，其中所述功能基因包括以下中的一种或多种：抗生素抗性基因、可移动遗传元件、毒力因子基因、和杀菌剂和金属抗性基因，优选为抗生素抗性基因；将组装获得的所有重叠群序列的开放阅读框与功能基因数据库进行比对，以获得第一重叠群序列，其中所述第一重叠群序列含有所述抗生素抗性基因、可移动遗传元件、毒力因子基因和杀菌剂和金属抗性基因中的任一种功能基因；将所述第一重叠群序列的开放阅读框与所述功能基因数据库进行比对，以获得第二重叠群序列，其中所述第二重叠群序列含有其它三种功能基因中的至少一种；统计所述第一重叠群序列和所述第二重叠群序列中含有的功能基因组合类型及其对应的重叠群序列数量，并对所述重叠群序列数量赋予权重；和基于所述权重，计算得到所述环境样品的环境耐药组风险指数。在一些实施例中，对所述权重进行加和以得到所述环境样品的环境耐药组风险指数，其中所述功能基因数据库选自由以下中的一个或多个：SARG、MGE database、ACLAME、VFDB和BacMet。
在一些实施例中，所述第一重叠群序列和所述第二重叠群序列中含有的功能基因组合类型及其所赋权重如下表所示：
所述风险指数的计算公式为：
（1），
其中：
（2）
（3）
（4）
nContig为所述组装获得的所有重叠群序列的总数。
本申请实施例提出的基于宏基因组的环境耐药组分析方法及其相关应用实现了如下有益效果：
（1）本发明围绕基因注释和定量、宿主信息识别、水平基因转移分析等研究内容，选用性能较为稳定和高效的软件，建立了“短序列（读长序列，reads）-重叠群序列（contigs）-基因组序列（MAGs）”多层次的生物信息学数据分析流程，根据不同需求将复杂的分析流程进行模块化，提高了应用的便捷性，同时输出更多样的分析结果，为后续个性化统计分析提供可能。
（2）本发明覆盖了更多与抗生素耐药性密切相关的功能基因，分别提供了基于非拼接序列或组装序列的基因比对和定量方法，既可以进行环境样本的大规模快速筛查，又可以严格判别目标基因的赋存状况，降低基因鉴定的假阳性，为使用者提供多样化的选择。
（3）在ARGs宿主鉴别方面，本发明提出了三种不同的宏基因组分析策略，特别是一种基于短序列直接识别环境中ARGs宿主的策略，与宏基因组组装和分箱策略相比，具有简化计算流程、获得的宿主多样性高等优势。
（4）在ARGs的水平基因转移筛查方面，本发明分别针对组装序列和基因组序列提出了相应的分析方法，构建了微生物组水平基因转移网络，为从基因和群落层面预测ARGs的传播过程提供了解决方案。
（5）本发明基于多种风险指示基因在宏基因组序列上的共存关系，提供了一种综合风险指数的计算方法，可以用于评价不同环境样本中耐药组的潜在风险水平，具有高通量、快速简便、可标准化等特点，节约了时间、人力和相关实验成本，适合大规模的环境调查研究，为建立环境耐药组风险的评价标准提供依据。
附图说明
为了更清楚地说明本申请实施例中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为根据本申请实施例1-3的基于宏基因组的环境耐药组分析方法流程图；
图2为根据本申请实施例4的环境耐药组风险评价方法流程图；
图3为根据本申请实施例1的基于reads分析的不同基因组成和相对丰度，(a)ARG; (b) BMG; (c) MGE; (d) VFG；
图4为根据本申请实施例1的基于contigs分析的不同基因组成和相对丰度. (a)ARG; (b) BMG; (c) MGE; (d) VFG；
图5为根据本申请实施例1的基于MAGs分析的不同基因组成和数量. (a) ARG;(b) BMG; (c) MGE; (d) VFG；
图6为根据本申请实施例2的基于reads分析的ARGs宿主组成；
图7为根据本申请实施例2的基于contigs分析的ARGs宿主组成；
图8为根据本申请实施例2的基于MAGs分析的ARGs宿主组成；
图9为根据本申请实施例3的不同ARGs的遗传背景分析；
图10为根据本申请实施例3的基于contigs的水平基因转移事件识别；
图11为根据本申请实施例3的基于MAGs的水平基因转移网络. (a) LH1; (b)LH2；
图12为根据本申请实施例4的BMG权重为0、0.5和1时计算所得的耐药组风险指数展示图；
图13为根据本申请实施例4的环境耐药组风险评价方法与对比例中现有方法的评价结果比较。
具体实施方式
下面结合具体实施方式对本发明进行进一步的详细描述，给出的实施例仅为了阐明本发明，并非限制本发明的范围。以下提供的实施例可作为本技术领域普通技术人员进行进一步改进的指南，并不以任何方式构成对本发明的限制。
实施例中未注明具体技术或条件的，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品。
本申请是发明人基于以下认识做出的：
相关技术中，有多种方法被用以检测ARGs，包括传统的细菌培养方法、基于PCR扩增的分子生物学方法和宏基因组学方法等。然而其中细菌培养方法既耗时又容易受培养条件的影响；基于PCR扩增的方法会受到引物序列可得性的限制，无法有效追踪环境中极为丰富多样的ARGs，因此很难满足环境耐药组（基因组或微生物群落中所有的ARGs及其前体）的相关研究。基于高通量测序技术的宏基因组学方法克服了这些传统方法的局限性，包括提取环境样本中的总DNA并对其进行高通量测序，主要是通过BLAST算法将获得的宏基因组序列数据与特定的抗生素抗性基因数据库（如ARDB、CARD等）进行比对，从而获知环境样本中ARGs的多样性组成。然而，现有宏基因组学方法和工具较为分散，缺乏合理、高效、***的数据分析流程，无法满足环境耐药组研究内容的多目标需求。
为此，本申请实施例提供了一整套基于宏基因组短序列、重叠群序列和基因组序列的生物信息学分析流程，涵盖四类风险指示基因（即ARG、可移动遗传元件（MGE）、毒力因子（VFG）、杀菌剂和金属抗性基因（BMG））的注释和定量、耐药基因宿主识别以及水平基因转移事件筛查等分析。本申请实施例提出的该分析流程围绕基因注释和定量、宿主信息识别、水平基因转移分析等研究内容，选用性能较为稳定和高效的软件，建立了“短序列-重叠群序列-基因组序列”多层次的生物信息学数据分析流程，根据不同需求将复杂的分析流程进行模块化，提高了应用的便捷性，同时输出更多样的分析结果，为后续个性化统计分析提供可能；同时，该分析流程覆盖了更多与抗生素耐药性密切相关的功能基因，分别提供了基于非拼接序列或组装序列的基因比对和定量方法，既可以进行环境样本的大规模快速筛查，又可以严格判别目标基因的赋存状况，降低基因鉴定的假阳性，为使用者提供多样化的选择。此外，本申请实施例的分析流程还包括ARGs宿主鉴别及其水平基因转移筛查方法，其中在ARGs宿主鉴别方面，本申请实施例的分析流程提出了三种不同的宏基因组分析策略，特别是一种基于短序列直接识别环境中ARGs宿主的策略，其与宏基因组组装和分箱策略相比，具有简化计算流程、获得的宿主多样性高等优势；同时，在ARGs的水平基因转移筛查方面，分别针对组装序列和基因组序列提出了相应的分析方法，构建了微生物组水平基因转移网络，为从基因和群落层面预测ARGs的传播过程提供了解决方案。
本申请实施例中，“环境耐药组”是指宏基因组或微生物群落中所有的ARGs及其前体。可以理解的是，通过对宏基因组或微生物群落中的更多的功能基因（MGE、VFG、BMG等）进行分析，可以对该环境的耐药性基因的扩散传播进行有效分析和评估。
本申请实施例中，“环境样品”为所分析的宏基因组数据的DNA的提取来源。在本公开实施例中，环境样品可以为肠道微生物、水源微生物、土壤微生物等，样本来源对应可以为肠道、水源、土壤等。
本申请实施例中，“功能基因在环境样品中的相对丰度”是指含有该功能基因的序列在样品中的数量值，其中含有该功能基因的序列可以为含有该功能基因片段或全长的reads、contigs和MAGs。可以理解的是，通过在reads、contigs和MAGs多个层次统计计算功能基因的丰度，可以对功能基因在该环境中的赋存情况进行整体性评价，且有效降低了基因鉴定的假阳性。
本申请实施例中，“功能基因的组成信息”是指该功能基因和可选的其相似序列在环境样品中的所有存在形式和占比。
本申请第一方面实施例提出了一种基于宏基因组的环境耐药组分析方法，参考图1，该方法可以包括以下步骤：
（1）采集环境样品，使用标准方法提取和纯化样品中的总DNA样品，进行宏基因组高通量测序，并对原始测序数据（raw reads）进行质量控制，获得高质量的清洁序列（cleanreads）；
（2）对clean reads进行与功能基因数据库的比对和注释，识别抗生素抗性基因（ARG）、可移动遗传元件（MGE）、毒力因子基因（VFG）和杀菌剂和金属抗性基因（BMG）等功能基因，并计算其在每个环境样品中的相对丰度（即第一丰度）；
（3）对每个环境样品的clean reads进行组装，获得重叠群序列（contigs）；
（4）预测contigs的开放阅读框（open reading frame，ORF），获得各功能基因的核酸序列及其翻译后的蛋白序列，根据序列相似性对各功能基因进行聚类，并计算各功能基因的相对丰度，获得非冗余基因集；
（5）将非冗余基因集与不同的功能数据库进行比对和注释，获得ARG、MGE、VFG和BMG等功能基因的组成信息和丰度信息（即第二丰度）；
（6）对每个环境样品的contigs进行分箱（binning），获得微生物基因组草图（MAGs），对MAGs进行去冗余，并根据完整度和污染率对MAGs进行质量评估，筛选出符合要求的非冗余高质量MAGs；
（7）对非冗余高质量MAGs的ORF进行预测，获得功能基因的核酸序列及其翻译后的蛋白序列，将蛋白序列与抗生素抗性基因数据库进行比对，识别出含有ARG序列的MAGs（ARG-carrying MAGs，ACMs），并统计ACMs携带的ARG种类及数量（即第三丰度）。
基于此，本申请实施例进一步提供了在reads、contigs和MAGs多个层次上的基于所述功能基因进行宿主分析策略（8），包括：
（8-1）识别clean reads中的ARG相似序列（ARG-like reads，ALRs），对ALRs进行物种注释，并统计物种序列数即为丰度（即第一宿主丰度）；
（8-2）根据步骤（4）和（5）所述的方法，识别含有ARG序列的contigs（ARG-carryingcontigs，ACCs），对ACCs进行物种注释，并计算每个ACCs的相对丰度（即第二宿主丰度）；
（8-3）对步骤（7）所述的ACMs进行物种注释，并计算每个ACMs的相对丰度（即第三宿主丰度）。
进一步地，本申请实施例还提出了在reads、contigs和MAGs多个层次上的基于所述功能基因进行侧向基因转移潜力分析方法（9），包括：
（9-1）从所有contigs中识别染色体序列和质粒序列，并根据步骤（4）和（5）所述的方法，分析两类序列所携带的ARGs组成；
（9-2）从所有contigs中识别可能发生水平基因转移的contigs，并根据步骤（4）和（5）所述的方法，分析该可能发生水平基因转移的contigs所携带的ARGs组成；
（9-3）基于步骤（6）和（7）所述的MAGs和ACMs，预测ACMs之间、以及ACMs与非ACMs（nACMs）之间潜在的水平基因转移过程，构建基于“供体-受体”配对关系的基因转移网络。在一些实施例中，可以进一步对其进行可视化和网络分析。
在一些实施例中，步骤（1）中，可以使用FastQC（v0.11.8）检查raw reads的整体质量，使用KneadData软件（https://bitbucket.org/biobakery/kneaddata）实现序列的质量控制，包括利用软件Trimmomactic （v0.39）去除序列的双端接头序列和低质量序列，通过软件Bowtie2（v2.3.5）去除宿主污染。
在一些实施例中，步骤（2）中，使用ARGs-OAP（v2.0）软件获得样品中ARGs、MGEs、VFGs和BMGs的组成，分别使用Ublast（v11.0.667）工具从clean reads中预先筛选目标基因相似序列，再使用Blastx（v2.2.28+）工具将上述相似序列与功能基因数据库进行比对。
在一些实施例中，步骤（2）中，所述的功能基因数据库具体可以为SARG（v2.2，https://github.com/biofuture/Ublastx_stageone）、MGE database (https://github.com/KatariinaParnanen/MobileGeneticElementDatabase)、VFDB （VFDB_setA_prot.fas，http://www.mgc.ac.cn/VFs/main.htm）和BacMet（v2.0, http://bacmet.biomedicine.gu.se/）。可以理解的是，可以选择其它功能基因数据库进行比对和注释分析，只要能够实现样品中ARGs、MGEs、VFGs和BMGs的有效检出并获得其组成信息和丰度信息即可，本申请对所使用的功能基因数据库不做限制。
在一些实施例中，步骤（2）中，功能基因丰度的计算公式为（1）-（3）：
（1）
（2）
（3）
式中：Abundance16S表示用16S rRNA基因序列进行标准化的基因丰度，Abundancecell表示用细胞数标准化的基因丰度，Cnumber表示每个样品的细胞数量；Ni(mapped reads)是宏基因组clean reads中比对到某条特定基因参考序列的序列条数；Lreads是高通量测序获得的序列长度；Li(reference reads)是某条特定基因参考序列的长度；N16S是宏基因组clean reads中被鉴定为16S rRNA基因的序列条数；L16S是所有16S rRNA基因序列的平均长度；N16S copy number是微生物细胞中16S rRNA基因序列的平均拷贝数；n 是样品中被注释到属于某个基因的参考序列个数。
在一些实施例中，步骤（3）中，用MEGAHIT（v1.1.3）软件对clean reads进行denovo拼接，获得重叠群序列 (contigs)，筛选拼接长度在500 bp以上的contigs。
在一些实施例中，步骤（4）中，使用Prodigal（v2.6.3) 软件对contigs进行ORF预测，使用CD-Hit (v4.8.1) 软件对获得的ORF进行聚类，获得非冗余基因集；使用Salmon(v0.13.1) 软件计算出每个基因的reads count，并进一步标准化，获得TPM (transcriptsper kilobase million) 值，用以表征基因在各个样品中的丰度信息。
在一些实施例中，步骤（5）中，使用Blastp（v2.6.0）工具将ORF对应的蛋白序列与功能基因数据库进行比对。
在一些实施例中，步骤（6）中，使用metaWRAP (v1.2.1)软件提取（binning模块）、纯化（bin_refinement模块）各个样品中的微生物基因组草图 (MAGs)，设置MAGs的筛选阈值为污染度≤10%、完整度≥50%，使用dRep (v2.6.2) 软件对获得的MAGs进行去冗余。
在一些实施例中，步骤（7）中，使用metaWRAP (v1.2.1)软件annotate_bins模块预测MAGs的基因序列，使用Blastp（v2.6.0）工具将基因序列与所述抗生素抗性基因数据库进行比对。
在一些实施例中，步骤（8-1）中，使用Ublast（v11.0.667）工具从clean reads中筛选出ARG相似序列（ALRs），使用Kraken2（v2.0.8-beta）软件和GTDB（r89）数据库进行物种注释。
在一些实施例中，步骤（8-2）中，使用Kraken2（v2.0.8-beta）软件和GTDB（r89）数据库对所有ACCs进行物种注释，使用CoverM（v0.6.1）软件计算每条ACCs的相对丰度。
在一些实施例中，步骤（8-3）中，使用GTDB-Tk（v1.0.2）软件 和GTDB（r89）数据库对MAGs进行物种注释，使用metaWRAP（v1.2.1）软件的quant_bins模块计算每个MAGs的相对丰度。
在一些实施例中，步骤（9-1）中，使用Plasflow （v1.1）软件对contigs进行分型（将序列长度小于1000bp的contigs排除在外），识别出染色体序列和质粒序列。
在一些实施例中，步骤（9-2）中，使用WAAFLE软件（https://github.com/biobakery/waafle）识别可能发生水平基因转移事件的contigs。
在一些实施例中，步骤（9-3）中，使用MetaCHIP（v1.10.10）软件预测不同MAGs之间潜在的水平基因转移事件，筛选“ACM-ACM”和“ACM-nACM”的配对关系，构建有向基因转移网络。
在一些实施例中，步骤（9-3）中，使用Gephi（v0.9.2）软件进行网络可视化。
由此，本申请实施例中所提出的基于宏基因组的环境耐药组分析方法涵盖了reads、contigs、genome（MAGs）多个层次的多个与耐药性相关的功能基因的分析，并在其中提供了有关宿主分析和水平基因转移的分析流程模块，解决了现有技术中仅针对ARGs基因的分析单薄难以有效支撑该环境中的耐药基因的整体分析的问题；同时本申请实施例提出的方法对繁多的分析工具和数据库进行了流程性整合，提高了应用的便捷性，同时输出更多样的分析结果，有效实现了针对环境样品的多层次多基因的整体性的高效数据分析。
本申请第二方面实施例还提出了一种基于宏基因组的环境耐药组风险评价方法，其中该方法基于上述第一方面任一实施例中提出的宏基因组的环境耐药组分析方法所获得的分析结果，对该分析结果进行了进一步整合。图2为根据本申请实施例的环境耐药组风险评价方法流程图，参考图2，在一些实施例中，该风险评价方法具体可以包括：
（1）从环境样品中提取总DNA进行宏基因组测序，对原始数据（raw reads）进行质量控制并去除宿主（如人类）序列，得到高质量的清洁序列（clean reads）；
（2）对每个样品的clean reads进行组装，获得重叠群序列（contigs），并统计每个样品的contigs数量；
（3）对所有contigs预测开放阅读框（ORF），将所有ORF与ARG数据库进行比对，识别出含有ARG序列的contigs（ACCs），并统计每个样品中ACCs的数量；
（4）将所有ACCs的ORF与MGE、VFG、BMG等相关功能数据库进行比对，识别含有MGE、VFG和BMG序列的ACCs，统计含有不同基因组合的ACCs数量；
（5）根据所含基因的类型，对不同的ACCs数量赋予相应的权重，并进行加和，结合每个样品中所有contigs的数量进行标准化，计算得到每个样品的综合风险指数，即可代表环境耐药组风险水平。
在一些实施例中，步骤（1）-（4）的具体操作流程、使用软件和参数及功能基因数据库可参照本申请第一方面实施例中的相应流程，在此不再赘述。
在一些实施例中，步骤（5）中，筛选出8种类型的ACCs，根据含有目标基因的类型数对不同的ACCs数量进行赋权，其中ARG、MGE和VFG表征了环境耐药组的直接风险，权重为1，而BMG则表征了环境耐药组的间接风险，权重为0.5。每个样品的综合风险指数根据公式（1）-（4）计算：
（1）
（2）
（3）
（4）
式中：nContig表示contig数据集中contigs的总数，α表示ACCs的数量，β1表示携带MGE的ACCs数量，β2表示携带VFG的ACCs数量，β3表示携带BMG的ACCs数量，γ1表示携带MGE和VFG的ACCs数量，γ2表示携带MGE和BMG的ACCs数量，γ3表示携带BMG和VFG的ACCs数量，δ1表示同时携带MGE、VFG和BMG的ACCs数量。
可以理解的是，本申请第二方面实施例中提出的基于宏基因组的环境耐药组风险评价方法，基于宏基因组的环境耐药组分析方法的分析结果，并通过对含有不同功能基因类型的contigs数量赋予合适的权重，为环境样品中的耐药基因的存在和转移提供了有效的评价指标（即环境耐药组风险指数）；同时，本申请实施例提出的风险评价方法并不仅仅根据ARGs的相对丰度、移动性和病原菌内的存在情况对ARGs进行风险分级再对环境样本的耐药性风险进行综合评价，而是基于ARGs及更多的相关功能基因的存在及其类型并对其进行合理赋权以直接进行风险评价，由此照顾到了抗生素耐药性与微生物其他抗性机制的共选择过程，实现了更为完整和准确的环境耐药组的风险评价。
下述实施例中的实验方法，如无特殊说明，均为常规方法，按照本领域内的文献所描述的技术或条件或者按照产品说明书进行。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
如无特殊说明，以下实施例中的定量分析试验，均设置三次重复实验，结果取平均值。如无特殊说明，以下实施例中的数据分析流程中所使用的参数及命令行均为对应软件说明书中的默认参数和默认流程，或为本领域的常规参数或流程。
实施例1
本实施例以污水处理厂的耐药组研究为例示出了本申请提出的耐药组整体性分析流程。
1.1 样品采集及预处理
选择某城镇污水处理厂作为研究对象，采集进水(inf)、二次沉淀池出水(sc)和最终出水样品(eff)，采用0.22 μm混合纤维素酯滤膜 (Millipore, 德国) 过滤截留水样中的DNA，每个样点获得3~4张滤膜后存放到-20℃冰箱直至处理。
1.2 DNA提取与宏基因组测序
将0.22 μm滤膜剪碎为2~3mm尺寸的方块，并基于这些方块，采用PowerSoil DNAIsolation Kit (Mo Bio, 美国) 试剂盒进行DNA提取，具体操作参见试剂盒的操作说明书。检测合格的DNA 样品用于文库构建，对于质检合格的文库将采用Illumina HiSeq 2500高通量测序平台进行PE150测序，获得原始测序序列 (raw reads)。
1.3 宏基因组数据分析
1.3.1 宏基因组序列质控、组装和分箱
（1）首先使用FastQC (v0.11.8)软件检查raw reads的整体质量，使用KneadData流程 (https://bitbucket.org/biobakery/kneaddata) 进行序列的质量控制，包括Trimmomactic (v0.39)软件去除序列的双端接头序列和低质量序列 (参数设置：SLIDINGWINDOW:4:20 MINLEN:50)，Bowtie2 (v2.3.5)软件去除人类宿主污染，得到质控过滤后的无接头高质量reads，即clean reads。
（2）使用软件MEGAHIT (v1.1.3)对clean reads进行de novo拼接（默认参数设置），获得重叠群序列 (contigs)，筛选拼接长度在500 bp以上的contigs进行后续分析。
（3）采用软件Prodigal (v2.6.3)对contigs进行ORF预测 (参数设置：-p meta)，利用软件CD-Hit(v4.8.1)对获得的ORF进行聚类 (参数设置：-aS 0.9 -c 0.95)，获得非冗余基因集；使用软件Salmon (v0.13.1)计算出每个基因的reads count，并进一步标准化，获得TPM(transcripts per million)值，用以表征基因在各个样品中的丰度信息。
（4）基于质控后的clean reads和组装后的contigs，利用软件metaWRAP (v1.2.1)提取、纯化各个样品中的微生物基因组草图 (MAGs)，并预测MAGs的基因序列，设置MAGs的筛选阈值为污染度≤10%、完整度≥50%，使用软件dRep (v2.6.2)对获得的MAGs进行去冗余(参数设置：-sa 0.95, -nc 0.30)。
1.3.2 不同风险指示基因的定性和定量识别
（a）基于clean reads：使用ARGs-OAP (v2.0) 流程（https://github.com/biofuture/Ublastx_stageone）获得样品中不同基因的组成。第一步是使用Ublast工具(v11.0.667) 从clean reads中预先筛选目标基因相似序列，第二步是使用Blastx工具(v2.2.28+) 按照默认的参数设置 (e-value≤10-7，序列相似度≥ 80%，比对长度≥25 aa)将上述相似序列与相应功能数据库进行比对，并进行基因的注释和分类，其中所使用的功能数据库包括SARG（v2.2，https://github.com/biofuture/Ublastx_stageone）、MGEdatabase (https://github.com/KatariinaParnanen/MobileGeneticElementDatabase)、VFDB （VFDB_setA_prot.fas，http://www.mgc.ac.cn/VFs/main.htm）和BacMet（v2.0,http://bacmet.biomedicine.gu.se/），分别用来注释ARG、MGE、VFG和BMG。根据下列公式，对目标基因的丰度进行标准化，由此可在不同样品或不同基因类型之间进行自由平行比较。
（1）
（2）
（3）
式中：Abundance16S表示用16S rRNA基因序列进行标准化的基因丰度，Abundancecell表示用细胞数标准化的基因丰度，Cnumber表示每个样品的细胞数量；Ni(mapped reads)是宏基因组clean reads中比对到某条特定基因参考序列的序列条数；Lreads是高通量测序获得的序列长度；Li(reference reads)是某条特定基因参考序列的长度；N16S是宏基因组clean reads中被鉴定为16S rRNA基因的序列条数；L16S是所有16S rRNA基因序列的平均长度；N16S copy number是微生物细胞中16S rRNA基因序列的平均拷贝数；n 是样品中被注释到属于某个基因的参考序列个数。
（b）基于contigs：利用Blastp工具 (v2.6.0) 将非冗余基因集的蛋白序列分别与不同的功能基因数据库进行比对，获得ARG、MGE、BMG和VFG等基因的多样性和丰度信息。
所使用的功能数据库与上方（a）中的数据库相同。
（c）基于基因组水平：对1.3.1的步骤（4）中获得的非冗余高质量MAGs的ORF进行预测，获得基因的核酸序列及其翻译后的蛋白序列，将蛋白序列与抗生素抗性基因数据库进行比对，识别出含有ARG序列的MAGs（ARG-carrying MAGs，ACMs），并统计ACMs携带的ARG种类及数量。
同理，使用相同流程，分别识别出含有MGE、BMG和VFG序列的MAGs，并统计其分别携带的MGE、BMG和VFG种类及数量。
结果：
下方表1汇总了本实施例中基于不同水平的针对废水样品中ARG、MGE、BMG和VFG的检出数量。如图3所示，对两个污水处理***不同废水样品的宏基因组的clean reads进行基因注释。图3（a）展示了ARGs在不同样品中的相对丰度和组成，结果表明进水样品中ARG的检出数量和相对丰度更高，污水处理过程对ARG的多样性有一定的削减作用。其中，多重耐药类 (multidrug)、磺胺类 (sulfonamide)、氨基糖苷类 (aminoglycoside)、大环内酯-林可酰胺-链霉素类 (macrolide-lincosamide-streptogramin, MLS)、β-内酰胺类 (beta-lactam) 和四环素类 (tetracycline) 抗性基因是样品中主要的ARG类型。图3（b）展示了BMGs在不同样品中的相对丰度和组成，LH2的最终出水样品中BMG的相对丰度显著高于其他样品，其中金属Hg的抗性基因占比最高。图3（c）和（d）分别展示了不同样品中MGE和VFG的相对丰度和组成，LH1_sc和LH2_eff的相对丰度高于其他样品，转座酶基因（transposase）在各个样品中占据主导地位，主要的VFG类别包括粘附力（adherence）、运动性（motility）和免疫调节（immune modulation）等。
图4展示了基于contigs水平识别的不同样品中ARG、BMG、MGE和VFG等基因的相对丰度和组成，从基因类型来看，不同样品的相对丰度和组成结果与基于reads的分析较为一致，说明本申请实施例的分析方法能够在各个水平上产出稳定且一致性高的结果。图5则展示了基于MAGs水平识别到不同基因的数量在各个样品中的分布，由于筛选条件更为严格，可以检出到的基因数量相比于contigs水平的结果少，但是这样识别到的基因更可能在环境样品中广泛存在，并通过某些特定的微生物发挥相应功能。
表1 三种分析水平识别到的不同基因的数量
由此，可以看出本申请实施例中提出的基于宏基因组的环境耐药组分析方法在多个层面（reads、contigs及MAGs）上均能产出较高数量的耐药相关基因，说明本申请实施例的方法能够灵敏且全面地捕获环境耐药组信息；同时，由图3-5可以看出，各个层面所产出的耐药相关基因的分布相似，说明本申请实施例提出的方法能够在多个层面上提供较为稳定且准确的耐药组分析结果，从而通过综合多个层面的耐药组信息，能够合理高效的进行***性的环境耐药组风险评估。
实施例2
本实施例基于实施例1中的数据进一步进行了多水平的ARGs的宿主信息识别，具体流程包括：
（1）基于clean reads：首先使用Ublast工具 (v11.0.667) 和Blastx工具(v2.2.28+)将clean reads与SARG (v2.2) 进行比对，从中筛选出ARG相似序列 (ARG-likereads, ALRs) (e-value≤ 10-7，序列相似度≥ 80%，比对长度≥25 aa)；之后使用软件Kraken2 (v2.0.8-beta) 和GTDB (r89) 数据库对提取到ALRs进行物种信息的匹配，统计ALRs在各样本中的分布情况，得到ARB群落在不同样本中的组成信息。
（2）基于contigs：首先利用Blastp工具 (v2.6.0)和SARG (v2.2)数据库筛选出携带ARG的contigs (ACCs)，然后使用软件Kraken2 (v2.0.8-beta) 和GTDB (r89) 数据库对所有ACCs进行物种信息的匹配，利用软件CoverM (v0.6.1, https://github.com/wwood/CoverM) 计算每条ACCs的相对丰度，用比对到ACCs不同位点的平均序列数 (mean) 来表示。
（3）基于MAGs：利用BLASTP工具 (v2.6.0) 将MAGs蛋白序列与SARG (v2.2) 数据库进行比对 (e-value ≤ 10-5, identity ≥ 80% , coverage ≥ 70%)，对MAGs中的ARGs进行注释，获得携带ARGs的MAGs（ACMs）；之后使用软件GTDB-Tk (v1.0.2) 和GTDB (r89)数据库对ACMs进行物种注释，软件metaWRAP (v1.2.1) 的quant_bins模块给出了每个ACMs的相对丰度信息，表示为每百万序列中基因组拷贝数 (genome copies/ppm reads)。
结果：
图6、图7和图8分别展示了ALR、ACC和ACM三种策略识别到的不同废水样品中ARG宿主组成，结果表明Gammaproteobacteria是最主要的ARG宿主，呈现出较高一致性。基于ALR策略能够识别到的ARG宿主多样性更高，包括低丰度物种，同时可以建立宿主丰度与ARGs丰度的直接关联，具有更简洁的分析流程，可显著节约计算资源和时间。而基于ACCs和ACMs的策略可以深入分析ARGs侧翼序列的功能信息，识别ARG宿主的基因组特征，为表征其环境风险并研究靶向性管控措施提供参考。
实施例3
本实施例基于实施例1和2中的数据进一步进行了多水平的ARGs的水平基因转移预测，具体流程包括：
（1）基于contigs：利用软件Plasflow (v1.1) 按照默认参数 (-c 0.7) 对contigs (根据软件开发人员的建议，将序列长度小于1000bp的contigs排除在外) 进行分型，识别出染色体序列和质粒序列，并分析两类序列所携带的ARGs组成。利用软件WAAFLE筛选不同宏基因组数据集中发生水平（侧向）基因转移事件的contigs (>500bp)。
（2）基于MAGs：在实施例2步骤（3）分析的基础上，区分MAGs数据集中的ACMs和非ACMs（nACMs），使用软件MetaCHIP提取ACMs之间、或者ACMs与nACMs之间潜在的水平基因转移事件，并可以预测基因转移的方向，MAGs核酸序列及其物种注释信息作为输入文件。最终获得“ACMs- ACMs”或“ACMs-nACMs”的配对关系，构建有向基因转移网络，并通过软件Gephi(v0.9.2) 进行可视化。
结果：
基于序列分型的结果（图9），样品中编码在质粒序列上的ARGs相对丰度高于其他两类序列，说明废水样品中的大部分ARGs普遍具有较高的移动性，其中氨基糖苷类、四环素类、磺胺类、MLS类和β-内酰胺类抗性基因更多地编码在质粒序列上。利用WAAFLE软件识别了不同样品contigs数据集中发生HGT事件的contigs数量，如图10所示，结果表明发生HGT的contigs数量均在进水样品中最高，且随废水处理过程逐级下降。如图11所示，基于ACMs的识别，针对两个污水处理***分别构建了水平基因转移网络模型，表2为与网络相关的特征参数数据。图中每个节点表示一个MAG，其中红色代表ACM，蓝色代表非ACM，节点大小表示携带ARG的数量，每条边表示两个节点之间可能发生HGT事件，边的粗细表示发生HGT的数量。通过比较两个网络结构及相关参数，LH2中ARGs发生HGT的可能性和频率更高。
表2 水平转移网络特征
实施例4
本实施例以我国杭州湾沿岸4个独立的污水处理***和杭州湾海域作为研究对象，其中污水处理***包括生活污水处理***（SD）、工业废水处理***（SI）和两个混合污水处理***（LH1和LH2），采集进水（inf）、二次沉淀池出水（sc）和最终出水样品（eff），在杭州湾海域布设10个采样点，采集表层沉积物样品（HB1-HB10）。
4.1 样品预处理
同实施例1的1.1步骤。
4.2 样品DNA提取与宏基因组测序
同实施例1的1.2步骤。
4.3宏基因组数据分析
（1）基于KneadData流程 (https://bitbucket.org/biobakery/kneaddata) 实现序列的质量控制，包括利用软件Trimmomactic (v0.39)去除序列的双端接头序列和低质量序列 (参数设置：SLIDINGWINDOW:4:20 MINLEN:50)，通过软件Bowtie2 (v2.3.5)去除宿主污染。
（2）采用软件MEGAHIT (v1.1.3)对clean reads进行de novo拼接 (参数设置：k-min 27 k-max 191 step 20)，获得重叠群序列 (contigs)，筛选拼接长度在500 bp以上的contigs进行后续分析。
（3）采用软件Prodigal (v2.6.3)对contigs进行ORF (open reading frame) 预测 (参数设置：-p meta)。利用BLASTP工具 (v2.6.0) 将ORF与SARG（v2.2，https://github.com/biofuture/Ublastx_stageone）进行比对，识别携带ARGs的contigs（ACCs）
（4）利用BLASTP工具 (v2.6.0)将ACCs的ORF分别与ACLAME (v0.4，http://aclame.ulb.ac.be/)、BacMet (v2.0, http://bacmet.biomedicine.gu.se/)和VFDB(VFDB_setA_prot.fas，http://www.mgc.ac.cn/VFs/main.htm) 等数据库进行比对，获得MGEs、BMGs和VFGs等基因在ACCs上的出现情况。
4.4耐药组风险指数计算
4.4.1 将每个样品的所有ACCs分成8种类型，根据含有目标基因的类型数对不同的ACCs数量进行赋权，其中ARG、MGE和VFG表征了环境耐药组的直接风险，权重为1，而BMG则表征了环境耐药组的间接风险，将其权重分别设置为0、0.5和1并按照下方式（1）计算综合风险指数。
（1）
4.4.2 结果
图12为BMG权重为0、0.5和1时计算所得的耐药组风险指数。如图12所示，污水处理***样品的风险水平普遍高于杭州湾海域的沉积物样品，特别是工业废水处理***的最终排水依然表现出较高的耐药组风险，这与先验知识相符，由此证明了本实施例风险评估的可信性。其次，在出水样本中，如果完全忽略BMG的风险(权重=0)，最终风险评分将显著降低。这表明，在废水环境中，由于各种化学污染物和重金属的存在，ARGs和BMGs更容易共同发生，并且BMGs的存在也增加了ARB在外部环境压力下更稳定地增殖和扩散的能力。因此，BMGs对抗生素耐药菌综合风险的贡献不容忽视，尤其是在污染严重的环境中。将BMG的权重值设定为0.5和1时，最终风险评分差别不大，将BMG权重设定为0.5能够更加谨慎地对待其直接的风险贡献。
由此，基于BMG权重为0.5，每个样品的综合风险指数根据公式（1）-（4）计算：
（1）/>
（2）
（3）
（4）
式中：nContig表示contig数据集中contigs的总数，α表示ACCs的数量，β1表示携带MGE的ACCs数量，β2表示携带VFG的ACCs数量，β3表示携带BMG的ACCs数量，γ1表示携带MGE和VFG的ACCs数量，γ2表示携带MGE和BMG的ACCs数量，γ3表示携带BMG和VFG的ACCs数量，δ1表示同时携带MGE、VFG和BMG的ACCs数量。
对比例
本对比例使用现有的风险评价方法(metacompare)流程(https://github.com/minoh0201/MetaCompare)，对实施例4中的样品进行了风险评价，其与实施例4提出的方法所得的结果的比较见图13。
如图13所示，本申请实施例提出的风险评价方法与现有方法在废水样品中识别出一致的风险变化模式，由此证明了本申请提出的方法的可行性和准确性。进一步地，与metacompare计算的风险评分相比，本申请实施例提出的方法能够识别较高风险的采样点，如HB5和HB6。与metacompare相比，经探索及验证，本申请实施例提出的方法为识别contigs上的风险基因设定了更严格的阈值，同时增加了对BMGs等相关功能基因的考虑，实现了更为全面和准确的综合风险水平评估。
为达上述目的，本公开实施例还提出了一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述计算机程序被处理器执行时实现如上任一实施例所述的基于宏基因组的环境耐药组分析方法或基于宏基因组的环境耐药组风险评价方法。
为达上述目的，本公开实施例还提出了一种非临时性计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如上任一实施例所述的基于宏基因组的环境耐药组分析方法或基于宏基因组的环境耐药组风险评价方法。
为达上述目的，本公开实施例还提出了一种计算机程序产品，包括计算机程序，其中所述计算机程序在被处理器执行时实现如上任一实施例所述的基于宏基因组的环境耐药组分析方法或基于宏基因组的环境耐药组风险评价方法。
为达上述目的，本公开实施例还提出了一种计算机程序，其中所述计算机程序包括计算机程序代码，当所述计算机程序代码在计算机上运行时，使得计算机执行如上任一实施例所述的基于宏基因组的环境耐药组分析方法或基于宏基因组的环境耐药组风险评价方法。
需要注意的是，前述对基于宏基因组的环境耐药组分析方法或基于宏基因组的环境耐药组风险评价方法实施例的解释说明也适用于本公开实施例的计算机设备、非临时性计算机可读存储介质、计算机程序产品和计算机程序，此处不再赘述。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (6)
1.一种基于宏基因组的环境耐药组分析方法，其特征在于，所述方法包括：
基于环境样品，获得所述环境样品的宏基因组数据；和
基于所述宏基因组数据，识别功能基因，并计算所述功能基因在所述环境样品中的相对丰度，包括：
a.基于所述宏基因组数据，对读长序列进行比对和注释以识别所述功能基因，并计算所述功能基因在所述环境样品中的第一丰度；
b.对所述读长序列进行组装以获得重叠群序列，并基于所述重叠群序列进行比对和注释以识别所述功能基因，并计算所述功能基因在所述环境样品中的第二丰度；和
c.对所述重叠群序列进行分箱以获得微生物基因组草图，并基于所述微生物基因组草图进行比对和注释以识别所述功能基因，并计算所述功能基因在所述环境样品中的第三丰度，
其中所述功能基因包括以下中的一种或多种：抗生素抗性基因、可移动遗传元件、毒力因子基因和杀菌剂和金属抗性基因，
其中步骤b中，基于所述重叠群序列进行比对和注释以识别所述功能基因具体包括：
预测所述重叠群序列的开放阅读框，以获得所述重叠群序列对应基因的核酸序列及其翻译后的蛋白序列；
基于所述重叠群序列对应基因的核酸序列及其翻译后的蛋白序列的序列相似度对所述重叠群序列对应基因进行聚类，计算所述重叠群序列对应基因的基因丰度并获得非冗余基因集；和
基于功能基因数据库，对所述非冗余基因集进行比对和注释，以识别所述功能基因，并获得所述功能基因的组成信息和丰度信息，其中所述比对中，所述非冗余基因集中e-value≤10-5、序列相似性≥80％、覆盖度≥70％的基因序列被识别为所述功能基因；并且
其中步骤c中，所述基于所述微生物基因组草图进行比对和注释以识别所述功能基因还包括：
对所述微生物基因组草图进行去冗余和质量控制，以获得非冗余高质量的微生物基因组草图；
预测所述非冗余高质量的微生物基因组草图的开放阅读框，以获得所述非冗余高质量的微生物基因组草图中各对应基因的核酸序列及其翻译后的蛋白序列；和
基于功能基因数据库，对所述各对应基因的翻译后的蛋白序列进行比对和注释，以识别所述功能基因，并获得所述功能基因的组成信息和丰度信息，
其中，基于完整度和污染率对所述微生物基因组草图进行所述质量控制，其中要求所述完整度≥50％并且所述污染率≤10％，
其中所述功能基因在所述环境样品中的第一丰度的计算基于16S rRNA基因序列数和细胞数量，计算公式为：
其中Abundance16S表示用16S rRNA基因序列进行标准化的功能基因丰度，Abundancecell表示用细胞数标准化的功能基因丰度，Cnumber表示每个环境样品的细胞数量；Ni(mapped reads)是宏基因组reads中比对到某条特定基因参考序列的序列条数；Lreads是高通量测序获得的序列长度；Li(reference reads)是某条特定基因参考序列的长度；N16S是宏基因组reads中被鉴定为16S rRNA基因的序列条数；L16S是所有16S rRNA基因序列的平均长度；N16S copy number是微生物细胞中16S rRNA基因序列的平均拷贝数；n是所述环境样品中被注释到属于某个基因的参考序列个数。
2.根据权利要求1所述的方法，其特征在于，所述方法还包括：
基于所述功能基因进行宿主分析，包括：
i.识别所述读长序列中的功能基因相似序列，对所述功能基因相似序列进行物种注释，以获得功能基因相似序列在各物种的序列数作为第一宿主丰度，其中所述功能基因相似序列与所述功能基因的序列相似度为≥80％，比对长度≥25aa，并且e值≤10-7；
ii.识别所述重叠群序列中含有功能基因序列的重叠群序列，对所述含有功能基因序列的重叠群序列进行物种注释，以获得所述含有功能基因序列的重叠群序列在各物种的序列数作为第二宿主丰度；和
iii.识别所述微生物基因组草图中含有功能基因序列的微生物基因组草图，对所述含有功能基因序列的微生物基因组草图进行物种注释，以获得所述功能基因序列的微生物基因组草图在各物种的序列数作为第三宿主丰度。
3.根据权利要求2所述的方法，其特征在于，所述方法还包括：
基于所述功能基因进行侧向基因转移潜力分析，包括：
x.识别所述重叠群序列中的染色体来源序列和质粒来源序列，并获得所述染色体来源序列和质粒来源序列所携带的所述功能基因的组成信息；
y.识别所述重叠群序列中可能发生侧向基因转移的重叠群序列，并获得所述可能发生侧向基因转移的重叠群序列所携带的所述功能基因的组成信息；和
z.基于所述微生物基因组草图和所述含有功能基因序列的微生物基因组草图，预测所述含有功能基因序列的微生物基因组草图之间和所述含有功能基因序列的微生物基因组草图与不含有所述功能基因序列的微生物基因组草图之间的侧向基因转移过程，构建“供体-受体”配对关系的侧向基因转移网络。
4.根据权利要求2所述的方法，其特征在于，步骤ii还包括：利用软件CoverM计算所述第二宿主丰度，其中取所述含有功能基因序列的重叠群序列在各物种的序列数的平均序列数为所述第二宿主丰度；并且
步骤iii还包括：使用1.2.1版本的metaWRAP的quant_bins模块计算所述第三宿主丰度，其中取每百万序列中基因组拷贝数表示所述第三宿主丰度。
5.根据权利要求3所述的方法，其特征在于，
步骤x还包括：对所述重叠群序列进行过滤，以去除序列长度短于1000bp的所述重叠群序列；和
使用1.1版本的Plasflow软件对过滤后的重叠群序列进行分型，以识别出所述过滤后的重叠群序列中的染色体来源序列和质粒来源序列，并且
所述步骤y还包括：使用WAAFLE软件识别所述重叠群序列中可能发生侧向基因转移的重叠群序列，其中基于序列长度大于500bp使用WAAFLE软件筛选出不同宏基因组数据集中发生侧向基因转移的所述重叠群序列。
6.一种基于宏基因组的环境耐药组风险评价方法，其特征在于，所述方法包括：
根据权利要求1至5中任一项所述的基于宏基因组的环境耐药组分析方法获得环境样品的功能基因分别在读长序列、重叠群序列和微生物基因组草图水平上的组成信息和丰度信息，其中所述功能基因包括以下中的一种或多种：抗生素抗性基因、可移动遗传元件、毒力因子基因、和杀菌剂和金属抗性基因；
将组装获得的所有重叠群序列的开放阅读框与功能基因数据库进行比对，以获得第一重叠群序列，其中所述第一重叠群序列含有所述抗生素抗性基因、可移动遗传元件、毒力因子基因和杀菌剂和金属抗性基因中的任一种功能基因；
将所述第一重叠群序列的开放阅读框与所述功能基因数据库进行比对，以获得第二重叠群序列，其中所述第二重叠群序列含有其它三种功能基因中的至少一种；
统计所述第一重叠群序列和所述第二重叠群序列中含有的功能基因组合类型及其对应的重叠群序列数量，并对所述重叠群序列数量赋予权重；和
基于所述权重，计算得到所述环境样品的环境耐药组风险指数，其中，对所述权重进行加和以得到所述环境样品的环境耐药组风险指数，
其中所述功能基因数据库选自由以下中的一个或多个：SARG、MGE database、ACLAME、VFDB和BacMet，
其中，所述第一重叠群序列和所述第二重叠群序列中含有的功能基因组合类型及其所赋权重如下表所示：
所述风险指数的计算公式为：
其中：
β＝2β1+2β2+1.5β3 (2)
γ＝3γ1+2.5γ2+2.5γ3 (3)
δ＝3.5δ1 (4)
nContig为所述组装获得的所有重叠群序列的总数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311389711.8A CN117174165B (zh) | 2023-10-25 | 2023-10-25 | 基于宏基因组的环境耐药组分析方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311389711.8A CN117174165B (zh) | 2023-10-25 | 2023-10-25 | 基于宏基因组的环境耐药组分析方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117174165A CN117174165A (zh) | 2023-12-05 |
| CN117174165B true CN117174165B (zh) | 2024-03-12 |
Family
ID=88945246
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311389711.8A Active CN117174165B (zh) | 2023-10-25 | 2023-10-25 | 基于宏基因组的环境耐药组分析方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117174165B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118230820A (zh) * | 2024-03-19 | 2024-06-21 | 浙江洛兮医学检验实验室有限公司 | 基于宏基因测序数据的耐药基因物种来源鉴定方法 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108334750A (zh) * | 2018-04-19 | 2018-07-27 | 江苏先声医学诊断有限公司 | 一种宏基因组数据分析方法及*** |
| CN110819699A (zh) * | 2019-11-27 | 2020-02-21 | 辽宁大学 | 一种水环境中人类粪便指示物的定量检测方法 |
| CN111944914A (zh) * | 2020-07-16 | 2020-11-17 | 中国科学院生态环境研究中心 | 一种基于抗性基因及毒力因子基因评价水体健康风险的方法 |
| CN113689912A (zh) * | 2020-12-14 | 2021-11-23 | 广东美格基因科技有限公司 | 基于宏基因组测序的微生物对比结果校正的方法和*** |
| KR20220069626A (ko) * | 2020-11-20 | 2022-05-27 | 대한민국(농촌진흥청장) | 가금류 맹장내 마이크로바이옴의 메타게놈 자원을 이용한 가금류 사육환경의 정보 제공방법 |
| CN114822697A (zh) * | 2022-03-28 | 2022-07-29 | 南京农业大学 | 一种利用宏基因组分析溯源土壤耐药基因污染的方法 |
-
2023
- 2023-10-25 CN CN202311389711.8A patent/CN117174165B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108334750A (zh) * | 2018-04-19 | 2018-07-27 | 江苏先声医学诊断有限公司 | 一种宏基因组数据分析方法及*** |
| CN110819699A (zh) * | 2019-11-27 | 2020-02-21 | 辽宁大学 | 一种水环境中人类粪便指示物的定量检测方法 |
| CN111944914A (zh) * | 2020-07-16 | 2020-11-17 | 中国科学院生态环境研究中心 | 一种基于抗性基因及毒力因子基因评价水体健康风险的方法 |
| KR20220069626A (ko) * | 2020-11-20 | 2022-05-27 | 대한민국(농촌진흥청장) | 가금류 맹장내 마이크로바이옴의 메타게놈 자원을 이용한 가금류 사육환경의 정보 제공방법 |
| CN113689912A (zh) * | 2020-12-14 | 2021-11-23 | 广东美格基因科技有限公司 | 基于宏基因组测序的微生物对比结果校正的方法和*** |
| CN114822697A (zh) * | 2022-03-28 | 2022-07-29 | 南京农业大学 | 一种利用宏基因组分析溯源土壤耐药基因污染的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117174165A (zh) | 2023-12-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gruber-Vodicka et al. | phyloFlash: rapid small-subunit rRNA profiling and targeted assembly from metagenomes | |
| Ren et al. | Identifying viruses from metagenomic data using deep learning | |
| Yue et al. | Evaluating metagenomics tools for genome binning with real metagenomic datasets and CAMI datasets | |
| CN110349629B (zh) | 一种利用宏基因组或宏转录组检测微生物的分析方法 | |
| Nayfach et al. | Toward accurate and quantitative comparative metagenomics | |
| Miller et al. | Autometa: automated extraction of microbial genomes from individual shotgun metagenomes | |
| Bickhart et al. | Assignment of virus and antimicrobial resistance genes to microbial hosts in a complex microbial community by combined long-read assembly and proximity ligation | |
| Freitas et al. | Accurate read-based metagenome characterization using a hierarchical suite of unique signatures | |
| Zielińska et al. | The choice of the DNA extraction method may influence the outcome of the soil microbial community structure analysis | |
| Nishimura et al. | The OceanDNA MAG catalog contains over 50,000 prokaryotic genomes originated from various marine environments | |
| CN117174165B (zh) | 基于宏基因组的环境耐药组分析方法 | |
| Louvel et al. | meta BIT, an integrative and automated metagenomic pipeline for analysing microbial profiles from high‐throughput sequencing shotgun data | |
| Saheb Kashaf et al. | Recovering prokaryotic genomes from host-associated, short-read shotgun metagenomic sequencing data | |
| US20150376697A1 (en) | Method and system to determine biomarkers related to abnormal condition | |
| Gostel et al. | Microfluidic Enrichment Barcoding (MEBarcoding): A new method for high throughput plant DNA barcoding | |
| Bidovec-Stojkovič et al. | Prospective genotyping of Mycobacterium tuberculosis from fresh clinical samples | |
| Haryono et al. | Recovery of high quality metagenome-assembled genomes from full-scale activated sludge microbial communities in a tropical climate using longitudinal metagenome sampling | |
| Jurburg et al. | The community ecology perspective of omics data | |
| Basset et al. | Comparison of traditional and DNA metabarcoding samples for monitoring tropical soil arthropods (Formicidae, Collembola and Isoptera) | |
| Hickl et al. | binny: an automated binning algorithm to recover high-quality genomes from complex metagenomic datasets | |
| Chang et al. | Circulating cell-free RNA in blood as a host response biomarker for detection of tuberculosis | |
| Kumar et al. | A comprehensive overview of microbiome data in the light of machine learning applications: categorization, accessibility, and future directions | |
| CN115938491B (zh) | 一种用于临床病原诊断的高质量细菌基因组数据库构建方法及*** | |
| Xi et al. | SiftCell: A robust framework to detect and isolate cell-containing droplets from single-cell RNA sequence reads | |
| Fremin et al. | Comparative genomics identifies thousands of candidate structured RNAs in human microbiomes |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |