シンプソンズと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を使ったターゲットエンリッチメントの例があります。

続く