カバレージで精度が上がることを確率で証明（？）

分子生物学会も終わり、今年も残すところあとわずかになりましたね。
年末といえばジャンボ宝くじ！　宝くじといえば当選確率、ということで、無理やりですが今日は確率の話です。


PacBioのリードには、エラーがランダムに入ります。
その割合は前にも書きましたが13～15%くらい。　
ということでリード1本あたりの精度がおよそ85%、というふうになるわけです。

で、この精度を上げるために、カバレージを増やします。

この図は、カタログなどにも載せてますが、これは実測値だそうです。
でも精度85%のリードを重ねただけで、そんなに正確性が上がるのはなぜ？
という疑問を持たれる方もいるかもしれませんね。
高校の確率・統計で苦しんだ私がご説明致します。

話を単純にするために、正しい塩基を多数決で決めることとします。
また、PacBioのエラー率は実際より少し高めの、16.7%とします。

つまり、16.7% = 1/6 の確率でランダムに間違える。
リードが1本の場合、間違える確率は1/6 です。　ここまではOK?

次にリードが2本の場合、
1本が間違える確率　P(1 of 2) + 2本とも間違える確率 P(2 of 2) が、間違える確率です。
ちょっと格好良く言うと、1/6の確率で起こる事象があり、それを2回試したとき、１回以上起こる確率です。
サイコロを2回振って、1が1回以上出る確率と同じですね。

リードが3本の場合、
2本が間違える確率 P(2 of 3) +3本とも間違える確率 P(3 of 3) が、間違える確率です。
同じく格好良く言うと、1/6の確率で起こる事象があり、それを3回試したとき、2回以上起こる確率です。
サイコロを3回振って、1が2回以上出る確率と同じです。

絵にするとこんな感じで、上から1本のとき、2本のとき、3本のとき、の「間違い」リードの組み合わせです。間違いリードは×印です。

さて、ではリード4本のとき、
2本が間違える確率　P(2 of 4) ＋ 3本が間違える確率 P(3 of 4) + 4本とも間違える確率 P(4 of 4)
=＞1/6の確率で起こる事象があり、それを4回試したとき、2回以上起こる確率
（サイコロを4回振って1が2回以上出る確率）

リード5本のときは？

3本が間違える確率　P(3 of 5) + 4本が間違える確率 P(4 of 5) + 5本とも間違える確率 P(5 of 5)
=＞1/6の確率で起こる事象があり、それを5回試したとき、3回以上起こる確率
（サイコロを5回振って1が3回以上出る確率）

さてさて、これらの確率を実際に計算してみましょう！
1回の試行で特定の事象が起きる確率がpのとき、n回の試行で特定の事象がm回起きる確率Pは、

 

です。　これ、関数電卓なら頑張れば計算できますけど、今は便利なサイトがあるもので、

こちら、その名もke!san http://keisan.casio.jp/#menu



ホーム/ 計算応用集/ 確率・統計/ サイコロの目の出る確率 で、上の式の計算が簡単にできます。

そうすると、先ほどのリードの場合は以下のような数字になります。

ここまで 0.167 -> 0.306 -> 0.074 -> 0.132 -> 0.035 というふうに間違い確率「エラー率」は上がり下がりしながら全体としては下がってきました。

これを、n - > 50 までやってみるとどうなるか？

11カバレージでQV23、19カバレージでQV34、30カバレージでQV46、40カバレージでQV59

このグラフ、上がったり下がったりしながらも、0に近似されていく、・・・　美的ですね。

 

これを、1000bpのシュミレーションデータで検証した方がおります。

Sergey Koren et al., “Hybrid error correction and de novo assembly of single-molecule sequencing reads” Nature Biotechnology  (2012)

 

先のグラフと非常に似ていますね。

 

エラーが「ランダムに起こる」ということがこの大前提にあります。
ランダムであれば、リードを重ねれば全体としてのエラー率は確実に下がっていくのです。
どうでしょうか？
85%精度のリードでも、重ねると100%に近くなる、というのは本当なんですよ。

と、今日はデータの精度を、確率論で説明を試みました。
数学者からはツッコミが入りそうですけれど。
どうでしょうね。

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は、レポート機能が充実している、というか、手軽にデータを確認できるという点で優れていると思います。

SMRT Cellに使用期限がある理由

NHKの、「探検バクモン」という番組で、東京国立近代美術館フィルムセンターの潜入をやっていました。　

名画から普通の映画、珍しい映画や成人映画まで、およそ6万5000本の国内映画のマスターフィルムを、国が保存しているそうです。こんな施設があったとは知らなかった！
「仕分け」されなくて良かった。

目からウロコだったのは、フィルムというのは原料にセルロースを含んでいるので、高温多湿に弱く、ほっておくと酸化してしまってすぐにダメになってしまうそうです。
ここでは温度5度、湿度40%の秘密の倉庫で管理しています。
フィルムって、有機物だったんだ、って思いました。

有機物って言えば、SMRT Cellもそうなんです。
PacBioでシーケンスに使われる、あのタテヨコ1cmくらいのチップです。

SMRT Cellにある15万もの穴は、ZMWと呼ばれ、その底でシーケンス反応が進みます。
穴の底にはビオチン分子があらかじめ接着されていて、反応溶液＋ライブラリーを入れたときに、ライブラリーに結合しているポリメラーゼがビオチン分子に結合し、底に固定化されるわけです。

Korlach et al., Selective aluminum passivation for targeted immobilization of single DNA polymerase molecules in zero-mode waveguide nanostructures. より

ビオチン分子をZMWの底に接着させる作業は、Cell作製の比較的早い段階で行われます。

かなりシークレットな情報なので大雑把に言うと、半導体のウエハーってわかりますか？　丸い大きな薄い板みたいなやつです。
あれ一枚から数百個のCellができるのですが、最初にZMWという穴をあけて、ビオチン分子を接着させます。　
その後、穴以外を薄いガラスでコーティングしたり、洗浄したり、１個１個にしたり、60以上の工程を経て、晴れてSMRT Cellが完成します。
以上はダストフリーの、本当に半導体工場みたいなところで行われています。

というわけで、SMRT Cellには、出荷段階で各ZMWに、すでにビオチン分子がくっついています。
これは紛れもなく有機物なので、温度と湿度に弱いのです。
使用期限があるのはそのためです。



ところで最近見た映画は何？と聞かれて、
「アベンジャー」と答えるのはダメですか？

ちゃんとフィルム映画をみなきゃいけませんね。

LSC Long-Readのエラーコレクション

秋は学会・展示会シーズンですね。　
私は10月10日のBioJapan、10月15日のCBI、に行きました。　
来週の日本遺伝学会は残念ながら行けないのですが、来月サンフランシスコで行われるアメリカ人類遺伝学会には行く予定です。
今から楽しみです！

10月10日のBioJapanでは、前回のブログでも書きましたが、沖縄県のセミナーでPacBioの発表がありました。
そこでちょっとサプライズゲストとして、PacBioの創始者で現CTO、Steve Turner氏の挨拶を頂きました。
その場にいたひとは聞いたと思いますが、PacBioはこれからも進化します。　
リード長ももっと長くなります。
新酵素の開発、新ケミストリーの開発、Movie開始時期の改良、ソフトウェアの改良・・・。

そこで良く聞かれるのが、「リードは長くなっても、精度は低いままじゃあちょっと・・・」という声。

リード単位の精度はもっと上がるのか？　
リファレンスにアラインしたとき、塩基がリファレンスと異なる（InDelも含め）率が約15%
これが劇的に改善されるか？
本音はというと、ちょっと・・・

そこで、バイオインフォマティクスでエラー率を低くする方法がたくさん開発されているのです。

7月5日のブログにも書きましたが、エラーコレクションといえば、Nature Biotecで発表されたpacBioToCA が有名？でしょうか。　最近のセミナーでも、さかんにこれを使った発表がされています。

そんな中、先週別のアルゴリズムが論文化されました。　これはRNA-Seq用にも適したエラーコレクションアルゴリズムだそうです。

Au et al.,  Improving PacBio Long Read Accuracy by Short Read Alignment
PLoS One 7(10), e46679.

Homopolymer Compression という方法で例えば、GCGAAAATA　＝＞　GCGATA
に情報を圧縮します。　Pylosequencerなどでのエラーコレクションに使われる手法で、これをPacBioにも応用しています。
情報量は意外と失われないそうで、私も近いうち試してみたくなりました。

 

pacBioToCAと同様、ショートリードでPacのロングリードを「修正」するアルゴリズムですが、Mammalianのトランスクリプトーム用に開発されただけあって、メモリー消費やラン時間に工夫がされているようです。　pacBioToCAは処理時間が長い、LSCよりメモリを食う、と。
LSCというのがこのアルゴリズムの名前です。

軽いなら是非試してみたい！
ツールも公開されていますので。
そのうちこのツールを使用した発表も、学会などで多く見られるようになるでしょう。

