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

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分のムービーなので十分ありえます。

続く・・・

以前のコメント

以前のMovieってという記事にて、コメントがありました。　ありがとうございます。
気づくの遅すぎですね。　すみません。　今度からコメントもチェックしようっと。


①何故1度の撮影でZMWは75,000ウェルしか観測できないのか、照明系かCMOSサイズの問題だろうか。
A:　なぜでしょうね。　もっと大きなSMRT CELLをつくればいいのに。　
そういう仕様だと聞いて今まで納得していました。　
カメラサイズの問題だと思いますが、今度Pacのエンジニアに聞いてみます。

②ポリメラーゼの合成速度について、ポリメラーゼ改変体を使っているとのことだが、ポリメラーゼの基質として蛍光ラベルdNTPを使っているため、通常のdNTPより合成速度が遅くなると思うだが、その速度への影響はどの程度なのか、また蛍光ラベルdNTPを用いたことによってDNA合成の正確性は変化するのか。
A:　蛍光ラベルdNTPも確かに速度に影響しますが、私が聞いた説明ですと、カメラの性能が追い付いていないため、わざとポリメラーゼを遅くしている、ということでした。　ここも確認が必要ですね。
蛍光dNTPがDNA合成の正確性に影響を及ぼすか、というのは特に聞いていません。　ゼロではないでしょうが。

一つ前の記事について
③下記表現について
>「ポリメラーゼによって合成される1本のリードに含まれる塩基が、平均して85%は正しく読まれる」　
この15%の高いミスリーディングは何に由来するものか、また欠失・重複等はどの程度含まれるのか。
A:　Hybrid Assemblyの記事　にも書きましたが、挿入エラーが多いのがPacの特徴です。　パルスデータがはっきりと正しい数の塩基に変換されないのでしょう。　
どこに原因があるのか、ムービーからトレースデータに変わるところでそもそもエラーが起きるのか？　
トレースデータからパルスデータに変わるときにエラーが起きるのか？
何か特定の配列が読みづらいのか？
これは私も前から興味を持っていたのでいつか書きます。　

④
>確率的には1/3の割合（厳密にはポアソン分布）で、正しく合成が行われる状態のZMWができるわけです。
場合の数が3通りあるだけで、その３通りが等確率とは限らないと思うだが、この値は理論値か、それとも実験値か。
A:　理論値です。　仰る通り、実験すると33%よりずっと低い値になることもあります。　

今日はもう遅いので。ここまで。

New Chemistry ! - C2 とは

PacBioのシーケンスアプリケーションについて報告です！

現在、インストールしているシーケンサーは、C2という新しいバージョンです。
これからインストールされるものは全て、C2バージョンになります。
世界的にも、すべてのマシンが次々にアップグレードされていく計画です。

今までのC1バージョンとどこが変わったのかというと、

1. より長く読める（10Kb のライブラリーを読むことができる）
2. より多く読める（ZMWの穴に1つちゃんと入る「DNA‐ポリメラーゼ」が多くなった）
3. 解析ソフトウェアがちょっと良くなった
4. ソフト上でデータのグループ管理ができるようになった

ポリメラーゼ酵素がより安定するように改善されましたので、より長く、より安定して読めるようになったのです。
それに伴って、今までStrobeという、長いライブラリーを飛び飛びに読むプロトコールが無くなりました。　
ショートリードでいう、メイトペア、のかわりにStrobeがあったのですが、もうそのまま10キロ読んじゃえ！ってことでこれは無くなったのです。
実際どれくらい長く読めるのか？　というのが気になるところでしょう。
10Kbのライブラリーを作って読んでも、すべてのリードが10Kbというわけではありません。
平均は数Kb、何割かが10Kbいくかどうか。
正確な数値を伏せているのはまだオフィシャルになっていないかもしれないので。
近くお知らせします。

解析ソフトについては、細かい機能がちょこちょこアップデートしたり、バグが直ったりと、いまいち地味ですが、使い勝手はいいと思います。
各セルごとに、塩基ごとのシグナル:ノイズ比がわかるようになりました。
使ったことはまだありませんが、GATKとのインテグレーションができるようになったそうです。

もう少ししたら、カタログもプレゼンも、すべてC2用に変わります。
今までC1で覚えたことが、C2ではいろいろ変わっていて大変ですが、すぐ慣れるでしょう。

