百度360必应搜狗淘宝本站头条
0xalinux查看python版本信息python怎么改成中文python安装opencv库1到10的阶乘之和是多少
当前位置：网站首页 > 技术文章 > 正文

宏基因组菌群注释之MetaPhlAn3（细菌宏基因组测序）

itomcoil 2025-03-14 18:06 30 浏览

微信搜一搜：百合木瓜，更多生信分析工具分享。

MetaPhlAn3用于分析宏基因组鸟枪法测序数据（非16S）中微生物群落组成的软件，基于物种marker gene进行菌群分类，降低假阳性，但只能检出数据库内的物种，对于测序量较少的宏基因组样本可能会检测不到部分真实存在的物种。
例如：使用MetaPhlAn3 对E. coli基因组原始测序数据进行分析，准确分类到E. coli，该基因组数据约2.6G，共包含14562039条序列，双端共比对86863条序列，有效reads利用率0.59%。

MetaPhlAn3使用bowtie2将宏基因组样本中的原始reads比对到物种marker数据，支持单端或双端的fastq测序数据。MetaPhlAn3估计每个marker的覆盖率，并计算同一进化枝的marker覆盖率的稳健平均值作为进化枝覆盖率；然后将将进化枝覆盖范围在所有检测到的进化枝上进行归一化，以获得每个分类单元的相对丰度。

MetaPhlAn3优化：

  • 参考数据库扩增，marker genes是由~100,000参考基因组（~99，500细菌、~500真核生物）确定，约1.1M；
  • 进化枝稳健平均值参数('stat_q')，不包含marker丰度最高和最低部分，默认0.2(i.e. excludes the 20% of markers with the highest abundance as well as the 20% of markers with the lowest abundance)；
  • 质量控制：丢弃低质量序列(短于70bp）和比对结果（MAPQ值<5)；
  • 增加当前数据库中不存在的分类群相应的'unknown'部分，通过从总reads数中减去每个分类单元的平均reads深度（taxon-specific average genome length）来计算的；
  • 支持CAMI输出格式；

（B）MetaPhlAn 3 was applied to a set of 113 total evaluation datasets provided by CAMIrepresenting diverse human-associated microbiomes and five datasets of non-human-associated microbiomes . MetaPhlAn 3 showed increased performance compared with the previous version MetaPhlAn 2, mOTUs2, and Bracken 2.5.(C) MetaPhlAn 3 better recapitulates relative abundance profiles both from human and murine gastrointestinal metagenomes as well from non-human-associated communities compared to the other currently available tools . Bracken is reported both using its original estimates based on the fraction of reads assigned to each taxon and after re-normalizing them using the genome lengths of the taxa in the gold standard to match the taxa abundance estimate of the other tools.

安装

conda
conda install -c bioconda python=3.7 metaphlan

docker
docker pull biobakery/metaphlan

Pypi
pip install metaphlan

数据库

Metaphlan 3是基于ChocoPhlAn3数据库的物种marker序列进行菌群注释，最新版数据库是基于17,000 reference genomes (~13,500 bacterial and archaeal, ~3,500 viral, and ~110 eukaryotic)筛选而来。ChocoPhlAn3数据库marker序列筛选过程请参考：宏基因组数据库ChocoPhlAn3的巧妙设计

#使用Pypi方式安装metaphlan软件，未安装bowtie2。bowtie2无需安装，下载解压后即可运行，下载地址：https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.5

#ChocoPhlAn3数据库下载及建库
metaphlan --install

#添加自己数据，或使用自己数据库。
#可按mpa_v30_CHOCOPhlAn_201901.fna和mpa_v30_CHOCOPhlAn_201901.pkl 件格式整理自己数据；
#ChocoPhlAn3数据下载链接，5分钟左右即可下载完成：http://cmprod1.cibio.unitn.it/biobakery3/metaphlan_databases/mpa_v30_CHOCOPhlAn_201901.tar

