ログイン 概要 よくある質問

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's gravatar image

touran
1524


適当な解答の罪滅ぼしに、自分でも同じことしてみました。結果、ffは機能しているように見えます。 何か、別の原因があるのではないでしょうか?

以下の処理で実行しているのは、この場合はRNA-Seqの公開ペアデータを使用していますが、 本質的にはtouran様のデータと変わらないように思います。

やったことは、以下の3データセットそれぞれに対し

  1. "リード1そのまま"-"リード2そのまま"のペア
  2. "リード1そのまま"-"リード2の逆鎖"のペア
  3. "リード1逆鎖"-"リード2そのまま"のペア

それぞれ、fr/rf/ffの各オプションを試しています。

# 頭から10,000エントリー抜き出す。
head -n 40000 Data/SRR315301_1.fastq > Dataext/SRR315301_1.ext.fastq
head -n 40000 Data/SRR315301_2.fastq > Dataext/SRR315301_2.ext.fastq

# リード1とリード2の逆鎖をそれぞれ変換して別のファイルにします。
/toolDir/fastx_toolkit-0.0.13/bin/fastx_reverse_complement -i Dataext/SRR315301_2.ext.fastq -o Dataext/SRR315301_2.ext.revComp.fastq -Q 33
/toolDir/fastx_toolkit-0.0.13/bin/fastx_reverse_complement -i Dataext/SRR315301_1.ext.fastq -o Dataext/SRR315301_1.ext.revComp.fastq -Q 33

# ffの場合
bowtie --ff /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.fastq -2 Dataext/SRR315301_2.ext.fastq > ff.12.bowtie
bowtie --ff /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.fastq -2 Dataext/SRR315301_2.ext.revComp.fastq > ff.12rev.bowtie
bowtie --ff /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.revComp.fastq -2 Dataext/SRR315301_2.ext.fastq > ff.1rev2.bowtie

# frの場合
bowtie --fr /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.fastq -2 Dataext/SRR315301_2.ext.fastq > fr.12.bowtie
bowtie --fr /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.fastq -2 Dataext/SRR315301_2.ext.revComp.fastq > fr.12rev.bowtie
bowtie --fr /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.revComp.fastq -2 Dataext/SRR315301_2.ext.fastq > fr.1rev2.bowtie

# rfの場合
bowtie --rf /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.fastq -2 Dataext/SRR315301_2.ext.fastq > rf.12.bowtie
bowtie --rf /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.fastq -2 Dataext/SRR315301_2.ext.revComp.fastq > rf.12rev.bowtie
bowtie --rf /dataDir/hg19multi/bowtie-0.12.7/hg19.fa -1 Dataext/SRR315301_1.ext.revComp.fastq -2 Dataext/SRR315301_2.ext.fastq > rf.1rev2.bowtie

#結果のサマリ
$ wc -l *.bowtie
     6 ff.12.bowtie
  1494 ff.12rev.bowtie
    14 ff.1rev2.bowtie
  1584 fr.12.bowtie
     2 fr.12rev.bowtie
     0 fr.1rev2.bowtie
    24 rf.12.bowtie
     2 rf.12rev.bowtie
     4 rf.1rev2.bowtie

上記の結果、まともにヒットしたペアリードのある組み合わせは以下の2つのみです。

  1. ff.12rev.bowtie
  2. fr.12.bowtie

この内ff.12revと書いてある条件は、リード1はそのままに、リード2は逆鎖をとったデータに対して"ff"オプションを使用した際の 結果であるはずなので、touranさんの当初の意図が達成されていることになるかと思います。 以下、実際ヒットしたものを抜き出していますが、ffのペアは++で当っていることが分かるかと思います。

