yokaのblog

琵琶湖で微生物の研究してます

USEARCHを使った16S rRNAアンプリコンシーケンス析②

前回よりもさらにマニアックなんだけど、備忘で書いておく。
 USEARCHでNGSデータを分析するとき、454のMIDタグ付きのリードファイルであれば、前回の方法でそのまま処理できるんだけど、illuminaのリードの場合は、

①ペアエンドのリードがそれぞれ別のfastqファイルで提供されている
②しかもフォワードとリバースが混ざっていて向きがそろっていないことがある
シーケンサーからデータが出てくる段階でMIDタグごと(サンプルごと)に別のファイルに分割されて出てくることがある

という違いがあるので、そのままでは"fastq_filter"コマンドで読むことができない。そこで以下のようなステップで、前加工をして454のリードと同じファイル形式に揃える必要がある。
(ちなみに454で読んだfastqであっても、投稿者がデータベース登録前にご丁寧に③をしてくれている場合が結構あるので、その場合は同じ作業をやる必要がある)

f:id:yokazaki:20150728212251j:plain



まず、ペアエンドについては、USEARCHのfastq_mergepairsというコマンドで結合することができる。このコマンドの良いところは、それぞれのリードのQ値を引き継ぐだけでなく、2つのペアでリードが被る部分については、読みが一致しているかどうかを判断して、自動でQ値を付け直してくれるところだ(詳細は本家のHPを参照)。なので、illuminaのアウトプットをそのまま放り込めば自動的にコンセンサスのfastqが出てくる。


つぎに、リードの向きがそろっていない場合。USEARCHのコマンドは、シーケンスを前から順番に読んで処理する(逆相補は見てくれない)ので、全てのリードの向きをそろえる必要がある。ここは正しい方法なのかわからないけど、図のように、"fastq_revcomp"で逆相補にした後に、"fastq_strip_barcode_relabel2.py"をかけるという方法を使っている。なぜなら、"fastq_strip_barcode_relabel2.py"は指定したタグ・プライマーが最初に出てこない配列は全てミスマッチとして切り捨ててくれるからだ。これを元配列と逆相補のそれぞれでかけて、ファイルをマージ(シェルのcatコマンドを使う)すれば、正しい方向で正しくタグ・プライマーがマッチしたシーケンスだけのファイルが手に入る。ちなみに、③のようにタグが無い配列の場合は、タグの配列にプライマーの1文字目(全てのシーケンスで1文字目に出てくる塩基)の文字を入れておいて、サンプル名は"sample"のように仮でおいておけば、(本来の使い方ではないけど)きちんと読んできてくれる。


最後に、サンプルごとにリードが分かれている場合。これは、USERACH自身には専用のコマンドはないようで、シェルのsedとcatコマンドを使う方法が本家では紹介されている。この方法に従って、sedコマンドで各fastaの各シーケンス1行目に"barcodelabel=xxx"の形でサンプル名を示すラベルを入れておき、そうして作った各fastaファイルをcatコマンドでくっつければ、454のタグ付きリードを"fastq_strip_barcode_relabel2.py"で処理した時と同じアウトプットが手に入るので、これでそのままfastq_filterにかけられるようになる。ちなみに②でタグが無くて"sample"として仮ラベルを付けた場合は、"sample"の文字列だけを本当のサンプル名に置換すればOK。


 これでやっとSRAに転がっているデータを自由に分析できるようになりました!とはいっても、実際にやってみると、投稿者ごとに色んなフォーマットのファイルが登録されているし、プライマーやタグの条件も論文に書いてあったりなかったりで、かなり大変・・・自分が登録するときは、ちゃんとメタデータを揃えて登録しようと思いました。