macでインフォマティクス

macでインフォマティクス

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

MetaPhlAn2によるメタゲノムデータの解析

 

 MetaPhlAn2は、メタゲノムシーケンスデータから、どのような生物がどのくらいの割合でいるのか評価するツールである。種の同定が可能なのは、著者らが要した100万以上のマーカー遺伝子が生物と紐付けされていて、そのデータベースの配列にアライメントを行うからである。マーカー遺伝子のみを対象とすることで、ゲノムサイズに対するノーマライズなしに生物種の割合を出すことを可能にしている(DESeq2による正規化は行われている)。

 無論このツールで検出されたからといってメタゲノムが完全に定量的に解釈できると鵜呑みにしてはならないが(データにないゲノムの生物がいるかもしれない)、推定値が出ればqPCRで再現性を調べることもできるし、培養実験の結果と比較することも可能になる。動作が高速なMetaPhlAn2はファーストスクリーニングかつファストスクリーニングにも適したツールと言える。

 

 2016年には、同じ研究チームからPanPhlAnも発表されている。PanPhlAnは種内の多様性を遺伝子の有無で評価するツールである。MetaPhlAn2でメタゲノムにどんな種のゲノムがあるか調べ、PanPhlAnで種内の多様性を調べることが可能になった。

やや稚拙な例えであるが、例えば院内感染で抗生剤の効かない強毒性の菌が増殖したとき、メタゲノムシーケンス→MetaPhlAn2で主要な菌を同定(E.coliとする)→PanPhlAnでゲノム構成を分析、ある遺伝子クラスターに欠損があり→そのクラスターをblast解析すると、そのうちの1つが多剤排出系のMarRによく似たドメインを持つタンパク質だった。という風に、これらのツールを用いることで、菌を同定し、そのgenotypeまで掴むことが可能になる。上記の例で言えば、アウトブレイクが起きた背景が瞬時に掴めるということである。

MetaPhlAn2の使い方を見ていく。

 

 

PanPhlAnは以前紹介しました。

  

マニュアル

https://bitbucket.org/biobakery/biobakery/wiki/metaphlan2

オーサーらのツール一覧

http://segatalab.cibio.unitn.it/tools/ 

 

 

インストール

brewで導入できる。

brew install biobakery/biobakery/graphlan
brew install metaphlan2

 ただしutil/が入らないので別途hg cloneして本体をダウンロードしておく。

hg clone https://bitbucket.org/biobakery/biobakery

 

 

入力ファイル

fasta、fastq、tar.bz2、gzなど。

 

 

ラン 

テストデータダウンロード。ここでは1MB程度のgzがダウンロードされる。

curl -O https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/input/SRS014459-Stool.fasta.gz 
curl -O https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/input/SRS014464-Anterior_nares.fasta.gz
curl -O https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/input/SRS014470-Tongue_dorsum.fasta.gz
curl -O https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/input/SRS014472-Buccal_mucosa.fasta.gz
curl -O https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/input/SRS014476-Supragingival_plaque.fasta.gz
curl -O https://bitbucket.org/biobakery/biobakery/raw/tip/demos/biobakery_demos/data/metaphlan2/input/SRS014494-Posterior_fornix.fasta.gz

 ダウンロードしたデータを移動させる。

mkdir metaphlan2_analysis #新しくディレクトリを作る。
mv *fasta.gz metaphlan2_analysis/ #ダウンロードしたファイルを移動させる
cd metaphlan2_analysis/ #移動

 

fastaデータを1つだけ解析するには以下のように打つ。

metaphlan2.py SRS014476-Supragingival_plaque.fasta.gz  --input_type fasta --nproc 4 > SRS014476-Supragingival_plaque_profile.txt
  •  --nproc 
  • --input_type 

終わると2つのファイルが出来る。

  1. SRS014476-Supragingival_plaque.fasta.gz.bowtie2out.txt 1行1リード形式でマーカー遺伝子とのマッチング結果がまとめられている。中間ファイルだが、どのような生物のマーカー遺伝子とhitしたのか調べる時に使える。
  2. SRS014476-Supragingival_plaque_profile.txt 解析結果。同定できたtaxisonomyがプリントされている。

中身を確認する。

 

今回は4データあるので、それぞれランする。

metaphlan2.py SRS014464-Anterior_nares.fasta.gz --input_type fasta --nproc 4 > SRS014464-Anterior_nares_profile.txt 
metaphlan2.py SRS014470-Tongue_dorsum.fasta.gz --input_type fasta --nproc 4 > SRS014470-Tongue_dorsum_profile.txt
metaphlan2.py SRS014472-Buccal_mucosa.fasta.gz --input_type fasta --nproc 4 > SRS014472-Buccal_mucosa_profile.txt
metaphlan2.py SRS014494-Posterior_fornix.fasta.gz --input_type fasta --nproc 4 > SRS014494-Posterior_fornix_profile.txt

 公式サイトではforで回す例も記載されている。

 

 

4つの結果をマージする。

metaphlan2/utils/merge_metaphlan_tables.py *_profile.txt > merged_abundance_table.txt

ここではワイルドカードを使っているが、順番にファイルを指定しても良い。処理は1旬で終わる。

出力のmerged_abundance_table.txtを確認。

 

less -S merged_abundance_table.txt |head -10

ID SRS014459-Stool_profile SRS014464-Anterior_nares_profile SRS014470-Tongue_dorsum_profile SRS014472-Buccal_mucosa_profile SRS014476-Supragingival_plaque_profile SRS014494-Posterior_fornix_profile

