macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

MetaPhlAn3

2022/02/24 kronaのコマンド追記

 

 微生物群集の培養によらない解析は、特にショットガン・メタゲノミクスによる生物学的プロファイリングの手法の進歩により、この10年で劇的に進歩した。マルチオミクス、微生物参照ゲノム、株レベルの多様性へのアクセスがより容易になり、改善のチャンスは加速度的に増え続けている。このような背景から、メタゲノムの分類学的、株レベル、機能的、系統的プロファイリングのための統合的かつ改良された手法であるbioBakery 3を発表する。MetaPhlAn 3は分類学的プロファイリングの精度を高め、HUMAn 3は機能的潜在性と活性の精度を向上させた。これらの手法は、CRC(1262メタゲノム)およびIBD(1635メタゲノムおよび817メタトランススクリプトーム)への応用において、新しい疾患とマイクロバイオームとの関連を検出した。StrainPhlAn 3とPanPhlAn 3を用いた4077個のメタゲノムの系統レベルプロファイリングにより、これまで15個の分離ゲノムしか報告されていなかった一般的な腸内細菌Ruminococcus bromiiの系統および機能構造が明らかになった。オープンソースの実装とクラウド展開可能な再現性の高いワークフローにより、bioBakery 3プラットフォームは、微生物群研究のためのマルチオミクスプロファイリングの解像度、規模、精度を深めるために研究者を支援することができる。

 

ここではMetaPhlAn 3について紹介します。

 

HP

https://huttenhower.sph.harvard.edu/metaphlan/

Tutorial

https://github.com/biobakery/biobakery/wiki/metaphlan3

 

バージョン3の新機能(MetaPhlAnのレポジトリより)

  • UniRefをベースにしたChocoPhlAnの新バージョンで抽出したMetaPhlAnマーカー遺伝子を新たに追加
  • パラメータ --unknown_estimation による未知の微生物が構成するメタゲノムの推定
  • パラメータ --index latest による最新のMetaPhlAnデータベースの自動取得とインストール
  • ウィルスプロファイリング --add_viruses
  • 特定のクレードにマッピングされたリードの推定を向上させるためにメタゲノムのサイズを計算
  • NCBI分類IDの出力ファイルへの取り込み
  • CAMI(分類学)プロファイリング出力形式を含む
  • MAPQ値の低いリードの除去

 

インストール

https://github.com/biobakery/MetaPhlAn/wiki/MetaPhlAn-3.0

#conda (link)
mamba create python=3.7 -n metaphlan3 -y
conda activate metaphlan3
mamba install -c bioconda metaphlan -y

#pip
pip install metaphlan

#docker
docker pull biobakery/metaphlan

> metaphlan -h

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.

 

AUTHORS: Francesco Beghini (francesco.beghini@unitn.it),Nicola Segata (nicola.segata@unitn.it), Duy Tin Truong, Francesco Asnicar (f.asnicar@unitn.it), Aitor Blanco Miguez (aitor.blancomiguez@unitn.it)

 

COMMON COMMANDS

 

 We assume here that MetaPhlAn is installed using the several options available (pip, conda, PyPi)

 Also BowTie2 should be in the system path with execution and read permissions, and Perl should be installed)

 

========== MetaPhlAn clade-abundance estimation ================= 

 

The basic usage of MetaPhlAn consists in the identification of the clades (from phyla to species ) 

present in the metagenome obtained from a microbiome sample and their 

relative abundance. This correspond to the default analysis type (-t rel_ab).

 

*  Profiling a metagenome from raw reads:

$ metaphlan metagenome.fastq --input_type fastq -o profiled_metagenome.txt

 

*  You can take advantage of multiple CPUs and save the intermediate BowTie2 output for re-running

   MetaPhlAn extremely quickly:

$ metaphlan metagenome.fastq --bowtie2out metagenome.bowtie2.bz2 --nproc 5 --input_type fastq -o profiled_metagenome.txt

 

*  If you already mapped your metagenome against the marker DB (using a previous MetaPhlAn run), you

   can obtain the results in few seconds by using the previously saved --bowtie2out file and 

   specifying the input (--input_type bowtie2out):

$ metaphlan metagenome.bowtie2.bz2 --nproc 5 --input_type bowtie2out -o profiled_metagenome.txt

 

*  bowtie2out files generated with MetaPhlAn versions below 3 are not compatibile.

   Starting from MetaPhlAn 3.0, the BowTie2 ouput now includes the size of the profiled metagenome and the average read length.

   If you want to re-run MetaPhlAn using these file you should provide the metagenome size via --nreads:

$ metaphlan metagenome.bowtie2.bz2 --nproc 5 --input_type bowtie2out --nreads 520000 -o profiled_metagenome.txt

 

*  You can also provide an externally BowTie2-mapped SAM if you specify this format with 

   --input_type. Two steps: first apply BowTie2 and then feed MetaPhlAn with the obtained sam:

$ bowtie2 --sam-no-hd --sam-no-sq --no-unal --very-sensitive -S metagenome.sam -x ${mpa_dir}/metaphlan_databases/mpa_v30_CHOCOPhlAn_201901 -U metagenome.fastq

$ metaphlan metagenome.sam --input_type sam -o profiled_metagenome.txt

 

*  We can also natively handle paired-end metagenomes, and, more generally, metagenomes stored in 

  multiple files (but you need to specify the --bowtie2out parameter):

$ metaphlan metagenome_1.fastq,metagenome_2.fastq --bowtie2out metagenome.bowtie2.bz2 --nproc 5 --input_type fastq

 

------------------------------------------------------------------- 

 

 

========== Marker level analysis ============================ 

 

MetaPhlAn introduces the capability of characterizing organisms at the strain level using non

aggregated marker information. Such capability comes with several slightly different flavours and 

are a way to perform strain tracking and comparison across multiple samples.

Usually, MetaPhlAn is first ran with the default -t to profile the species present in

the community, and then a strain-level profiling can be performed to zoom-in into specific species

of interest. This operation can be performed quickly as it exploits the --bowtie2out intermediate 

file saved during the execution of the default analysis type.

 

*  The following command will output the abundance of each marker with a RPK (reads per kilo-base) 

   higher 0.0. (we are assuming that metagenome_outfmt.bz2 has been generated before as 

   shown above).

$ metaphlan -t marker_ab_table metagenome_outfmt.bz2 --input_type bowtie2out -o marker_abundance_table.txt

   The obtained RPK can be optionally normalized by the total number of reads in the metagenome 

   to guarantee fair comparisons of abundances across samples. The number of reads in the metagenome

   needs to be passed with the '--nreads' argument

 

*  The list of markers present in the sample can be obtained with '-t marker_pres_table'

$ metaphlan -t marker_pres_table metagenome_outfmt.bz2 --input_type bowtie2out -o marker_abundance_table.txt

   The --pres_th argument (default 1.0) set the minimum RPK value to consider a marker present

 

*  The list '-t clade_profiles' analysis type reports the same information of '-t marker_ab_table'

   but the markers are reported on a clade-by-clade basis.

$ metaphlan -t clade_profiles metagenome_outfmt.bz2 --input_type bowtie2out -o marker_abundance_table.txt

 

*  Finally, to obtain all markers present for a specific clade and all its subclades, the 

   '-t clade_specific_strain_tracker' should be used. For example, the following command

   is reporting the presence/absence of the markers for the B. fragilis species and its strains

   the optional argument --min_ab specifies the minimum clade abundance for reporting the markers

 

$ metaphlan -t clade_specific_strain_tracker --clade s__Bacteroides_fragilis metagenome_outfmt.bz2 --input_type bowtie2out -o marker_abundance_table.txt

 

------------------------------------------------------------------- 

 

positional arguments:

  INPUT_FILE            the input file can be:

                        * a fastq file containing metagenomic reads

                        OR

                        * a BowTie2 produced SAM file. 

                        OR

                        * an intermediary mapping file of the metagenome generated by a previous MetaPhlAn run 

                        If the input file is missing, the script assumes that the input is provided using the standard 

                        input, or named pipes.

                        IMPORTANT: the type of input needs to be specified with --input_type

  OUTPUT_FILE           the tab-separated output file of the predicted taxon relative abundances 

                        [stdout if not present]

 

Required arguments:

  --input_type {fastq,fasta,bowtie2out,sam}

                        set whether the input is the FASTA file of metagenomic reads or 

                        the SAM file of the mapping of the reads against the MetaPhlAn db.

 