そんなところで、明日も頑張ります。

de novo Hybrid Assembly

参照配列が無い生物の配列決定や、転座などの構造変化を見つけるために欠かせないのが、de novo Assembly です。
NGSでde novo Assembly といえば 定番中の定番かもしれませんね。

異なるプラットフォームのリードデータを合わせてアセンブることを Hybrid Assembly と言います。
例えばRoche FLXのロングリードと、イルミナGAIIxのショートリードを一緒にするAssembly のことです。
なぜこんなことをするかというと、ロングリードとショートリードにはそれぞれ一長一短があって、合わせることでそれぞれの短所を補いつつ、長所を生かしてより長い、より正確なContig を作れるからです。


先日、The Plant and Animal Genomes Meeting というのがサンディエゴでありました。　実際行ってはいませんが、そこで Cold Spring Harbor Laboratory の Michael Schatz 博士が、PacBioを使ったHybrid Assembly の講演をしましたので紹介します。
プレゼンのスライドは、Schatz博士のサイトにアップされていますので興味のあるかたはどうぞ（ここをクリック！）。

Schatz博士は、PacBio だけで酵母ゲノムを読んだところ、エラーに特徴があることがわかりました。
それは、
ミスマッチエラー（1.4%）よりも挿入（11.5%）や欠損（3.4%）が多いということでした。
全体的なエラーは83.7%で、Pacとしてはまあまあの値です。
ちなみに65個のSMRT Cellを使用して、734,151本のリード（フィルタリング後）を得ています。
12メガbpのゲノムで、です。

Schatz博士のプレゼンより引用

Pacのリードの特徴は、ご存じの通り、長く読めることです。　最新バージョンでは運が良ければ最高10kbまで読める!? かも。
長く読めるということは、大きなギャップやリピート配列も読める、という利点があります。
GCバイアスが無いということも相まって、読める場所と読めない場所に偏りが少ない、という利点もあります。

しかーし、出力される本数が少ないですし、エラー率がちょっと高め? なので、そこはハイスループットのショートリードシーケンサーたちの出番となるわけです。

そこでHybrid Assembly が必要になります。

では、Schatz博士の話に戻ります。
昨年、他のポスターにも発表していましたが、Schatz博士は以下の手順でパイプラインを組んでいます。

1. ショートリード配列をトリムし、フィルタリングする
2. きれいになったショートリードをPacBioのロングリードにマッピングする
3. カバレージとギャップに基づいて、明らかに読み間違いと思われる場所はトリムする
4. それぞれのロングリードw/ショートリードマッピングをもとにコンセンサス配列を作成し
5. エラーが無くなったリードを使ってアセンブルする

近いうちに詳細は論文で明らかになりますが、面白いのは、アセンブラーに、Celera Assembler というクラシカルなもの（の改良版）を使っていることです。
Sanger からNGS に行き、また戻ってきた懐かしい感がしました。
また、上記、エラーを除去しているステップは、pacBioToCAというオープンソースで公開されています。

Schatz博士は講演の最後に、高等生物のde novo Assembly についても触れています。
ゲノムサイズ1.23Gbpのオウムです。
3パターンのアセンブルを試していました。

1. Illumina の194X （インサート長が220、500、800bpのペアエンドと、同じく2K、5K、10Kbpのメイトペア）
2. Roche 454 の15.4X（FLXとFLX Plus、3K、8K、10Kbpのペアエンド）
3. Roche 454 の15.4X とPacBio 3.75X

IlluminaのみだとContig数が24,181で、N50が47,383bp、最長Contigは1,050,202bp
RocheのみだとContig数が16,574で、N50が75,178bp、最長Contigは75,729bp
PacとRocheだとContig数が15,081で、N50が99,573bp、最長Contigは1,238,843bp
（個人的にはPacとIlluminaがほしかった・・・）

PacとRocheの組み合わせで、他と比べて、より長いContigが多く作られたことがわかります。
とはいっても、ゲノムサイズ1.23Gのうち、1.07Gしか読めていませんでした。
それにまだ、15,000以上のContigがつながらずにできてしまうのですね。
de novo Assembly は大変です。
他の生物（バクテリアとか）でも試されているようですので、詳しくは先程のリンクからSchatz博士のプレゼンＰＤＦの最後の方のテーブルを参照のほど。

これはその一部