そういえば今年の7月くらいにPacBio本社で、cDNA解析のPacでの可能性についてディスカッションをしたとき、マッピングの方法と合わせて、エラーコレクションの話をしたことを思い出しました。　
ちなみにPac社内では、「error correcton」という言葉を嫌うひともいます。
correction = 修正です。　修正するってことは、もとが間違っているということ。　なるほど。
ですので、私は気を利かせて「improve accuracy = 精度の向上」と言うようにしています。
ま、単なる言葉ですけど。

ショートリードでロングリードを「修正」する方法は目新しくなくなりましたが、「修正」すること無しにAccuracyを良くする方法も出てきています。
それがQuiverというツールです。　
これに関してはまた今度。

最後に、精度には2つの意味があります。
リード単位の精度とマッピング精度です！

リード単位の精度ではショートリードの方がずっと良い。　これは認める。　
では、リファレンスやゲノムにマッピングさせたときの、真の場所に正しくマッピングされるという意味での精度は、どうか？
Yes!　ロングリードであればある程、良い、でしょう。

本当はリード一本で精度がphred 40 以上くらいあって、かつ数十キロ以上読めれば（＋低コストで）、向かうところ敵無しなんでしょうが。
まだ無理かなあ。

PacBioを使った国内での研究

東京大学農学部2号館で開催された、インフォマティクスオープンセミナーに行ってきました。
"Challenge to de novo sequence of relatively large genomes with new sequence technologies"
BGI、東京大学 大学院農学生命科学研究科
アグリバイオインフォマティクス教育研究ユニット
新学術領域「複合適応形質進化の遺伝子基盤解明」共催
http://www.iu.a.u-tokyo.ac.jp/main_event.html

演題は5つあり、そのうち3つにて、PacBio RSのデータを使った話があるとのこと、これは "Must Lissen to" だろうと思っていましたが、まさにその期待通りの内容でした。

de novo Assembly がテーマなので、非モデル生物のゲノムプロジェクトの話、どうやって数ギガのサイズのアセンブリを遂行させるか、PacBioを使ったときの良いところ、改善が必要なところ、エラーコレクションの方法、リピートの処理について　…etc.

改めて言うまでもないかもしれませんが、高等生物のゲノムアセンブリの場合、リピートの扱いがとても厄介なんです。　複数の繰り返し配列、segmental duplication などは、ショートリードではちゃんと正確に読むことが難しいかほぼ不可能。

そのためロングリードが救世主となるわけなんですが、Pacといえども短所があります。
長さが一定では無い（例えば出力リードの長さが "平均" 3Kbであること）
シーケンスエラーが他のテクノロジーと比べると高い（15%くらい）
ランコストがまだまだ高価（というと怒られるかもしれませんが）
などです。
そこを考慮しても、他のシーケンサーでは得られない「超ロングリード」という長所が、de novoのアセンブリにはとても魅力的なんですね。

話のほとんどがまだ論文途中のものでした。
ですが皆さん、結構オープンに話していて、会場も質疑応答が絶えませんでした。

知っての通り、PacBioはまだ日本に来て1年も経っていません（2012年9月現在）。　ですから、実際どんなもの何だろう？　本当に使えるのか？と疑問に思っている研究者は多いと思います。

アメリカでは、Early Accessと言って2年以上前から導入されていたところもありますし、現在も多くのユーザーがいて、ノウハウもだいぶ蓄積されつつありますが、なかなか日本には情報が入って来ないと思います。
そんな中、日本のPacBioユーザー、あるいはBGIのような日本のPacBioユーザーの共同研究先が、本音で語ってくれるこのような会はとても貴重でした。
我々メーカーが宣伝するより何倍もインパクトがありますからね。
オーガナイザーの長谷部先生、西山先生、門田先生には感謝です!!


さて、今回の他にも、Pacを使った日本での研究例をもっと知りたい！という方にお勧めのセミナーがあります。

10月10日（水）、パシフィコ横浜で開かれるBio Japan 2012で、「先端シーケンサーが拓く沖縄生物資源」と題したスポンサーセミナーがあります。　詳細プログラム
沖縄工業技術センターの照屋先生の講演がPacBioを使ったゲノム解析についてです。
当日、ちょっとしたサプライズも用意していますので、時間と興味のある方は是非参加したらいかがでしょう？
無料ですが、参加にはBioJapanの公式ホームページから登録する必要があるそうです。

論文紹介：SNP検出のバリデーションに使用！

最近、PacBioを使っての論文がたくさん出ています。
一番最近のは、PacBioの超ロングリードに対応したMappingツール、BLASRのアルゴリズム
Mapping single molecule sequencing reads using Basic Local Alignment with Successive Refinement (BLASR): Theory and Application (PMID: 22988817)
このBLASRというMappingツールは、PacBioの二次解析ソフトウェアについてきます。
もちろんフリーなので、PacBioのユーザー以外でも使うことができます。

このBLASRがまだ前のバージョンだった頃、Broad Instituteでは、既存の有名なMappingツール、BWA-SWを使ってPacBioで読んだターゲットリードをゲノムにMappingし、その後GATKでSNPを検出、既にショートリードシーケンサーで検出したSNPを検証しています。
検証にはPacBioの他にFluidigmも使用しています。

Medulloblastoma exome sequencing uncovers subtype-specific somatic mutations. (PMID: 22820256)

92人の髄芽腫と健常人のサンプルを比較しての変異解析です。
Whole-Exomeキャプチャー法で33Mbの範囲を、HiSeq2000 2x76bpでディープシーケンスし、統計的有意な12の遺伝子変異を見つけています。
この変異遺伝子のコーディング領域をFluidigm Access ArrayとPacBio RS で読んで、検証しています。

髄芽腫の、Whole Genome Amplification サンプル（PCR増幅）と、Native DNAサンプルそれぞれ48ずつ、バーコードシーケンスしています。　ちなみにPacでバーコードはまだ製品になってません。

20ラン行い、フィルタリング後のリードが389,215本、サブリードは 4,367,852本、そのうちサンプル由来だとわかったのは64%の2,834,170本のサブリード

このサブリードを、hg19に対してBWA-SWマッピング、GATKを使ってBase-quality recalibrationとSNP 検出。
検出した変異候補は、IGV上でマニュアルキュレーションし、カバレージが十分にあった19の変異を検証しています。

結果ですが、リードごとのエラー率は高いけれど、カバレージを稼いでコンセンサス配列にするとSNPバリデーションに十分使え、実際に2bp detetion をCTDNEP1遺伝子で確認しています。
他のショートリードテクノロジーは、エラーが起こる場所に配列的特異性があり、ランダムでは無いので、コンセンサスにしてもエラーの高い場所と低い場所ができます。
PacBioのエラーは全くランダムに起こる（と言われている）ので、そのぶん、ショートリードで見つけた変異の検証に向いているということです。
しかし、PCRで十分増幅されなかった他の変異箇所、またGenomc DNAを十分な量得られなかった変異箇所は検証に使用できませんでした。

これの補助的論文として、
Pacific biosciences sequencing technology for genotyping and variation discovery in human data. (PMID: 22863213)

があります。
BWA-SWのパラメータは、

• Missmatch penalty を3から5へ
• Gap open penaltyを5から2へ
• Gap extension penaltyを2から1へ
• Heuristics （信頼性・精度）を1から20へ　（Ampliconは精度が高い）

と工夫しています。こちらの論文はAmpliconを読んでいます。　同じパラメータを髄芽腫論文でも使用しているそうです。

 

この論文の途中、次のような一文があります。

“…Pacific Biosciences does not provide a simple measure of base quality score. Instead the software uses the instrument’s estimated insertion error probability as the base’s quality score.”

疑問に思って著者に聞いてみました。

この論文が書かれた時は、PacBioのSMRT Analysisという二次解析ソフトでBLASRマッピングをしてBAMを出すときに、Qualityがinsertionのエラーだけを基本に算出されていたので、実際の質とかけ離れていると問題だったそうです。

それで彼らもBLASRを使わずにBWA‐SWを使ったわけです。

今のPacBioのソフトは彼らが考えた、"PacBio Processing Pipeline" を取り入れ、Insertion以外のエラーも基本に「ある程度正しい」塩基クオリティを出すようになったそうです。

 

今度は誰かが、そのクオリティはどれくらい確かか、を研究対象にした論文を出しそうな予感が・・・。

 

 

 

DNAはたくさん必要

もうすぐ9月だというのに、まだ暑いですねえ。
私は蚊が嫌いです。　夜中でも殺すまでは寝られません。
書いている今も、どこかに蚊がいるんですよね。　「虫よけ当番」を枕元に置かなきゃ。

