今年最後は･･･Quiver と Arrowの話 我ながら地味

2016年、大晦日です。
紅白見ながら書いてます。PPAP、ゴジラ、真田丸、ブラタモリ。
今年は面白いな。

さて今年最後の投稿は、地味なネタです。
QuiverとArrowの違いであまり知られていないこと

GitHubのサイトでもあまり詳しくは説明されていないので、少しお話しします。

【はじめに】
QuiverとArrowは、どちらも、CCS（Circular Consensus Sequence）を作るアルゴリズムに由来します。
QuiverはRSIIのデータ用、ArrowはSequelのデータ用の、コンセンサス配列作製アルゴリズムです。
SMRT Analysis v2.3 まではQuiverが標準搭載されていて、v3.x からはArrowになりました。
v3.xは、SMRT Linkというパッケージソフトになったのですが、引き続きSMRT Analysisという名前も解析用ソフトとして使われています。
SMRT LinkでもP6-C4のRSIIのデータは解析できるよう、Arrowも対応されました。

【QuiverよりArrowの方が精度が高い】
CCSを作るとき、Quiverベースのアルゴリズムでは、何回パスを重ねても平均QVフレッドスコア３０くらいで精度に限界がありました。
また、ホモポリマーもうまく認識できないという欠点もあったそうです。
これは様々なクオリティ値を使ったエラーモデルが複雑すぎて余計ノイズになったとか。
そこで開発されたのがArrowというアルゴリズム。
トレーニングモデルやエラーモデルをより単純にし、またZMWごとの違いを考慮するようにした。（2016年アジアユーザグループミーティングより）

Arrowを使ったCCSは、パスを増やすとQVをフレッドスコア５０、６０に高めることが可能になりました。

Old = Quiver, New = Arrow

【QuiverよりArrowの方が使っているQV値は少ない】
普通、シークエンサーから出力されるリードの塩基には、その塩基をその塩基であろうとしたクオリティ値（QV）があります。
つまりAである塩基をAとコールしたとき、Cではなく、Gでもなく、Tでもなく、Aであったという「確からしさ」を表現する方法として、QVがあります。

RSIIは、ベースコールの結果、各塩基ごとのQVのほかに詳細QV (DeletionQV, InsertionQV, SubstitutionQV, SubstitutionTag, DeletionTag)という複数のQVを出力し、Quiverはこれらの情報も利用してコンセンサス配列を補正します。

Sequelは、これら詳細QVは出力せず、各塩基ごとの全体のQV（QUAL/FASTQ qv）のみを出力します。
（実データを見たことがあるかたは、現在のSequelリード BAMファイルではQVの情報が全て「！」になっているのに驚くかもしれません。この問題は次バージョンで改善される予定です）
ですがArrowではそもそも、塩基ごとのQVを使用しません！

ではArrowは何のデータをもとにエラー補正&コンセンサス作製をするのでしょう？
Arrowは、

1. 各リードのシグナルノイズ比
2. 各塩基のパルス幅の情報

この２つを使用して、マッピング後のコンセンサス配列を補正して作ります。
これらの情報は今もリードBAMファイルに記載が有ります。

P6-C4のRSIIデータも、Arrowを使えます。
このときは、RSIIデータ（bax.h5）をSequelのBamファイルに変換し、リードごとのSN比と各塩基のパルス幅情報を使って、Arrowが使われます。
結局、補正に使うデータを単純にしたほうが精度が上がったということでしょう。

今後、Arrowで塩基ごとのQV（詳細QVではなく、全体のQV）が使われるようになるかは、現在のところ何ともいえません。
精度が向上するようであれば、使われるようになる可能性はあります。

ああー、もうすぐ2016年も終わります。

しかし笑ったのは、白組の方がダントツで票が入ったのに何故か赤組が勝った、というアメリカ大統領選挙よりも摩訶不思議な紅白歌合戦。

明日の新聞で説明あるかな？

また来年も宜しくお願いしまーす！
２月７日の「PacBioセミナー＠秋葉原」も宜しく

Rocheとの別れ

年末になってまたまたビッグニュースが舞い込んで来ました。
2013年から続いていたRocheとの提携が、解消されることになったという。

いろんな噂が出る前に、ちゃんとした情報を取っておこうと思い、投資家向けテレカンを聞きました。
細かい内容は明らかにしてはいけないのでここでは書けませんが、これはプレスリリースにもあるとおり、Roche側から提携の延長をしないという判断がされたそうです。

