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 については、発現量の解析をすることになるでしょう。
新規候補の解析は、少しやっかいだが、アルゴリズムが無くはないのです。

つづく