RNA-Seqのカウントって....どうする？

“ There are known knowns; there are things we know we know.
We also know there are known unknowns; that is to say we know there are some things we do not know.
But there are also unknown unknowns -- the ones we don't know we don't know."

これは2002年 2月12日、記者から、イラクにあるとされる大量破壊兵器がテロリストの手に渡ったという確かな証拠が無いことを指摘されて、これに答えたドナルド・ラムズフェルド国防長官（当時）の言葉です。

当時、あまりに意味難解な返答に、たくさんのコメディアンがジョークのネタにしていたことを覚えています。

既知なこと、知っていること、わかっていること。　これらが存在することを知っていること、わかっていること。

わからないことがあるということを知っていること。　あるいは知らないこと。

・・・哲学的ですね。　

これとRNA-Seqとが何の関係があるかというと、ノーマライズの方法について、知っていることと知っていると仮定していること、知らないこと、をちゃんと知っておく必要があるからです。

タグの数から発現量を算出する時の、RPKMという方法を以前書きました。

RPKMでは、タグカウントを、サンプル内の全カウントで割って、さらに転写産物の長さで割って、定数を乗じて、ノーマライズするのでした。
くわしくはWikiなどを参照。

ここで、
タグ（リード）数はKnownですか？
サンプル内総タグ数はKnownですか？
転写産物の長さはKnownですか？

Isoformの存在、未知Exonの存在などを加味すると、転写産物の長さは、厳密にはUnknownですね。　しかしこれをKnownと仮定して、計算します。
タグ数も、この仮定の転写産物の長さにマップされた数としてカウントされるのですから、厳密にはUnknownです。
こうなってくると、サンプル内総タグ数は、転写産物にMappableなタグという意味では、Knownとするべきでしょうか？

とは言っても、転写産物の長さを厳密に認識するのは、ショートリードでは不可能ですから、みなさんRefSeqなどで公開されている既知のmRNAの長さを用いるのです。　
でも頭のどこかに、これは仮定の長さかもしれない、と覚えておくのは意味があるでしょう。

さて、
複数サンプル間で発現量を比較したいとき、遺伝子の長さがとても影響するということがわかっています。
Oshlack et al., Transcript length bias in RNA-seq data confounds systems biology. Biology Direct 2009, 4:14.
 
転写産物が長い遺伝子ほど、異なるサンプル間で発現量の変動が大きく計算される傾向があるのです。
この論文では遺伝子の長さで補正していません。　発現変動の大きさと、遺伝子の長さをそれぞれ軸にプロットすると、発現変動の大きい遺伝子ほどmRNAが長い遺伝子だということがわかったそうです。

aとcがRNA-Seqのデータで、縦軸があるしきい値で発現変動ありとされた転写産物の割合（％）、横軸がその転写産物の長さ（bp）です。　bは同じくマイクロアレイのデータです。

RNA-Seqの「発現変動ありとする」やりかたに問題があるのは、長い遺伝子ほど変動ありとされている遺伝子の割合が高いことから、明らかです。

そこで、RPKMのように転写産物の長さで補正することが必要になるのですが、これまた、RPKMも万能ではありません。

つづく

ExomeとSNVについてのメモ２

Exomeというのは、ゲノム配列からExon、つまり転写領域のみに注目して、その部分の配列のみをシーケンスする方法です。
ゲノム全体を読むよりも、はるかに効率的に、遺伝子領域を高いカバレージで読むことができます。

ゲノムから興味のある配列のみを選択的にセレクトする技術を、ターゲットキャプチャーと呼びますが、いくつかのメーカーからキットが販売されています。
イルミナ社のTrueSeq、ライフテック社のTargetSeq、ロッシュ社のSeqCap、アジレント社のSureSelectがあります。
ここまでは前にも書きました。

では、初めてのひとは、どのキットを使ったら良いのでしょう？
選択にあたっての基準はあるのでしょうか？

当然ながら、アジレント社を除く3社は、自分たちのシーケンサーにもっとも適したプロトコールを用意してキットを設計しているはずです。
しかし、メーカーに問い合わせてみると、「他社のシーケンサー用にも使えますよ」、という答えが返ってきます。　（「でも保証はしませんが」、と念を押されるかもしれませんけど）
考えてみればキャプチャーは、シーケンサーにかけるずっと前の作業なので、キャプチャーに使用した試薬等がきちんとWash outされれば、ライブラリー作成などに問題はおこさないはずです。