#序列文件
>1000373__GeneID:11604944
ATGGATTCCACAGACAACATAGAGCACGACACATACAGAGAGGACATCCGGT
>100053__V6HSA9__LEP1GSC062_4244 UniRef90_V6HSA9;k__Bacteria|p__Spirochaetes|c__Spirochaetia|o__Spirochaetia_unclassified|f__Leptospiraceae|g__Leptospira|s__Leptospira_alexanderi;GCA_000243815
GACGCGGACTACGGCAAAGACGCTCGTAGAGTTGAACCACAT

#bowtie2建库
bowtie2-build --quiet --threads 4 -f mpa_v30_CHOCOPhlAn_201901.fna mpa_v30_CHOCOPhlAn_201901

#注释文件pkl
#使用python的pickle包查看，插入注释信息

import pickle
import bz2
db=pickle.load(bz2.open('mpa_v30_CHOCOPhlAn_201901.pkl','r'))

#db.keys()
#dict_keys(['markers', 'taxonomy', 'merged_taxon']) 
#共包含三个部分
#list(db['markers'].items())[0]
#('244590__GeneID:2658371', {'ext': [], 'score': 0.0, 'clade':'s__Sulfolobus_spindle_shaped_virus_2', 'len': 147, 'taxon':'k__Viruses|p__Viruses_unclassified|c__Viruses_unclassified|o__Viruses_unclassified|f__Fuselloviridae|g__Alphafusellovirus|s__Sulfolobus_spindle_shaped_virus_2'})
#list(db['taxonomy'].items())[0]
#('k__Viruses|p__Negarnaviricota|c__Monjiviricetes|o__Mononegavirales|f__Rhabdoviridae|g__Nucleorhabdovirus|s__Maize_mosaic_nucleorhabdovirus|t__PRJNA14920',('10239|2497569|2497574|11157|11270|11306|279896', 12133))
#list(db['merged_taxon'].items())[0]
#(('k__Archaea|p__Euryarchaeota|c__Euryarchaeota_unclassified|o__Euryarchaeota_unclassified|f__Euryarchaeota_unclassified|g__Euryarchaeota_unclassified|s__candidate_divison_MSBL1_archaeon_SCGC_AAA833K04', '2157|28890|||||1698258'), [('k__Archaea|p__Euryarchaeota|c__Euryarchaeota_unclassified|o__Euryarchaeota_unclassified|f__Euryarchaeota_unclassified|g__Euryarchaeota_unclassified|s__candidate_division_MSBL1_archaeon_SCGC_AAA833F18', '2157|28890|||||1698257', 1),('k__Archaea|p__Euryarchaeota|c__Euryarchaeota_unclassified|o__Euryarchaeota_unclassified|f__Euryarchaeota_unclassified|g__Euryarchaeota_unclassified|s__candidate_division_MSBL1_archaeon_SCGC_AAA833K04', '2157|28890|||||1698258', 1)])

#Add the information of the new marker as the other markers
db['markers'][new_marker_name] = {
                                   'clade': the clade that the marker belongs to,
                                   'ext': {the GCA of the first external genome where the marker appears,
                                           the GCA of the second external genome where the marker appears,
                                          },
                                   'len': length of the marker,
                                   'taxon': the taxon of the marker
                                }

#Add the taxonomy of the new genomes
db['taxonomy']['7-levels taxonomy with clade names of genome1'] = ('7-levels NCBI taxonomy id of genome1', length of genome1)
db['taxonomy']['7-levels taxonomy with clade names of genome2'] = ('7-levels NCBI taxonomy id of genome1', length of genome2)

#保存新pkl文件
with bz2.BZ2File('mpa_v30_CHOCOPhlAn_New.pkl','w')as ofile:
   pickle.dump(db,ofile,pickle.HIGHEST_PROTOCOL)

使用

