2021年12月16日木曜日

環境 DNA メタバーコーディングとユニバーサルプライマー

20211216_Blogger0017

環境 DNA メタバーコーディングでは自分たちが興味ある特定の分類群に属する生物の DNA をまとめて増幅するための ｢ユニバーサルプライマー｣ が使用されます。調査･研究の目的に合ったユニバーサルプライマーを選ぶことは環境 DNA メタバーコーディングにおいて最も重要なステップの一つです。

この ｢適したユニバーサルプライマーを選ぶ｣ という作業、簡単なようでその難易度は研究目的や解析対象の分類群によりかなりばらつきがあります。今回はユニバーサルプライマーに関連して、以下の項目の記事を書いてみようと思います。

1. ユニバーサルプライマーを選ぶ際に考慮する条件
2. 環境 DNA メタバーコーディングで使用されるユニバーサルプライマー ･ 生物分類群 ･ シーケンサー
3. 自分の研究対象の分類群/種はユニバーサルプライマーで増えるのか ?

ユニバーサルプライマーを選ぶ際に考慮する条件

ユニバーサルプライマーの選択の際に考慮すべき条件は概ね以下のとおりです。

  1. 生物分類群 : 自分たちが興味ある分類群の DNA を増幅できること。これは必須の条件です。また、自分たちが興味ない分類群の DNA を増幅しないこと/増幅しすぎないこと、も重要です。
  2. (種) 同定の解像度 : 解読された DNA 配列を用いて ｢どの程度の階層 (種･属･科など) まで生物の分類群を同定できるか｣ も重要な指標です。
  3. 増幅産物の長さ : 環境 DNA メタバーコーディングでは高スループットなシーケンサーを用いて配列解読を行いますが、シーケンサー毎に解読可能な配列長が異なります。したがって、どのシーケンサーで配列を解読予定か、予め考慮しなければなりません。

ユニバーサルプライマーの開発論文では、コンピュータ上での計算や実際にサンプルを用いた実験により、上記に関する様々な指標の評価が行われています。なお、本記事ではユニバーサルプライマーの開発方法そのものについては解説しません。

環境 DNA メタバーコーディングで使用されるユニバーサルプライマー ･ 生物分類群 ･ シーケンサー

環境 DNA メタバーコーディングで使用可能なユニバーサルプライマー･対応する生物分類群とシーケンサーをいくつか列挙します。当然 (?) 列挙しきることは不可能なので下の方で参考文献･ウェブサイトも紹介します (それでも足りないですが)。

生物分類群とユニバーサルプライマー

その他、哺乳類 (Ushio et al. 2017)、両生類 (Sakata et al. 2021)、水生昆虫 (Leese et al. 2020)、エビ･カニ類 (Komai et al. 2019) などについてのユニバーサルプライマーも報告されています。真核生物の 18S rRNA 遺伝子領域については PR2 Primer Database で多くの情報が手に入ります。また、Environmental DNA: For Biodiversity Research and Monitoring の巻末にあるユニバーサルプライマーリストも非常に便利です。Earch Microbiome Project のウェブサイト や 東北大の田邉さんのウェブサイト にも記載があります。

対応するシーケンサー

よく使われるプライマーとシーケンサーの組み合わせは以下のとおりです。一般的に解析する DNA が短い方がシーケンサーの選択肢が多く、シーケンス試薬も安い傾向があります。一方で、短い DNA は解析する領域が同じ場合は、種同定の分解能 (科まで同定可能か、種まで同定可能か、など) は低い傾向にあります。

  • 150 bp x 2 Paired End (iSeq 100, MiSeq, HiSeq, NovaSeq など)

    • MiFish-U-F/R
    • 1391F - EukBr

  • 250 bp x 2 Paired End (MiSeq, NovaSeq など)

    • mlCOIintF - HCO2198
    • 515F - 806R
    • 515F-Y - 926R

  • 300 bp x 2 Paired End (MiSeq のみ)

    • ITS1_F_KYO1 - ITS2_KYO2

また、一般的な組み合わせは以上のとおりですが、515F - 806R や ITS1_F_KYO1 - ITS2_KYO2 を iSeq 100 のシングルエンドで 300 bp 程度連続で読んだりすることもできます。このあたりは、使える予算･かけられる時間･求める種同定の分解能などを考慮して決めることになります。

また、近年は Oxford Nanopore MinION などロングリードシーケンサーの登場により環境 DNA でもより長い領域を解析するような研究も報告されつつあります。

ユニバーサルプライマーで自分の興味ある種の DNA は増えるのか?

さて、調査対象の生物分類群全般を増幅してくれるユニバーサルプライマーは非常に頼もしいですが、当然ながら対象分類群に含まれる生物種の DNA を 100% 増幅してくれるわけではありません。例えば、魚類全般を対象とした MiFish-U-F/R を用いても DNA が増えない or 増えにくい魚種が存在します。 例えば ヤツメウナギの DNA は増えないそうですし、エイ･サメの仲間の DNA も増えにくいようです (なので MiFish-E-F/R がデザインされている; Miya et al. 2015)。

このような状況下で、例えば次のような疑問がよく生まれます:
“昆虫を mlCOIintF - HCO2198 のプライマーを用いて環境サンプル中の DNA を網羅的に解析予定だが、このプライマーは自分が特に興味あるハチの仲間の DNA も (もし含まれていたら) 増幅してくれるのだろうか？”

最近、特にこういった質問をよくされます。そこで、いくつか典型的なやり方をまとめてみました。 (調べましたが、以下の 3 つ以外の方法はよく分かりませんでした。もし他にいい方法があったら教えて下さい)

1. テキストファイルで検索 (眺める)

最もシンプルなのはこの方法です。興味ある生物の DNA 配列 (FASTA ファイル) がすでにある場合は、そのファイルをテキストエディタで開き、プライマーの配列を検索する、ということが可能です。

例えば、以下のようにクサフグのミトコンドリア DNA の配列が手元にあるとします。

クサフグのミトコンドリア DNA の配列

このテキストファイルに MiFish-U-F の配列 (GTCGGTAAAACTCGTGCCAGC) で検索をかけます。

MiFish-U-F の配列で検索

MiFish-U-F の配列が 1 箇所だけにヒットしました。どうやら一致する領域があるようです。さらに、図は出していませんが、MiFish-U-R の配列の逆相補配列 (“CAAACTGGGATTAGATACCCCACTATG”) に一致する場所も見つかります。したがって、クサフグの DNA は MiFish-U-F/R のプライマーセットで増幅できそう、ということになります。

似たようなことは FASTA ファイルの表示に特化した他のソフトウェアでも可能でしょう。

この方法は、特別な技術やソフトウェアが不要で直感的にわかりやすいのでお手軽かもしれません。ただし、配列が手元にないと調べられませんし、調べたい配列数やプライマー数が増えると大変です。また、上記 515F-Y - 926R などに見られるような縮重塩基を含むようなプライマーと一致する領域を検索で探すのも大変です。さらに、プライマー - 鋳型 DNA のミスマッチを許容しつつ探す、ということはできません。

2. NCBI + seqkit amplicon

2 番目はもう少し網羅的な方法で、少し配列解析っぽくなります。

配列データをダウンロード