問題は、どれだけキャプチャーの効率と、キャプチャーに用いるプローブの設計、何塩基の間隔でプローブが設計されているか、アレイと使うか液体中のビーズを使うか、などによる違いが大きいでしょう。

これらを比較した論文があります。　

Performance comparison of exome DNA sequencing technologies.
Nat Biotechnol. 2011 Sep 25;29(10):908-14.　PMID：21947028

時間が無いのでサマリーだけ読みたい方はこちら
ヒトのExonキャプチャーキットの比較です。
これによると、プローブの設計にだいぶ特徴があるみたいです。

プローブ長もさることながら、設計場所の「くせ」もあるようです。

Nimblegenはもっとも多くのExon領域を、micro RNA を含め、たくさんカバーしているそうです。　
AgilentはExonバリアントの検出に強いそうです。　
Illuminaは、UTRの配列もキャプチャーできるそうです。　
また、Illuminaはリードデータをもっとも多く必要とするそうです。　最新のHiSeq2000やSOLiD5500ならリードは膨大に出力されるので問題なさそうですが、バーコードを付けて多くのサンプルを読むときは注意が必要ですね。

あくまでヒトExonキャプチャーの比較ですが、バリアントもキャプチャーしたいならAgilent、UTRに注目したいならIllumina、micro RNA も含めてもっとも高いカバレージでキャプチャーしたいならNimblegen、と言えるでしょうか。　かなり大雑把ですが。

私はキャプチャーキットに関しては詳しく知りませんでした。　違いは実験の差くらいに考えていました。
ですが、そのような私にもわかるように、キットの特徴を第三者の目で比較してくれる、こういう論文はうれしいです。

P.S.　SNP検出に関するキャプチャーキットの差はあまり無いそうですが、検出されるSNPの差はマッピングやSNP検出ツールのパラメータによるものが大きいと思います。

Exome と SNVについてのメモ 1

最近ブログを更新していませんでした。
先月末から何ていうのでしょう、一種のスランプみたいな、そんな気分でしたので。
もう脱しつつあるので、大丈夫です。

昔からGWASなどでSNPを解析しているひとには常識かもしれませんが、NGSが出た最近やっとSNP（SNV）を始めてみようかな、というひとには意外と「目からうろこ」なことがあります。

私もSNPを昔からやっていたわけではないので、たまに論文を読んでいて気づかされることがあります。
で、そういうことを他の研究者と話すと、意外とその人も知らなかったり。

具体的には、NGSを使ってWhole Exomeをやっている解析で、SNVを見つける、そのときのワークフローです。
当たり前のように行っていたフィルタリングのプロセスで、「なぜそのフィルタリングを行うのか」「基本となる考え方・仮定は何か」という基本的なところがおろそかになっていたと。
そこについてメモ書きのようですが、まとめました。

Exome解析の最終目的は、おおよそ、「疾患原因となるSNVの候補を絞り込む」ことになると思います。
そこで、

1. ヒトならヒトのExon領域をカバーしたターゲットシーケンスを行い、その領域で十分な厚みを持ってマップされた場所から、SNVを見つけます。　（例：SureSelectなどでターゲットキャプチャーしたあとIlluminaマシンで大量に読み、BWAでヒトゲノムにマッピング、キャプチャー領域にてSamtoolsでSNVを検出）
2. 検出されたSNVの場所とリストをもとに、dbSNPなどのデータベースに無いものを抽出
3. さらにその中からNon-synonymousのSNV（アミノ酸置換を伴うSNV）のみを抽出
4. そして残ったnsSNVのうち、患者複数のサンプルで共通するものを選び出す

ここから先はその疾患原因SNV候補のバリデーションになるわけですが、上記の１～４の処理はそれぞれ理由、というか前提や仮定があります。
それを理解していないでなんとなくパイプライン的に解析を処理していると、後から他の研究者に突っ込まれてタジタジ・・・となってしまう！

以下はSNPをやってたひとには当たり前すぎることを書いていると思います。　どうぞお許しを。