#单端数据
metaphlan metagenome.fastq --input_type fastq --bowtie2out metagenome.bowtie2.bz2 --nproc 5 -o profiled_metagenome.txt

#双端数据（并不使用双端信息）
metaphlan metagenome_1.fastq,metagenome_2.fastq --bowtie2out metagenome.bowtie2.bz2 --nproc 5 --input_type fastq -o profiled_metagenome.txt

#先bowtie2进行数据库比对，然后使用metaphlan进行菌群注释
bowtie2 --sam-no-hd --sam-no-sq --no-unal --very-sensitive -S metagenome.sam -x metaphlan_databases/mpa_v30_CHOCOPhlAn_201901 -U metagenome.fastq
metaphlan metagenome.sam --input_type sam -o profiled_metagenome.txt

#估计宏基因组中未知微生物（数据库外的物种）
metaphlan metagenome.fastq --bowtie2out metagenome.bowtie2.bz2 -nproc 5 --input_type fastq --unknown_estimation -o profiled_metagenome.txt

#合并多个宏基因组样本菌群注释结果,仅当使用相同版本的 MetaPhlAn 数据库执行分析时，才能合并输出文件。
merge_metaphlan_tables.py metaphlan_output1.txt metaphlan_output2.txt metaphlan_output3.txt > output/merged_abundance_table.txt
merge_metaphlan_tables.py metaphlan_output*.txt > output/merged_abundance_table.txt

#计算UniFrac距离矩阵，相关脚本在/usr/local/lib/python3.8/dist-packages/metaphlan/utils
#安装R
#sudo apt install r-base-core
#安装R包
#install.packages('vegan')
#install.packages('rbiom')
#install.packages('ape')
#加权
Rscript calculate_unifrac.R merged_abundance_table.txt mpa_v30_CHOCOPhlAn_201901_species_tree.nwk unifrac_merged_mpa3_profiles.tsv
#非加权
Rscript calculate_unifrac.R merged_abundance_table.txt mpa_v30_CHOCOPhlAn_201901_species_tree.nwk unifrac_merged_mpa3_profiles.tsv unweighted

参数

usage: metaphlan --input_type {fastq,fasta,bowtie2out,sam} [--force]
                 [--bowtie2db METAPHLAN_BOWTIE2_DB] [-x INDEX]
                 [--bt2_ps BowTie2 presets] [--bowtie2_exe BOWTIE2_EXE]
                 [--bowtie2_build BOWTIE2_BUILD] [--bowtie2out FILE_NAME]
                 [--min_mapq_val MIN_MAPQ_VAL] [--no_map] [-tmp_dir]
                 [--tax_lev TAXONOMIC_LEVEL] [--min_cu_len]
                 [--min_alignment_len] [--add_viruses] [-ignore_eukaryotes]
                 [--ignore_bacteria] [--ignore_archaea] [--stat_q]
                 [--perc_nonzero] [--ignore_markers IGNORE_MARKERS]
                 [--avoid_disqm] [--stat] [-t ANALYSIS TYPE]
                 [--nreads NUMBER_OF_READS] [--pres_th PRESENCE_THRESHOLD]
                 [--clade] [--min_ab] [-o output file] [-sample_id_key name]
                 [--use_group_representative] [--sample_id value]
                 [-s sam_output_file] [--legacy-output] [--CAMI_format_output]
                 [--unknown_estimation] [--biom biom_output] [--mdelim mdelim]
                 [--nproc N] [--install] [--force_download]
                 [--read_min_len READ_MIN_LEN] [-v] [-h]
                 [INPUT_FILE] [OUTPUT_FILE]


DESCRIPTION
 MetaPhlAn version 3.0.14 (19 Jan 2022): 
 METAgenomic PHyLogenetic ANalysis for metagenomic taxonomic profiling.