Sequelの開発は、Rocheの資金協力のおかげもあったので、長い目でみたらお互いWin-Winの関係を築けていただろうに。そう思うと残念です。

とは言ってもなぜRocheが提携をやめる判断をしたのか、は、Roche側のビジネスの問題です。
3年前と今とではインビトロ診断領域のビジネス環境が変わった、というのがRocheの判断につながったのでしょうか。
装置を売るよりも、どのメーカーの装置にも対応したアッセイキットを作って売った方が儲かる！と考えたのなら理にかなう、かな？

でも、RocheのおかげでSequelができて、PacBioのターゲットが大きく広がったのは事実なので、ありがとうございました。

つい先日発表されたPacBioのrevenueも、大きく増えたし、Sequel装置もパフォーマンスの向上でだいぶ良くなった。

今後はどうなるのか？　
株は売られましたが、まだ大丈夫でしょう。

どこかと提携する可能性はあります。
個人的には、本社すぐ隣のFacebookとかが興味示してくれたらなー、と思っています。
シリコンバレーという場所柄、お金持ちが多いですからねえ。

2017年もやります! PacBioセミナー＠秋葉原

12月も中旬になると、東京も急に寒くなってきました。
毎年、分生が終わったこの時期は、忙しいはずなんだけど・・・「今年も良く働いた！」というへんな達成感がある。あとは飲み会？Year End Party

さて、来年ですが、2月7日（火）に、また秋葉原UDXで「PacBioセミナー」やります。
今回は午後からの半日セミナー。
なぜ午後からかというと、朝から5時まででは長い！という意見が前回あったから。
そこで今回は1時から5時までで、英語での発表はひとつ、PacBio本社からのみ。

NGS現場の会のメーリスにも流しましたが、

今回の名前はシンプルに「PacBioセミナー」。

なぜ「PacBio現場の会」にしなかったかというと

1. 同じ年にNGS現場の会があるので、社内で「現場の会の準備」と言っとき、「え？どっちの？」とややこしいから
2. お客さんに紹介したとき、「え？どっちの？」とややこしいから

いう理由でした。

PacBioセミナー　2017

ロングリードが世界を変えた！

【日時】

2017年2月7日（火）

午後1時スタート（受付開始は12時30分～）

午後5時半頃終了予定（その後懇親会あり）

参加料金は無料！ 

【場所】

秋葉原UDX

【事前登録制】

ここからかならず登録してくださいね（そのページの一番下に入り口あります）

日本にPacBioが導入されてからまもなく6年目を迎えます。（ということは私も今の仕事してから6年目ということか･･･）

バクテリアゲノムから大型真核生物ゲノムまで、様々なところでSMRTテクノロジーのスーパーロングリードは活躍してきました。今回は、ヒトゲノム、植物ゲノム、そして新型機Sequelにフォーカスしたセミナーを行ないます。もちろん、そのほかのアプリケーションやバクテリアゲノムも忘れません！

今回も、豪華な演者の先生方をお招きします。

• 芦澤哲夫先生（Houston Methodist Neurology, Director）；ショートタンデム反復の拡張と神経遺伝疾患との関連性の研究がご専門。米国NIHの組織するSCAコンソーシアムのリーダー
• 安田純先生（東北大学　東北メディカル・メガバンク機構　教授）；日本人基準ゲノム配列（JRG）をPacBioを用いて決定。このプロジェクトの中心メンバー
• 榊原康文先生（慶應義塾大学理工学部生命情報学科　教授）；11月にNature Communicationsに発表され注目を浴びたアサガオゲノム。その解析チームリーダーでバイオインフォマティクスがご専門
• 豊田敦先生（情報・システム研究機構　国立遺伝学研究所 特任教授）；ゲノム支援・先進ゲノム支援の枠組みの中でPacBio RSIIを様々なプロジェクトに活用され、日本で最もPacBioを使われているユーザのおひとり
• ほか交渉中演者数名

PacBio本社からは、アメリカでのPacBio応用例をご紹介します（Luke Hickey, PacBio）。

さらに、PacBioアプリケーションアップデートの発表や、今話題のDovetail Genomics（ゲノムアセンブリ専門の受託サービス）、Swift Biosciences（少量DNAからのショートリードライブラリ作製キット）の紹介、NGS関連商品の紹介などもする予定です。

今回もPacBioに興味のある現場の方や、ロングリードの可能性に期待する先生方が一同に集まる貴重な場になると思います。

是非この機会にご参加下さい。

しつこいですが事前登録制です。

