初心者の菌叢解析 Qiime2で解析(8) ノイズの除去~DADA2~
今回は前回確認したクオリティを元にノイズ除去（デノイズ）を行います。ここではノイズ除去にDADA2というプラグインを用います。forward配列とreverse配列のマージやクオリティの低い配列（Read）の除去、プライマー配列の除去等を行います。
詳しい内容については概念図を用いた良いサイトがありましたので、参考にしていただければと思います。
1．DADA2の解析
まずは目的のフォルダ（Qiime2_test）に移動し、Qiime2を起動します。
cd /Users/ユーザー名/Desktop/Qiime2_test
conda activate qiime2-2021.2先日作成した「demux.qzv」をQiime2 Viewにドラッグ＆ドロップし、「Interactive Quality Plot」のタブを開きます。ここで、クオリティの著しく減少する部分を確認します。
私は赤枠内の「Middle of Box」の値を参考にしています。ここの値が30を下回り始めたところ以上の配列を全て切っています。今回はforward配列の場合250 bp以上から段々とクオリティが減少しており、reverse配列は240 bpあたりからクオリティが減少しています。
また、場合によっては配列の先頭にプライマー配列を含む場合があります。その時はプライマー配列の分だけカットする必要があります。今回は上流の30 bp程度をカットします。
もしくは「cutadapt」といったプラグインでプライマー配列を切り取ることも出来ます。
以上を踏まえて、以下のコマンドを入れます。
qiime dada2 denoise-paired \
--i-demultiplexed-seqs paired-end-demux.qza \
--p-trim-left-f 30 \
--p-trim-left-r 30 \
--p-trunc-len-f 250 \
--p-trunc-len-r 240 \
--o-representative-sequences rep-seqs-dada2.qza \
--o-table table-dada2.qza \
--o-denoising-stats stats-dada2.qza \
--verbose1行目の「qiime dada2 denoise-paired」はqiimeから始まるコマンドですので、今回行う内容を指示するコマンドになります。DADA2のペアエンド配列用のデノイズを行うということが読み取れると思います。
2行目は「--i-demultiplexed-seqs paired-end-demux.qza 」ですが、今回の動作でインプットするファイルを指示するコマンドです。Qiime2の場合、「--i-」で始まるコマンドはインプットを指示するコマンドです。
3行目は「--p-trim-left-f 30」です。「--p-」はおそらくパラメータの「p」だと思います。forward配列の上流（左側）を切り取り（トリム）を指示するコマンドです。今回は「30」と入れていますので、上流30 bpを除去します。
4行目は「--p-trim-left-r 30 」です。3行目と同じで、reverse配列の上流30 bpの切り取りを指示しています。
5行目は「--p-trunc-len-f 250」です。こちらはforward配列の250 bpより上流を切り取ることを指示するコマンドです。
6行目は「--p-trunc-len-r 240 」です。こちらも5行目と同様で、reverse配列の240 bp以上の切り取りを指示しています。
7~8行目は「--o-」から始まるコマンドで、アウトプットを指示するコマンドです。後ろの.qzaが出力されるファイルの名前になります。好きな名前をつけることも出来ます。例えば解析し直しの2回目ならば「_2」を.qzaの前に追加すれば、区別が出来ると思います。
9行目の「--verbose」は解析の進行具合を表示してくれるコマンドです。「verbose」自体は「詳細、おしゃべり」といった意味なので、このコマンドを入れると進み具合が分かります。
2．stats-dada2.qzvの作成
DADA2の解析は非常に時間がかかります。うまくいけば以下の3つのファイルが生成されているはずです。
1．table-dada2.qza
2．rep-seqs-dada2.qza
3．stats-dada2.qza
次に、DADA2の結果を確認する為に、「stats-dada2.qza」の中身を確認します。.qzaファイルは中身が確認できないので、以下のコマンドを入れます。
qiime metadata tabulate \
--m-input-file stats-dada2.qza \
--o-visualization stats-dada2.qzvこちらはすぐに終わります。出力された「stats-dada2.qzv」を確認してみましょう。Qiime2 Viewにドラッグ＆ドロップしてください。
上の画像のように各サンプルのDADA2での結果が見られます。
「input」はfastqファイルの中のRead数になります。まずフィルターがかけられ、マウスのサンプルでは85%程度、大阪湾では60%程度がパスしていることが分かります。
その後ノイズを除去し、forward配列とreverse配列を結合（マージ）しています。マージ出来た配列はマウスの方で80%程度、大阪湾の方で60%程度となっています。
最後にキメラ配列を除去しており、最終的にはそれぞれ70%程度と50%程度の配列がDADA2をクリアしたようです。
3．DADA2の注意点
3-1．時間について
DADA2は今回の解析の中で最も時間のかかる部分になります。サンプル数が多ければ1日以上かかったこともありました。私は帰る直前にDADA2を実行しそのまま帰宅することもあります。
解析時間はPCのスペックの問題とサンプルの量に依存すると思います。PCのスペックによっては解析できない場合もあります。
3-2．トリミングについて
demux.qzvにてクオリティをチェックしDADA2を実施しますが、きれいな配列を使えばそれだけ最終的な菌叢解析の精度が上昇します。ただ、配列の右側を切り取ることになりますので、切り取りすぎるとforward配列とreverse配列のマージが出来なくなってしまいます。マージの為に20 bp程は残しておく必要があります。
今回では220 bp程度でDADA2を行うとマージが出来ない割合が増えると思います。その時は「stats-dada2.qzv」の「merged」の部分で、Read数が急激に減少すると思います。
ですので、場合によってはクオリティが低くてもある程度は諦めてDADA2を実施する必要があります。
3-3．プライマー除去について
プライマー除去の部分を行わないで今回のDADA2を実行した場合の「stats-dada2.qzv」の結果を以下に示します。
「non-chimeric」の部分で大幅にRead数が減少していることが分かります。つまり、プライマー配列が残っている配列はキメラ配列と認識され除去されているということになります。
プライマー配列は「cutadapt」のプラグインかDADA2で除去する必要があります。また、「non-chimeric」の部分で大幅にRead数が減少した場合はプライマー配列の残存を疑ってみてください。
4．DADA2を予備試験的に試しにやってみる
DADA2は時間がかかりますので、私はサンプル数を極端に減らしたmanifestファイルを作成し、試しにここまでの解析を行ってみることが多々あります。
サンプルが1つだけのmanifestファイルを作成し、「stats-dada2.qzv」の出力まで進めてみます。様々なエラーが出ないことを確認してから全てのサンプルを解析する場合もあります。
トリミングの長さで最終的に残る配列のパーセンテージが大きく変わりますので、是非色々と試してみてください。
ここまで来れば大きな壁は越えた感じがします。皆様が無事進んでいることを願っています。
読んでいただきありがとうございました。
この記事が気に入ったらサポートをしてみませんか？