$ head -n 6 ff.12rev.bowtie
SRR315301.653 JOHNLENNON_0014:1:1:4723:1141 length=76/1 +       chr13   39587153        GGCCGTTGAGCTAGGATGAGGGAGAGATGTGGAGGGGGCAGCGGTAGATGTGGTGGCGATTAATGACCAGGGGGGG      GGDGGDDGGGDDGGG>AGCEECAC@###################################################      0       0:T>G,35:T>G,50:A>T,57:T>G,68:C>A,69:T>G,73:A>G,75:A>G
SRR315301.653 JOHNLENNON_0014:1:1:4723:1141 length=76/2 +       chr13   39587320        TGGGAGTNGAGGAACCATGTGGTGGAAGGNNAATAGAGGAAAGAGGATTTGAACCATTTAAAGGAGTACTCAAGCT      >??A@?=#A?AIEEDFHHIIH?==?=59@##GGGGGIIDHIGGGBEDGEGGEEGGGI=GIFIDGIIEGD>DIGGGI      0       7:G>N,29:G>N,30:A>N
SRR315301.803 JOHNLENNON_0014:1:1:1949:1141 length=76/1 +       chr5    75537629        GTTTTCTCCACATCAATCCAGCTTTGGAGACATTCTATTAGTGACATATGCCCCTTCCCCCAAAAACAACAATGAA      GGBGGHHHB@GGGEGHH<HBDAGGDEC>>A?=A??GDEBCC<C>@BEBCED:DCDEEEDAB@1:@;9325249823      0       41:C>T
SRR315301.803 JOHNLENNON_0014:1:1:1949:1141 length=76/2 +       chr5    75537797        TTCAAACNGGACAGTCACCAGAAGTACACNNTTATCAAAAATGCACACACTTCACTTGGCATCTCCAGCACCTTCA      ###########B:>+B?CC??78353<79##?:B?B?AAEA<;?;B9@.A5AF;F@CCCFC=B;@BEFE<FIGIII      0       7:T>N,9:T>G,29:A>N,30:G>N
SRR315301.817 JOHNLENNON_0014:1:1:4676:1152 length=76/1 +       chr6    7287801 TGCGAGATCCATGCTTTATGAACACATTATTCTCACACATCCGGGGATAGGTTTGATACACAAGATCCTTAAAGTC      HHHHBHHGHHHHHHHHHHHHHHH<HHHGHHHHHHHHHHHHHHGHHHHGHF@EHHHHHFGHHGDGDBDB<B===A=A      0       0:A>T
SRR315301.817 JOHNLENNON_0014:1:1:4676:1152 length=76/2 +       chr6    7287964 AAACTTCATCCAAATATTTACATTTTAAAAAGACATTTAAAAAATAGGCTACTTATTAAGTGTGCTTGAACAAACC      EBBDB>IIDGIHHIFFEEFDIIEECIIHIHIIGIIHIIIIFIIIFEIHIIIIIEGGGHIIHIIIIHIIIIIIIIII      0
$ head -n 6 fr.12.bowtie
SRR315301.653 JOHNLENNON_0014:1:1:4723:1141 length=76/1 +       chr13   39587153        GGCCGTTGAGCTAGGATGAGGGAGAGATGTGGAGGGGGCAGCGGTAGATGTGGTGGCGATTAATGACCAGGGGGGG      GGDGGDDGGGDDGGG>AGCEECAC@###################################################      0       0:T>G,35:T>G,50:A>T,57:T>G,68:C>A,69:T>G,73:A>G,75:A>G
SRR315301.653 JOHNLENNON_0014:1:1:4723:1141 length=76/2 -       chr13   39587320        TGGGAGTNGAGGAACCATGTGGTGGAAGGNNAATAGAGGAAAGAGGATTTGAACCATTTAAAGGAGTACTCAAGCT      >??A@?=#A?AIEEDFHHIIH?==?=59@##GGGGGIIDHIGGGBEDGEGGEEGGGI=GIFIDGIIEGD>DIGGGI      0       45:A>N,46:G>N,68:G>N
SRR315301.719 JOHNLENNON_0014:1:1:11792:1141 length=76/2        +       chr14   32625445  CACAGGTAAATTTCACTAGGTTATATTTTGTGTAGTAAAGAAAAANNNTATTTGGTCAATGTTATCNNNNTTCATA    IDFGIIF>IGGGGGGEDDGDDGGGGGGGDGDGEFGEGFGGA>>==###:68888;BIEEIIDFEE###########      0       45:A>N,46:A>N,47:T>N,66:T>N,67:T>N,68:A>N,69:A>N
SRR315301.719 JOHNLENNON_0014:1:1:11792:1141 length=76/1        -       chr14   32625609  TCCTAAAATAAATATTGTCTCTCCCAACTGTTAAGTTCTAGGTATTGTACTTCCAATTTTAACTTCAGAACCAAGA    CCGCFCADD@DA2=8=<BD8DBEBEE@GBDDHHBDHFG>FGAECCCDEDGGDHHHHBBGGGGG@DGDGGGGHGHHD      0
SRR315301.803 JOHNLENNON_0014:1:1:1949:1141 length=76/1 +       chr5    75537629        GTTTTCTCCACATCAATCCAGCTTTGGAGACATTCTATTAGTGACATATGCCCCTTCCCCCAAAAACAACAATGAA      GGBGGHHHB@GGGEGHH<HBDAGGDEC>>A?=A??GDEBCC<C>@BEBCED:DCDEEEDAB@1:@;9325249823      0       41:C>T
SRR315301.803 JOHNLENNON_0014:1:1:1949:1141 length=76/2 -       chr5    75537797        TTCAAACNGGACAGTCACCAGAAGTACACNNTTATCAAAAATGCACACACTTCACTTGGCATCTCCAGCACCTTCA      ###########B:>+B?CC??78353<79##?:B?B?AAEA<;?;B9@.A5AF;F@CCCFC=B;@BEFE<FIGIII      0       45:G>N,46:A>N,66:T>G,68:T>N

