QIIME2の"Moving Pictures" TutorialをArtifact APIでやってみた

注意（2022/2/7追記）
公式チュートリアルの2021.11版 (https://docs.qiime2.org/2021.11/tutorials/moving-pictures-usage/) で、ついに（？）Artifact APIでもチュートリアルを進められるようになったようです。なのでこの記事は参考にせず、公式チュートリアルを参考にしましょう。

はじめに

このエントリはwet出身のプログラミングとかよくわからん人が、細菌叢解析に使われるポピュラーなパッケージ、QIIME2のチュートリアルをjupyter notebook上で実施するものです。QIIME2は2010年に発表された細菌叢解析パッケージQIIMEの後継で（ https://en.wikipedia.org/wiki/QIIME ）、2019年にNatureの姉妹誌で公表されています。前身であるQIIMEは万単位の論文で引用されており、非常にポピュラーなことがわかります。QIIME2も今後多数の使用が見込まれ、実際qiitaや外部ブログでも複数の方がQIIME2に関する記事を書いていらっしゃいます ¹²³。QIIME2自体についてはわざわざこれ以上説明する必要はないと思いますので詳細は他の記事などを参照ください。

QIIME2には3つのインターフェースが実装されていますが（公式には4つだが、うち一つは解析の共有と可視化にしか使えない）、解析においてはコマンドラインインターフェースを用いる方が多いと思います。ここではpythonによるデータ分析に用いられるjupyter notebookに最適化されたArtifact APIでQIIME2のチュートリアル、“Moving Pictures” tutorial を実施します。

Artifact APIでQIIME2を動かすメリットは次の通りです。

• jupyter notebookはコードとその実行結果を記録します。そのため細かい条件を変えた再解析が必要な時も、再現性の高い解析が可能です。もちろんコマンドラインでもシェルスクリプトを書けばいい話ですが、逐次実行したい自分にはjupyter notebookの方が好都合でした。
• 些細な話ではありますが、コマンドラインで動かす場合はコマンド毎にファイルの読み込みと書き出しが必要になります。Artifact APIで動かす場合は一々ファイルの書き出しをしませんので、その分実行速度が速くなるはずです。（確認はしていませんが）
• 自分の場合、pythonの理解度>>シェルの理解度なので、自分用にスクリプトを書く際にpythonでかけた方が助かります。（個人の感想です）
• QIIME2は解析結果としてVisualizationと呼ばれるオブジェクト（コマンドラインで実行する場合は.qzvファイル）を返しますが、Visualizationを表示するためには一々コマンドを打つか、ブラウザを立ち上げてGUIによる操作を行う必要があります。jupyter notebook上で実行すれば、セル内に結果の表示までできます。

このような動機からArtifact APIでQIIME2を動かそうとした私ですが、公式ドキュメントのArtifact APIの説明はきわめて雑です。また日本語のドキュメントも皆無でしたので、実際に使ってみるとエラーの連発で苦労しました。例えばコマンドラインから一部のメソッドを実行する場合、メタデータとその列を別々に引数として与えてやりますが、Artifact APIで同じようなことをする場合はMetadataオブジェクトのget_columnメソッドを使用してメタデータの列を指定します。このような細かい違いに難儀して中々作業が進まなかった経験がありますので、自分の備忘録も兼ねてここに記録しておきます。どなたかの参考になれば幸いですが、いかんせんdryに弱いので（wetの知識も大してないですが）、細かい用語の間違いやコードの間違いがあれば お許しください。

コード

Qiime2 2019.7をVirtual Box上で実行しています。また、元となっているのは公式チュートリアルの2019.10版（ https://docs.qiime2.org/2019.10/ ）です。
それぞれの工程の意味や、プラグインの詳しい説明は公式ドキュメントや他の記事を参照してください。

必要なライブラリのimport

import os
import pandas as pd
from qiime2 import Artifact
from qiime2 import Metadata
from qiime2.plugins import demux
from qiime2.plugins.dada2.methods import denoise_single
from qiime2.plugins import feature_table
from qiime2.plugins.phylogeny.pipelines import align_to_tree_mafft_fasttree
from qiime2.plugins import diversity
from qiime2.plugins import emperor
from qiime2.plugins.feature_classifier.methods import classify_sklearn
from qiime2.plugins import taxa
from qiime2.plugins import composition

作業ディレクトリの作製と移動

os.makedirs('qiime2-moving-pictures-tutorial', exist_ok=True) # >= Python 3.2
os.chdir('./qiime2-moving-pictures-tutorial')

データのダウンロード

正直この辺の工程はコマンドラインでやったほうが手っ取り早い

# download metadata
os.system("wget -O sample_metadata.tsv https://data.qiime2.org/2019.10/tutorials/moving-pictures/sample_metadata.tsv")

os.makedirs('emp-single-end-sequences', exist_ok=True) # >= Python 3.2

# download barcodes
os.system("wget -O emp-single-end-sequences/barcodes.fastq.gz https://data.qiime2.org/2019.10/tutorials/moving-pictures/emp-single-end-sequences/barcodes.fastq.gz")

# download sequences
os.system("wget -O emp-single-end-sequences/sequences.fastq.gz https://data.qiime2.org/2019.10/tutorials/moving-pictures/emp-single-end-sequences/sequences.fastq.gz")

データ（とメタデータ）のインポート

• Artfact APIインターフェースの場合、メタデータもArtifactととしてインポートする必要がある（CLIの場合も内部的にはそうしていると思われる）のであらかじめインポートしておく
• CLIと違い、Artifactは自動的に保存されないので、保存する場合はsave()メソッドを用いる
• データの取り込みやノイズ除去の過程は時間がかかるので、出力を保存しておいたほうが再解析時に便利

rawseq = Artifact.import_data(type='EMPSingleEndSequences', view='emp-single-end-sequences')
# rawseq.save('emp-single-end-sequences.qza') # save an Artifact as a .qza file
# rawseq = Artifact.load('emp-single-end-sequences.qza') # load an Artifact from a .qza file

metadata = Metadata.load('sample_metadata.tsv')
metadata.to_dataframe() # visualization of loaded metadata

Demultiplexing

CLIの場合、作成した.qzvファイルをいちいちコマンドで開かないと可視化できないが、jupyter notebook上ならノートブック上でそのまま可視化できる

demux_seq, demux_details = demux.methods.emp_single(rawseq, metadata.get_column('barcode-sequence'))
"""
demux_seq.save('demux.qza')
demux_details.save('demux_details.qza')
"""

結果の要約

demux.visualizers.summarize(demux_seq).visualization

Quality ControlとFeature Tableの作製

今回はDADA2（というかDeblurを使ったことがない）

table, rep_seq, stats = denoise_single(demux_seq, trim_left=0, trunc_len=120)
table.save('table_dada2.qza')
rep_seq.save('rep_seq_dada2.qza')
stats.save('stats_dada2.qza')
"""
ここでチュートリアルだとmetadataプラグインのtabulateを利用してstatsを可視化しているが、
Artifact APIでは型が違うと怒られて動かない（し、エラーの回避法もよくわからないので誰か教えてください）
extractメソッドを用いればノートブックで中身を見ること自体は可能だが、ここでは省略
"""

FeatureTableとFeatureDataの可視化

feature_table.visualizers.summarize(table, sample_metadata=metadata).visualization

feature_table.visualizers.tabulate_seqs(rep_seq).visualization

系統樹の作製

aligned_rep_seq, masked_aligned_rep_seq, unrooted_tree, rooted_tree = align_to_tree_mafft_fasttree(rep_seq)

α多様性解析とβ多様性解析

core_metrics_results = diversity.pipelines.core_metrics_phylogenetic(table, rooted_tree, sampling_depth=1103, metadata=metadata)

rarefied_table, faith_pd_vector, observed_otus_vector, shannon_vector, evenness_vector, \
unweighted_unifrac_distance_matrix, weighted_unifrac_distance_matrix, jaccard_distance_matrix, \
bray_curtis_distance_matrix, unweighted_unifrac_pcoa_results, weighted_unifrac_pcoa_results, \
jaccard_pcoa_results, bray_curtis_pcoa_results, unweighted_unifrac_emperor, weighted_unifrac_emperor, \
jaccard_emperor, bray_curtis_emperor = core_metrics_results # unpack core_metrics_results

unweighted_unifrac_emperor # display PCoA result

diversity.visualizers.alpha_group_significance(faith_pd_vector, metadata=metadata).visualization

diversity.visualizers.alpha_group_significance(evenness_vector, metadata=metadata).visualization

diversity.visualizers.beta_group_significance(unweighted_unifrac_distance_matrix, metadata=metadata.get_column('body-site'),                                               pairwise=True).visualization

diversity.visualizers.beta_group_significance(unweighted_unifrac_distance_matrix, metadata=metadata.get_column('subject'), 
                                              pairwise=True).visualization

emperor.visualizers.plot(unweighted_unifrac_pcoa_results, metadata=metadata, 
                         custom_axes=['days-since-experiment-start']).visualization

emperor.visualizers.plot(bray_curtis_pcoa_results, metadata=metadata, 
                         custom_axes=['days-since-experiment-start']).visualization

希薄化曲線

alpha_rarefaction_curve = diversity.visualizers.alpha_rarefaction(table, max_depth=4000, phylogeny=rooted_tree, 
                                                                  metadata=metadata)

# alpha_rarefaction_curve.visualization.save('alpha-rarefaction.qzv')
alpha_rarefaction_curve.visualization

系統解析

# download classifier
os.system("wget -O gg-13-8-99-515-806-nb-classifier.qza https://data.qiime2.org/2019.10/common/gg-13-8-99-515-806-nb-classifier.qza")

classifier = Artifact.load('gg-13-8-99-515-806-nb-classifier.qza')
taxonomy, = classify_sklearn(rep_seq, classifier) 
# ここもmetadataプラグインのtabulateを利用して可視化しているが、ノートブック上での実行方法がわからない

taxonomy.save('taxonomy.qza')
taxonomy.view(pd.DataFrame) # 代わりにpadasのデータフレームとして可視化

taxa.visualizers.barplot(table, taxonomy, metadata).visualization

ANCOM

# creat a new feature table containing only the gut samples
gut_table, = feature_table.methods.filter_samples(table, metadata=metadata, where="[body-site]='gut'")
# add pseudocount
comp_gut_table, = composition.methods.add_pseudocount(gut_table)

composition.visualizers.ancom(comp_gut_table, metadata=metadata.get_column('subject')).visualization

gut_table_l6, = taxa.methods.collapse(gut_table, taxonomy, 6)
comp_gut_table_l6, = composition.methods.add_pseudocount(gut_table_l6)
composition.visualizers.ancom(comp_gut_table_l6, metadata=metadata.get_column('subject')).visualization

1. https://qiita.com/HiguchiRyusuke/items/90ca05c60ab8210c2d7a ↩

2. https://qiita.com/kohei-108/items/5d29bd8bbe1d19fdeec5 ↩

3. https://qiita.com/danryo_official/items/e45afa5aaa0f270fbea8 ↩