NGSのdata解析の全くの初心者です。 一本鎖DNA由来のライブラリーをpair-endで読んだ後のBowtieによるMapping時のoptionについて質問です。 MeDIP-seqの要領で、Illumina Truseqのライブラリーから一本鎖を選択的に精製したライブラリーをPCRで増幅後pair-endのsequenceingを行いました。BowtieでMappingする際、精製したライブラリーのstrandの情報を維持した状態でmappingの結果を得るためには、一つ目(-1)を"Read1"、二つ目 (-2)を ”Read2のreverse complement”として、pair-endのoptionを"--ff"とすれば、Read1とRead2が同じstrandにmapされた結果が得られるのでしょうか? Galaxyでは、--ffはSOLID用のoptionとして扱われているような気がしますが、IlluminaのR1とR2rev-compを用いても結果は問題ないのでしょうか? ご存知の方で、アドバイス頂けましたら幸甚です。 よろしくお願い致します。 質問日 Apr 14 '12 at 15:22 touran |
適当な解答の罪滅ぼしに、自分でも同じことしてみました。結果、ffは機能しているように見えます。 何か、別の原因があるのではないでしょうか? 以下の処理で実行しているのは、この場合はRNA-Seqの公開ペアデータを使用していますが、 本質的にはtouran様のデータと変わらないように思います。 やったことは、以下の3データセットそれぞれに対し
それぞれ、fr/rf/ffの各オプションを試しています。
上記の結果、まともにヒットしたペアリードのある組み合わせは以下の2つのみです。
この内ff.12revと書いてある条件は、リード1はそのままに、リード2は逆鎖をとったデータに対して"ff"オプションを使用した際の 結果であるはずなので、touranさんの当初の意図が達成されていることになるかと思います。 以下、実際ヒットしたものを抜き出していますが、ffのペアは++で当っていることが分かるかと思います。
可能性としては、ヒットの確認使っているビューワ等が、ffリードのペアに対応していないとか。 回答日 Apr 21 '12 at 02:16 nob_fj ♦ nob_fjさま 回答ありがとうございます。 --ffのpair-end mapping自体うまくいくようになりました。公開されているdataを使ってやってみる限りではうまくいきましたが、自分のdataではうまくmappingされていなかったのでdataとbowtieとの相性を調べていました。お恥ずかしい限りですが、自分のdataを使う限りbowtieでpair-endのmapping自体がうまくいっていないように見えるdataセットがあることにようやく気づきまして、原因がpair-endのライブラリーのインサートの長さであることがわかりました。IPしたライブラリーには長鎖の断片が優先的に濃縮されてくる性質があることを考慮してませんでした。以下のサイトでも議論されていますが、Bowtieのデフォルトのインサートサイズが250b?と短めなので、長めのライブラリーにはoptionが必要なのですね。 http://www.biostars.org/post/show/9090/bowtie-pair-end-broken/ -X 800などのインサートのサイズを大きめに設定してやるとR1&R2の--frでもR1&R2rvcompの--ffでもmappingがうまくいくようになりました。ご迷惑をお掛けいたしましたが、いろいろ親身になって回答頂き大変助かりました。ありがとうございました。
(Apr 23 '12 at 13:30)
touran
こちらこそ、ChIPのノウハウは乏しいため、貴重な情報提供ありがとうございました。また、お役に立てて何よりです。
(Apr 23 '12 at 15:12)
nob_fj ♦
一応誤解のないように。IP(免疫沈降)といってもChIPではなくMeDIPやBrdU-IPのように抗体の抗原がDNA自体になる場合は、長いDNAの方が抗原となるものが多く含まれる率が高くなるため、飽和していない条件では長鎖の断片の方が優先的にIPされてくる性質があります。
(Apr 24 '12 at 13:23)
touran
誤解しておりました。DNA自体を抗原にした場合のIP独自の問題だったのですね。大変参考になります。
(Apr 24 '12 at 16:41)
nob_fj ♦
|
直接の解答ではないですが、ご参考までにですが、 入力ファイルの先頭10000エントリーとかを取り出して、小さいサブデータセットを作り、 考えうる選択肢でマップしてみて、どれがうまくいくか試すのは、パラメータの検討や、確認でよくやります。 たまに、マニュアルの記載の方が間違っている場合もあるらしいですし。コマンドでやるなら以下のようなかんじでしょうか。 比較はSAMからBAMにして、IGVやら、samtoolsのtviewでやるなり、いろいろ考えられます。 Galaxyで同じことやるやり方は良く知りません。
回答日 Apr 16 '12 at 14:22 nob_fj ♦ nob_fjさま アドバイスありがとうございました。小さいサブデータで試してみましたが、R_1とR_2のreverse complementの組み合わせの-ffはうまくいかないようでした(ほとんどがマップされない)。R_1とR_2のreverse complementを別々にsingle readでmapしたものはきちんとマップされますし、サブデータの中ではR_2とR2_revはきちんと同じ場所の相補鎖にほとんどすべてがmapされてますので、とりあえずread countを稼ぐためにsingle readとしてwetの実験自体がうまくいったかを検証しています。mappingのoptionの問題は時間をかけて解決したいと思います。貴重なアドバイスありがとうございました。
(Apr 20 '12 at 21:55)
touran
こちらこそ、貴重な情報報告いただいて、大変参考になりました。上記のスクリプトを貼りつけた後に、よくよく考えると、40000では、少ないこともあるかもと思い、若干心配しておりましたが、わずかながらでもお役に立てて何よりです。
(Apr 20 '12 at 22:02)
nob_fj ♦
|
ちなみに、今更、質問内容を何となく理解しましたが、
という理解であっておりますでしょうか。 試していないので、確証はありませんが、リード2だけ予め配列のcomplementを取った上で"-ff"オプションでマップすれば、上記の状況と異なり
となるのでないでしょうか。 逆鎖の取得はfastx tool kitなどでもできますし、
Illuminaのfastqであれば、以下のようなperlのワンライナーでも(多分問題なく)変換できます。
いちおう、上記の結果は私が試した限りでは、結果は一致しました。
上記のいずれか、または他の適当な手段でリード2のファイルのみ相補鎖に変換した上で"-ff"でマップとかならうまくいくかもしれません。
※余談ですが、fastx tool kitを使う場合は、最近の配列は-Q 33が必要だと思いますが、そのオプションをつけないとエラーが出るときだけつけて下さい。tool kitのバージョンや、配列のバージョンによって、いらない場合もあるかもしれません。 回答日 Apr 20 '12 at 22:51 nob_fj ♦ nob_fjさま 回答ありがとうございます。これまで質問の意図がうまく伝わっていなかったようですみません。”・・・となるのでないでしょうか。”とおっしゃるとおりの質問でした。Fw-Rvの組み合わせのreadなら-fr、Rvのreadをrv-compにしてやってFwに変えてやりFw-Fw'の組み合わせにすれば-ffでmappingされると思っていたのですが、なぜかうまくいっていません。 当方toolのインストールなど不慣れですので、FASTXなどはGalaxyなどを使ってR_2のreverse complementをおっしゃるとおりつくり、fwのR_1とfwの向きになおしたR_2(R_2rv-comp)を使って-ffのmappingをやってみましたがほとんどmappingされませんでした(GalaxyのBowtieでも私のMacのBowtieでも)。なにか途中でミスをしてるのか、他のオプションの設定があるのか、何かがまちがっているのか理解不能な状況に陥っています。 P.S. -Q 33とは何のことでしょうか?(初心者質問ですみません)
(Apr 20 '12 at 23:31)
touran
理解しました。理解できておらず、こちらこそお手数おかけしました。何某かのBowtieの仕様なのでしょうね。気になるようなら開発者に連絡してみるのも一つの手かもしれませんね。 -Q 33は、fastx toolkitの使用時のオプションで、galaxy経由で実行する場合は、関係ない気がします。上のスクリプトの隠れた部分に書いてあります。質問自体には直接関係ないので、忘れて下さい。
(Apr 20 '12 at 23:41)
nob_fj ♦
|