2012年11月27日火曜日

CLC Genomics WorkbenchでPacBioデータ

11月はあっという間に過ぎていきました。
サンフランシスコで人類遺伝学会があり、本当はその報告もしたいのですが、日常の仕事に追われ・・・。
PacBioのセミナーは、既に論文やプレゼンで知っている内容だったので、私にとっては驚きは少なかったのも影響してるのかな? 

アカデミックな報告ではありませんが、Linked-Inでつながっていてまだあったことの無かったひとと、会場で会い、ランチをご馳走になりました。
前ブログ「ショートリードの憂鬱」での、PacBioの最初の紹介が、Expression Analysisというボストンの受託会社のウェブセミナーを見たあとだったのですが、当時からその会社の担当者とはLinkedInでやり取りしていて、今回やっと会えました。 CEOにも会ってきました。
そこでの会話の内容はお伝えできませんが、同じPacBioユーザ同士、共感できるものがありました。

また展示会場のCLC-Bio社のブースで、前職にいたとき知り合った、CLCデンマーク本社のひとと1年ぶりくらいに会いました。 懐かしかった。
というわけで今日はCLCの話。

Genomics Workbenchは、使い勝手の良いソフトですが、まだPacBio用の入口がありません。
しかし、諦めるのはまだ早い! 
まず、Pacのリード(filtered_subreads.fastq)は、Import>Illumina ( !! )からインポートします。
Quality value はNCBI/Sanger を選びます。 Paired-endとかのチェックは全て外します。
PacBioのデータはHiSeqとかに比べると桁違いに小さいので、Windows版のDesktop Editionでも十分使えるでしょう。

インポートが無事終わると、このようになります。 シーケンスの塩基数に注目! 10000bp 越えもある!
 
今のはfilterd_subreads、この意味は、SMRT PortalというPac付属のソフトで、
  1. Read長が5000 base 以上
  2. Sub-read長が5000 base以上
  3. Read Qualityが75以上
などという基準でフィルタリングした後のsubreadということです。
まだエラーコレクションも何もしていません。
 
Pacのエラーコレクションで、pacBioToCAというものがあります。 これを流すと、CCSなどのHighクオリティショートリードでCLRの精度が高まります。
こいつのOutputが、fasta + qual ファイルです。
fasta + qual形式ファイルの「エラー修正subread」を取り込むには、Import > Roche 454 ( !! ) 
fasta と qual ファイルを選んだらNextボタンでインポート
というふうにするとちゃんと入ります。 

まあ、fasta と qual があれば、fastqファイルが作れますけどね。
既に誰かが作ったfasta+qual --> fastq converterはいろいろあると思いますが、私はSEQanswersのこのスレッドに出ていた、このスクリプトはいいのでは?と思います。Perlで動きます。
メモリに全部展開するのが嫌?なひともいるでしょうけど。

さて、リードを入れて何が面白いかというと、CLCのGenomics Workbenchはレポート機能が充実しているからそれをやってみると、Pacのデータの特徴が簡単に見えるのです。GCの割合や、リード長の分布、なども1つのレポートにまとまっているから便利。
下図の左はエラーコレクション前、右はエラーコレクション後の、クオリティ(Phredスコア)の分布図です。 例えばこんなグラフを見て、確認しながら、LinuxのPacサーバに戻ってアセンブる。

 
Pacサーバでアセンブるわけは、Genomics Workbenchのアセンブラーが、de Bruijn Graph 式アセンブラーだからです。
ハイスループットのショートリードには向いていますが、長さがバラバラの超ロングリードには向いていないかも(あくまで「かも」ですが)しれません。
現バージョン5.5.1では、k-mer を64以上には増やせないようですので、これ以上のk-mer では試して無いのですが、どうでしょう?
k-mer=64 ではアセンブルの成績は良くありませんでした。


Genomics Workbenchでリードが取り込めれば、リファレンスがあればマッピングも、できます。
上側がエラーコレクション後、下がコレクション前です。E.coliのデータです。
マッピングパラメータはデフォルトでやったのですが、この例では、エラーコレクション前はSubread全体の37%しかマッピングできませんでした。それがエラーコレクションをした後は、99%以上がマッピングされました! 


SMRT Portal/Analysis というPacBio付属のアセンブラーやマッピングツールも勿論よくできています。 Pacデータをガッツリ解析するにはこっちの方が向いているでしょう。
GUI以外にコマンドラインも充実しています。
また、最初のsubreadフィルタリングなどには必須です。

でもCLCのGenomics Workbenchは、レポート機能が充実している、というか、手軽にデータを確認できるという点で優れていると思います。