INPUT_FILE：输入文件；
--input_type TYPE：指定输入文件类型，包括fastq,fasta,bowtie2out,sam；
OUTPUT_FILE：输出结果文件，菌群分类相对丰度，tab制表符分隔的文本文件；
--bowtie2db METAPHLAN_BOWTIE2_DB：参考数据库目录，默认/usr/local/lib/python3.8/dist-packages/metaphlan/metaphlan_databases；
-x INDEX, --index INDEX：参考数据索引；
--bt2_ps BowTie2 presets：bowtie2预设参数，值包括：sensitive、very sensitive、sensitive-local、very-sensitive-local。默认：very-sensitive；
--bowtie2_exe BOWTIE2_EXE：bowtie2安装路径；
--bowtie2out FILE_NAME：保存bowtie2结果文件，结合参数‘-t’可用于后续不同分析内容分析，文件内容例如：EAS18:1:1:1:16:530:0/1__1.743    562__P76656__CD803_05075
--min_mapq_val MIN_MAPQ_VAL：最小比对质量值（MAPQ），默认5；


Post-mapping arguments:
  --tax_lev TAXONOMIC_LEVEL  The taxonomic level for the relative abundance output，默认'a':
                        'a' : all taxonomic levels
                        'k' : kingdoms
                        'p' : phyla only
                        'c' : classes only
                        'o' : orders only
                        'f' : families only
                        'g' : genera only
                        's' : species only
--min_cu_len ：进化枝中markers最小总核苷酸长度，用于在不考虑进化枝丰度情况下估计丰度。默认：2000
--min_alignment_len：最小比对长度，比对长度低于指定阈值，丢弃。默认None；
--add_viruses         Allow the profiling of viral organisms
--ignore_eukaryotes   Do not profile eukaryotic organisms
--ignore_bacteria     Do not profile bacterial organisms
--ignore_archaea      Do not profile archeal organisms
--stat_q              Quantile value for the robust average 。默认 0.2
--perc_nonzero：错误鉴定物种markers的相对丰度，默认：0.33
--ignore_markers IGNORE_MARKERS File containing a list of markers to ignore.
--avoid_disqm: 停止消除准marker歧义的过程。 通常建议保留消歧程序以最大程度地减少误报
--stat   将marker丰度转换为进化枝的丰度的统计方法，默认：tavg_g；
                        'avg_g'  : clade global (i.e. normalizing all markers together) average
                        'avg_l'  : average of length-normalized marker counts
                        'tavg_g' : truncated clade global average at --stat_q quantile
                        'tavg_l' : truncated average of length normalized marker counts (at --stat_q)
                        'wavg_g' : winsorized clade global average (at --stat_q)
                        'wavg_l' : winsorized average of length normalized marker counts (at --stat_q)
                        'med'    : median of length-normalized marker counts

Additional analysis types and arguments:
-t ANALYSIS TYPE 要执行的分析类型，默认rel_ab：
                         * rel_ab：根据相对丰度分析宏基因组
                         * rel_ab_w_read_stats：根据相对丰度分析宏因组，并估计来自每个进化枝的reads数。
                         * reads_map：从reads到进化枝的映射（仅read命中marker）
                         * clade_profiles：具有至少一个非空marker的化枝，markers是按进化枝获得丰度；
                         * marker_ab_table：归一化marker计数（仅当 >0.0 并且指定 --nreads 则按宏基因组大小归一化）
                         * marker_counts：非标准化的marker计数[谨慎用]
                         * marker_pres_table：样本中存在的marker列（如果没有使用--pres_th 指定，则阈值为1.0
                         * clade_specific_strain_tracker：特定进化的marker列表，用--clade指定进化枝，及其所有子进化枝
--nreads NUMBER_OF_READS：原始宏基因组中的读取总数。 仅当指定 -t marker_table 以使用宏基因组大小对长度归一化计数才使用它；
--pres_th PRESENCE_THRESHOLD：Threshold for calling a marker present by the -t marker_pres_table option
--clade：指定进化枝；
--min_ab：当clade_specific_strain_tracker 分析时进化枝的最小丰度百分比；