1. Exon領域をキャプチャーする大前提は、テーマとしている疾患の変異を、遺伝子のコーディング領域から探そうとしているわけです。　つまりコーディングされない場所にいくら疾患原因SNVがあったとしても、それは見ないことにしよう、としています。　
   もちろんコーディング領域は重要です。　しかしそれ以外の転写制御領域のSNVも大変重要ですが、これらは一般的なExomeでは解析対象から外されます。　
2. dbSNPに無い、SNVのみを探すというのは、一見、新規性を探索するのに合理的な方法です。　この疾患原因SNVは非常にレア、珍しいという前提があります。　しかし、データベースに登録されるSNVは、現在爆発的に増え続けていますので、既知であっても自分の研究疾患では新規ということは十分あり得ます。
3. Non-synonymousのSNVに絞るということは、アミノ酸置換、あるいはフレームシフト、ミッセンス等を伴うような変異は特に重要だという前提に立っています。　「この疾患の原因は、アミノ酸変異を伴うSNVだ」　という仮定がそこにあります。　そうでない仮定の場合にNon-synonymous SNVに絞るのは適当でないでしょう。
   また気をつけるべきは、Alternative Exon（スプライスバリアント）の存在です。　これを無視してSNVを見つけているケースが多いのではないでしょうか。　ゲノムのSNVではなく遺伝子コーディング領域中のSNVとなると、スプライスバリアントが異なればそのSNVの有る無しの意味は大きいはずですよね。　（言葉足らず）
4. 複数の患者のサンプルで同じSNVを見つけるというのは、そのフェノタイプに共通するSNVを見つけ出すということです。　SNV（原因）があれば必ずそのフェノタイプ（結果）を引き起こすことをComplete Penetranceと言います。　一方、そのフェノタイプ（結果）には必ずそのSNV（原因）が見られるときをComplete Detectance　と言います。
   つまりこのフィルタリングは、全ての患者に共通するSNVがあるはずだ、という前提に立っているわけです。

まとめると、上記のワークフローで解析するための前提・仮定は、
「全ての患者に共通するSNVがあって、かつそのSNVを持っていればかなりの確率でその疾患にかかり、そのSNVはとても珍しく、そしてタンパク質のアミノ酸コードを変えることで疾患フェノタイプを引き起こす」　
ということになりますでしょうか。

今日は字ばっかりですみませんでした。

PacBio データの特徴 2

PacBioの実験ワークフローで、とてもユニークなのがテンプレートの形です。
ダンベル型のSMRTBellというそうです。

Travers et al., A flexible and efficient template format for circular consensus sequencing and SNP detection (2010). Nucleic Acids Research 38, e159.

DNAをフラグメント化して、断片の両端をエンドリペア後、輪っか配列（プライマーのため）を付けてAのダンベル型テンプレートができます。
左端の、輪っかの部分がシングル鎖になっていますね。ここにプライマーがくっつきます。
ポリメラーゼは、DNAの二重鎖のところを押しのけて、プライマーの先に、相補鎖に正しい塩基を合成していきます。
テンプレートは、ずーっといくと一周しますね。
そうしたらまた、Bのように、今合成した通りにもう一度合成をするのです。
つまり、インサートの長さに依存することなく、何度でもシーケンスが可能です。

これがダンベル型テンプレート。
このままの形がスタンダードです。

インサートの長さがものすごーく長い（数10kbとか）とき、全部読むことはできないので輪っかの近いところだけを読むことになります。　これはMait-Pairのような感じです。
Strobeと呼ぶそうです。

もうひとつはインサートの長さは適度で、テンプレートを何周も読む方法。
何度も何度も同じ配列を読むことになるので、精度は上がります。
Circular Consensusと呼ぶそうです。

ここから先は、スタンダードのテンプレートの話です。
SMRTでは、インサートを全部読むので、Read長＝インサート長、になります。（SMRTについては前回のブログを参照）
しかし、全部読み切れなかったインサートもあるでしょう。
そんな、「途中まで読めたインサート」配列が混じっている状態のReadのことを、Subreadと呼ぶそうです。
Subreadがいくつか集まってRead（=インサート配列）を形成する、ということです。