ま、それは置いておいて、PacBioは1分子シーケンスを売りにしています。
基本的にPCR増幅はしません。
想像はつくと思いますが、そうするとスタート時にDNAはたくさん必要になるのです。
ちなみに今は、2本鎖のDNAだけ、からスタートできます。

スタート時に必要なDNA量ですが、サイズによっても違いますが、大体、10Kbのライブラリーで読むときは、9μgくらいは必要でしょうか？
と、これはアメリカの受託会社、Expression Analysis （EA）のドキュメントが公開している数値です。
Webで公開しているので書いても良いでしょう。
ドキュメントはここ（http://www.expressionanalysis.com/platforms/category/pacific_biosciences/）
のページの、右上の、Comprehensive PacBio RS Overview　というところからダウンロードできます。

実は、この会社との出会いが、私をPacBioにのめり込ませたのですが、この会社は、アメリカでマイクロアレイやシーケンサーの実験・解析受託をしています。
昨年の6月に、PacBioのサービスを始めて、安定してビジネスを展開しているそうです。

さて、DNA必要量の話に戻りますが、9 μg という数値はかなり多いでしょう。
シェアリングして、End-repairして、Ampure Washして、…と進めていくと、ライブラリーにするころには、750 ~ 1250 ng が予想される回収量になるそうです。
（なるそうです、と書きました。いえ、私も予想回収率くらい知っているんですよ。でもPacBioのプロトコールはWeb公開されていないので、ブログで書くのはまずいからあくまで第三者のふりしてます。詳しく聞きたい方は別途）
EA社はこの場合、20 ~ 35のSMRT Cellが使用できる、と書いています。

で、データはどれくらい出てくるか、ですが、10Kbで90分Movieで読んだとき、生データで最長20Kb、クオリティフィルタリングした後で14.6Kb
平均リード長は、生データで3.7Kb、フィルタリング後で2.6Kb
フィルタリング後のリード数は44,232本、塩基数は116Mb

このデータは、比較的現実に沿った数値です。
チャンピオンデータではありませんが、いろいろなサンプルをやった経験値でしょう。

フィルタリングの条件は、リードにはそれぞれRead Qualityという精度の指標があるのですが、これが 75/100 以上であることと、クオリティの良い部分のリードの長さが50bp以下ではないこと、です。
これにひっかかるリードは除去されます。
そうして残ったリードの、平均長が 2,600 bp、Maxが14,600 bp ということなんですね。

Read Quality 75以上、リード長50bp以上、というのはデフォルトの値です。
50bpでは短すぎるな、と思うときは私はこれを上げていいと思います。


サンプル調整自体は、私もトレーニングを受けましたが、それほど難しいステップはありません。
シェアリングもCovarisの機械任せだし、精製もマグネットビーズでピペット処理だし、アダプターダイマー（アダプター同士の接着）をこれまたビーズで取り除くだけ。
基本プロトコールは順調です。
でも、裏プロトコール？みたいなトラブルシュート的なプロトコールでやったときはほとんどDNAをロスってしまいました。
ラボを離れて約7年・・・の私には、裏プロトコールは難しかった。
とは言っても、基本プロトコールで作った10Kbのライブラリーの方は、ライブラリーまでの回収率は33%で、Good！
ランの結果も、上記EA社の数値と同じくらいでした。


ところで、より少ないスタートDNA量でランできないものでしょうか。
たくさんDNAがとれないサンプルでも、長く読んでみたいという要望はあります。
この要望に応えるための方法は、・・・
公式にPacのWebSiteで公開されたらここでも書こうかと思います！
まだ書けないーっ！

Base Modification 塩基修飾 （新

気づいたら8月ですね。　7月はブログ、サボっていました。
「最近更新してませんね」って言われるようになり・・・。

全然関係ないですが、西日本では、クマゼミが鳴いていますね。
沖縄からずっと愛知県くらいまではクマゼミ、見ますね。
でも東京（少なくとも私の住んでる板橋区）では見たことありません。
子供の頃から、図鑑でしか見たことのないクマゼミを、大人になって初めて見たのは京都に出張した時でした。　3年前くらいのこと。
今でも出張先で見つけると、「ああ、捕まえたいなあ。東京に持って帰りたいなあ」って思います。


さてさて、今日はBase Modificationの話です。
塩基修飾とも言います。　
Pacの機械はシーケンサーなのに、DNA修飾を検出することができる！というのが売りのひとつです。　
営業的には正しいのですが、過度な期待は禁物、ということで、ここではもう少し正確に書こうと思います。


今まで、「ショートリードの憂鬱」というブログの中で書いたことがあるんですが、
「メチレーションをダイレクトに検出するシーケンサー」
この中に書かれていることをちょっとだけ修正します。　まあ、こちらのブログを書いたときは、私は別の会社にいましたので、情報があまり無かった、という言い訳も成り立ちますよね。

どこを修正するかというと、後半部分に書いた、何百回も読まなくても同じ場所を5回くらい読めばメチル化と非メチル化のDNAを区別できる、という所。　これは人工配列でできた、というレベルのことで、実際にE.coliなどの配列で区別するには、もっともっと深いカバレージが必要です。

ここからが本題！

PacBioに付属する解析専用ソフト、SMRT Analysisの最近のバージョンに、新しくBase Modification 検出機能がつきました！
って言うと格好いいんですが、正確には、Base Modが起こったかも知れない場所の検出機能がつきました！　というのが正しい。

どんな種類の塩基修飾なの？　どれくらいの信頼性で？　というところは次バージョンに期待してください。

先のブログにも書いていますが、メチル化サンプルと非メチル化サンプルとを両方持っていて、それぞれPacBioで読んだ場合、これら2つのサンプルの塩基のIPDを比較します。
IPDというのは、Inter Pulse Duration の略で、ある塩基が読まれたあと次の塩基が読まれるまでの間の時間と同じ意味です。
メチル化された塩基がテンプレートにあると、Polymeraseがポージング（ちょっと止まる・スピードが落ちる）を起こし、塩基取り込みまでに時間がかかります。

挿絵は　Flusberg et al. (2010) Nature Methods 7: 461-465

例えば6mAの場合、テンプレートがmAのとき、Tが結合するまでの時間は、テンプレートがただのAの時と比べて平均５～６倍長い。その比（IPD-Ratio）を見ることで、ここに修飾が起きていそうだな、と予測できるのです。

これも同じ論文から。　6mAの場合です。　IPD-Ratioが縦軸で、横軸は塩基の場所です。

このピークの高さと並び方をシグナチャーと呼んでいます。

必ずしもメチル化されている箇所でPolymeraseがポージングするわけではなく、メチル化箇所の周囲の塩基でもポージングされることがわかっています。

そしてそのパターンが、塩基修飾の種類によって、様々に異なることもわかってきました。

PacBioでは、色々なタイプの塩基修飾、DNAダメージのIPD-Ratioシグナチャーを集めています。

（ブログでは明らかにできませんが、プレゼンでなら紹介できますので興味あるかたはお知らせ下さい）

この情報を使えば、配列を読んだだけで、IPD-Ratioのパターンから、どんなタイプの塩基修飾が起こっているのか、予測できる、というわけです。

と言っても、現バージョン（2012年8月現在）ではそこまではできません。　

では、今は何ができるか！

1. 読まれた部分の配列の全IPD-Ratioの情報を使って、バックグラウンドから有意なIPD-Ratioを見つける　（ここに塩基修飾が起こってたかも？と期待させる）
2. IPD-Ratioのコントロールは無くても良い。コンピュータで求めたin silico コントロールを使うことが可能　（これは約12塩基の長さで、全塩基組み合わせのIPDをあらかじめシュミレートしたもので、生物種を問わず使うことができる）
3. 意味ありげなIPD-Ratioの場所をリストしてくれる
4. そのリストは、Viewerでも表示できるし、前後20塩基の配列とともにGFFファイル、塩基単位のCSVファイルで出力
5. GFFからはモチーフを見つけたりできる

真ん中の青と赤の棒がIPD-Ratioです。　ここに塩基修飾があったかもしれない。

他の場所と比べて抜きん出ていますね。

実際このIPD‐Rが周りと比べて有意かどうかは、IPD‐Rの高さだけでなく、Coverageの深さも問題になります。

そこでCoverageについて

IPD-Rシグナチャーの種類によって大きく異なります。

はっきりしているシグナチャーを持つ塩基修飾タイプ（先の6mAなど、５－６倍のRatioを持つもの）については、15－20xが最低必要で、Pacでは50xを勧めています。

IPD-Rが2くらいの、ゆるーいシグナチャーを持つタイプは、最低200xは必要です。

ではでは、はっきりしたシグナチャーを持つ6mAを読み取ろうとした場合、どれくらいのSMRT Cellが必要なのか？

ターゲットサイズが5Mbだとします。

＋/－両方のストランドを読みますから　5Mb x 2 = 10 Mb

50カバレージ必要だとすると、10Mb x 50 = 500Mb