このページの一番下から、Googleフォームへのリンクへ行って、レジストしてください。

入り口がわからない方は、ここから直接どうぞ

懇親会は、会場近場で18:00から予定しています。

是非こちらもご参加下さい。.

タイトル、アブストなどが決まり次第、またお知らせします。

来年は、手のひらサイズシークエンサーが話題になりそうですが、PacBioはまだまだ負けてませんよ！

ランニングコスト、データ量、精度、供給安定性、サイエンスとしてのデータ、いずれも客観的に見てPacBioの方が上回っている、ということをこれからもお示ししていきます！

なーんて、大きなこと言ってますが、アレですね、市場が活性化するのはいいことです。

同じロングリード同士、競いながらもうまく住み分け、するのかな。

これについてはまたの機会にしっかり書きます。

では

IGVでの格好いい見せ方の、もとデータ

前回、IGVでPacBioを格好良く見せる方法を紹介しました。
その時例に使ったデータは何かと言うと、Sequelでヒトゲノム NA12878 を10xくらいの深度で読んだものです。
具体的には、ここのデータ

ライブラリの長さは25kb、Blue Pippinを使って15kbにサイズセレクション

• 使用したSMRT Cell 1M の数：10
• トータルラン時間：60時間
• 出力塩基数：32.8 Gb
• リード数：340万本
• リード長のN50 ：11.823 bp

このとき使用した試薬は、旧バージョン、v.1.2 のもの
なので今ならセルあたりの出力はもっと多いはず。

とにかくこれでヒトゲノムの10倍のデータが出た。

このデータをヒトゲノムにマップするのですが、ここで使ったツールは、NGM-LR　＋　PBHoney
PBHoneyは構造変異を検出するツールです。

NGM-LR って何？　という方、これはロングリード用のマッピングツールです。
Next-Gen Mapping tool for Long Read、だったかな？何かそんな名前。
Githubにもあるので、興味のある方はここからどうぞ。

PacBioリードは1本が長いので、例えば 1kb 程度の挿入・欠損をまたいで読むことが可能。
しかし通常のマッパーでは 1kb の変異を考慮してほかの配列を綺麗にゲノムにマップすることができなかった。NGM-LRは、二箇所に分かれてマップするような、ロングリード独特な性質をフルに発揮できるマッピングツール。

BWAとNGM-LRのマッピング結果　Aaronのスライドより

さて、こうしてマッピングした結果は、もちろん参照できます。

先ず、DNAnexusのデモアカウントと作りましょう！
いえいえ、決して私はDNAnexusの手先ではありません。
仕方無いんです。ここにアクセスした方が、データ参照が楽だから。

私はアカウント持っているのですが、ロングインするとこんな感じです。

左下の、"PacBio Sequel Data" というところをクリックします。
これが例の10カバレッジのSequelデータ

"Sequel Data" を開きます

NA12878.reads.ngm.bamというファイルが、マッピングファイルですが、20Gbもあって大きいです。
そこで、indexsession ファイルをダウンロードします。
これは、IGVに取り込むと、DNAnexusのサーバにアクセスしてデータを表示してくれるインデックスファイルです。
IGVはこのように、必ずしもローカルに大きなサイズのマッピングファイルを持っておく必要がありません。

さ、IGV を開きましょう。前回のあれ、ですよ。わからないひとはちょうどこの前の記事をチェック！

ゲノム配列が "Human hg19" であることをチェックして（違ったら Human hg19 を選ぶ）

File > Open Sessions 

今ダウンロードした　indexsession ファイルを選択。

何か聞かれるけどOKをする。

ゲノムのポジションを入力する場所に、試しに、

chrX:116453100-116453795

と入れてGo！

いぇーーーーい !!

ほかにもチラッと見てみたいポジションは、こんなところかな？

chrX:116454160-116454859
chr10:92213800-92216245

InDelの変異箇所は、先の　DNAnexusのアカウントから、

NA12978_Output を選び、ここから ... del.bed  ins.bed ファイルをダウンロードしてきて、取り込むと面白いかも！

いかがだったでしょうか？

IGVでの格好いい見せ方

Integrative Genomics Viewer (IGV) といえば、NGSやっているひとなら一度は聞いたことのある、ゲノムビューワーですね。
フリーで使えてさくっと見るのにはとても便利。
先日、本社とのウェブミーティングで、このIGVを使って、PacBioのデータを格好良く見せる方法を教わりました。
皆さんも知っていると便利なツールだったのでシェアします！