ですから、SMRTで読んだ解析結果には必ず、
Subread の平均の長さ
Read の平均の長さ、95パーセンタイルの長さ、・・・
など、Subread とReadの両方の情報があるのが普通でしょう。

前回紹介した、Expression Analysis社のWebセミナーで紹介していたのは、SMRTで大腸菌のゲノムを読んだときのデータ解析についてです。
深夜2時過ぎのことでしたので私の記憶が正しく無い箇所もあるかもしれません。

2kbと6kbの2種類のインサートを読んでいました。
SMRTのひとつのセル（フローセルのような単位）で、34～108 MbのMappableなデータが出力されました。
平均1500塩基長、プロトコルを改良すると2700塩基長を一度に読めていました。

特徴のひとつは、大腸菌ゲノムに対するカバレージが一定だったことです。

普通、GC含量のばらつきなどにより、読みにくかったりする箇所がどうしても出てくるものですが、
それが比較的無くカバレージが一定で、
2kbインサートが6.4～17.1カバレージ
6kbインサートが6.7～11.2カバレージをゲノム全体で保持していたそうです。

精度ですが、Subreadのレベルの精度は85%程だったそうです。
これがSMRTの精度だと一般的に言われると、15%は間違い？　結構高いね、と思われそうですが、SubreadではなくRead、あるいは何度も繰り返して読んだ時のConsensusともなると、96.46～99.39%まで向上していたそうです。　

何度も読むというのは、賢い方法です。
ダンベル型テンプレートならでは？

今はまだ、バクテリアサイズのゲノムシーケンスに用いられているPacBioのSMRTです。
ランニングコストがどれくらいなのかはわかりませんが、もっと汎用的に使われる時はもうすぐ来るのでしょうか。
個人的には、機械が大きすぎ！　って思います。

もうちょい、コンパクトにならないのか・・・

PacBio データの特徴 1

PacBioのシーケンサーから出てくるデータはどんなものなんだろう？

その前に、「4分でわかるPacBio-SMRTテクノロジー」をYouTubeで見つけたので、SMRTって何だ？ってひとはまずはここから。

テクノロジーについてもっと詳しく知りたい方は、PacBioのホームページか、この論文がお勧めです。
Schadt et. al., "A window into third-generation sequencing". Hum. Mol. Genet. (2010) 19 (R2): R227-R240
SMRTの特徴は、DNAを1分子単位で、PCRによる増幅無しで、結構長く読めるということです。
Single Molecule Real Time の略です。

ガラススライド上に無数の小さなポケット（ZMWs (zero-mode waveguides) ）があり、その穴の底には、DNAポリメラーゼが固定されています。 この穴は非常に小さく、波長600nmまでの可視光線は通り抜けできない構造になっています。

シーケンスは、プライマーがアニールされたDNA鋳型がポリメラーゼに取り込まれ、反応準備が完了します。
4種類の異なる蛍光の付いた塩基が、ZMWsの外から入り込み、DNA鋳型に取り込まれると蛍光が切り離されてZMWsの穴の底で光ります。その光は底から30ナノメートルまでしか届きません。
スキャナーはその光を検出します。　底にはポリメラーゼが固定されているので、その近くで観測される光は、取り込まれた塩基による光しかない、というわけです。

でもランダムにZMWsの外からも蛍光付きの塩基が入り込むので、偶然ポリメラーゼ近辺に来て、蛍光が検出されてしまうことはないのでしょうか。　
IlluminaやSOLiDのような、蛍光検出一回につき余分な蛍光を洗い流す、Wash-and-Detectはしていません。　
あくまでもReal-Timeにポリメラーゼが自然に塩基がとりこまれるままで読んでいるのです。

トリックは、ポリメラーゼに塩基がとりこまれる速度はミリ秒単位、それ以外のZMWs内の塩基の出入りはマイクロ秒単位、という時間の差にあります。
ミリ秒単位で光らないと、DNA鎖の塩基として認識しないのです。

先日、Expression Analysis（http://www.expressionanalysis.com/）というアメリカの会社がPacBioのデータ解析Webセミナーをやっていたので、深夜2時だったけど参加してみました。
PacBioのデータって、どんなんだろう？　すごく興味があったせいか、3時までずっと起きていました。　次の日会社・・・

