*2021.5.25 追記: Claident のマニュアルの置き場所が変わったため、記事内の Claident マニュアルへのリンクは機能していません。マニュアルの新しい置き場所はこちら。metabarcodingtextbook2.ja の pdf 版もしくは html 版が見やすいかなと思います。
環境微生物群集の解析方法
環境サンプル中の微生物群集の分析方法には多くの種類があります。例えば、単離培養・顕微鏡観察・核酸情報に基づいた解析などです。この中でも核酸情報に基づいた分析は近年の核酸シーケンス技術の急速な発展に伴って非常にポピュラーなものになっています。
核酸情報に基づいた分析にも多くの種類があります。例えば、2018年に発表された論文では微生物嚢 (microbiome; 要するに微生物群集) 研究の方法が概観されています (Knight et al. 2018 Nature Review Microbiology)。その中では ｢Marker gene｣ ｢Metagenomics｣ ｢Metatranscritomics｣ の 3 つが紹介されており、Table 1 でそれぞれの方法の長所短所が列挙されています。この 3 つの方法の中で、｢Marker gene (analysis)｣ は “Quick, simple and inexpensive sapmle preparation and analysis” とされていて、要するに環境微生物群集の特徴を網羅的に、さほどコストをかけずに知りたいときには有効な方法です。
この記事では、比較的「手軽に･安く」微生物群集を解析する方法である Marker gene analysis について、特にデータ解析に焦点を当てて解説します。配列データを得るまでの過程、つまりサンプル取得・保存・DNA抽出・ライブラリ調整・シーケンス、も重要かつ手間のかかる作業ですが、それについては別の機会に解説したいと思います。
｢Metagenomics｣｢Metatranscriptomics｣ に関しては今の所それほど自分自身の経験がないこともあり、この記事では解説しないこととします。
Marker gene analysis
Marker gene analysis はゲノム中の特定の領域を PCR により増幅したのちにシーケンスを行ってサンプル中の微生物群集の組成を解析するものです。対象となる分類群によってよく利用される領域が異なりますが、例えば原核生物だと 16S rRNA の配列を増幅して解析するのが一般的です。16S rRNA のどの部分を増やすのかは研究によってまちまちですが、例えば Earth Microbiome Project では 16S rRNA の 515番目-806番目の領域が分析されています。自分も原核生物を解析するときはよくこの領域のシーケンスを行っています (例えば、Ushio 2019 Methods in Ecology and Evolution)。
16S rRNA は全長 1.5 kbp 程度で、Illumina MiSeq などのシーケンサーを利用するときはその性能上、16S rRNA の部分配列しか解読できません。通常 300 bp 程度の断片を読みます。このような短い marker 配列を解析することを short-read amplicon sequence と言ったりもします (amplicon は ｢PCR で増幅されたもの = amplify されたもの｣ という意味です)。最近ではより長い DNA 断片を読めるシーケンサーが登場しており (PacBio Sequel や Nanopore MinION)、16S rRNA 全長の解析も可能です。そのような解析を long-read amplicon sequence と言ったりもします。
以下では自分のところで行っている short-read amplicon sequencing で得られた 16S rRNA 配列のデータ解析について詳しく解説します。
配列データ解析の流れ
現在のところ、自分は以下の流れに沿って配列解析を行っています。
1. Illumina MiSeq から生データを吸い出す
2. bcl2fastq というプログラムを使って fastq ファイルを生成
3. Claident と呼ばれるパイプラインの機能を使って fastq ファイルをサンプルごとに分割 (demultiplex と呼ばれる操作)
4. DADA2 と呼ばれる R のパッケージを利用してエラー配列の除去やペアエンド配列のマージを実行し、サンプル × 配列の (しばしば巨大な) 行列データに変換
5. 再び Claident を利用して配列の分類群を推定
以下、細かい解説を述べていきます。
1. Illumina MiSeq から生データを吸い出す
MiSeq によるシーケンスが無事完了したら、まずデータをマシンから吸い出さなければなりません。Illumina BaseSpace に無料アカウントを作っておけば、シーケンス時にランをモニタリングすることができ、かつラン終了後に BaseSpace から demultiplex 済みの fastq ファイルをダウンロードできます。自分も当初はこの方法で fastq ファイルを入手していましたが、ここ何年かは行っていません。その理由については ｢3. Claident による demultiplex｣ で解説します。
現在は、シーケンスが終了したら、ハードディスクを MiSeq がおいてある部屋に持っていき、直接データをコピーしています。MiSeq のマシンの中に MiSeq_Output というフォルダが有り、内部に個別のランのデータが格納されています (BaseCall というフォルダが含まれます)。それを丸ごとハードディスクにコピーします。
もしくは、OS が Linux であれば BaseMount というプログラムを使うことで遠隔で生データをダウンロードすることも可能です。ただし、生データは結構重いので (〜 数十 GB)、マシンから直接データをコピーできるときはその方がいいように思います。
2. bcl2fastq を利用した fastq ファイルの生成
MiSeq から生データをコピーしたら、次は bcl2fastq という Illumina が提供しているプログラムを利用して fastq ファイルを生成します。以下のコマンドで v2 系の bcl2fastq が入手・インストールできます (Linux の端末上から操作してください; バージョンについては適宜変更してください)。現在のところ bcl2fastq は Linux のみで動くようです。下記は Claident のマニュアルの付録 B に記載されているものです。v1系と v2系では若干コマンドが異なるので注意が必要です。
# Install rpm2cpio
sudo apt-get install rpm2cpio cpio
# Get bcl2fastq2
wget -O bcl2fastq2-v2-18-0-12-linux-x86-64.zip http://bit.ly/2drPH3W
# Unzip bcl2fastq2
unzip -qq bcl2fastq2-v2-18-0-12-linux-x86-64.zip
# Exstracting command
rpm2cpio bcl2fastq2-v2.18.0.12-Linux-x86_64.rpm | cpio -id
sudo mv usr/local/bin/bcl2fastq /usr/local/bin/
sudo cp -R usr/local/share/css /usr/local/share/
sudo cp -R usr/local/share/xsl /usr/local/share/
bcl2fastq が入手できたら以下のコマンドで fastq ファイルを生成できます。一点注意が必要なのは、Claident のマニュアルにもある通り、SampleSheet.csv というランフォルダ以内にあるファイルの名前を変えておく必要があることです (書き換える名前は何でも良いです)。
bcl2fastq --processing-threads CPUコア数 --use-bases-mask Y250n,I8,I8,Y250n --create-fastq-for-index-reads --runfolder-dir ランフォルダ --output-dir 出力フォルダ
v1系ではコマンドが異なるので上記のコマンドでは動きません。--processing-threads オプションにはCPUのコア数を、--use-bases-mask オプションにはサイクル数・インデックス配列長を記載しています。上記は V2 500 cycle 向けのコードで、仮に V3 600 cycle の試薬キットを利用した場合は --use-bases-mask Y300n,I8,I8,Y300n となります。処理が順調にすすむと全てのサンプルのデータがひとまとめになった巨大な fastq ファイルが (上記の例では) 4 つ作成されるはずです。
3. Claident による demultiplex
bcl2fastq による fastq ファイルの生成が終了したら、Claident (Tanabe & Toju 2013 PLoS ONE) を利用して fastq ファイルの demultiplex を行います。この処理を採用する理由は Claident のマニュアルの 2.3.1 に記載されているとおりです: “Illumina社のシーケンサーでは、標準でdemultiplexする機能がありますが、タグ配列部分の品質値を一切見ていないため、タグ配列の品質が低い配列がファイルに含まれています”。
Claident で処理することによりサンプル識別用のタグ配列の品質 (Q-score) が低いものを除去して demultiplex ができます。すなわち、以後の解析でそのサンプルに由来するかどうか怪しい配列を除去して解析できるため、解析の信頼性が向上します。また、Claident による demultiplex では (配列の前半部分の) プライマー領域を同時に除去してくれます。
Claident による demultiplex の処理については本家のマニュアルの 2.3.2 から丁寧に解説されているのでそちらを参照してください。
4. DADA2 (Amplicon Sequence Variant 法) による配列解析
Claident により demultiplex ができたら次は DADA2 と呼ばれるパイプラインによる処理を行っています。
2016年頃までは amplicon sequence で得られた配列は Operational Taxonomic Unit (OTU) というユニットを用いて解析されていました。これは、例えば配列の相同性が 97% 以上である配列を操作上 ｢種｣ のようなものとみなして解析を進めていく方法です。これによって、例えば 1 塩基の読み間違いから生じてしまったような配列を誤って ｢種｣ として扱ってしまうことを防ぐことができます。しかし、一方で本当に 1 塩基しか違わない ｢種内変異｣ のような違いを扱うことはできませんでした。
そんな中 2016 年に DADA2 (Divisive Amplicon Denoising Algorithm 2) と呼ばれる手法が提案されます (Callahan et al. 2016 Nature Methods)。DADA2 ではエラーを明示的な統計モデルに基づいて予測することで、 1 塩基程度しか違いのない 2 つの配列をどちらかの配列がPCRエラーやシーケンスエラーによって生成されたものなのか、それとも 2 つとも生物学的に意味のある種内変異のようなものなのかを区別します。このように 1 塩基違いを識別するための配列解析方法は Amplicon Sequence Variant (ASV) 法、もしくは Exact Sequence Variant (ESV) 法と呼ばれており (Callahan et al. 2017 The ISME Journal)、OTU に基づいた方法に変わって微生物生態学で標準的な方法となりつつあります (もう標準になっている気もします)。また、DADA2 以外にも ASV法を実装しているアルゴリズムが複数あります (Deblur や UNOISE3)。理論的な背景は原著論文に譲るとして、以下では R と dada2 パッケージを用いたデータ解析方法について概説します。dada2 パッケージは Windows, macOS, Linux に対応しています。
DADA2 には丁寧なチュートリアルが用意されています。基本的にはそのページに書かれてある情報に従って、R の dada2 パッケージのインストールを行い、配列解析を進めることで ASV 法が実行できます。ただし、チュートリアルでは BaseSpace が出力した fastq ファイルを解析するという前提でコードが書かれているため、コードを多少変更する必要があります。以下では、変更部分のみコードを記載します。その他の部分については本家のチュートリアルを参照してください。
注意: V2 500 cycle のキットで 16S rRNA の 515F-806R 領域を読んだ場合はあまり問題になりませんが、もし真菌の ITS のような配列長が変異に富んだ領域を解析する場合は注意が必要です。すなわち、MiSeq によるシーケンスで増幅した配列を読み切って、リバースプライマーやその後ろのアダプターの配列まで読んでしまうことがあります。これらは生物由来の配列ではないので、このまま解析するとよくありません。この ｢反対側のプライマー｣ は Claident による demultiplex でも除去されないので、DADA2 ITS Pipeline に従って除去する必要があります。
変更点: フィルタリングのプロセスでは Claident による demultiplex を採用したことによるコードの改変が必要です。具体的には以下のコードを利用しています。
out <- filterAndTrim(フォワード側入力配列,
フォワード側フィルタリング配列,
リバース側入力配列,
リバース側フィルタリング配列,
#trimLeft = c(0,0),
maxN = 0, maxEE = c(2,2), truncQ = 2, rm.phix = F,
compress = TRUE, multithread = TRUE)
# On Windows set multithread=FALSE
Claident が予めプライマー配列を除去してくれているので、DADA2 のチュートリアルでは使用されている trimLeft のオプションを使用していません。(また、こちらではプライマー配列に NNNNNN を入れてクラスターの識別を改良しており、シーケンス時に PhiX を使用していないため、rm.phix = F となっています)
考慮点: ペアリードを結合する場合の最小の重なり (overlap) は 12 となっていますが、これは解析する領域のよっては変更を検討しても良いと思います (以下では minOverlap = 20 となっています)。また、trimOverhang = TRUE のオプションを指定することで ｢はみ出た配列｣ を切り取ることができるようです。
mergers <- mergePairs(dadaFs, derepFs,
dadaRs, derepRs,
trimOverhang = TRUE, minOverlap = 20,
verbose=TRUE)
変更点: DADA2 チュートリアルではこのあとに assignTaxonomy()を用いて分類群を割り当てています (dada2 による分類群割り当て)。しかし、自分は DADA2 の機能は利用せず、ここで再び Claident のお世話になっています。詳しくは Claident の原著論文 (Tanabe & Toju 2013 PLoS ONE) に譲りますが、分類群の割り当て機能は Claident が最も注意深くかつ網羅的にできるように設計していると考えているためです。
DADA2 による配列解析の結果得られた ｢代表配列｣ を R で出力し、それを Claident に投げる、という戦略を取っています。例えば、以下のようにして代表配列を出力しています。
seqs <- colnames(配列表) # 列名に配列 (ATGC...) が記載されているはずです
seqs_out <- as.matrix(c(rbind(sprintf(">Taxa%05d", 1:length(seqs)), seqs)), ncol = 1)
write.table(seqs_out, "代表配列ファイル名", col.names = FALSE, row.names = FALSE, quote = FALSE)
5. Claident による分類群の割り当て
代表配列のファイルが得られれば Linux の端末 (ターミナル) 上から以下のコードを実行して Claident による分類群の割り当てを行っています。
# 1. overall_class データベースによる同定
clmakecachedb --blastdb=overall_class --numthreads=CPUコア数 代表配列ファイル名 overall_classの出力キャッシュファイル名
clidentseq --blastdb=overall_classの出力キャッシュファイル名 --numthreads=CPUコア数 代表配列ファイル名 overall_classによる近隣配列の取得結果
classigntax --taxdb=overall_class --maxpopposer=0.05 --minsoratio=19 overall_classによる近隣配列の取得結果 overall_classによる同定結果のファイル名
# 2. overall_genus データベースによる同定
clmakecachedb --blastdb=overall_genus --numthreads=CPUコア数 代表配列ファイル名 overall_genusの出力キャッシュファイル名
clidentseq --blastdb=overall_genusの出力キャッシュファイル名 --numthreads=CPUコア数 代表配列ファイル名 overall_classによる同定結果のファイル名
classigntax --taxdb=overall_genus --maxpopposer=0.05 --minsoratio=19 overall_genusによる近隣配列の取得結果 overall_genusによる同定結果のファイル名
# 3. 同定結果の結合 (overall_class + overall_genus)
clmergeassign --priority=descend overall_classによる同定結果のファイル名 overall_genusによる同定結果のファイル名 最終の同定結果のファイル名
Claident には分類群や解像度の異なる分類群割り当て用のデータベースがいくつか用意されています 。 16S rRNA のみが分析対象の場合は ｢prokaryota_16S_genus｣ などのデータベースを使用するだけで十分かもしれませんが、自分のプロジェクト上、対象の生物以外が検出されることもあり、それらも重要になってくるので ｢overall_class｣｢overall_genus｣ のデータベースを利用しています。ただし、これらは解析のためにメモリを大量に消費するので、通常の 16S rRNA の解析であれば ｢prokaryota_16S_genus｣ など、分類群を絞ったデータベースを使用するほうが良いかもしれません。
｢1｣｢2｣ では clmakecachedb で解析高速化のためのキャッシュを作成しています。これは特に overall_XXXX を利用するときには有効です。その後、clidentseq によりターゲットの配列の近隣配列群を取得し、classigntax により分類群を割り当てています。--maxpopposer=0.05 --minsoratio=19 のオプションを指定することで多少の誤同定の可能性を許容し、できるだけ下位の分類階層まで同定しようとしています。
｢3. 同定結果の結合｣ の部分では clmergeassign を用いて overall_class + overall_genus の同定結果を統合しています。
その後の解析など
DADA2 による配列解析と Claident による分類群の割り当てが終了した後は得られた結果を用いてより生態学的な解析を行います。この際に、よく phyloseq (McMurdle & Holmes 2013 PLoS ONE) というパッケージを利用しています。
また、上記で解説した一連の解析では基本的にはサンプル中の微生物群集の組成しか出せません。つまり、A という分類群が X %、B という分類群が Y %、という具合です。一方で、微生物生態学では絶対量も重要な情報です (例えば、A という分類群が X g/ml water とか)。そこで、絶対量を見積もりつつ、網羅的に (微) 生物群集の DNA を解析するために定量 MiSeq という方法を考案しました (Ushio et al. 2018 Metabarcoding & Metagenomics; Ushio et al. 2019 Methods in Ecology and Evolution)。厳密に絶対量を測定できるわけではないのですが、それでも微生物量の推定ができるので、以後に適用できる統計解析の幅が広がり、議論できることが増えます。
機会を見つけて phyloseq や定量 MiSeq についても解説ができればと思っています。
16S rRNA を利用した微生物群集の解析のコードを通しで見てみたい方は Ushio et al. 2019 Methods in Ecology and Evolution で公開したコード (Zenodo) が参考になるかと思います。
Written with StackEdit.