まず、普通にIGV（ここではv.2.3.88）を使ってPacBioのアライメントファイルを開くとこう見えます。

ぜんぜん綺麗じゃない！　
SubreadをアラインしているのでInDelエラーが目立ってるんです。

これをある方法を使うと、このように見ることができます。

とてもすっきりしてますねえ。

これ、Development Snapshot版を使っています。

先ず、バイナリを落としてきます（上図カーソルの場所から）。
Zipを解凍して、中に作られる、igv.batファイルをダブルクリックして起動させます（テスト環境はWindows 7&10）

最初は、これまでIGVで表示していた、デフォルトで、PacBioアライメントが表示されるかもしれません。ここで、View > Preferences　を選択

Alignmentsタブを開き、とりあえず以下のように数字を入れてみて下さい。

特に、

• Label indels > 1：　2塩基以上のindelに表示が出る
• Hide indels < 20：　19塩基以下のindelは表示されない

としてチェックを入れて、OKします。

ちなみに、普通のIGV(v.2.3.88)でも、これに似た画面はありますが、上記オプションはありませんね。

これはv.2.3.88のIGV

さてこれで、先ほどのような綺麗な表示ができるようになりました！

こんな感じに。

しかし例えば、下図のような場所があったとします。

大きなDeletionがある領域です。

そのほかにも、小さなInDelが点在するようですね。

では、Allele Frequencyが小さいものはエラーの可能性が高いので、除くことにします。

例えば、35％以上のアレルだけを残したい、そんな場合は、再び View > Preferences > Alignments から、Coverage allele-fraction threshold: を、0.35としてみる

そうするとこんなにすっきりします！

大きなDeletionがあるところの上流側に、ハプロタイプ（ここでは黄緑のA、リファレンスはT）があるのに気がつきます。

この黄緑色Aを右クリックして、Group alignments by > base at [塩基の場所]　を選択

するとこうなります

大きなDeletionと上流SNV（T -> A）がリンクしているのがわかると思います。

IGVはあくまでビューワーですので、これ自体が解析するソフトではありませんが、データを見せる方法としてはとても使えるツールだと思います。

尚、この機能は、PacBioのプログラマー、Aaronさんによって作られました。

さて、次回は、この例に使用したヒトゲノムデータについてお話ししたいと思います。

多分誰でも使えるデモデータ、のはず。。

10X Genomics のデータだけに頼ってはいけない理由

 10X Genomics （10XG）という会社をご存知ですか？
Synthetic Long Read、またの名をLinked Read と呼ばれるデータを出力し、Scaffoldを作ります。
Synthetic Long Readというのは訳すと、合成ロングリード？擬似的ロングリード？でしょうか。
でも、これは本当の意味でロングリードじゃない！　という意見があったのか、最近は Linked Read、リンクしたリード？、という呼び方がされるようです。

http://www.10xgenomics.com/products/より

バンクーバーのASHGでも、ランチョンやその他のワークショップでも存在感を出していた10XGですが、この会社の対象マーケットは、大きく、1)デノボアセンブリ、2)シングルセル発現解析、の2つがメインだと言えます。
シングルセル発現解析のほうは既存の技術に比べて何十倍ものハイスループット解析ができる、といった優位性があると思います。
ですが、もうひとつのデノボアセンブリは、どうかなあ・・・
10XGに頼って解析を進めていると、大きな情報を見落としていることに気づかないことがあるでしょう。今日はそういう話。

ここから先は、PacBioのJasonさん（Falconなどの開発者）が解析してくれた結果。
彼はバイアスをかけない解析をしますので、私は、とても客観的なデータだと思います。


まず、NA19240サンプルを、PacBioでアセンブルした結果と、10XG v.1.3 でアセンブルした結果があります。
Contiguityを比較するのは大人気ない、というかそもそもContigとScaffoldを比較するのはあまり意味が無い。
そこで、アセンブルした後の結果の、構造変異を検出した場所（Segmental Duplication以外の場所で）を見てみます。

• 10xGアセンブリ、PacBioアセンブリ共に、数千～1万強の、挿入、欠損部位を検出した
• 10xGアセンブリの方が多くの欠損変異を検出し、5,573個は、PacBioアセンブリでは検出していなかった
• 反対に、2,693個の欠損変異は、PacBioアセンブリでは検出され、10xGアセンブリでは検出されなかった

【欠損：10XGで検出されたがPacBioで検出されなかった5,573のうちの１例】

