Answers to: Bowtieのpair-end mappingのオプションについて

Answer by nob_fj

nob_fj — Sat, 21 Apr 2012 02:16:36 +0900

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

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

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

"リード1そのまま"-"リード2そのまま"のペア
"リード1そのまま"-"リード2の逆鎖"のペア
"リード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つのみです。

ff.12rev.bowtie
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リードのペアに対応していないとか。

Answer by nob_fj

nob_fj — Fri, 20 Apr 2012 22:51:11 +0900

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

入力ファイルはリード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のバージョンや、配列のバージョンによって、いらない場合もあるかもしれません。

Answer by nob_fj

nob_fj — Mon, 16 Apr 2012 14:22:38 +0900

直接の解答ではないですが、ご参考までにですが、入力ファイルの先頭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