Mapping arguments:

  --force               Force profiling of the input file by removing the bowtie2out file

  --bowtie2db METAPHLAN_BOWTIE2_DB

                        Folder containing the MetaPhlAn database. You can specify the location by exporting the DEFAULT_DB_FOLDER variable in the shell.[default /home/kazu/mambaforge/envs/metaphlan3/lib/python3.7/site-packages/metaphlan/metaphlan_databases]

  -x INDEX, --index INDEX

                        Specify the id of the database version to use. If "latest", MetaPhlAn will get the latest version.

                        If an index name is provided, MetaPhlAn will try to use it, if available, and skip the online check.

                        If the database files are not found on the local MetaPhlAn installation they

                        will be automatically downloaded [default latest]

  --bt2_ps BowTie2 presets

                        Presets options for BowTie2 (applied only when a FASTA file is provided)

                        The choices enabled in MetaPhlAn are:

                         * sensitive

                         * very-sensitive

                         * sensitive-local

                         * very-sensitive-local

                        [default very-sensitive]

  --bowtie2_exe BOWTIE2_EXE

                        Full path and name of the BowTie2 executable. This option allowsMetaPhlAn to reach the executable even when it is not in the system PATH or the system PATH is unreachable

  --bowtie2_build BOWTIE2_BUILD

                        Full path to the bowtie2-build command to use, deafult assumes that 'bowtie2-build is present in the system path

  --bowtie2out FILE_NAME

                        The file for saving the output of BowTie2

  --min_mapq_val MIN_MAPQ_VAL

                        Minimum mapping quality value (MAPQ) [default 5]

  --no_map              Avoid storing the --bowtie2out map file

  --tmp_dir             The folder used to store temporary files [default is the OS dependent tmp dir]

 