まず、NCBI (https://www.ncbi.nlm.nih.gov/nuccore) などの公共のデータベースで興味ある種や分類群の DNA 配列を検索します。｢1｣ と同じく Takifugu niphobles (もしくは Takifugu alboplumbeus?) を検索します。ミトコンドリア DNA の配列が欲しいので、｢Takifugu niphobles AND mitochondrion｣ で検索してみます。

NCBI Nucleotide で検索

3 件ヒットしました。これらの配列をダウンロードしたい場合は、チェックボックスにチェックを入れます。

図4. 検索結果

その後、ページの上の方にある ｢Send to｣ のボタンをクリックして、｢Complete Record｣ → ｢File｣ → ｢FASTA｣ → ｢Create File｣ を選択すると FASTA ファイルとして DNA 配列がダウンロードできます。

配列のダウンロード

seqkit を使ってプライマーと一致する領域を抽出

さて、ダウンロードした配列ファイルは seqkit というツールで調べます。

seqkit のウェブサイト

配列ファイルを操作する様々な機能を備えており、かつ軽くて高速、使い勝手もいいです。非常におすすめのツールです。インストール方法はこちら。

macOS の場合は上記のウェブサイトで seqkit_darwin_amd64.tar.gz を選択してダウンロードします。その後、ターミナルで以下をタイプします。

  1. tar -zxvf seqkit_darwin_amd64.tar.gz
  2. sudo cp seqkit /usr/local/bin/
  3. seqkit version とタイプしてバージョンが表示されればインストールされています。

Windows の場合は上記のウェブサイトで seqkit_windows_amd64.exe.tar.gz を選択します。その後、コマンドプロンプトで以下をタイプします。

  1. tar -zxvf seqkit_windows_amd64.exe.tar.gz
  2. seqkit.exe を C:\WINDOWS\system32 にコピー。
  3. seqkit version とタイプしてバージョンが表示されればインストールされています。

インストールができたら、NCBI からダウンロードした配列 (sequence.fasta というファイル) があるフォルダに移動して、ターミナルかコマンドプロンプトで以下をタイプして、seqkit amplicon の機能を使用します。

seqkit amplicon -F GTCGGTAAAACTCGTGCCAGC -R CATAGTGGGGTATCTAATCCCAGTTTG -w 0 -m 0 sequence.fasta

-F にフォーワードプライマー配列、-R にリバースプライマー配列を指定します (ここでは MiFish-U-F/R)。-w 0 は出力する DNA 配列の改行なしのオプション、-m 0 はミスマッチ数 0 の仮定のもとでの抽出を行うオプションです。

すると次のように画面上にプライマーと一致する領域およびその間の領域、つまりは PCR で増幅されうる領域が切り出されます。

seqkit amplicon でプライマーと一致する領域を抽出

もしくは、下記のように入力することで出力結果をファイルに保存できます。

seqkit amplicon -F GTCGGTAAAACTCGTGCCAGC -R CATAGTGGGGTATCTAATCCCAGTTTG -w 0 -m 0 sequence.fasta > amplified.fasta

この方法を使えば手持ちの配列がなくても、興味ある生物の DNA が特定のユニバーサルプライマーで増えそうかどうかを調べることができます。また、seqkit amplicon の -m オプションを指定すればミスマッチの数を指定して調べることもできます。さらに、プライマー配列に縮重塩基が含まれていても大丈夫です。

3. カスタムスクリプトを作る

多少のコーディング作業がありますが、｢2｣ のやり方でほぼ完璧なように思えます。しかし、いくつもの生物分類群やプライマーを調べる必要がある場合は、上記の検索作業を何度も繰り返すことになります。検索して、ファイルをダウンロードして、ターミナルを開いて seqkit amplicon を実行、を繰り返すのはなかなか面倒なものです。

基本的には ｢配列のダウンロード｣ と ｢seqkit amplicon｣ が使えれば全て解決するはずなので、全ての作業を R で記述し、一つの関数にすることにしました。少し探したのですが、幸い rentrez という良さそうなパッケージがあったのでそれを利用することにしました。

R から配列のダウンロード

まず、検索ワードを指定して配列の ID などを収集します。

# rentrez のインストールと読み込み
devtools::install_github("ropensci/rentrez")
library("rentrez"); packageVersion("rentrez") # v1.2.3

# 検索ワードを指定 (Trachurus かつ mitochondrion かつ 配列長 1000~20000)
SEARCH_QUERY <- "Trachurus AND mitochondrion AND 1000:20000[SLEN]"

# 検索
rentrez_search <- entrez_search(db = "nucleotide", term = SEARCH_QUERY)

配列をダウンロードして、一度書き出して保存します。

# 配列をダウンロード
entrez_fasta <- entrez_fetch(db = "nucleotide", id = rentrez_search$ids, rettype="fasta", parsed = FALSE)

# FASTA ファイルとして書き出し
write.table(entrez_fasta, file = "entrez_fasta.fa", quote = FALSE, col.names = FALSE, row.names = FALSE)

R の system 関数を使って seqkit amplicon を R から実行して、プライマー (ここでは MiFish-U-F/R) と一致する領域を持つ配列を抽出します。

system("seqkit amplicon -F GTCGGTAAAACTCGTGCCAGC -R CATAGTGGGGTATCTAATCCCAGTTTG -w 0 -m 1 entrez_fasta.fa > amplified.fa")

これで amplified.fa という FASTA ファイルに MiFish-U-F/R で増幅されうる領域が書き出されます。以上の操作を関数としてまとめれば、検索ワードを変更するだけでどんどん ｢自分の研究対象種/分類群はユニバーサルプライマーで増えるのか ?｣ を調べることができます。

おまけ！

実は seqkit amplicon は ｢縮重塩基を含むプライマー + -m オプションを使ったミスマッチ数の同時指定｣ ができません (これ、昔はまって seqkit の開発者に質問したことがあります; 正確に言うと seqkit 内で解決できる問題ですが、手作業でやると少々ややこしいです)。-F もしくは -R に縮重塩基を含むプライマーを指定すると -m オプションは無視されます。

そのへんの問題を自作の関数で解決しつつ、上記のカスタムスクリプトを関数としてまとめてみました。ついでに以前から興味があった R package の作成に挑戦するのにいい機会 (?) ということで、単純な関数ではありますがパッケージ化を試みました。成果物が以下です。

もし ｢自分の研究対象種/分類群はユニバーサルプライマーで増えるのか ?｣ を調べる必要があったら使ってあげてください。

rDoAMP (DOes my primer set AMPlify my targets?)
https://github.com/ong8181/rDoAMP

詳しくは Github のページを読んでもらうといいのですが、devtools::install_github("ong8181/rDoAMP") でインストール可能です。seqkit が正しくインストールされていれば macOS と Windows で動作します。Linux でも動くはずですが未確認です。

縮重塩基を含むプライマーの問題は expand_degenerate_primer() という関数を用意して、あり得るプライマーの組み合わせを tsv ファイルとして書き出し、それを用いることで解決しています。例えば以下のようなコマンドを使っています。-p プライマーリスト のオプションと -m のオプションは同時使用可能です。

seqkit amplicon -p primer_list.tsv -m 3 download.fa | seqkit rmdup -w 0

また、entrez_search はヒットした ID を上から順番に取ってくるので、似通った配列の ID が連続して取得されることが多いです。そこで、デフォルトでは予めダウンロード予定の配列数の 10 倍の ID を取得しておいて、その中からランダムにダウンロードする配列を選ぶ、という風にしています。こうすることで取得する配列の多様性が挙がります。

(いろいろ小細工しましたが、ほぼほぼ seqkit amplicon です。seqkit すごいです)

R package の作成は以下のウェブサイトを参考にさせていただきました (どうもありがとうございます)。楽しかったです (もっと勉強します)。

参考文献･情報

Written with StackEdit.

時刻: 3 件のコメント:

2021年8月14日土曜日

定量的な環境 DNA メタバーコーディング

20210814_Blogger0016

近年、環境 DNA メタバーコーディングは ｢野外や実験室でのプロトコルの確立｣ ｢シーケンス技術の発達｣ ｢データ解析パイプラインの充実｣ などのおかげで、野外における生物相調査の方法として非常に popular な技術になりつつあります (もうなっている？)。しかし、まだまだ手法として検討･改善の余地のある部分も多くあります。環境 DNA をどのように定量するか、はその一つであり、現在でも様々な改善が試みられています。環境 DNA 分析を行う人、行いたい人にとっても定量は興味があるところのようで、最近とくに定量的環境 DNA メタバーコーディングに関して多く質問をもらいます。

そこで今回は、定量的な環境 DNA メタバーコーディングについての記事を書こうと思います。

なお、本記事内では ｢環境 DNA｣ に微生物 DNA も含んで話を進めます。また、環境 DNA = Environmental DNA = eDNA と記載します。eDNA そのものについての基礎的な部分はすでに Wikipedia にも日本語の解説があります。また、最近、『環境 DNA ―生態系の真の姿を読み解く―』という日本語の本も出版されましたので、興味ある方はご一読いただければと思います (自分も一部執筆しています)。

eDNA メタバーコーディングにおける ｢定量｣ とは

まず、｢eDNA を定量的に分析する｣ と言った場合に、その意味することに２つのパターンがありえます。

1. サンプル中の eDNA の濃度を定量する

｢eDNA を定量的に分析する｣ を文字通りに解釈するとこちらの意味となります。例えば、水サンプル中に調査対象となる生物種の DNA がどの程度含まれているか、その濃度を測ります。最終的な結果は “xx copies/ml water” とか、“xx ng/ml water” といった単位で表されるでしょう。 対象種が 1 種であれば、定量 PCR を用いることで定量が達成できます。

2. eDNA の濃度から調査対象の種の個体数やバイオマスを推定する

文字通りの意味とはやや異なりますが、しばしば勘違いされるのがこちらの意味です。実際の生物調査では、調査する人が最終的に欲しいのは eDNA の濃度そのものではなく、調査対象種の個体数やバイオマスの情報であることがほとんどです。現状では、調査対象種の個体数やバイオマスを知るために、｢eDNA の濃度が調査対象種の個体数やバイオマスを反映している｣ と仮定して eDNA を定量していることが多いです。この仮定は、魚由来の eDNA を水サンプルから抽出して定量する場合はどうやら概ね正しいようです (Takahara et al. 2012; Yamamoto et al. 2016; Stoeckle et al. 2017)。また、個体 (細胞) そのものをまるごと分析することになる微生物が対象の場合も概ね正しいでしょう。ただし、他の生物や他のサンプルではどうか、は注意深く考える必要があります。例えば、森林の池から陸生哺乳類の eDNA が検出可能であることが分かっていますが (Ushio et al. 2017)、その濃度が森林全体の陸生哺乳類の量とどのように関係しているかはわかりません。また、魚や微生物の場合でも、eDNA の濃度のみから厳密な個体数の推定は容易ではありません。魚であれば、種･行動･成長段階によって eDNA の放出速度が変わるでしょうし、オープンな水域であれば eDNA の拡散過程や分解速度を考慮する必要もあるでしょう (労力をかけてこれを行ったのが Fukaya et al. 2021)。細菌の 16S rRNA 遺伝子を対象とした解析の場合は、各細菌種がもつ 16S rRNA 遺伝子コピーの数が異なることも問題になるでしょう。

いずれにしても、｢2｣ の意味での定量はかなりハードルが高く、未だに革新的な技術は生み出されていません。今回はその前段階である ｢サンプル中の eDNA 濃度を定量する｣ について主に解説します。

定量 PCR と eDNA メタバーコーディング: なぜ通常の eDNA メタバーコーディングは定量的でないのか

さて、eDNA を定量したいとして、対象種が 1 種であれば種特異的なプライマーをデザインして、定量 PCR を行えば問題は解決します。プライマーがきちんとデザインできて、検量線がきちんと引けていれば、PCR 阻害などの問題は起こるかもしれませんが、感度の高い定量が可能でしょう。

一方で、調査対象種が複数いる場合、または調査する系に何種類の生物 (例えば魚) がいるかは分からないが、とにかく網羅的にそれらを検出して解析したい場合は Illumina MiSeq などの High-throughput sequencer を用いた eDNA メタバーコーディングを行うことになります。

なお、定量 PCR と eDNA メタバーコーディングは eDNA 分析で用いられる最も代表的な手法でそれぞれメリット･デメリットがあります (図 1 参照)。

図 1. 定量 PCR と eDNA メタバーコーディング

eDNA メタバーコーディングを行えば、サンプルから網羅的に調査対象種を検出することが可能ですが、定量性は失われます (“xx copies/ml water” などという値で結果は出ない)。なぜでしょうか？

eDNA メタバーコーディングでは配列解析をした後、｢種 A 由来 (と思われる) 配列があるサンプルから XX 本｣ というデータが得られます。この ｢配列 XX 本｣ というのは、最終的にシーケンサーに投入したあるサンプル中の種 A 由来の DNA 濃度に概ね比例します。しかし、eDNA メタバーコーディングでは、一般的にはシーケンサーに投入するまでに、増幅･精製･切り出し･濃度調整などのサンプル処理プロセスが必要です (｢ライブラリ調整｣ などと呼ばれます)。これらのプロセスの内、特に定量性を損なう原因となるのは増幅と濃度調整のプロセスです。 (ただし、実験手法にもよりますし、異論もあるとは思います。個人的見解が入っています。)

増幅のプロセスでは解析のターゲットとなるゲノム領域が増幅されますが、環境サンプルでは、サンプルごとに PCR の効率がかなりばらつきます。抽出 DNA に含まれる PCR 阻害物質の量がサンプルごとに違うことが要因の一つと考えています (が、他にも理由があるかもしれません)。これにより、同じ量の DNA が各サンプルに入っていたとしても、増幅後の PCR 産物の濃度はサンプルごとに変わってしまいます。

また、eDNA メタバーコーディングでは、しばしば各サンプルから得られる配列数を同程度にするために、各サンプルの増幅･精製後の PCR 産物の濃度をあえて均一に揃える、という操作をすることがあります。これはあえて定量性を無くしているとも言えますが、各サンプルからの取得配列を同程度にすることで、各サンプルに対する ｢調査努力量｣ を揃える効果があります。

その他、PCR 産物の精製に AMPure を用いたビーズ精製を利用する場合は、PCR 産物の濃度が濃すぎる場合にビーズに吸着しきらない DNA が出てくるため、このプロセスも定量性に影響がありそうです。これら様々なプロセスを経て、｢最終的にシーケンサーに投入したあるサンプル中の種 A 由来の DNA 濃度｣ が決まります。そのため、基本的には｢あるサンプルから得られる種 A 由来の配列 XX 本｣ ≒ ｢最終的にシーケンサーに投入したあるサンプル中の種 A 由来の DNA 濃度｣ は ｢抽出 DNA に含まれていたあるサンプル中の種 A 由来の DNA 濃度｣ を反映しません。

定量的 eDNA メタバーコーディング: 方法の概要

eDNA メタバーコーディングは網羅的に種を検出できる、という点で非常に強力ですが、上記で述べたように定量性があまりありません。eDNA メタバーコーディングに定量性を付加することができれば、定量性と網羅性を同時に達成する eDNA 分析ができることになります。定量性を持った多種の eDNA モニタリングデータは適用できる解析の幅も広がり、定量性を持たない eDNA データに比べてより多くの情報を提供してくれます。

そこで、Ushio et al. (2018) Metabarcoding & Metagenomics では定量的 eDNA メタバーコーディングの手法を開発し、その性能を検証しました。

方法

Ushio et al. (2018) で採用した手法は極めてシンプルで、｢PCR による増幅の前にサンプルに濃度既知の標準 DNA を添加する｣ というものです (図 2)。この方法自体は完全に新規なものというわけではなく、類似の方法は RNA-seq や土壌微生物の研究で使用されていました (例. Smets et al. 2016; 魚の eDNA に適用したのはおそらく Ushio et al. 2018 が初めて)。

この手法では、標準 DNA を複数種類 (3-5 種類) 用意し、濃度を調整して ｢PCR 前｣ にサンプル (実際には PCR の反応溶液) に添加し、その後通常のライブラリ調整を行います。図 2 では、赤･青･緑の 3 つの標準 DNA を濃度 1, 2, 3 で PCR 前に加え、ライブラリ調整をしています。

図 2. 標準 DNA を用いた定量的 eDNA メタバーコーディング

途中、増幅や精製を経ても、最終的に ｢投入した標準 DNA の濃度｣ と ｢得られた標準 DNA の配列数｣ に線形の対応が見られれば (この線を検量線、と呼んでいます)、ライブラリ調整やシーケンスのプロセスをすべて含めた結果、｢サンプル中の単位 DNA 量あたり何本の配列が取得されるか｣ が分かることになり、標準 DNA 以外にもともとサンプル中に含まれていた DNA の濃度も推定できます。

この手法では、実際に ｢投入した標準 DNA の濃度｣ と ｢得られた標準 DNA の配列数｣ の間に線形な関係が見られるかどうかが問題になりますが、これまでに実施したほぼすべてのケースで、良好な直線関係が得られています。

例えば、下の図 3 は Ushio (2020) bioRxiv で公開中の査読前プレプリントの Fig. S3A, B です。

図 3. 標準 DNA の濃度と取得配列数の直線的関係 (Ushio 2020 の Fig. S3A, B).

Ushio (2020) では、人工水田から採取した 610 の水サンプルから DNA を抽出し、標準 DNA を添加した定量 eDNA メタバーコーディングを実施しました。標準 DNA には人工合成 DNA を利用しました。610 サンプルの 4 つの DNA 領域 (16S rRNA, 18S rRNA, COI, ITS) を分析しており、すなわち合計 2,440 サンプルのシーケンスが行われたことになります。

図 3 A (左のパネル) は 16S rRNA を分析した結果で、横軸に投入した標準 DNA の濃度、縦軸にシーケンス後に得られた配列数をプロットしています。3 本の直線が引かれていますが、標準DNAの配列数 = 傾き × 標準DNAの濃度 として回帰した時の回帰係数 (= 傾き) が最も高かったサンプル、中間のサンプル、最も低かったサンプル、の結果を例として示しています。この図 3 A からは２つのことがわかります。１つ目は、投入した標準 DNA の濃度と取得配列数の間には非常に綺麗な直線関係があるということです。つまり、DNA 1 copy あたり取得されうる配列数がサンプルごとに推定可能ということです。２つ目は、サンプルによって直線の傾きが大きく違う、ということです。つまり、あるサンプルでは 50,000 copies/µl の DNA を投入すると 5,000 配列程度の配列が取得されるのに対して、違うサンプルでは 1,000 配列程度しか取得されない、ということが起こっています。これを見ても、｢配列数そのものは DNA 濃度の指標とならない｣ ことがわかります。

図 3 B はこの直線の R2R^2 (決定係数) の分布を解析対象した 4 つの領域ごとに示しています。概ね、ほとんどのサンプル･領域で R2>0.9R^2 > 0.9 です。つまり、ほとんどのサンプル･領域について、投入した標準 DNA の濃度と取得配列数の間には綺麗な直線関係があるということがわかります。

定量的 eDNA メタバーコーディング: 利点･不利点

標準 DNA を用いた定量的 eDNA メタバーコーディングは多種に対して同時に eDNA 量の推定を可能にします。これは大きな利点ですが、同時にいくつかの不利点もあります。ここではこの手法の利点･不利点をリストアップしてみます。また、利点とも不利点とも言えないが、実際の運用や結果の解釈において注意深く考慮すべき点もありますので、それらもリストアップします。

ちなみに微生物の研究に関する情報がメインですが、Harrison et al. (2020) Molecular Ecology Resources “The quest for absolute abundance: The use of internal standards for DNA-based community ecology” においても、以下に上げた項目のうちいくつかと共通する内容が議論されています。

利点

  • 網羅的な種検出と同時に DNA 量の推定ができる
    言うまでもなく一番の利点は網羅的な種検出と同時に DNA 量の推定ができることです。手前味噌ではありますが、標準 DNA を用いたメタバーコーディングが他の DNA 定量手法と比べてどのような測定値を出すかについて、これまでいくつかの検討を行っています (図4, 5)。

    図 4. 定量 eDNA メタバーコーディングによる定量と他の手法との比較 (Ushio 2020 の Fig. S3C, D, E)。© 蛍光により定量した全 DNA 量 (Quant-IT を使用) と定量メタバーコーディングにより推定した全 DNA コピー数の関係。(D) 定量 PCR と定量メタバーコーディングにより定量した 16S rRNA 遺伝子の全コピー数の関係。(E) 定量 PCR で測定した 16S rRNA 遺伝子の全コピー数と 定量メタバーコーディングの際に取得した 16S rRNA 遺伝子の全配列数の関係。点線は Ⅱ 型回帰による回帰直線。点の色はサンプルを取得した水田プロット。

    図 4C では、蛍光法で定量した全 DNA 量と定量メタバーコーディングで定量した全 DNA コピー数の関係を示しています (手法の細かい部分は Ushio 2020 bioRxiv の Supplementary Information を参照)。軸は両対数ですが、それにしてもきれいな対応関係が見られます。図 4D は 16S rRNA 遺伝子のコピー数を定量 PCR と定量メタバーコーディングで測っていますが、こちらもきちんと対応関係がみられます (ただし、1:1 のラインからは少しずれています)。一方、シーケンサーから出力される配列数と定量 PCR による 16S rRNA 遺伝子の定量結果の間にはきれいな対応関係は見られません (図 4E)。
    図 5. 定量 eDNA メタバーコーディングによるカタクチイワシの定量と種特異的な定量 PCR との比較 (Ushio 2018 の Fig. 4b)。点線は 1:1 のライン。実線は線形回帰の直線。点の色は定量メタバーコーディングの際に得られた各サンプルの検量線の傾きを示す。詳しくは Ushio et al. (2018) を参照。

    また、魚の eDNA に関しても図 5 で示すような検討を行っています。舞鶴湾のカタクチイワシ eDNA を定量メタバーコーディングと定量 PCR (TaqMan) で測定した結果ですが、きれいな対応関係が見られます。もちろんここまできれいな対応関係が見られない場合もあって、その例や考えられる理由は Ushio et al. (2018) で考察しています。とは言うものの、定量 eDNA メタバーコーディングは、特に優占的な種については既存の手法とよい対応関係が得られるような結果を出してくれるようです。

  • 異なる研究の間で結果の比較ができる
    eDNA の定量ができる、ということはその値を異なる研究間で比較できる、ということになります。これはメタ解析などを行う上ではかなり有益と考えています。一方で、配列数しか得られない場合は、配列数は実験手順や使用するシーケンサーによって大きく変わるので、各種の割合 (群集組成) に注目した (メタ) 解析しかできない、ということになります。

  • 適用できる統計解析が広がる
    eDNA の濃度が分かると割合データには適用できなかったデータ解析が適用可能になります。例えば、何度かこのブログで紹介した Empirical Dynamic Modeling (EDM) は基本的には割合データには適用できません (この記事の最後の方に書いてます)。そもそも、一連の eDNA の網羅的定量化の方法は eDNA データに EDM を適用するために開発しました。

不利点

  • 標準 DNA を準備しなければならない
    定量 eDNA メタバーコーディングを行うためには標準 DNA を準備しなくてはなりません。標準 DNA は ｢分析予定のサンプルには標準 DNA の配列が含まれていない｣ という条件を満たさなければなりません。舞鶴湾の魚類 eDNA を分析した Ushio et al. (2018) では東南アジア･アフリカの淡水魚の抽出 DNA を標準 DNA として利用しました (舞鶴湾には決して居ないはずの魚)。Ushio (2020) bioRxiv では人工合成 DNA を利用しています。人工合成 DNA は (i) 使用予定のプライマーと完全一致するプライマー結合部位を持ち、(ii) インサート領域には解析対象の生物群類群 (例えば、原核生物) と同じ保存領域を持ち、(iii) しかし可変領域は現在登録されているどの解析対象生物とも一致しない、という 3 つの特徴を持つようにデザインしました。

  • 投入する標準 DNA の濃度を決めなければならない
    定量 eDNA メタバーコーディングを実行する上でこれが一番厄介だと感じています。つまり、投入する標準 DNA の濃度が薄すぎると、標準 DNA の配列が検出されないサンプルが出てきてしまい、そのようなサンプルは検量線が引けず、取得した配列数を DNA 濃度に変換できないことになります。一方、投入する標準 DNA の濃度が濃すぎると、検出される配列のほとんどが標準 DNA になってしまい、肝心のサンプル中にもともと含まれていた生物の DNA が検出されない、ということになってしまいます。投入する標準 DNA の濃度は以下のような方法で見積ることができます。

    1. 定量 PCR で解析予定サンプルの解析対象領域を定量する。
    2. 解析予定サンプルと標準 DNA の解析予定領域をメタバーコーディングのライブラリ調整と同様に PCR で増幅し、増幅産物の濃度を定量する。
    3. 少数の代表サンプルについて、実際にいくつかの濃度で標準 DNA を混ぜて、ライブラリ調整をしてシーケンスまで行う。MiSeq であれば MiSeq V2 試薬キット Nano/Micro といったシーケンスのスケールが小さく価格も安い試薬キットがあるので、それらを使って予備シーケンス。

    時間とお金に余裕がある場合は ｢3｣ の選択肢が一番確実です。Ushio (2020) bioRxiv では一度のシーケンスで 700 近くのサンプルをシーケンスする必要があったので、MiSeq V2 試薬キット Nano を使って 50 サンプルほどの代表サンプルを予備シーケンスして、投入する標準 DNA の濃度を決めていました。

  • 出力される配列の少なくない割合が標準 DNA に "喰われる"
    上記の問題と関連しますが、適切に投入する標準 DNA の濃度を決めたとしても、それなりの割合の配列が標準 DNA に割り振られます (いつも “喰われる” と言ってます)。これまでのシーケンスでは、だいたい 30% 程度の配列が標準 DNA に喰われていますが、ちょっともったいない感じがするのは否めません。経験上、検量線はうまく引ける場合がほとんどなので、Ushio et al. (2018) や Ushio (2020) で使用した 5 点で検量線を書く方法ではなく、3 点で検量線を描く方法でもいいのかな、と思い始めています。実際、Tsuji et al. (2020) Molecular Ecology Resources では 3 点の検量線を用いてアユの種内多型を定量的に解析していますが、きれいな結果を得ることができています。

  • 実験室が汚染される恐れがある
    定量 eDNA メタバーコーディングは標準 DNA として同じ配列を何度も大量に増幅します。結果として、実験操作のミスや廃棄物の処理を適切に行なっていないと、実験室環境を同じ配列の増幅産物で汚染してしまうことになりかねません (しばしばその危険性を指摘されます)。現在のところ、自分の実験環境では深刻な汚染が起こったことはありませんが、これから導入を考えている場合は注意が必要でしょう。

考慮点

  • DNA 抽出効率は考慮されない
    ここで紹介した方法は、DNA を抽出した後に標準 DNA を添加するので、標準 DNA 添加以前の過程は考慮できません。例えば、サンプルによって DNA 抽出効率に差がある場合でも、｢各サンプル間で DNA 抽出効率に差はない｣という仮定をおいて解析を進めていることと同じになります。環境サンプルの場合はサンプル間で抽出効率が多少なりとも違うということは起こりそうです。土壌微生物を解析した Smets et al. (2016) では DNA 抽出を始める前に土壌に検量線を引くのに使う微生物 (土壌にはまずいないであろう微生物) を添加しており、サンプルごとの DNA の違いを考慮したものとなっています。ただし、微生物種間での DNA 抽出効率の違いはこの方法でも考慮できません。

  • プライマーの match/mismatch に起因する生物種間での増幅効率の差は考慮されない
    標準 DNA はしばしばプライマー領域が使用予定のプライマーと完全一致するようにデザインされています。しかし、実際の解析対象生物の中にはプライマーとミスマッチがあるプライマー領域を持つものもいるでしょう。そのような場合、PCR における増幅効率が標準 DNA と異なることが想定されますが、このような違いも今回紹介した方法では考慮できません。研究上重要な特定の解析対象生物種がいる場合は、あえてプライマー領域にミスマッチのある標準 DNA をデザインして検量線を引くのもいいかもしれません。

標準 DNA を用いる以外の定量的メタバーコーディングの方法

今回紹介した標準 DNA を用いる以外にも定量的･網羅的な eDNA データを得る方法があります。ここでは２つ、概要だけ紹介します。

Unique Molecular Identifier (UMI) を利用する方法

UMI はランダムな塩基を複数 (例えば、12個) つなげてつくる ｢分子特異的なタグ｣ です。DNA の塩基は 4 種類なので、例えば 12 個のランダム塩基 (N で表す) をつなげるとすると、NNNNNNNNNNNN となりますが、その中身は 412=16,777,2164^{12} = 16,777,216 種類の塩基配列となります。この UMI とプライマー配列 (例えば MiFish-F) をくっつけると NNNNNNNNNNNNGTCGGTAAAACTCGTGCCAGC となります。1677万コピー以下しか MiFish 領域が存在しないサンプルに対して、このようなプライマーを用いて PCR を 1 回分だけ行うと、理想的には MiFish プライマーと 1677万種類のうちどれか一つの ｢タグ｣ がついた産物ができます (確率的な挙動を考えると間違っていますが、話を単純にしています)。このタグ配列が以後置き換わらないようにライブラリ調整を行い、シーケンスをします。すると、理想的には NNNNNNNNNNNN はかぶることなくテンプレートの DNA と結合するので、シーケンス後に NNNNNNNNNNNN の種類数を数えればそれがそのままテンプレート DNA 中の解析対象 DNA の濃度になります。

この技術も RNA-seq などで頻繁に使われ、細菌の 16S rRNA 解析にも使用されています (Hoshino & Inagaki 2017 PLoS ONE)。また、最近、水槽水の魚の eDNA 定量にも応用されました (Hoshino et al. 2021 Scientific Reports)。

この方法は標準 DNA を必要としないため、標準 DNA の用意･投入する標準 DNA の濃度の決定･配列が標準 DNA に喰われる、などの問題を回避できます。一方で、1 回目の反応でタグ付きのプライマーが狙った DNA と結合するかどうかが鍵のようで、解析対象の eDNA の濃度が低い場合 (= 魚類などのマクロ生物を野外の水で検出する場合など) について、どれほどの感度で定量が可能なのか、慎重に検討する必要があります。また、アンプリコンシーケンスにおいて UMI を効率的に処理するパイプラインもまだあまり充実していないようです。探した限りでは Clements et al. (2018) Bioinformatics による AmpUMI というツールが使えるかもしれません (RNA-seq の場合は他にも色々ツールがあるようです)。

他の手法で補正する方法

もう一つの手法は、通常の eDNA メタバーコーディングを行った後に他の手法で解析対象の生物全体の絶対量を測定し、その絶対量に eDNA メタバーコーディングで分かる割合データ (相対優占度) を掛けて、各生物種の絶対量とする、というやり方です。例えば、細菌が解析対象であれば、定量 PCR による 16S rRNA 全コピー数の定量や顕微鏡観察による細菌の全細胞数の計測などによって細菌の全体量を推定し、そこに eDNA メタバーコーディングによる割合データを掛ければよいでしょう。また、通常の eDNA メタバーコーディングを行った後、(ほぼ) すべてのサンプルから出てくる生物種について、後から定量 PCR で eDNA の絶対量を測る、というやり方も考えられます。例えば、魚であればクロダイやアジ、フグの仲間は日本の多くの沿岸域で見られると思いますが、それらの eDNA 量を定量すれば、通常の eDNA メタバーコーディングのデータと合わせて、サンプルごとに単位 eDNA 量あたり得られる配列数の推定が可能になり、その他の魚種についても eDNA の絶対量が推定できます。ただし、この定量 PCR での事後補正は、サンプルごとに異なるかもしれない PCR 阻害の効果の補正はできません。また、ほとんどのサンプルから出てくる生物種、というのが見つからない場合は使用できません。一方で標準 DNA を用意する必要がないため、上記で紹介した標準 DNA 関連の不利点は回避できます。

まとめ

さて、いろいろと紹介させてもらいましたが、現時点では標準 DNA を用いた eDNA メタバーコーディングは、網羅的かつ定量的に eDNA を分析したい場合、有力な候補となる方法です。特に、現状 UMI での正確な定量が難しそうな低濃度なマクロ生物の eDNA を分析する際は有効でしょう。一方で、標準 DNA の準備や投入する濃度の決定などが手間であることも事実です。これらの点をなんとか克服して、定量的･網羅的な eDNA 分析が簡便･迅速に実行できるようにする、というのが当面の課題です。

また、環境サンプル中の eDNA のより長い配列がより正確に分析できて、eDNA から ｢個体識別｣ ができるようになれば、個体数を数えることが可能になり、｢eDNA から調査対象の種の個体数やバイオマスを推定する｣ ことが一気に現実的になってきます。こちらもずいぶんチャレンジングですが、このようなことができれば、eDNA 分析が提供できる情報の価値がさらに増すでしょう。

参考文献

Written with StackEdit.

時刻: 0 件のコメント:

2021年4月28日水曜日

phyloseq による DADA2 処理後の統計解析

20210428_Blogger0015

以前のブログで DADA2 と Claident を用いたアンプリコンシーケンスのデータ解析について触れました。今回は DADA2 と Claident を用いて「分類群 × サンプルのテーブル」ができた後の解析について解説します。DADA2 は Amplicon Sequence Variant (ASV) を生成するので、ここではこの「分類群 × サンプルのテーブル」を ASV テーブルと呼ぶことにします。分類群 = Operational Taxonomic Unit (OTU) の場合は OTU テーブルと呼びます。

DADA2 後の統計解析

DADA2 で ASV テーブルを作成した後はいくつかお決まりの解析が待っています。ここでは以下の 5 点について、R と phyloseq パッケージを利用した実装とともに解説します。

  • 群集組成バターンを棒グラフで可視化
  • 種多様性パターンの可視化
  • Rarefaction カーブの作成
  • 次元圧縮によるパターンの可視化
  • ASV の再クラスタリング

phyloseq は McMurdie & Holmes (2013) PLoS ONE で発表されたアンプリコンシーケンスのデータの統計解析を行うための R のパッケージです。2021年4月時点で5000回以上引用されているようです。公式チュートリアルはこちら。

デモデータ

今回も Ushio (2019) Methods in Ecology and Evolution のデータを使います。全配列データは DRA (こことここ) に、全解析コードやメタデータは Zenodo (ここ) にアーカイブしてありますが、今回のデモでは一部のデータだけ使用します。

研究は DNA 抽出法の検討をするために行いました。複数のサイトから環境水をサンプリングし、16S メタバーコーディングを行っています。DNA 抽出法は 3 種類試しました。デモ用のデータは、２サイト (Site = “sea”, “pond”)、３抽出法 (Method = “w_beads”, “wo_beads”, “destruction”)、５反復 + ネガコン のデータを含みます。

DADA2 や Claident の部分は飛ばして、ASV table からスタートします。Github にデモデータやスクリプトを置いたので、自分で再現したい方はダウンロードしてください (https://github.com/ong8181/random-scripts/tree/master/05_DADA2Phyloseq)。

準備

まずは以下のコードで必要なパッケージをロードします。

# Load library and functions
library(tidyverse); packageVersion("tidyverse") # 1.3.0
library(phyloseq); packageVersion("phyloseq") # 1.34.0
library(cowplot); packageVersion("cowplot") # 1.1.1
library(ggsci); packageVersion("ggsci") # 2.9
theme_set(theme_cowplot())

tidyverse はデータ整形を効率的に行うため、cowplot と ggsci は可視化のための ggplot2 パッケージの補助パッケージです。ggplot2 は tidyverse をロードするときに同時にロードされています。theme_set(theme_cowplot()) で図のテーマを指定しています (これはただの好みです)。

パッケージのロードができたら以下のコードでサンプルメタデータ、ASV テーブル、 分類群情報を読み込みます。分類群情報は Claident で推定したものです。

sample_sheet <- read.csv("data_subset/sample_sheet.csv", row.names = 1)
asv_sheet <- read.csv("data_subset/asv_sheet.csv", row.names = 1)
tax_sheet <- read.csv("data_subset/tax_sheet.csv", row.names = 1)

phyloseq へのデータの読み込み

デモデータの読み込みができたら、まずはデータがきちんと phyloseq オブジェクトに変換できる形になっているかどうかを確認します。

dim(sample_sheet); dim(asv_sheet); dim(tax_sheet)
all(rownames(asv_sheet) == rownames(sample_sheet))
all(colnames(asv_sheet) == rownames(tax_sheet))

以下の条件を満たしている必要があります:

  • asv_sheet の行数 = サンプル数 = sample_sheet の行数
  • asv_sheet の行名 = サンプル ID = sample_sheet の行名
  • asv_sheet の列数 = 分類群数 = tax_sheet の行数
  • asv_sheet の列名 = 分類群 ID = tax_sheet の行名

今回は ASV テーブルの行 = サンプル、列 = 分類群、となっていますが、行と列を逆にしたテーブルも読み込み可能です。データが phyloseq オブジェクトに変換できる形になっていることを確認できたら、以下のコードで phyloseq オブジェクトを生成します。

ps_all <- phyloseq(otu_table(asv_sheet, taxa_are_rows = FALSE),
                   sample_data(sample_sheet),
                   tax_table(as.matrix(tax_sheet)))

ASV テーブルの行 = 分類群、列 = サンプル、となっている場合は taxa_are_rows = TRUE としてください。tax_sheet は as.matrix() をかまさないと読み込めないのですが、理由は分かりません。また、系統樹情報も同時に読み込むことができますが、今回は利用しません。

phyloseq にはいくつものデータハンドリングのための関数が用意されています。

  • サンプルを絞る
## Select "sea" samples
subset_samples(ps_all, Site == "sea")

## Select samples with total reads > 10^4
prune_samples(sample_sums(ps_all) > 10^4, ps_all)
  • 分類群を絞る･優占種を選ぶ
## Select "Proteobacteria" phylum
subset_taxa(ps_all, phylum == "Proteobacteria")

## Select top 50 taxa
top_taxa <- names(sort(taxa_sums(ps_all), decreasing = TRUE)[1:50])
prune_taxa(top_taxa, ps_all)

subset_xxx() と prune_xxx() が似ていますが、選びたいサンプル名や分類群名のリストがある場合は prune_xxx() が便利です。例えば、select_samples <- c("R001", "R002", "R003") とある場合は prune_samples(select_samples, ps_all) として選べます。ps_all を入れるべき位置が subset_xxx() と prune_xxx() で違う理由は分かりません。

  • 配列の値を変換
## Transform sample counts
transform_sample_counts(ps_all, function(x) 100* x / sum(x))
transform_sample_counts(ps_all, log)

上では配列を相対値に変換、log を取る、という操作をしています。

全体のパターンを棒グラフで可視化

群集組成の可視化

plot_bar() で全体のパターンを確認できます。また、 + を続けることで ggplot2 の機能が利用できます。scale_fill_igv() は ggsci パッケージ内のカラーパレットです。

# plot_bar
plot_bar(ps_all, x = "Sample", fill = "phylum") + scale_fill_igv()

facet_wrap() を使って、サイトと手法ごとに分けて結果を表示。

plot_bar(ps_all, x = "replicate", fill = "phylum") +
  facet_wrap(.~ Site + Method) + scale_fill_igv()

種多様性パターンの可視化

plot_richness() で多様性パターンを可視化できます。以下ではシンプルに検出 ASV の数をプロットしています。

# plot_richness
plot_richness(ps_all, measures = "Observed")

measures の引数を変更するとシンプソンやシャノンの多様度指数も図示することができます。

psmelt() と ggplot() による可視化

psmelt() を使うと元データをバラすことができ、より自由度の高い可視化が可能になります。例えば、以下は最も優占していた分類群の配列数を箱ひげ図 + ジッタープロットで表示したものです。

# Visualize patterns (scatterplot)
ps_m1 <- psmelt(ps_all) %>%
  select(OTU, Sample, Abundance, Site, Method) %>%
  filter(OTU == top_taxa[1], Site == "pond")
ggplot(ps_m1, aes(y = Abundance, x = Method)) +
  geom_boxplot(outlier.shape = NA, width = 0.2) +
  geom_jitter(height = 0, width = 0.2)

Rarefaction カーブ

｢自分の配列データが多様性をどの程度を捉えているのか｣ を配列数が増える毎にどのくらい検出 ASV の数が増加していくのかで把握します。この「検出 ASV 数 v.s. 配列数」の曲線は rarefaction カーブと呼ばれ、微生物群集関連の研究では頻繁に登場します。以下では、ここの記事を参考に rarefaction カーブを ggplot で書いています。

まず、群集生態学関連のパッケージ vegan でデータを rarefy しています。ちなみにこの時点でも図が出力されます。これ以後のコードは単に ggplot で描くためだけのものです。

## From 
# Rarefy samples
ps_rr <- vegan::rarecurve(as.matrix(otu_table(ps_all)), step = 50, label = TRUE)

rarecurve で得た出力 (ps_rr) を ggplot 向けに整形します。

rare <- lapply(ps_rr, function(x){
  b <- as.data.frame(x)
  b <- data.frame(ASV = b[,1], raw_read = rownames(b))
  b$raw_read <- as.numeric(gsub("N", "",  b$raw_read))
  return(b)
})
# Label list
names(rare) <- rownames(sample_data(ps_all))
# Convert to the data frame
rare_df <- map_dfr(rare, function(x) return(data.frame(x)), .id = "sample")

最後に ggplot で図示します。概ね 5,000 – 10,000 配列以上取得しても出現する ASV の数に変化はないようです。ただし、DADA2 による解析結果は原理上シングルトンが出現しないので解釈には注意が必要です。

# Visualize with ggplot
ggplot(rare_df, aes(x = raw_read, y = ASV, color = sample)) +
  geom_line() + scale_color_igv() +
  xlab("Reads") + ylab("The number of ASV")

次元圧縮によるパターンの可視化

次元圧縮 (次元削減) による高次元データの可視化は全体のパターンを把握する上で非常に有効です。以下では、NMDS、t-SNE、UMAP を適用して微生物群集データを 2 次元で可視化しています (各手法の原理の詳細はここでは解説しません)。

ここではデータのうち海から取得した配列データ (Site == "sea" のサンプル) に対して次元圧縮を適用してみます。また、ネガコンを除いたサンプルのみ (sample_nc == "sample") を解析します。

ps_sea <- subset_samples(ps_all, Site == "sea" & sample_nc == "sample")

NMDS

Nonmetric Multidimensional Scaling (NMDS; 非計量多次元尺度法) は次元圧縮の古典的な手法です。phyloseq では ordinate() と plot_ordination() で解析を実行できます。ランダムな過程を含むので、再現性確保のために乱数のシードを設定しています。

# Nonmetric multidimensional scaling
set.seed(1234)
ps_bray <- ordinate(ps_sea, "NMDS", "bray")
plot_ordination(ps_sea, ps_bray, color = "Method", shape = "Method") +
  stat_ellipse(geom = "polygon", alpha = 0.1, aes(fill=Method)) +
  geom_point(size = 2) + scale_color_startrek() + scale_fill_startrek() +
  xlab("Axis 1") + ylab("Axis 2") + ggtitle("NMDS")

上のコードでは stat_ellipse() を使って各グループの信頼区間 (範囲?) も表示しています。scale_xxxx_startrek() は ggsci パッケージの関数で、カラーパレットを指定しています。主成分分析 (PCA) なども ordinate() 関数を用いた同じコードで実行できます。

t-SNE

t-distributed Stochastic Neighbor Embedding (t-SNE) は確率分布に基づいて次元圧縮する手法で 2008年に提案されています (Maaten & Hinton 2008)。アンプリコンシーケンスの微生物群集データよりは RNA-seq データの解析などでよく見る印象です。phyloseq そのものには t-SNE は実装されていませんが、tsnemicroboita パッケージを使うことで phyloseq オブジェクトのまま t-SNE を実行できます (単に t-SNE を実行するだけのパッケージは他にも色々あります)。

# t-SNE
#library(tsnemicrobiota)
tsne_res <- tsnemicrobiota::tsne_phyloseq(ps_sea, distance='bray', perplexity = 5, rng_seed = 1234)
tsnemicrobiota::plot_tsne_phyloseq(ps_sea, tsne_res, color = "Method", shape = "Method") +
  stat_ellipse(geom = "polygon", alpha = 0.1, aes(fill=Method)) +
  geom_point(size = 2) + scale_color_startrek() + scale_fill_startrek() +
  xlab("Axis 1") + ylab("Axis 2") + ggtitle("t-SNE")

UMAP

Uniform Manifold Approximation (UMAP) は 2018 年に提案された手法です (Mcinnes et al. 2018 arXiv)。t-SNE よりもデータをきれいに分けることができ、かつ高速とされています。ただし、今回のデモデータでは分類群の数もサンプルの数も非常に少ないので、その効果は実感できません。アンプリコンシーケンスの微生物群集データに適用されているのはほとんど見たことがありません (自分が知らないだけかもですが)。phyloseq オブジェクトをそのまま扱って UMAP を実行するパッケージは見つけられませんでした。

# UMAP by uwot package
#library(uwot)
set.seed(1234)
umap_res <- uwot::tumap(otu_table(ps_sea)@.Data, n_neighbors = 5, n_components = 2)
umap_df <- cbind(data.frame(sample_data(ps_sea)), umap_res)
colnames(umap_df)[(ncol(umap_df)-1):ncol(umap_df)] <- c("UMAP1", "UMAP2")
ggplot(umap_df, aes(x = UMAP1, y = UMAP2, color = Method, shape = Method)) +
  stat_ellipse(geom = "polygon", alpha = 0.1, aes(fill=Method)) +
  geom_point(size = 2) + scale_color_startrek() + scale_fill_startrek() +
  xlab("Axis 1") + ylab("Axis 2") + ggtitle("UMAP")

ざっと並べてみると、サンプル数や ASV 数がそれほど多くないのであれば NMDS で十分な印象を受けます。データサイズが大きくなれば (RNA-seq のデータとか)、t-SNE や UMAP が良いかもしれません。いずれにしてもどちらも試してみることをおすすめします。

ASV の再クラスタリング

ASV はアンプリコン同士の違いを 1 塩基レベルで識別する非常に強力な方法ですが、種内多型のような変異も拾うため、多様性 (= ASV の数) が増大し、ときに見たいパターンが見えなくなるときがあります。そのような場合、DADA2 でエラー除去した ASV を OTU として再クラスタリングすることが有効な場合があります。

以下では、R を用いて DADA2 の出力を OTU に再クラスタリングする方法を紹介します。以下のコードは DADA2 の issues のページで議論されていたものを参考にしています。Biostrings と DECIPHER (http://www2.decipher.codes/) というパッケージを使用しています。また、R 内で解析を完結されることにこだわらないのであれば vsearch などのツールを使っても同様のことができるはずです。

クラスタリングプロセスは非常にシンプルです。以下では 97% 一致で OTU を作成しています。

# Preparation
n_cores <- parallel::detectCores() # set to number of cpus/processors to use for the clustering
dna <- Biostrings::DNAStringSet(tax_sheet$seq)
# Find clusters of ASVs to form the new OTUs
aln <- DECIPHER::AlignSeqs(dna, processors = n_cores)
aln_dist <- DECIPHER::DistanceMatrix(aln, processors = n_cores)
# use `cutoff = 0.03` for a 97% OTU 
clusters <- DECIPHER::IdClusters(aln_dist, method = "complete",
                                 processors = n_cores, type = "clusters",
                                 cutoff = 0.03, showPlot = F)

クラスタリングができたらクラスター情報を利用して ASV テーブルを変換します。

colnames(clusters) <- "cluster_id"
clusters$cluster_id <- factor(clusters$cluster_id, levels = unique(clusters$cluster_id))
length(unique(clusters$cluster_id))
# Use dplyr to merge the columns of the asv_sheet matrix for ASVs in the same OTU
merged_asv_sheet <- asv_sheet %>% t %>%
  rowsum(clusters$cluster_id) %>% t

また、phyloseq オブジェクトにクラスタリング結果を適用することもできます。speedyseq パッケージ内の merge_taxa_vec() という関数を使用しています。

# Make OTU-based phyloseq objects
# Import DADA2 ASV output to phyloseq
asv_sheet2 <- asv_sheet; colnames(asv_sheet2) <- tax_sheet$seq
ps_tmp <- phyloseq(otu_table(asv_sheet2, taxa_are_rows = FALSE))
# Merge taxa in "ps" using cluster information
ps_otu <- speedyseq::merge_taxa_vec(ps_tmp, group = clusters$cluster_id)

以下は解析そのものとは関係ないですが、ASV と OTU の対応を確認するためのサマリーを作成しています。

# Get OTU sequences and OTU table
otu_seqs <- colnames(otu_table(ps_otu))
otu_sheet <- otu_table(ps_otu)@.Data

# Sanity check (make clustering summary data frame)
otu_only <- data.frame(taxa_id = sprintf("OTU%05d", 1:length(otu_seqs)),
                       seq = otu_seqs)
clusters$seq_sum <- colSums(asv_sheet)
clusters$asv_seq <- tax_sheet$seq
asv_1st <- clusters %>% group_by(cluster_id) %>% summarize(otu_seqs = asv_seq[[1]])
clusters$otu_seq <- unlist(asv_1st[match(clusters$cluster_id, asv_1st$cluster_id), "otu_seqs"])

お役立ちパッケージ

その他、phyloseq にはいろいろな拡張パッケージがあります。一部を紹介します。

  • speedyseq (https://github.com/mikemc/speedyseq)
    phyloseq の ｢ちょっと遅いなぁ｣ と感じる関数を高速化したパッケージです。以下では psmelt() 関数を試しています。データサイズが大きいとかなり威力を発揮します。
#library(speedyseq)
system.time(psmelt(ps_all))             # 0.103 sec elapsed
system.time(speedyseq::psmelt(ps_all))  # 0.025 sec elapsed

#library(fantaxtic)
ps_top20 <- fantaxtic::get_top_taxa(ps_all, n = 20, relative = TRUE,
                                    discard_other = FALSE, other_label = "Other")
fantaxtic::fantaxtic_bar(ps_top20, color_by = "phylum",
                         label_by = "phylum", other_label = "Other", order_alg = "alph")

まとめ

微生物群集のアンプリコンシーケンス解析では、サンプルメタデータ･配列データ･ASV (or OTU) テーブル･分類群データ･系統樹など、複数のデータを組み合わせて解析することが求められます。phyloseq はそのようなデータを整然と扱えるようにデザインされていますので、うまく使って再現性のある解析を進めていただければと思います。

References

Written with StackEdit.

時刻: 0 件のコメント:
登録: 投稿 (Atom)

Empirical Dynamic Modeling with rEDM: (2) Near-future forecasting (Simplex projection)

20231015_Blogger0021 *This is an English version of my previous post (first translated by ChatGPT, checked and edited ...