1 SMRT Cellで100MbのMappableなシーケンスが出てくると仮定すると、500Mb / 100Mb = 5 個のSMRT Cell が必要という計算になります。

（大体1kb前後のライブラリーを作成して読むことを想定）

無論、Mappableなシーケンスの数は生物種によって様々なので、最初に1個か2個ランして確かめる必要はあります。

＋/ - それぞれのストランドで50x必要、というのが忘れがちですが大事ですね。

5個だとすると、SMRT Cellは1連8個入りですから一回のランで済みます。

45分x2回のムービーでシーケンスするとして、90分 x 5 = 7.5時間

機械にセットしてからのロボットによる試薬調整、ベースコールや転送に+2時間として、ランにかかる時間は、合計10時間弱

オーバーナイトで行う感じでしょうか。

解析には2、3時間くらいかかります。　

論文紹介： PacのエラーのCorrectionについて１

祝！　1万ビュー突破！！

　　

　　こんなマニアックなブロブなのに、始めてから約半年で、10,000 ビューを達成しました！

　　ありがとうございます！　

　
　　さて、以前、紹介しました、PacBioデータのエラー修正に関する論文がNature Bitechから発表されました！

“　“Hybrid error correction and de novo assembly of single-molecule sequencing reads.”

　　 　Sergey Koren et al.,　

·   　“A hybrid approach for the automated finishing of bacterial genomes.” 

　　　Ali Bashir et al.,

　　

　　

PacBioのロングリードは、そのままではエラーが多い（15%くらいと言われています）ので、それを修正するために、CCS（短いリードだけれども何回も同じ所を読んで精度を上げたコンセンサス配列）や、Roche454のリードや、Illumina GAIIやHiSeq2000などのショートリードで同じサンプルを読み、これらショート（といっても数百塩基まで）を、Pacのロングに対してマッピングして、Pacのロングの精度を飛躍的に上げる、という方法論。

ツール自体は今までもあったのですが、論文が出ましたので、少しホッとしています。

方法論については、かなり細かく書かれています。

Celera Assembler というアセンブルツールで最終的にアセンブるのですが、エラーを直さないままのPacロングリードでは全くダメ。

そこで、精度の高いCCSやショートリードをマッピングしてエラーを直してあげることで、ロングの精度を上げる。



http://seqanswers.com/forums/showthread.php?t=18478)を見ると、このツールを使うのに苦労しているひとはいるようです。　
結構時間がかかるそうです。


そこで実際に、Pacのサイトから誰でも自由にダウンロードできるデモデータを使ってみました。


大腸菌MG1655のSequenceデータです。
Long reads: 1 cell of 10 kb shears for CLR　（1X90 min movies）
Short reads: 16 cells of 2 kb shears for CCS　（2X45 min movies）


詳しくはまた色々試して集計してから書こうと思いますが、とりあえず時間だけ言うと、2時間弱で終わりました。　専用サーバですが、16コアCPUでメモリ32GBなものです。