さて、SMRTテクノロジーで検出される塩基は、ポリメラーゼがReal Timeに合成する速度と関係が深くなります。
ということは、塩基の検出される時間幅は、一定では無いということです。
Webセミナーの中でもそれを言っていました。

同じAGCCATと読んでいても、ポリメラーゼは機械ではありませんので、1塩基読むのに時間のばらつきが生じます。　これが画像からのシグナル変換、いわゆる一時解析を複雑にしている原因なのです。
とは言っても、PacBioの機械・ソフトは、これを何とか解決しているのでしょう。　でなければ商品化していませんよね。

次は、PacBioのデータを解釈するにあたって絶対必要な、リード、サブリード、ダンベル型のSMRTBellテンプレート、について説明します。

NGSでmicroRNA 2

前回はNGSでmicroRNAを読むときの一般的なワークフローで、マッピングまで書きました。

既知のmiRNA配列にマップしたときは、その発現量がマップされたリードの本数で表現できる、ということは直感的にわかるでしょう。
以前、RNA-Seqの発現量推定のところでも書きましたが、リードの数はサンプルごとのノーマライズ（正規化）が必要です。

（その既知miRNAにマップされたリードの本数）/ （そのサンプルランで取れた全リードのうちmiRNAにマップされたリードの合計数）

出力されたリードの数、特にマッパブル（Mappabl）なリードの総数は、サンプルによって異なります。　ですから、各miRNAにマップされたリードの本数を全体で割ってあげて、異なるサンプル間で比較可能にするのです。
分母を、そのサンプルのリード本数全体　として計算している方もおられますが、私は、マップされたリード本数全体　とした方が精度が上がると考えています。
だってマップされないリードまでノーマライズの分母に含めるのは不自然でしょう。

これだけでは値が非常に小さくなるので、100万倍して、RPM（Read Per Million）という表現にするのが良いかも知れません。

さて、そうして正規化された値は、当然ながらmiRNAの発現量そのものではありません。
あくまで比です。
ということで、この値は、ほかのサンプルの正規化miRNA発現値と比較して、何倍発現量が多い、少ない、というような解析に持っていきます。
マイクロアレイと同じですね。

Nが3以上あれば、正規化後の値を対数変換してから、通常の統計手法（T‐Testなど）で比較すると良いでしょう。

DESeqという発現量比較のツールも使えるかと思います。
いろいろ調べてみると結構miRNA関連論文で、このDESeqが使われています。　
元々は発現変動遺伝子を解析するツールです。　
遺伝子は全部で数万あり、その中でほとんどは変動していない、という仮定の下、変動有意遺伝子を見つけるわけです。
microRNAは全部でもたかだか数百、これで遺伝子発現と同じモデル（負の二項分布）式を使っても良いのか、という疑問は残りますが、ここではそれを置いておくとして。

Nが1しか無い場合、割り算して比較するしかありませんね。

