細菌の16S rRNAアンプリコンシーケンスの分析を一通り勉強したので備忘で。今回はUSEARCHというパイプラインを使って分析を行うことにした。全体像をまとめると以下の通り。
①・②がデータをきれいにするステップ、③～⑥がOTUをつくるステップ、⑦・⑧が作ったOTUにそれぞれのリードを割り当てるステップ、って感じ。オレンジで囲ったのがコマンドで、これをUSEARCHをインストールしたPCのコマンドラインから打ち込んでいけば、処理後のfastaやfastqが指定したファイル名ではき出されてくる仕組み。各コマンドにはオプションでいくつか変数を入れていく必要があるんだけど、特にこちらで考えて指定しなければならない変数を灰色で示した。実際のコマンドの例はUSEARCHのウェブサイトでも見ることができる。
指定すべき変数を見ていくと、まず①ではフォワードプライマー配列と、サンプルごとにつけたmultiplex identification tagの配列をFastaで読ませる必要がある。これをもとにして、tag情報をFastqの1行目に付加し、タグの配列を削除した状態の配列がはき出されてくる。
次の②では、2つの変数をフィルタリングの閾値として指定する必要がある。一つは、expected errors (E)という、fastqのQスコアを元に計算された配列のクオリティーを示す値だ。（expected errorsの説明は本家のページにかかれている。）デフォルトは１だけど、基準を厳しくするのであれば、値を小さくする。もう一つの変数は、取り出す配列の長さだ。ここは結果に結構影響するので、いくつかの値を試しながら慎重に選ぶ必要がある。USEARCHでは、全ての配列を同じ長さに揃えたうえで分析を進める必要があるため、たとえばここで"200"を指定すると、全ての配列がフォワードプライマーからみて200bpまでの配列に切られた状態のファイルがはき出されてくる。情報量を多くするためにはできるだけ長くとったほうがいいけれど、後ろのほうの塩基配列はクオリティが低いものが多いので、長くとればとるほど怪しい配列も増えるというトレードオフが存在することになる。FastQCを使って、あらかじめQ値の分布を見ておき、どの程度の長さまでなら信頼できそうか目星をつけたうえで、あとは実際にその後のステップでどの程度の配列が残るのか、結果がどう変わるかを見ながら、ベストな長さを選ぶ必要がある。このステップで、ファイルはfastqからfastaになる。
③～⑤のステップは、③重複配列を取り除く(dereplication)→④Singletonを取り除く→⑤OTUクラスタリング となっていて、とくにこちらが考えて指定する変数は無い。ちなみに⑤のOTUをつくるところが、UPARSEというアルゴリズムを使っていて、これがQiime等の他のパイプラインよりも精度が高い、という論文があったので、今回USERCHを使うことにしている。
⑥はキメラ配列を取り除くステップだ。ここは他のパイプラインでも使われている、UCHIMEというアルゴリズムを使う。このとき、リファレンスとなる配列集（キメラが混ざっていない、正確な配列）をこちらから指定する必要があるのだけど、これについては、本家のウェブサイトで「これ使えばOK」ってのが出ているので、今はこれを使わせてもらっている。
⑦では、⑥で作ったOTUをもとに、②ではき出されたDereplication前の配列を割り当てていく作業が行われる。このとき、どの程度の相同性があれば同一OTUとみなすか、という閾値を設定する必要があるのだけど、16Sの系統分析の場合、デファクトが97％になっているので、ひとまず"0.97"としておけば問題ない。ここでは.ucという形式のファイルが出てくるが、このままだと読みづらい。そこで⑧のステップで、このファイルをタブ区切りのテーブルに変換する。これをエクセルとかで読めば、「表側がOTU、表頭がサンプルID、中身が各OTUのリード数」という形式の表が手に入るので、あとはこれをRにかけて煮るなり焼くなり。
ちなみにここまでの分析では「どのOTUがどの分類か」という情報は一切手に入らないけど、これを簡単にやるのであれば、RDPのSeqmachに⑥のアウトプットのfastaファイルをそのまま放り込めばいい。結果を見れば、各OTUの分類と、これまでにGenbankに登録された配列で近いものの情報を出してくれる。
今のとこ、こんな理解です。間違っていたり、もっといい方法を知っている人は教えてください。
追記：パート②もあります
