bowtie2のmultiple hit時の挙動についてのメモ

bowtie2でマルチマップはどう処理されるのか という話。

背景

bowtie2のデフォルトの挙動では配列がmultiple hitしてもそのうち1箇所しか報告されない。ではその1箇所はどう選ぶのかというと、普通はアライメントスコアの最も高いものの中からランダムに、である。ドキュメントにもそう書いてある¹。
…のだが場合によってはランダムにならず固まってしまうらしい²のでその条件をはっきりさせたい。

といってもここ²にほとんど結論が書かれているのだが、この問題はbowtie2の乱数生成に起因している。Bowtie2の取扱説明書³にある--non-deterministicオプションの項に次のような説明がある。

Normally, Bowtie 2 re-initializes its pseudo-random generator for each read. It seeds the generator with a number derived from (a) the read name, (b) the nucleotide sequence, (c) the quality sequence, (d) the value of the --seed option. This means that if two reads are identical (same name, same nucleotides, same qualities) Bowtie 2 will find and report the same alignment(s) for both, even if there was ambiguity. When --non-deterministic is specified, Bowtie 2 re-initializes its pseudo-random generator for each read using the current time. This means that Bowtie 2 will not necessarily report the same alignment for two identical reads. This is counter-intuitive for some users, but might be more appropriate in situations where the input consists of many identical reads.

まとめると、bowtie2はリードごとに乱数生成器を初期化しており、そのシード値は次の4つに依存しているようだ。

• リードの名前
• 塩基配列
• クオリティ配列
• --seedオプションの値

要は、__名前も含めて同一の情報を持つリードは同じ場所にマップされる__ようにできている。
以下ではリードの名前・クオリティ配列・--seedオプションが異なればマッピング結果も異なることを確かめよう。

準備

マルチマップする塩基配列を基に、次のようなFASTQを用意する。（配列等はここ²の検証例を元に作成）

• 名前・塩基配列・クオリティ配列の全てが同一なリードでできたFASTQ
• 塩基配列は同じで、名前・クオリティ配列の{いずれか or 両方}がリードごとに異なるFASTQ

いちいち用意するのは面倒なので次の生成器で代用する。

gen_testseq.py

# !/usr/bin/env python
import argparse
import random

parser = argparse.ArgumentParser()
parser.add_argument("-r", "--reads", type=int, default=10000)
parser.add_argument("-n", "--name", action="store_true")
parser.add_argument("-q", "--quality", action="store_true")
args = parser.parse_args()

for i in range(args.reads):
    print "@test_sequence{}".format(random.random() if args.name else '')
    print "CAAAGTGTTTACTTTTGGATTTTTGCCAGTCTAACAGGTGAAGCCCTGGAGATTCTTATTAGTGATTTGGGCTGGGGCCTGGCCATGTGTATTTTTTTAAA"
    print '+'
    print ''.join([chr(random.randrange(33, 127)) for _ in range(101)]) if args.quality else 'I' * 101

print2行目の配列をhg19にマッピングすると9箇所にmulti hitし、うちアライメントスコアが最も高くマップ候補となるのは次の3箇所となる。

chrom	position	strand
chr1	11635	+
chr2	114359380	-
chr15	102519535	-

bowtie2はバージョン2.2.4を使用した。

検証

リードの情報が全て同一なFASTQの場合

こんな感じでマップ結果からポジションだけ取り出して個数をカウントする。

$ bowtie2 --no-hd -x hg19 -U <(./gen_testseq.py) | cut -f 3-4 | sort | uniq -c

結果

10000 reads; of these:
  10000 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    10000 (100.00%) aligned >1 times
100.00% overall alignment rate
  10000 chr15	102519435

すべて一箇所にマップされた。

--seedオプションを指定してみる

$ bowtie2 --no-hd --seed 17320508 -x hg19 -U <(./gen_testseq.py) | cut -f 3-4 | sort | uniq -c

10000 reads; of these:
  10000 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    10000 (100.00%) aligned >1 times
100.00% overall alignment rate
  10000 chr1	11635

値によっては場所が変わるものの、マップされるのは一箇所だけ。

--non-deterministicオプションを指定してみる

$ bowtie2 --no-hd --non-deterministic -x hg19 -U <(./gen_testseq.py) | cut -f 3-4 | sort | uniq -c

10000 reads; of these:
  10000 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    10000 (100.00%) aligned >1 times
100.00% overall alignment rate
   3302 chr1	11635
   3379 chr15	102519435
   3319 chr2	114359280

今度はランダムにマップされた。もちろん実行する毎に結果は異なる。

名前やクオリティ配列がリードごとに異なる場合

$ bowtie2 --no-hd -x hg19 -U <(./gen_testseq.py -n) | cut -f 3-4 | sort | uniq -c

10000 reads; of these:
  10000 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    10000 (100.00%) aligned >1 times
100.00% overall alignment rate
   3379 chr1	11635
   3294 chr15	102519435
   3327 chr2	114359280

今度はランダムにマップされた。ここで、

$ ./gen_testseq.py -n > test.fastq
$ bowtie2 --no-hd -x hg19 -U test.fastq | cut -f 3-4 | sort | uniq -c

上記のように一度ファイルに落とした場合は実行時毎に同じ結果が得られる。もちろん--seedオプションの値が異なる場合や--non-deterministicオプションを指定した場合はこの限りではない。

クオリティ配列が異なる場合もランダムにマップされる

$ bowtie2 --no-hd -x hg19 -U <(./gen_testseq.py -q) | cut -f 3-4 | sort | uniq -c

10000 reads; of these:
  10000 (100.00%) were unpaired; of these:
    0 (0.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    10000 (100.00%) aligned >1 times
100.00% overall alignment rate
   3351 chr1	11635
   3312 chr15	102519435
   3337 chr2	114359280

結論

• 普通はランダムにマップ先が選ばれると考えてよい（配列に加えてリード名まで重複する状況があまり考えられない）。
• 同一なリードが含まれる場合は--non-deterministicオプションを使用することでランダムにマップされる。
• ただし再現性は失われる。
• マルチマップ先が異なる結果が必要で再現性を確保したい場合は--seedオプションを使う。

参考

うまくばらつかない場合が以前あったらしい²けどバグだったのだろうか

1. -kオプションの項を参照のこと³ ↩

2. Biostars: ATTENTION: bowtie2 and multiple hits / (SEQanswers) ↩ ↩² ↩³ ↩⁴

3. Bowtie2 manual ↩ ↩²