Output arguments:
 -o output file, --output_file output file ：输出文件；
 --sample_id_key name  Specify the sample ID key for this analysis.Defaults to 'SampleID'.
--use_group_representative 使用一个物种作为物种组的代表。
--sample_id value     Specify the sample ID for this analysis.Defaults to 'Metaphlan_Analysis'.
-s sam_output_file, --samout sam_output_file ：输出bowtie2比对结果sam文件；
--legacy-output       Old MetaPhlAn2 two columns output
--CAMI_format_output  Report the profiling using the CAMI output format
--unknown_estimation  Scale relative abundances to the number of reads mapping to known clades in order to estimate unknowness
--biom biom_output, --biom_output_file biom_output：输出biom文件；
--mdelim mdelim, --metadata_delimiter_char mdelim    Delimiter for bug metadata: - defaults to pipe. e.g. the pipe in k__Bacteria|p__Proteobacteria 


--nproc N：进行数据库比对使用最小CPU数，默认：4；
--force_download：强制下载最新版本metaphlan数据库；
--read_min_len READ_MIN_LEN：指定'read_fastx.py' 脚本过滤最短reads长度，默认：7

输出文件示例

#bowtie2out

HWUSI-EAS1568_102539179:1:100:10001:7882/1     712117__F3PCC2__HMPREF9056_02717
HWUSI-EAS1568_102539179:1:100:10007:17628/1    712357__A0A0K2JD54__AMK43_09330
HWUSI-EAS1568_102539179:1:100:10017:5224/1     2047__E3H3C1__HMPREF0733_12099

#菌群结构
#mpa_v30_CHOCOPhlAn_201901
#/n/huttenhower_lab/tools/metaphlan3/bin/metaphlan SRS014476 Supragingival_plaque.fasta.gz --input_type fasta
#SampleID       Metaphlan_Analysis
#clade_name     NCBI_tax_id     relative_abundance     additional_species
k__Bacteria     2       100.0   
k__Bacteria|p__Actinobacteria   2|201174        100.0   
k__Bacteria|p__Actinobacteria|c__Actinobacteria 2|201174|1760  100.0   
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Corynebacterials    2|201174|1760|85007     65.25681        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Micrococcales       2|201174|1760|85006     34.74319        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Corynebacterials|f__Corynebacteriaceae      2|201174|1760|85007|1653        65.25681        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Micrococcales|f_Micrococcaceae      2|201174|1760|85006|1268        34.74319        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Corynebacteriales|f__Corynebacteriaceae|g__Corynebacterium   2|201174|1760|85007|1653|1716   65.25681        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Micrococcales|f__Micrococcaceae|g__Rothia    2|201174|1760|85006|1268|32207  34.74319        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Corynebacterials|f__Corynebacteriaceae|g__Corynebacterium|s__Corynebacterium_matruchotii    2|201174|1760|85007|1653|1716|43768     65.25681        
k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Micrococcales|f__Micrococcaceae|g__Rothia|s__Rothia_dentocariosa     2|201174|1760|85006|1268|32207|2047     34.74319        k__Bacteria|p__Actinobacteria|


参考：Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3

MetaPhlAn 3.0 · biobakery/MetaPhlAn Wiki · GitHub

metaphlan3 · biobakery/biobakery Wiki · GitHub

百合木瓜

相关推荐

编程学子看过来，竞赛刷题网站推荐

2022年编程竞赛已经公布，想要在今年取得竞赛成绩的学生，一定要把握寒假时间，学习知识的同时通过刷题，巩固所学知识，提升解题能力。小编为大家推荐几个刷题网站，想要竞赛的学生一定不要错过。USACO美国...

给大家推荐些好的c语言代码的网站

C语言，那就来推荐几个吧，部分含有C++：1、TheLinuxKernelArchives（kernel.org）Linux内核源码，仅限于C，但内核庞大，不太适合新手；2、redis（redi...