Post-mapping arguments:

  --tax_lev TAXONOMIC_LEVEL

                        The taxonomic level for the relative abundance output:

                        'a' : all taxonomic levels

                        'k' : kingdoms

                        'p' : phyla only

                        'c' : classes only

                        'o' : orders only

                        'f' : families only

                        'g' : genera only

                        's' : species only

                        [default 'a']

  --min_cu_len          minimum total nucleotide length for the markers in a clade for

                        estimating the abundance without considering sub-clade abundances

                        [default 2000]

  --min_alignment_len   The sam records for aligned reads with the longest subalignment

                        length smaller than this threshold will be discarded.

                        [default 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

                        [default 0.2]

  --perc_nonzero        Percentage of markers with a non zero relative abundance for misidentify a species

                        [default 0.33]

  --ignore_markers IGNORE_MARKERS

                        File containing a list of markers to ignore. 

  --avoid_disqm         Deactivate the procedure of disambiguating the quasi-markers based on the 

                        marker abundance pattern found in the sample. It is generally recommended 

                        to keep the disambiguation procedure in order to minimize false positives

  --stat                Statistical approach for converting marker abundances into clade abundances

                        '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

                        [default tavg_g]

 

Additional analysis types and arguments:

  -t ANALYSIS TYPE      Type of analysis to perform: 

                         * rel_ab: profiling a metagenomes in terms of relative abundances

                         * rel_ab_w_read_stats: profiling a metagenomes in terms of relative abundances and estimate the number of reads coming from each clade.

                         * reads_map: mapping from reads to clades (only reads hitting a marker)

                         * clade_profiles: normalized marker counts for clades with at least a non-null marker

                         * marker_ab_table: normalized marker counts (only when > 0.0 and normalized by metagenome size if --nreads is specified)

                         * marker_counts: non-normalized marker counts [use with extreme caution]

                         * marker_pres_table: list of markers present in the sample (threshold at 1.0 if not differently specified with --pres_th

                         * clade_specific_strain_tracker: list of markers present for a specific clade, specified with --clade, and all its subclades

                        [default 'rel_ab']

  --nreads NUMBER_OF_READS

                        The total number of reads in the original metagenome. It is used only when 

                        -t marker_table is specified for normalizing the length-normalized counts 

                        with the metagenome size as well. No normalization applied if --nreads is not 

                        specified

  --pres_th PRESENCE_THRESHOLD

                        Threshold for calling a marker present by the -t marker_pres_table option

  --clade               The clade for clade_specific_strain_tracker analysis

  --min_ab              The minimum percentage abundance for the clade in the clade_specific_strain_tracker analysis

 

Output arguments:

  -o output file, --output_file output file

                        The output file (if not specified as positional argument)

  --sample_id_key name  Specify the sample ID key for this analysis. Defaults to 'SampleID'.

  --use_group_representative

                        Use a species as representative for species groups.

  --sample_id value     Specify the sample ID for this analysis. Defaults to 'Metaphlan_Analysis'.

  -s sam_output_file, --samout sam_output_file

                        The sam output file

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

                        If requesting biom file output: The name of the output file in biom format 

  --mdelim mdelim, --metadata_delimiter_char mdelim

                        Delimiter for bug metadata: - defaults to pipe. e.g. the pipe in k__Bacteria|p__Proteobacteria 

 

Other arguments:

  --nproc N             The number of CPUs to use for parallelizing the mapping [default 4]

  --install             Only checks if the MetaPhlAn DB is installed and installs it if not. All other parameters are ignored.

  --force_download      Force the re-download of the latest MetaPhlAn database.

  --read_min_len READ_MIN_LEN

                        Specify the minimum length of the reads to be considered when parsing the input file with 'read_fastx.py' script, default value is 70

  -v, --version         Prints the current MetaPhlAn version and exit

  -h, --help            show this help message and exit

 

#version - MetaPhlAn version 3.0.14 (19 Jan 2022)

 

 

実行方法

fastqファイルを指定する。初回はデータベースがダウンロードされるので、余分に時間がかかる。20スレッド指定。

metaphlan input.fastq.gz --input_type fastq --nproc 20 > output.txt

 

複数のデータをランして結果をマージする。

#1 MetaPhlAn3のラン- for loop
#single-end
for i in illumina*.fq.gz; do metaphlan $i --input_type fastq --nproc 20 > ${i%.fq.gz}_profile.txt; done

#paired-end
#ペアエンドのデータは結合して指定する(ワイルドカードで指定するとR2側が消えたので注意)
for i in *_R1.fq.gz; do
zcat ${i%_R1.fq.gz}_R?.fq.gz > merged.fq
metaphlan merged.fq --input_type fastq --nproc 20 > ${i%_R1.fq.gz}_profile.txt
done

#2 merge
merge_metaphlan_tables.py *_profile.txt > merged_abundance_table.txt

#3 種のみの存在量のテーブルを作成
grep -E "s__|clade" merged_abundance_table.txt | sed 's/^.*s__//g'\
| cut -f1,3-8 | sed -e 's/clade_name/body_site/g' > merged_abundance_table_species.txt

#4 hclust (python2)でtop50をクラスタリングしてヒートマップで視覚化
mamba create python=2.7 -n hclust2 -c bioconda hclust2 -y
conda activate hclust2
hclust2.py -i merged_abundance_table_species.txt \
 -o abundance_heatmap_species.png \
--ftop 50 \
 --f_dist_f braycurtis \
 --s_dist_f braycurtis \
 --cell_aspect_ratio 0.5 -l \
--flabel_size 6 \
--slabel_size 6 \
 --max_flabel_len 100 \
 --max_slabel_len 100 \
 --minv 0.1 \
--dpi 300

クラスタリングで失敗する場合、--no_fclusteringと--no_sclusteringをつけてランする。うまくランできるなら、クラスタリングで失敗している。これは0が大半のデータを含んでると起きる可能性がある。対処方法として、クラスタリング方法を変えるか、非常に小さな値(0.01など)を全ての値に加える。

 

 

optional step,  kronaのラン

metaphlan2krona.pyを使う(condaでmetaphlan2を導入するとパスが通る)。

$ metaphlan2krona.py -h
Usage: metaphlan2krona.py [options]

Options:
  -h, --help            show this help message and exit
  -p PROFILE, --profile=PROFILE
                        The input file is the MetaPhlAn standard result file
  -k KRONA, --krona=KRONA
                        the Krons output file name

変換後、kronaをランする。テストした時はmetaphlan2krona.pyのランに失敗した。

#kronaのインストールとデータベースの準備
mamba install -c bioconda -y krona
rm -rf $HOME/miniconda3/envs/kraken2/opt/krona/taxonomy
mkdir -p krona/taxonomy
ktUpdateTaxonomy.sh krona/taxonomy/
cp -r krona/taxonomy/ $HOME/miniconda3/envs/kraken2/opt/krona/

#kronaのラン
ktImportTaxonomy -q 2 -t 3 report_bracken_species.txt -o krona.html

#全サンプル
ktImportTaxonomy -q 2 -t 3 *_bracken_species.txt -o krona.html

 

引用

Integrating taxonomic, functional, and strain-level profiling of diverse microbial communities with bioBakery 3
Francesco Beghini, Lauren J McIver, Aitor Blanco-Miguez, Leonard Dubois, Francesco Asnicar, Sagun Maharjan, Ana Mailyan, Paolo Manghi, Matthias Scholz, Andrew Maltez Thomas, Mireia Valles-Colomer, George Weingart, Yancong Zhang, Moreno Zolfo, Curtis Huttenhower, Eric A Franzosa, and Nicola Segata

Elife. 2021 May 4;10:e65088

 

関連


 2022/04/29

"あるメタゲノムサンプルでMetaPhlAnを実行すると、10種しか同定できないか、全く同定できないかもしれない一方で、同じサンプルでKrakenを実行すると、数万種まで同定できるかもしれないことに気づきました。"