これは、10XGのアセンブリでは176bpの欠損を示していた箇所。
リファレンスには確かに176bpあり、PacBioのアセンブリでも、アセンブリ前の補正後リード（p-reads）でも、確かに176bpがあった例。
この箇所は、GCの連続配列で、10xGが使用しているショートリードでは読めない。
なので間違ってコールされてしまった例。

ほかランダムに抽出した９例

最初の１例を合わせた10例のうち、8例は擬陽性だった。

【欠損：PacBioアセンブリでは検出されたが、10xGでは検出されなかった2,693例のうちの1例】

この配列も、GCが多い場所で、10xGでは読めていない。しかし、この場所には、リファレンスには無い58塩基の欠損領域が、PacBioアセンブリから検出できている。
しかもそれは、ヘテロ欠損。
研究者なら、このようなヘテロな場所の変異に、より興味が湧くのでは深いのでは？

ほかにランダムに抽出した9例

先の1例を除き、全てで10xGのアセンブリでは、読めていない箇所に、実際は欠損変異が存在した。

ということで結論、
読めない配列が原因でコールされた変異は擬陽性の可能性が高い
読めない場所から変異をコールするのは不可能
10xGのデータだけに頼って解析していては、重要な変異を取りこぼすことになる

今回は欠損変異だけに注目しました。（スライドはJasonさん作）

どんな技術もそうですが、それだけに頼って解析していては、見落とすものは必ずあります。　なのでバリデーションは大事。
特に、比較的新しい技術に関しては、一見素晴らしいように見えても、欠点もちゃんと意識して使わないと、レビューワーから突っ込みを受けることになりますね。

10xGユーザの方は、上記のような擬陽性の可能性もあることを考えて、一度、PacBioで読んでバリデートしましょう！！

アメリカ人類遺伝学会 ASHG PacBioネタ

ちょっと報告が遅れましたが、ASHG＠バンクーバーのPacBioネタです。

初日のマクロジェン社のランチョンでは、10Kアジア人ゲノムプロジェクトの話がありました。当然PacBioのにとも何人か参加していましたよ。もちろん私も。
韓国人ゲノムAK1は、最近論文になりましたよね。
これからはアジア人を10000人、読みまくる！という威勢のいい話でした。
もちろん全部をPacBioで読むわけではないでしょう。

二日目はPacBioのランチョンセミナー
これはもう、YouTubeで紹介されていますので、興味のあるひとはここを是非チェック！

「ヒトゲノム＋構造解析」は、これからのPacBioの大きなテーマです。
今までは、ゲノムアセンブリに注力してきたPacBioですが、アセンブリの分野ではもうグローバルスタンダード、になりました。
これ以上長く読める、技術的にも安定したシークエンサーは、今のところPacBio以外存在しない！　というのは過言ですが、客観的に言ってもPacBioのSMRTテクノロジーは素晴らしい技術です。
巷ではナノポア技術が競合とみなされているようですけれど、まだまだ実績が違う。
PacBioで読んだ結果は、ちゃんとこれまで1700本の論文になりましたから。
PacBio は、サイエンスが議論できるレベルのシークエンサーです！


さて、学会期間中に新商品やサービスをアナウンスするのは企業によくあるパターンですが、PacBioもやってくれました。

その名も、新試薬、バージョン 1.2.1 for Sequel
て、地味だ。地味すぎる・・・

いまのSequel 試薬が1.2で、それの新しい版が1.2.1
やっぱり地味だ。

でもこの試薬のすごいところは、Sequelの最初の必要DNA量は変わらないけれど、そこからランできるSMRT Cellの数が増えたこと。
ゲノムDNA20ugくらいからスタートして、今までSequel は数セルしかランできなかったけれど、10セル程度かそれ以上ランできるように、ライブラリ作製効率が上がった（注：この辺の数字は大きくサンプル依存）。
具体的には、プレート上にロードするDNA濃度は、RSIIと同じ、6.5～40 pMになった。
プライマーも新しくなった。
その他、試薬ではないけれどロボット系の動きも少し変わった。

ベータテストの結果、20kb以上の長いライブラリ（w/ サイズセレクション）を読んで、出力は1セルあたり5～7Gb
平均リード長は10kb、11kbくらい（平均ではなくてN50 ならもっと長いですけどね）
バクテリア、植物、ヒト細胞から同じような出力結果が出ているようで。

さらに来年春のアップデートではもっと大きく改善する予定。
と、その時に名前は 1.2.2 とかではなくて、もっとインパクトがあるものにして欲しい！