可能性としては、ヒットの確認使っているビューワ等が、ffリードのペアに対応していないとか。

回答日 Apr 21 '12 at 02:16

nob_fj's gravatar image

nob_fj ♦
50761328

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 touran's gravatar image

こちらこそ、ChIPのノウハウは乏しいため、貴重な情報提供ありがとうございました。また、お役に立てて何よりです。

(Apr 23 '12 at 15:12) nob_fj ♦ nob_fj's gravatar image

一応誤解のないように。IP(免疫沈降)といってもChIPではなくMeDIPやBrdU-IPのように抗体の抗原がDNA自体になる場合は、長いDNAの方が抗原となるものが多く含まれる率が高くなるため、飽和していない条件では長鎖の断片の方が優先的にIPされてくる性質があります。

(Apr 24 '12 at 13:23) touran touran's gravatar image

誤解しておりました。DNA自体を抗原にした場合のIP独自の問題だったのですね。大変参考になります。

(Apr 24 '12 at 16:41) nob_fj ♦ nob_fj's gravatar image

直接の解答ではないですが、ご参考までにですが、 入力ファイルの先頭10000エントリーとかを取り出して、小さいサブデータセットを作り、 考えうる選択肢でマップしてみて、どれがうまくいくか試すのは、パラメータの検討や、確認でよくやります。 たまに、マニュアルの記載の方が間違っている場合もあるらしいですし。コマンドでやるなら以下のようなかんじでしょうか。 比較はSAMからBAMにして、IGVやら、samtoolsのtviewでやるなり、いろいろ考えられます。 Galaxyで同じことやるやり方は良く知りません。

head -n 40000 file1.fastq > file1.ext.fastq
head -n 40000 file2.fastq > file2.ext.fastq
bowtie -fr index.fasta -1 file1.ext.fastq -2 file2.ext.fastq > out.fr.bowtie
bowtie -rf index.fasta -1 file1.ext.fastq -s file2.ext.fastq > out.rf.bowtie
bowtie -ff index.fasta -1 file1.ext.fastq -s file2.ext.fastq > out.ff.bowtie

回答日 Apr 16 '12 at 14:22

nob_fj's gravatar image

nob_fj ♦
50761328

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 touran's gravatar image

こちらこそ、貴重な情報報告いただいて、大変参考になりました。上記のスクリプトを貼りつけた後に、よくよく考えると、40000では、少ないこともあるかもと思い、若干心配しておりましたが、わずかながらでもお役に立てて何よりです。

(Apr 20 '12 at 22:02) nob_fj ♦ nob_fj's gravatar image

ちなみに、今更、質問内容を何となく理解しましたが、

  • 入力ファイルはリード1がforward, リード2がreverse
  • シングルエンドでリード1をマップすると、forwardに貼りつく領域では、シングルエンドでリード2をマップするとその近傍の下流にreverseで貼りつく。
  • -ffオプションを使用すると、リード1がforward、リード2もforwardに貼りつくペアしか有効なマップとみなされないため、ほとんどマップ結果が無い。
  • -frオプションでは、リード1がforward、リード2がreverseに貼りつくペアが有効なマップとみなされるため、多数がマップされる。

という理解であっておりますでしょうか。 試していないので、確証はありませんが、リード2だけ予め配列のcomplementを取った上で"-ff"オプションでマップすれば、上記の状況と異なり

  • 入力ファイルはリード1がforward, リード2がforward
  • シングルエンドでリード1をマップすると、forwardに貼りつく領域では、シングルエンドでリード2をマップするとその近傍の下流にforwardで貼りつく。
  • -ffオプションを使用すると、リード1がforward、リード2もforwardに貼りつくペアが有効なマップとみなされるため、多数がマップされる。

となるのでないでしょうか。

逆鎖の取得はfastx tool kitなどでもできますし、

$ /toolDir/fastx_toolkit-0.0.13/bin/fastx_reverse_complement -i file2.ext.fastq -o file2.ext.comp.fastq -Q 33

Illuminaのfastqであれば、以下のようなperlのワンライナーでも(多分問題なく)変換できます。

$ perl -lne 'if($. % 4==2){$_=reverse;tr/ATGC/TACG/;}elsif($. % 4 == 0){$_=reverse;};print;' file2.ext.fastq > file2.ext.comp.fastq_

いちおう、上記の結果は私が試した限りでは、結果は一致しました。

$ diff file2.ext.comp.fastq file2.ext.comp.fastq_
$

上記のいずれか、または他の適当な手段でリード2のファイルのみ相補鎖に変換した上で"-ff"でマップとかならうまくいくかもしれません。

bowtie -ff index.fasta -1 file1.ext.fastq -s file2.ext.comp.fastq > out.ff.2comp.bowtie

※余談ですが、fastx tool kitを使う場合は、最近の配列は-Q 33が必要だと思いますが、そのオプションをつけないとエラーが出るときだけつけて下さい。tool kitのバージョンや、配列のバージョンによって、いらない場合もあるかもしれません。

回答日 Apr 20 '12 at 22:51

nob_fj's gravatar image

nob_fj ♦
50761328

edited Apr 20 '12 at 22:57

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 touran's gravatar image

理解しました。理解できておらず、こちらこそお手数おかけしました。何某かのBowtieの仕様なのでしょうね。気になるようなら開発者に連絡してみるのも一つの手かもしれませんね。 -Q 33は、fastx toolkitの使用時のオプションで、galaxy経由で実行する場合は、関係ない気がします。上のスクリプトの隠れた部分に書いてあります。質問自体には直接関係ないので、忘れて下さい。

(Apr 20 '12 at 23:41) nob_fj ♦ nob_fj's gravatar image
あなたの回答
プレビューをトグルする

この質問をフォローする

By Email:

Once you sign in you will be able to subscribe for any updates here

By RSS:

回答

回答とコメント

タグ:

×6
×3
×1
×1

質問日: Apr 14 '12 at 15:22

閲覧数: 4,770 回

最終更新日: Apr 24 '12 at 16:41

powered by OSQA