ということなんですが、SEQanswers の書き込み(

論文：Capturing Single Cell Genome

学生時代を思い出してください。　有機生物学、これは覚えることが多くて大変でした。
教科書に出てくる亀の甲羅をノートいっぱいに描いて、ひたすら暗記したものです。
私はそのときアメリカの大学で勉強していて、先生はイギリス人で発音がいわゆる正統派のEnglishでかつ早口。　聞き取るのにとても苦労したのを思い出します。

さて、その中でも好きなのはPolyssacharide、多糖類です。　イギリス人教授の「ポリサッカライド」という発音が、とても印象に残ったのです。　書くときはよく綴りを間違えましたけど。

そんなPolyssacharide、あるバクテリアによって分解されることが知られています。
Verrucomicrobiaとう門のバクテリアです。
このバクテリアのゲノムを読んで、中の遺伝子を解析したのがこの論文です。

Martinez-Garcia et al., (2012) Capturing single cell genomes of active polysaccharide degraders: an unexpected contribution of Verrucomicrobia. PLoS.One. 7, e35314.

Laminarin

(C6H10O5)x

という多糖類の仲間があります。　海藻に多く含まれ、特にコンブの質量の４０％くらいを占めるそうです。　別名βグルカン。
このLaminarinを分解するバクテリアをとるため、4μMの蛍光付きLaminarinを食わせて、12 ~ 120 分後に蛍光がついた細胞（バクテリア）を回収したそうです（fluorescence-activated cell sorting）。
次に、細胞内のゲノムを、whole genome displacement amplificationという方法で増幅して、1細胞ごとのゲノムを読めるまでに量を増やし、同時に菌種を分類するために16S配列（rRNA遺伝子）をPCRで増やしました。

16S配列でゲノムの種類（バクテリアの種類）を門のレベルで分類すると、Laminarin分解に寄与したバクテリアは、ほとんど（90%以上）が Verrucomicrobiaでした（n=121）。

では、このバクテリアのゲノムの中はどうなっているのか？
このゲノム解読に、Illumina GAIIx＋PacBio RS のシーケンサーが使用されています！
サンプルDNAは、このバクテリアゲノムです。　1細胞単位のゲノムを読むために、ゲノムDNAを先の方法で増幅 して、Single-Cell Amplified Genomic （SAG）DNAを読んでいます。

M&Mを読んでみると、Illumina GAIIxの方はインサートサイズ320－540bpのPaired-endで、リード長は110bp、1レーンに3ライブラリーずつ、5レーン流しています。

PacBio RSの方は、先ず、SAG DNAをCovaris を使って2 kbに断片化しています。
精製した1 μgの2 kb DNA断片をもとにライブラリーを作成し、その後はPacのプロトコル通りに進めています。

Pacの場合、SMRT Cell ごとに出力されるリード数にばらつきがあります。　これはどうして？　って思う方もいるでしょう。　
ばらつきの原因はいろいろあるのです。　これだけで長い話になるので別の時にします。
ま、とにかく、SMRT Cellから出てくるリード数は結構ばらつきがあるので、サンプルあたり、10Xのカバレージを得るために、6個から20個のCellを使ったそうです。

Movieは40分x2回
フィルタリング条件として、リード長 > 100bp かつ Read Quality > 0.85　（リード長です。サブリードではありません。ちなみにRQ > 0.85 は少し厳しめですがこのくらいが良いのかなあ、って値です。）

こうして得た配列をもとにアセンブルには、

1. Velvetを使ってKmerを変えながらIlluminaデータからContigを作り、
2. PacBioのハイブリッドアセンブリツール「AHA」を使って、ContigどうしをPacのロングリードでScaffolding
3. Sequencher v4.10.1（Gene Codes社：http://genecodes.com/）というソフトを使って、さらにアセンブル
4. Contig、Scaffoldの末端を精査

で、結果

1. 一応、5サンプルを読んで、32～88%のゲノムリカバリー　（ばらつき多い？）
2. Contigは700～1400本　（んー、まだ多いかな。数本を期待。でもゲノムフィニッシュが目的では無いんですよね。ゲノムに含まれる遺伝子の種類と数を測るのが本研究の目的だそうです。）
3. ちなみに最長Contigは、116,000 bp

ドラフト配列はGenbankにSubmitされているそうです。

で、これらVerrucomicrobiaのゲノム解析でわかったことは、laminarinase, xylanaseなどのpolysaccharide分解酵素の遺伝子、Putativeなものも含めて様々なBiopolymerを分解する酵素遺伝子の配列を、他のバクテリアより多く含んでいた、ということ。


正直言って、「Pac使って良かったー！」的な論文ではありません。
多分、IlluminaのリードだけではScaffoldingが難しくて、そこにPacの出る幕があったのでしょう。
裏方で頑張った、って感じです。

今はC2ケミストリー、で出力も上がっていますから、もっと長いContigができるかも。

PacBio Bio-Infomatics Workshop in Tokyo

先週は大忙しでした。
火～木曜の3日間、PacBio Bio-Informatics Workshop と称して、日本と韓国のユーザー、台湾、香港、韓国の代理店、PacBio Asia Office (シンガポール) のメンバー、それにPacBio米国本社からのトレーナー、計30名余りを招いて、勉強会を開きました。

PacBioにとっても、初めての試みでしたが、大成功？に終わって良かったです。
どんなワークショップにするか、テーマはどうするか、時期はいつが良いか等々、
2か月前から東京と、Menlo Park（Pac本社があるところ）にて打ち合わせ、あとはWeb会議などで詰めていきました。
結構アメリカ人的なトレーナーの進め方を、「ゆっくり話す、ひとの話を遮らない」など、日本人向けに修正しました。
私もアメリカでトレーニングを受けたのですが、結構みんな質問すると途中から議論っぽくなってくるんですよね。　で、時間が押してくる。　日本ではありえない、って言いました。　言葉の壁もあるので。

やった内容は
1日目：PacBioのデータとソフトウェアの概要・特徴、コマンドラインツールの使い方、Hands-on
2日目：de novo Assembly とエラーコレクション、塩基修飾、Hands-on
3日目：Target Sequencing, cDNA Sequencing,

このような感じで進めました。　
英語で行われるので、連続45分以上やらないようにし、各セッション後に休憩時間を15分くらい設けて質問やディスカッションタイムを多く取りました。　
これが良かった！
皆さん、Pacのメンバーと直接話せる機会をフルに生かしていました。

ランチやディナー、飲み物の用意、ネットワーク環境の準備など、裏方で働いてくださった方にも感謝しています！

ワークショップで行った内容も、OKなものはここで公開していこうと思います。

リード長とインサート長の関係

先週の「現場の会」で私の喋りの後の質問の中で、ひとつ、深く考えさせられるものがありました。
「リード長を伸ばすために常に酵素を開発しているということだが、酵素をより良いものに変えるのと、インサートの長さと、どちらがクリティカルな問題なのか？」

つまり、酵素を良くすれば長いリードが読める。90分どころか120分、240分レーザー当てても死なないで塩基を読める酵素も出てくるかもしれない。
しかし、一方で、インサートの長さを今の10kbから20kb、30kbと伸ばしていって、これがちゃんと読めるようになるのか？
このせまい小さなZMWの穴に入ることができるのか？

化学的な条件と、物理的な条件と、どちらがよりクリティカルなのか？　ということですね。

今のC2ケミストリーでは、10kbまでのインサートを読むことができて、90分Movieで読んだとき、リード長は平均3.5kb、95パーセンタイルでは8kb、最長は15kbくらい、ロングリードでの精度は約85%です。

-----------------ここから先は仮定の話を前提にしていますので注意！--------------------

これが将来、もし、30kbまで読めるようになったとしましょう。　
でも、30kbまで読める、ということと、30kbのインサートが読める、ということは別ですね。
10kbのインサートを1周半（sense + antisense + sense）読めば30kbですから、30kb読める、ということにウソは無いですね。うん、正しい。
でもちゃんと説明しないと、普通、「すごい！それって30kbのインサートが読めるってことでしょ？」って誤解するでしょうね。　うんうん。
（そのときもっと長いインサートも読めたら最高ですが）

で、30kbのインサートを読めるか？ってことになると、今度は物理的な問題ですね。　ZMWの穴の中に効率的に入らなければいけませんから。
短いインサートの方がZMWに入りやすいってことは知られていますので、マグネットビーズでロングの（今は10kbですが）ライブラリーを効率的にZMWに誘導する方法は、開発されています。　学会のポスターでも発表されていますのでご存じの方もいるでしょう。



・・・・・　と、ここまでぎりぎり公開OKな話でした。

リード長が伸びたとしても、それは読める範囲（インサート長）が伸びた、ということとイコールではないので、誰かに説明するときはこれから注意します。
初めて聞いたひとは絶対誤解すると思いますから。

このブログを読んだあなたは、もう間違えませんよね。




これは私が「現場の会」で喋った、という証拠写真（現場の会・アルバムより拝借）

PacBio RS に付属するソフトの概要

この間「NGS現場の会」で、参加者のひとたちと話していて、気づきました。
PacBio RS に付属するソフトウェア、について今までブログで全然触れていなかった、と。

"PacBio RS Software and Data Analysis" という、最近出来立てのチラシがあるんですが、その中にこんな絵があります。

中央下のシーケンサーから時計回りに、

• RS Remote
• RS Touch
• SMRT Portal
• SMRT View
• DevNet

とあります。このうちDevNetはPacBioのサイトで、開発中のツールなどを提供しているもので、今日は含めません。

では、順番に説明します。

1. RS Remote：　ユーザの用意したWindows PC にインストールします。　ここでランの設定や、実行中のランの様子を確認します。　終わったランの、時間軸でのジョブの様子を見ることもできます（例えば何時に何番目のセルのランが行われていたか、など）
2. RS Touch：　シーケンサー本体にインストールされていて、タッチパネルから操作します。　主に、RS Remoteで設定し保存されたランを、実行する時に使います。　実行中のランの様子をここからも確認できます。　装置のエラー（温度やレーザーの異常、データの転送エラー）などがあれば、この画面に表示されます。
3. SMRT Portal:　ユーザが2次解析を行うときのウェブブラウザアプリです。　ユーザは、用意したLinuxサーバに、SMRT Analysis という2次解析ソフトをインストールします。　SMRT Portalは、このSMRT AnalysisソフトのGUIで、Internet Explore, Fire Fox, Google Chrome のブラウザに対応しています。　SMRT Analysisについては後ほど。
4. SMRT View：　いわゆるゲノムブラウザです。　Java で動き、Mapping の結果などを参照できます。

これらのソフトは似たような名前なので、最初は良くごっちゃになりました。　

特に、RS RemoteとRS Touchは機能が似ている上、名前も似ている。

そしてこれはデータの流れをまとめた図。 

ランの設定をするのがRS Remote（左下）

シーケンサー本体のRS Touch上でランを実行したら1次解析のデータ（ベースコール結果）が作られます。　この1次解析自体は、シーケンサー横の黒いベースコールサーバ（Blade Center）内で行われます。

1次解析データは、HDF5という階層型のデータフォーマットと、生のFastq、Fastaデータが含まれ、これらは自動的にユーザのストレージサーバに転送されます。　

次に、WebアプリであるSMRT Portal（中央下）から、ユーザは、SMRT Analysisにアクセスします。

SMRT Portal上で、ユーザは、ストレージサーバに転送された1次解析データをクエリに、2次解析（MappingとかAssemblyとか）の設定をします。

SMRT Analysisはストレージ内の1次解析データを呼び出し、計算をして、BAMやSAMなどの結果を指定した場所に書き込みます。

ちなみに2次解析のコマンド群をSMRT Pipeと呼んだりします。　

GUIで行いたい方向けに、SMRT Portalがあるんですが、CUIがいい！ってひとはこのSMRT Analysisをコマンドべースで使うのももちろんOK

SMRT Analysisの解析メニューは、

1. BLASR（ブレイザーと発音）というPacのロングリードに向いたアライメントツール
2. GATKのGenotyperを使用したSNP検出ツール
3. ALLORA（アローラと発音）という名のアセンブラー
4. ギャップフィルタリングやScaffoldingをするAHA（アハッ）
5. そしてMethylation 検出

の5つが柱。　5のMethylation検出は今度のバージョン1.3.1から新たに加わります。　

これらについては、そのうち、例を出しながら書く予定。

さて、SMRT Analysisはユーザがサーバにインストールするわけですが、どんなサーバを用意すれば良いのでしょうか？

OSはUbuntu 10.0.4以降、CentOS 5.6以降

MySQL、bash shell、Perl v5.8.8以降、Perl XML parserがインストールされていることが必要です。

ハードウェアの最低条件

ヘッドノード：　16GB～32GBのメモリ、250GBのスペース

子ノード（5つ）：　コアあたり2GBメモリで、ノードあたり8コア、250GBのスペース

ストレージ：　約10TB

これに満たなくても動くことは動きます。

ただ、遅いかも。

メモリは多いに越したことはないけれど、ストレージのディスクを速いのにしたほうが、実行速度は高まるらしいです。　

私は、感覚的にはそうかなあ、って感じですが。

ということで、ユーザがインストールしなくてはならないのはRS RemoteとSMRT Analysisの2つ。

RS Remoteは普通のWindows PCで動くで問題無いが、SMRT AnalysisはそこそこなスペックのLinuxサーバが必要。

特に、大きなゲノムサイズに挑戦するには、メモリは多い方が良いです。

現場の会

先週、「NGS現場の会・第二回研究会 in 阪急エキスポパーク」に行っていました。
私はその前日チュートリアル、というところで1時間時間を頂き、喋ってきました。
PacBioについて、その原理と特徴、PacBioを使ったアプリケーション例などなど。

プレゼンの前の日は大阪入りして、ホテルで夜3時までスライド作って練習していました。
当日は、最初は緊張してレーザーポインタが震え・・・

私の前はRoche、Illumina、LifeTechと続き、そしてPacBio
それぞれのシーケンステクノロジーの特徴が一日でわかる！って感じですばらしい。

終わってみて、話したいことは大体伝わったかな、という感じはします。

あと、 ブログを知っているひとはかなりいました。　最初、みんなに質問したのですが、300名弱ほどオーディエンスがいたそうですが、3分の1くらい手を上げてくれました。

こんなにマニアックなトピックなんですが。

うれしいことです！

いろんな参加者と話ができて、夜はお酒も入って、ゆるーい感じの学会？です。

ゲノム、メタゲノム、リシーケンス、計算技術、IT、メーカー・・・などなど、その筋の方々、特に現場の方々がそろってわいわいやる。

こういうの良いですね。

この業界は意外と横のつながり、というかメーカーや企業の壁を越えての知り合いは多いのです。

競合の会社でも、社員（現場）のレベルでは飲みに行ったりとか。

でも、これだけたくさんの研究者が一同に集まると、もっともっと輪が広がり、新しい出会いがありました。

前のブログ「ショートリードの憂鬱」を知っているひとから声をかけられ、一気に50人くらい新しい知り合いが増えました。

現場の会に感謝です！　

P.S.

私のプレゼンだけ、白黒コピーでした。　他の企業はカラーできれいだったのに。

これには理由があるんです。

23日の午前2時からKinko'sのサイトがメンテナンスになり、ウェブから印刷を注文できなくなったので、3時に書き上げた私は焦り、とりあえず4時前に寝て、朝8時ごろ北梅田のKinko'sまで20分かけて歩き、印刷を注文したんです。

白黒プリンターは2台あるので印刷が速いと言われたので、時間に追われて白黒印刷にしたわけです。

しばらくブログ更新をしていませんでしたが、これから行きますよ！

シンプソンズとPacBioの深い関係

アメリカでは結構人気のあるシンプソンズ。
コメディとしてはちょっと過激、というイメージ。
日本では、知らないひとが多いでしょうね。　私もPacBio にたずさわるまでほとんど知りませんでした。

何でPacBioとシンプソンズが関係あるかって？

それは、PacBioのシーケンサーを操作するプログラムの名前が、全部シンプソンズの登場人物名になっているからなんです!!

例えば、シンプソンズ一家の、ホーマー（父）、マギー（赤ちゃん）、バート（男の子）は、Pacエンジニアが頻繁に使う、あるいは良く耳にするプログラムの名前です。
シーケンサーを前に私たちがそんな言葉を発していても、「何言ってんだ？」と思わないで下さい。
ちゃんと仕事をしているので。

ロボットを操作したり、ベースコールサーバーにアクセスしたりと、ユーザが触ることはないプログラムの名前です。
でも、操作するプログラムの名前が「波平」、「カツオ」、とかだったらすごく恰好悪いのに、アメリカ人にはそういう感覚無いのかなあ。
それともシンプソンズがよっぽど親しみがあるアニメなのだろうか。

ちなみに、ユーザが触る2次解析サーバでも、少し登場します。　
インストールディレクトリのことを、SEYMOUR_HOME　と呼び、このSEYMOURはシンプソンズでは市長のこと。
自動で2次解析が始まるデーモンをKodosと呼び、これはエイリアンの名前。

最初、これらのキャラクター名の意味が分からず、Pacに聞いてました。
相手は、私がシンプソンズを知らないとは知らず、しばらーく無駄な会話のやりとりを続けたのでした。

今日は軽い話でした。




-------------------------------　おまけ　-------------------------------------

ところで、シンプソン一家が住んでいる街の名前は、Springfield です。
この名前も、Pacのあるものの名前なんです。
ここでクエスチョン！
その「あるもの」とは、何でしょう？

Pacの一次解析 その２

パルスデータまでできたら、次はベースデータです。
どうやってパルスからベースコールを出しているのか？

これは、「観測されたパルスがOの時に、そのテンプレート配列がTである確率が、最も高い、Tを、そのテンプレート配列の塩基とする」　ということにしています。
回りくどい言い方ですみません。　数学が得意な方はこちら

あ～あ、そういうことね。　ってわかった方はすごい！
つまり、パルスがこれくらいだったらこれはAである確率が一番高い、とかTである確率が一番高い、とかそういうパターンアルゴリズムがあって、それに当てはめているのです。

そこんとこもっと詳しく！　って思った方、しかしながら、トレースファイルからパルスデータ、その後の一連のベースコールアルゴリズムは、公開されていません。
ですので概念としてしかお伝えできません。　あしからず・・・



さて、PacBioのデータの特徴、エラーの特徴にはどんなものがあるのか。
トレース、パルスデータにヒントがあります。

以前、PacのデータはInsertionエラーがある、と書きましたが、それを含め、Pacのエラーはどんなときに起こるのか？

エンジニアに聞けばまた別の答え（カメラやレーザーだとか、窒素濃度だとか）があるのですが、分子生物学の立場からだと、その答えはポリメラーゼにあります。
（レーザーによるエラーはまた別の特徴があり、窒素濃度によるエラーもまた別の特徴があるので、InDelエラーはほぼポリメラーゼによる、と考えて良いでしょう）
DNAポリメラーゼは、人工的に変異を施しています。　
合成スピードをカメラの技術に合わせて、ダウンさせています。
CellのZMWの底にBiotinでくっつくようにしています。
さらに試薬のpHに耐えられるようにしています。
また、酸素は有害です。
そして何より、自身が蛍光を放つたびに、光子が出すエネルギーに絶えずさらされています。　これによるダメージも大きいのです。

以上の環境下では、ポリメラーゼ君は時に疲れて休んでしまうこともあるでしょう。
ちょっと休んでまた働くこともあれば、そのまま動かなくなってしまうことも。
されにサボって、合成しないときも。

1. サボってDNAを合成しないで、そのまま次に進んでしまうと　・・・・　その塩基ではパルスは検出されないので見かけ上Deletionになります。
2. 合成はしても、疲れて休んで動かないと、蛍光を外さないので2個＋分読まれてしまうことも　・・・　その塩基では2個分読まれてしまい、見かけ上のInsertion
3. さらに2個分の蛍光が、テンプレートと同じときと違うとき(間違いパルス)の2種類あるのでInsertionにも2種類あるのです。
4. 合成はしたけど間違ったパルスが検出された　・・・　これはMiscall

というわけで、ポリメラーゼの働きがInDelエラーに大きく影響を与えるのです。
もちろん、同じ蛍光が2つ並んだ時が、ポリメラーゼ君のお休みによるものなのか、真にホモポリマーによるものなのか、その見分けは、簡単ではありません。
モデル式を用意して、当てはめて、区別しているのです。
単純に偶数倍のパルスの幅があるからそれで割る、みたいなことよりもずっと賢い計算をしているそうです。

Pacの一次解析 その１

今日はPacの一次解析について書きます。
NGSではおなじみですが、ご存じない方のためにちょっと説明すると、一次解析、二次解析、三次解析というのあって、それぞれ、

1. 一次解析：　シーケンサーから出てきた生データから、塩基配列を求める（ベースコールする）
2. 二次解析：　塩基配列になったリードデータを、アセンブルまたはリファレンスにマッピングする
3. 三次解析：　マッピングされたデータをもとに、Contigをつなげたり、配列から何か意味のある現象（例えば発現量や変異など）を得る

という意味です。

一次解析は、最初の基本的なところなので、データの精度や結果に大きく影響しますね。
以下は全て、シーケンサーの横にあるBlade Center、というベースコールサーバーで行われます。

PacBioのシーケンサー生データは、Movieです。
Movieはとてつもなく大きなデータで、1秒間に3Gbも出てきます。
これはメモリ上に展開されるだけで、すぐに圧縮されてトレースデータ、という形で保存されます。
この、Movieからトレースデータへ、「Movie-to-Trace」 という部分が、Real Timeで行われているのです！
SMRT Cell のSingle Molecule "Real-Time" ですね。

Movieから変換されたトレースデータは、それでも1Cellあたり、50～150Gbもあります。　（Movieの時間によってサイズが異なります）。
トレースデータは、各ZMWから放出された、計測したPhotonの量（光量子束、photon flux）を数値化したものです。
各カメラから計測されたPhotonの推定量と、蛍光ごとによるPhotonの推定量とを使って、何やらノイズを取り除くノーマライズをしています・・・。
これ以上説明せよと言われてもできません。　ごめんなさい！

まあ、ざっくり言うと、動画から蛍光ごとに区別した波形データを作っているんですよ！

で、このトレースデータは、まだ塩基ごとの波形がばらばらなので、もうちょっと見やすく、さらにノーマライズしたのがパルスデータです。

トレースからパルスへ、「Trace-to-Pulse」変換が行われます。
この変換の最終目的は、ポリメラーゼがDNAを読んだ時の、パルスの正確な検出と、他との区別をすることにあります。

トレースデータはまだ、塩基ごとの区別があいまいなんです。　それを、閾値を設けてノイズをできるだけ取り除き、ほかの塩基としっかりと区別できるように変換したのがパルスデータなんです。

パルスデータには、重要なKinetics情報が含まれます。例えば、

1. パルスの高さ　Intensity, or Pulse Height
2. パルスの幅　Duration, or Pulse Width  (PW)
3. パルス間の時間　Spacing, or Interpulse distance  (IPD)
4. 塩基の種類　Content, or Base identity

などです。

次にようやく、パルスデータからベースデータへの変換、「Pulse-to-Base」が行われます。

Menlo Parkにて

久々の投稿です。
ちょっと最近、忙しくて書く暇が見つからなかった・・・。　
ネタはいっぱいあるんですけれど。

そういえば先週、Menlo Parkへ行ってきました。

カリフォルニア州のベイエリアのほとり、ちょうどFacebookの本社があるすぐそばに、PacBioの会社はあります。
社員およそ300人。
開発に携わるグループは、分子生物学、光学、計算科学、など各分野からの天才がそろっています。
R&Dセクションでは、常に新しい酵素のテストが繰り返され、今まで数千を超える酵素がランされたそうです。
工場では、パーツごとに外注で作られた部品が、慎重に組み立てられていました。

私が行ったところは、主にバイオインフォのセクションと、R&D、それにITシステムのセクションです。
ポリメラーゼの開発は、特に興味深いものでした。　
Publicには書けないのが残念です。

Pacには、フィールド・アプリケーション・スペシャリスト（別名FAS）という仕事があります。
私も、実際、Pacの分け方ではFASの一部で、FASのバイオインフォ担当、ということになるそうです。
そんなこんなで、全米各地から20～30人のFASが集まるトレーニングに参加してきました。

日本ではまだPacのシーケンサーが入ったばかりなので、正直言って、どんなサポートが必要なのか、未知のことが多いです。
アメリカではもう1年以上経験があるので、ユーザーがどんな質問をするのか、どんな要求をするのか、そこは参考になりました。

IT、ソフトウェアに関しては、私も前職で経験がありますが、改良すればバグがつきものです。
バグは見つけたら直す、ということの繰り返しです。
なので、私は最近、バグを見つけることにやりがいを感じています。
幸い致命的なバグは見つかっていませんが。

バイオインフォ関係は、以前もちょこっと書きましたが、アセンブリにCelera Assembler　を使うことができます。
Celeraをやってた技術者がそのままPacに転職してきたから、使うのに慣れていた、ということもありますが、ロングリードはCeleraが向いているそうです。
そのほか、マッピングやアセンブラは、まだまだ開発の余地が残っています。
開発というよりも、既存のマッパーやアセンブラーを使うか、少し改良して使う方が、精度を上げることができるように感じました。
実際、BWAでもちゃんとマッピングできます。
CLCでもできるかな？

さてさて、Menlo Parkというところは、Facebook効果もあってか、最近物価が上がってきているそうです。　
ミリオネア（つまり大金持ち）がたくさん引っ越してきているそうです。
近隣の町のダウンタウンも、いい感じのレストランが多くありました。
そして、ホームレスがいない。
ちょっと先のサンフランシスコとは大違いです。

Facebook本社前の「いいね」看板

ゲノムアセンブリ メソッドの種類

3月10日～12日まで立教大学で開会されていた日本ゲノム微生物学会に行ってきました。
11日のランチョンは、PacBioのプレゼンということで、満席でした。　良かった良かった。

 

バクテリアシークエンスについて講演して頂いたPacBioのPaulさんは、今でもピペットマンを持つ現役のサイエンティストです。　
講演を聞いた方はご存じでしょうが、PacBioのロングリードは平均15%のエラーがあります。　そこでショートリードでマッピングして、エラーをできるだけ取り除き、精度を増したロングリードでアセンブルすると、それ以外のどの組み合わせよりもN50、最長Contigなどで成績が良かったとのことです。

このグラフ↓は、その時のプレゼンでも紹介されました。

CSSというのはcircular consensus sequence のこと（CCSのこと）です。
PBcRは、PacBio correction Read で、エラーコレクションがされたロングリードです。
左から3つは、PacのCCSのみ50X、IlluminaショートリードPaired Endの100X、454の50Xです。
リード長はどれも数百塩基です。
これらはどれもN50が100Kb、Max 300Kb程度ですね。

次いでCCSでエラーコレクションしたロングリード（どちらも25X）、454でエラーコレクションしたロングリード（どちらも25X）。　この2つはほぼ同じような成績で、N50が200－300Kb、Max 600Kbでした。

Illuminaでエラーコレクションしたロングリード（どちらも50X）はN50が300Kb、Max 1Mbで、一番良かったのがCCSでエラーコレクションしたロングリードの50Xで、N50が500Kb、Max 1.2Mb

細かい数値やメソッドは、近く論文として紹介できると思います。


Pacのデータを使ってできるアセンブリの種類ですが、大きく3つ、あります。

1. SMRT Hybrid
2. SMRT Scaffolding
3. SMRT de novo

SMRT Hybridは、PacBioのロングリードをCCSやショートリードでマッピングして、エラーを取り除き、精度の高いContig/ロングリードを作っていくプロセスのこと。　
SMRT Scaffoldingは、精度が高くなったContigやロングリードを使って、Contig間をつなぐ（Scaffoldingする）こと。
SMRT de novo は、Pacのロングリードだけをそのまま使ってアセンブルすること。

です。

SMRT Hybridは一番クオリティが高いアセンブル手法です。　Illuminaや454データを持っている場合やCCSを組み合わせて実験したデータを持っている必要があります。
ゲノムサイズに制限の無いCelera Assemblerがお勧め、とのことです。
ALLPATHS-LGというアセンブラーもありますが、こちらは10Mbまでという制限があります。
制限といえば、Celera Assemblerも、インプットするロングリードの個々のリード長が30Kbを超えてはダメという制限がありますが、今のところPacのロングリードも30Kbは読めませんので問題無いでしょう。
将来的には・・・　30Kb読めてしまうかも

SMRT Scaffoldingは、ショートリードで作ったContigがある場合や、PacのCCSなどで作ったContig、Pacのロングリードをエラーコレクションしたリードがある場合、それらのContigの精度が高いと見なされる場合に、Contig同士をブリッジするのに有効です。
（Scaffoldingのちょうど良い日本語訳が見つからない！）

SMRT de novo は、ゲノムサイズが1Mb未満の小さいゲノムを、CCSを使わずに、ロングリードだけで読んで、クオリティは若干低くてもいいからドラフトゲノムを作りたいときなどに向いています。　
ウイルスなどは少ないCellでもカバレージも十分とれるので良いかもしれません。

以上のアルゴリズムについては、後々、書いていくつもりです。　っていうか、
書かざるを得ないでしょうねえ。

ターゲットリシークエンスと変異 3 - SureSelect

Agilent Technologies社のSure SelectというターゲットエンリッチメントシステムについてはショートリードのNGSでもおなじみ。
さて、これをPacBioのロングリードでやってみるとどうなるのでしょうか？

このポスターでは、2つのHapMapサンプルと、Tumor/Normalサンプルを使って、SureSelectでターゲットをキャプチャーし、それをPacBioで読んでいます。　
（直リンクを貼るわけにはいかないので、"Development of SureSelect Target Capture Methods for Sequencing on the PacBio RS"　をGoogleで検索してAgilent社のPDFをゲットして下さい。PacBioのサイトからも取れます。）

M&Mを簡単に言うと、1マイクログラムのgDNAをCovarisを使って250bpまたは2kbに切断した後、SureSelectシステムの中でベイトRNAとミックスし、0.5Mb Chr10＋X（Contig + Exons）領域をキャプチャーします。　その後はPacBioのライブラリー作成プロトコールに従い、250bp、または2kbの2種類のライブラリーを作っています。

シークエンスは45分ムービーを2回、読んでいます。　つまり、250bpのライブラリーはCCS(Circular Consensus Sequence)で何回もぐるぐる読み、2kbはCLR(Continuous Long Read)で読んでいることと同じです。
実際何個のSMRT Cellを使ったのか、という情報が欠けているのですが、CLRで200X～500XのSub-readのカバレージを得たそうです。（http://www.pacificbiosciences.com/applications/target　の、Technical Note: Targeted Sequencing on the PacBio RS using Agilent Technologies SureSelect Target Enrichment）
しかし実際は、ここまでカバレージを取らなくても良いでしょう。　その例はまた今度お見せします。


さて、250bpと2kbの2種類の塩基長でキャプチャーした場合、結果はどう違うでしょうか？
まず、マッピングの結果、2kbのほうが、250bpよりも均一にゲノムをカバーしました。
隣合うキャプチャーベイトのプローブをまたいでマッピングできたのも2kbのほうです。
DNA断片が2kbだと、ターゲット領域の上流・下流の配列も長くとれてくるので、のちのデータ処理では、ターゲット領域上に取れてきた配列のみにフィルタリングしたい場合、注意が必要でしょう。
話はずれますが、ショートリードの時は、

1. リードをゲノムにマッピングし　（BWA）
2. 冗長性のあるリード・Duplicateを除去し　（SamtoolsやPicard）
3. キャプチャーしたExon領域だけを取り出し　（Bedtools）

という流れが来ると思いますが、Pacのデータの場合、データ量が（HiSeqとかに比べると）少ないので、キャプチャー領域だけを取り出すのはもったいない気がします。個人的にですけど。

さて、それからSNPを見てみると、2kbでは、250kbでは見つからなかった既知SNPも見つかったそうです。　
上が2kbで下が250bpで、山の形はカバレージを表しています。　某遺伝子の3’-UTRの部分が拡大されていますね。　
青い囲みは2kbと250bpの両方で見つかった既知SNP、赤い囲みは2kbでのみ見つかったSNPだそうで。

こういうのがたくさん見つかるとしたら、2kbのライブラリーを作ってみたくもなるでしょうね。

あと、長いライブラリーのほうがGCバイアスにも強いし、ギャップも埋められる、という利点があるので、ショートでは見つからなかった変異も見つけることができる！　って、言われますが、まだインパクトのある実例が無いのが残念!!

で、結局、SMRT Cellは何個必要なの？

この答は、あくまで計算上でしかわかりません。
実際に何個使って、どれくらいのカバレージで、新規SNPを見ーぃつけた！　という論文やポスターがあれば参考になるのですが。
今年を期待しています。

さて、計算上はどうでしょうか。

例えば、10Mbの、ターゲットキャプチャーをして、その場所の20%の変異を見つけたい場合。
2kbのライブラリーを使うとして、最低必要カバレージは150Xと想定（前回のブログ参照）。
ターゲットの数は、単純に10Mbを2kbで割った5,000というわけではありません。
プローブがいくつ設定されているかは、ケースバイケースなので、ターゲットの数はとりあえずN としましょうか。
2kbライブラリーで「使える」CRL subreadは、1SMRT Cellあたり50,000とします（前回のブログ参照）。
サンプルバイアスを3とします。（3倍ぐらい多めに読んだ方が安全ということ。上記SureSelectの例でも200~500Xのカバレージを得たとありますので、最低カバレージ150とすると、3倍くらい多めに読むつもりでやっています）
するとSNP検出に必要リード数は、3x150xN= 450N   
必要なSMRT Cellの数は、これを50,000で割って、0.009N

ターゲット数が1,000ならSMRT Cell 9個
5,000なら45個
10,000なら90個

・・・・・・

結構な数ですよ、これは

ターゲットリシークエンスと変異 2

どれくらいのカバレージがあればSNPを検出できるのか、というのは興味のあるところだと思います。
以下はPacBioテクニカルノートに記載している方法です。　（うまく説明できていないかも）
同じゲノム位置の、ベースコールの集合（どれも同じQV）を2つ考えてみます。　1つは100%ワイルドタイプで、もう1つの集合は（50:50または80:20の）変異を持っているとします。　ランダムに2つの集合からベースコールをピックアップするとき、その数はそのゲノム位置でのカバレージと同じ意味です。　
2つのピックアップしたサンプルを、Fisherの正確度検定にかけて、どれだけ有意に変異を検出できるか（帰無仮説は、ワイルドタイプの集合からのサンプルも変異がある集合からのサンプルも、差はない）検定をします。

p-valueが0.05未満のとき、有意に検出できるとし、ピックアップする数を1から50に増やしながら（つまりカバレージが1から50のときまで）1000回繰り返します。
このようにしてシュミレーションした結果、(p-valueが0.05未満になるときの数)/1000 = 検出率とすると、0.95以上になるには、
50:50のサンプルでは、CCS-3パス（Q17）のとき17X、CLRフルパス（Q8.24）のとき29X必要で、
80:20のサンプルでは、CCS3-パスのとき55X、CLRフルパスのとき135X以上必要であることがわかりました。　
この数値はもちろん理論値です！

さて、以上をふまえて、ターゲットサイズ（アンプリコンのサイズ、インサートのサイズ）ごとに、SNP Detectionに必要なSMRT Cellの数を計算した表があります。
計算式はこちら

で、tはターゲットの数、cは最低必要なカバレージ（先の理論値より若干高めに設定）、βはサンプルバイアス（N倍くらい多めに読まないと期待される変異は取れないんじゃないか、という値）で3に設定。　tcβで必要リード数が計算されますので、それを前回の書き込みで紹介した、Cellごとの、ある程度クオリティが期待される「使える」リード数=rで割ると、必要Cell数が計算されます。

CLRフルパスの場合、500bp、1000bp、2000bp、の3パターンで、ターゲットの数はそれぞれ500、2000、1000です。　CCS-3パスの場合は、500bpと1000bpの2パターンでターゲットの数は500と2000です。　両方とも変異の割合を50:50と80:20で出しています。
（Table4の20% variantの500bpの、Estimated # SMRT Cellは2ではなくて3が正しい）
この2つの表からわかることは、500bpのターゲットでは、ターゲットの数が500のときCLRとCCSはどちらも1つくらいしかCellを必要としない、ということ。
でも、この表には無いですが、ターゲットの数を500から1000、1500、2000と増やしていったらどうなるのか、上記の式に代入して計算してみました。

ターゲットサイズ500bpのときと1000bpのときの2パターンでCLRとCCSのどちらがコスパがいいか？
ターゲット数が増えるほど、CLRでやったほうが、変異割合が50％のときも20%のときも少ないCell数で済む、と言えそうです。 
Again, あくまで理論値です!

では実験値ではどうなのか。
Agilent SureSelectを使ったターゲットエンリッチメントの例があります。

続く

ターゲットリシークエンスと変異 1

ゲノムの変異、広い意味での変異を見つけることは、次世代シーケンサーならではのテーマでしょう。

オーソドックスなのものでは一塩基変異（Single Nucleotide Variant）を見つけること。
ショートリードの世界では例えば、BWAでリードをゲノムにマッピングしてからSAMToolsでSNVを検出する、というようなパイプラインがほぼ確立されているので、多くの研究者の方がチャレンジしているのではないでしょうか。

SNVがコーディング領域のアミノ酸を変えればタンパク質の機能を、遺伝子の転写制御領域にあればmRNAの転写を、直接変化させます。

他にも、遺伝子や重要配列をごっそりDeleteしてしまう変異もありますし、何コピーもAmplifyしてしまう変異もあります。

あとは、染色体単位でのTranslocation、転写・遺伝子レベルではキメラRNAやFusion Gene、新規のAlternative Splicingなんかも変異になるのでしょうか。

ショートリードでは、SNVはそのまま検出可能ですが、大きなDeletion / Insersion、Amplificationは難しいでしょう。　
リードは長くても100－150bp（HiSeq、MiSeqなど）です。
マッピングアルゴリズムにもよりますが、10塩基以上のInDelを検出できるツールを私は知りません。（あったら教えて）　

そのSNVですが、ショートリードの世界では良く、カバレージは最低何X必要などと言われます。
これについてはいろいろな意見があると思いますが、以前学会や展示会で聞いて回ったところ、普通のDiploidのジェノタイプ（つまり50:50）には20X、ちょっとレアなSNV（80:20）には100～120Xは必要だと言われた記憶があります。
根拠が薄くてすみません。　でも今のHuman Exomeなどでは、HiSeqやSOLiDを使えば普通にこれくらいのカバレージは「平均で」出ますね。


では、PacBioのシーケンスではどうでしょうか？
SNVを見つけるには、どれくらいのカバレージが必要なのでしょう？

このテクニカルノートは最近出されたものです。
面白かったのでサマリーを紹介します。

まずは最初に単語の説明から
Continuous Long Read (CLR): 長いインサートを一本で読んだ時のリード
Circular Consensus Sequencing Read (CCS): 短いインサートを数回読んだ時のリード

CLRとCCSにはどんな特徴があるのでしょうか？

フルパスで一回読んだだけのCLRは、精度85－90%、塩基あたりのQVは8.25－10で、平均リード長は2kb、これを10XのカバレージにするとQVは30以上（99.9%以上）になるそうです。　

一方CCSは、2回パスを読んだものでリードの精度が97%、3回パスで98%、5回パスで99%ということになっています。　5回読めば1kbあたりのエラーは10塩基未満に抑えられるということ。
塩基のQVも、2パスで13、3パスで17、5パスで20以上となり、なんとか耐えられます。
10XのカバレージにするとQVは50以上になるそうです。

これと、重要なのが、SMRT Cell 1つで出力される「使える」リード数です。
「使える」リードとは、CLRでは1回のフルパスサブリード、CCSでは3回以上のパスで出力されるリード、と定義してみます。
45分のムービーを2回撮ったとき（C2ケミストリー）の設定です。
CLRは、500bpインサートで250,000本、1kbインサートで100,000本、2kbインサートで50,000本。
CCRは、500未満bpインサートで40,000本、500bpで30,000本、1kbで15,000本。

ここで問題：　あれあれCLRの500bpインサートが25万本を出力しています！　ZMWは15万/Cellしかないのに。

答え：　サブリードですので、1 SMRT Cell あたり2以上を出力している（かもしれない）のです。
CLRで500bpのインサートを完璧に1周したら2本のサブリードが出ます。
45分のムービーなので十分ありえます。

続く・・・