推荐几个编程入门学习网站_比较好的编程自学网站

有一些刚上大学的朋友和想对编程感兴趣的朋友经常会让我推荐学习网站，下面几个是我认为零基础学编程比较好的网站，希望大家都有收获！1.W3schoolhttp://www.w3school.com.c...

10个最值得收藏的编程学习网站_有什么学编程的网站

程序员是一个需要不断学习的职业。幸运的是，在这个互联网时代，知识就在那里，等着我们去获取。以下我列举一些免费的编程学习网站包含多个开发语言Java、php、html、javascript等多个。1、h...

6个超酷的练习算法，学习编程的网站

在不了解算法的前提下，您无法通过Google或Facebook的采访。那么为什么不现在学习。我是一位拥有15年以上经验的程序员。从高中开始的第一年，我在算法上学习和工作很多。在我毕业之前，我一直...

在线 python 编程的网站_python3在线编程,python3在线编译器,在线编辑器

以下是一些提供在线Python编程环境的网站：1.Repl.it：Repl.it提供了一个多语言在线编程平台，您可以使用它在任何地方编写、运行、共享代码。Repl.it支持多种编程语言，包括Pyth...

推荐 7 个能过招全球程序员的编程挑战网站，欢迎挑战!

作为程序员的你，是不是经常估不准自己的编程水平？下面推荐7个能过招全球程序员的编程挑战网站，助你磨练技巧，提升技能，最终问鼎代码江湖！1.HackerRank你可以参加各种编码竞赛，比如算法、数学...

盘点 20 个编程学习教程网站，建议收藏

欢迎关注@程序员柠檬橙私信回复「1024」获取海量编程学习资源！如果你想学习编程，现在互联网这么方便，不用着急报名培训班，有很多高质量的编程学习资源网站可供你学习，程序员日常浏览的技术教程网站有哪些...

Flask 数据可视化_flourish数据可视化

数据可视化是数据处理中的重要部分，前面我们了解了Flask的开发和部署，如何用Flask做数据可视化呢？今天我们来了解一下。Python语言极富表达力，并且拥有众多的数据分析库和框架，是数据...

【python 工具】selenium 浏览器操作

selenium的安装步骤：1.安装selenium，打开cmd控制台pipinstallselenium2.安装驱动程序(我这里安装的是chromedriver)，用来启动chrome浏览器...

可视化爬虫工具，EasySpider软件体验

现在提起爬虫，大家可能会联想到Python语言，然后就是各种使用无头浏览器去网页上爬取数据，使用Python的过程相较于使用其他语言来说，简单了不少。但毕竟是编程语言，也需要去学习来适配各种网...

cursor+mcp+playwright，让AI给你推荐五一旅游胜地

阅读本文前提当你已了解mcp是什么，若不知，猛击：https://github.com/modelcontextprotocol/servers。最近有个小需求，根据用户输入内容，使用大模型来理解用户...

Cursor+Claude+Playwright:AI 让自动化测试效率暴涨，快到飞起!

一、引言随着AI时代的到来，软件测试变得越来越复杂，如何高效、准确地进行自动化测试成了每一个开发团队必须面对的问题。在日常工作中，测试工作常常面临各种挑战，比如功能复杂、需求频繁变更、时间紧迫等。传统...

推荐一个检测 JS 内存泄漏的神器_js内存泄漏的几种情况

大家好，我是Echa哥。作为一名Web应用程序开发者，排查和修复JavaScript代码的内存泄漏一直是最困扰我的问题之一。最近，Meta开源了一款检测JavaScript代码内存泄漏...

Python+Playwright自动化实战:高效爬虫全攻略

一、为什么选择Playwright？在信息爆炸的时代，数据获取能力直接决定内容生产效率。Playwright作为微软开源的新型自动化工具，凭借以下优势成为技术创作者的新宠：支持Chromium/Web...

一周热门
最近发表
标签列表