有意な発現変動を示したmiRNAが見つかったら、そのターゲット遺伝子を探してみましょうか。
データベースとしてはまず、miRBase (http://www.mirbase.org)がいいでしょう。

ここからはTargetScan や Pictar　といった、ターゲット予測プログラムのサイトへもリンクできます。
ターゲット遺伝子を見つけたら、その遺伝子のGO（Gene Ontology）を検索して、どんな機能の遺伝子がターゲットになっていたのか、を調べることも良く行われます。

さて、これら以外に、NGSならではの解析とは何でしょう？

そう、既知ではなく、未知のmiRNAを見つけることです。

それには先ず、マッピングのデータから、既知miRNA以外の場所にマップしたリードが作っている、ある程度の長さのクラスターを見つけてこなければいけません。
ある程度の長さって、どれくらいだ？
クラスター配列は、2次構造を取り得る、いわゆるpalindromicの配列でなくてはなりません。
1次元で言うと、（mature miRNA配列）＋（ループ構造）＋（mature* miRNA配列）のパターンである必要があります。
そのような配列を取ってきたら、エネルギーが最小になるようなときにちゃんとヘアピンのようにフォールドするかどうか、を確かめます。
さらに、Droshaなどのタンパク質に切断されるサイトがあるか、調べます。

これを全部プログラムを作って行うのは大変ですが、幸い、アカデミックの方ならフリーで使えるツールがあります。　miRDeep というものです。

このドラえもんの手、みたいなのがmature miRNA-Loop-star の2次構造で、薄い線がマップされているリードです。　
matureの方により多くのリードがマップされているのは、より本物らしいmicroRNAだったことを裏付けます。

残念ながら、企業はフリーじゃないんですよね。　いくらか払うのかな。

2次構造だけの予測なら、できるソフトもありますが、miRDeep ほど特化しているものは少ないでしょう。　チャンスがあれば使ってみたいです。


microRNAをNGSで解析するというテーマで、よくまとめられているレビューがあるので紹介します。
Motameny et al. (2010) Next Generation Sequencing of miRNAs - Strategies, Resources and Methods. Genes. 1, 70.

こちらはmiRanalyzerというツールのワークフローです。フローは参考にしています。
Stark et al. (2010) Characterization of the Melanoma miRNAome by Deep Sequencing. PLoS One. 5, e9685.

発現量比較から新規miRNA予測、ターゲットのGO解析まで行った例です。　こちらも良い。
Dhahbi et al. (2011) Deep Sequencing Reveals Novel MicroRNAs and Regulation of MicroRNA Expression during Cell Senescence. PLoSOne. 6, e20509.

NGSでmicroRNA解析 1

microRNAというのは、タンパク質にはならない、いわゆる non-coding RNAの仲間です。　
それ自身で標的のmRNAの3’-UTRなどに結合し、標的mRNAを分解します。
microRNAは発現制御機能を持った、RNA分子なのです。

と、ごくごく簡潔にまとめましたが、microRNAについては、多くのレビューや教科書に記載がありますので、詳しくは書かないとして。　Wikipediaでもかなり情報が得られますよ。

さて、microRNAの解析は、長らく、TaqManアレイやその他のマイクロアレイ、PCRなどが主流でした。
近年、NGSでもこれを解析している例が見られます。
実験プロトコールとして一般的なのは、

1. トータルRNAを抽出して
2. ゲルに流した後18－30塩基のあたりを切り出し
3. Small RNAを精製して
4. ライブラリーキットにある通りアダプターを付け
5. cDNA合成し
6. ライブラリーをシーケンサーで短めに読む

というものです。

短めに読む、というのはフラグメントの長さが　18-30塩基を予想しているからで、35、6塩基も読めば十分全体をカバーするからです。
ゲノムや転写産物を読む場合は、できるだけ長く読めたほうが良いですが、microRNAは短くてもOK、むしろちゃんと正確に読めることが大切です。

さて、シーケンサーが無事データを出力しました。
ここからデータ解析です。

microRNAの解析では、リードのアダプタートリミングは大切です。　（microRNAに限ったことではありませんが）
読んでいるフラグメントの長さが短い分、5’側のアダプターから読んだ時に、フラグメントを読み切り、3’側のアダプターまで読んでしまうことがあるからです。
そのため、3’側のアダプターをトリミング（除去）することが必須です。

次に、トリミングした後のリードが、短くなりすぎると、後でマッピングする際にNon Specificにマップされる恐れがあるので、これも除きます。　15‐18塩基以下の短いリードは除去するといいかもしれません。

リードがきれいになったところで、いよいよマッピングです。　でも何に対して？

一例を示します。　今度参考文献も示しますね。

1. 先ずゲノムにマッピングし、マップされなかった配列はごみとして除去する
2. ゲノムに当たったリードを回収、miRBase の配列（precursor, mature, mature*）にマッピング
3. miRBaseの配列にぴったり当たったリードは、既知のmicroRNAとして保存
4. miRBaseの配列に当たらなかったリードは、次に、microRNA以外の、既知のnon-coding RNA （piwiRNA, snRNA, snoRNAなどの） 配列に対してマッピング
5. non-coding RNAにも当たらなかったリードは、念のため、RefSeqのmRNA配列に対してマッピング
6. それでも当たらなかったリードは、新規microRNA「候補」として保存

これにより、①既知のmicroRNA、②それ以外のnon-coding RNA、③新規microRNA 候補、の3種類のデータセットができる。

既知のmicroRNA、non-coding RNA については、発現量の解析をすることになるでしょう。
新規候補の解析は、少しやっかいだが、アルゴリズムが無くはないのです。

つづく