#SampleID Metaphlan2_Analysis Metaphlan2_Analysis Metaphlan2_Analysis Metaphlan2_Analysis Metaphlan2_Analysis Metaphlan2_Analysis

k__Bacteria 100.0 16.82175 100.0 100.0 100.0 100.0

k__Bacteria|p__Actinobacteria 0.0 14.07759 0.85102 0.0 90.32561 0.0

k__Bacteria|p__Actinobacteria|c__Actinobacteria 0.0 14.07759 0.85102 0.0 90.32561 0.0

k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Actinomycetales 0.0 14.07759 0.85102 0.0 90.32561 0.0

k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Actinomycetales|f__Actinomycetaceae 0.0 0.0 0.85102 0.0 0.0 0.0

k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Actinomycetales|f__Actinomycetaceae|g__Actinomyces 0.0 0.0 0.85102 0.0 0.0 0.0

k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Actinomycetales|f__Actinomycetaceae|g__Actinomyces|s__Actinomyces_graevenitzii 0.0 0.0 0.85102 0.0 0.0 0.0

k__Bacteria|p__Actinobacteria|c__Actinobacteria|o__Actinomycetales|f__Actinomycetaceae|g__Actinomyces|s__Actinomyces_graevenitzii|t__Actinomyces_graevenitzii_unclassified 0.0 0.0 0.85102 0.0 0.0 0.0 

 

species情報だけ抜き出す。

grep -E "(s__)|(^ID)" merged_abundance_table.txt | grep -v "t__" | sed 's/^.*s__//g' > merged_abundance_table_species.txt 

出力を確認。

head merged_abundance_table_species.txt

ID SRS014459-Stool_profile SRS014464-Anterior_nares_profile SRS014470-Tongue_dorsum_profile SRS014472-Buccal_mucosa_profile SRS014476-Supragingival_plaque_profile SRS014494-Posterior_fornix_profile

Actinomyces_graevenitzii 0.0 0.0 0.85102 0.0 0.0 0.0

Corynebacterium_matruchotii 0.0 0.0 0.0 0.0 58.21595 0.0

Corynebacterium_pseudodiphtheriticum 0.0 14.07759 0.0 0.0 0.0 0.0

Rothia_dentocariosa 0.0 0.0 0.0 0.0 32.10966 0.0

Bacteroides_cellulosilyticus 3.82206 0.0 0.0 0.0 0.0 0.0

Bacteroides_massiliensis 10.61295 0.0 0.0 0.0 0.0 0.0

Bacteroides_ovatus 4.08051 0.0 0.0 0.0 0.0 0.0

Bacteroides_stercoris 12.82765 0.0 0.0 0.0 0.0 0.0

Parabacteroides_distasonis 3.7393 0.0 0.0 0.0 0.0 0.0

 

hclust2を使いヒートマップを描く。ここではhclust2をダウンロードして、相対パスで実行している。

hg clone https://bitbucket.org/nsegata/hclust2

python2.7 hclust2/hclust2.py -i merged_abundance_table_species.txt -o abundance_heatmap_species.png --ftop 25 --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

hclust2はbrewで導入できる(今回はエラーになったため直接本体を叩いた)。

f:id:kazumaxneo:20170807130405j:plain

 

 ヒートマップが描かれた。黒はその菌がゼロに近いことを意味し、赤-->黄色になると非常に多い(優占種である)ことを意味する。今回の6サンプルは、データによって菌種に共通するものはないことがわかる。ダウンサンプリングされたデータなので、レアなゲノムのデータが閾値以下のゼロに近い値になってしまったのかもしれない。

 

サーバーで解析したなら、 図をscpでローカルに落として確認してください。

 

 

系統樹を描く。

 

GraPhlAnの入力ファイルを作る。

hg clone https://bitbucket.org/nsegata/graphlan

python2.7 graphlan/export2graphlan/export2graphlan.py --skip_rows 1,2 -i merged_abundance_table.txt --tree merged_abundance.tree.txt --annotation merged_abundance.annot.txt --most_abundant 100 --abundance_threshold 1 --least_biomarkers 10 --annotations 5,6 --external_annotations 7 --min_clade_size 1

graphlanもbrewで導入できる(brew install biobakery/biobakery/graphlan )。ここではエラーが出たので、直接ダウンロードして使用している。

 

cladogramを作成。

python2.7 graphlan/graphlan_annotate.py --annot merged_abundance.annot.txt merged_abundance.tree.txt merged_abundance.xml 
python2.7 graphlan/graphlan.py --dpi 300 merged_abundance.xml merged_abundance.png --external_legends

f:id:kazumaxneo:20170807132518j:plain

 

 赤い点が今回の4サンプルから検出された種となる。

 

公式サイトでは他にもチュートリアルがある(リンク)。

 

 

 

 

引用

MetaPhlAn2 for enhanced metagenomic taxonomic profiling

Duy Tin Truong, Eric A Franzosa, Timothy L Tickle, Matthias Scholz, George Weingart, Edoardo Pasolli, Adrian Tett, Curtis Huttenhower & Nicola Segata. Nature Methods 12, 902–903 (2015)

MetaPhlAn2 for enhanced metagenomic taxonomic profiling. - PubMed - NCBI

 

Compact graphical representation of phylogenetic data and metadata with GraPhlAn.

Asnicar F1, Weingart G2, Tickle TL3, Huttenhower C4, Segata N1.

PeerJ. 2015 Jun 18;3:e1029. doi: 10.7717/peerj.1029. eCollection 2015.