samtools rmdup の挙動に関する実験

samtools rmdup について

samtools rmdupは，potential PCR duplicates を除くのによく使われる便利なツールです．

samtoolsの公式では以下のように説明してあります． > samtools rmdup [-sS] > > Remove potential PCR duplicates: if multiple read pairs have identical external coordinates, only retain the pair with highest mapping quality. In the paired-end mode, this command ONLY works with FR orientation and requires ISIZE is correctly set. It does not work for unpaired reads (e.g. two ends mapped to different chromosomes or orphan reads). > > OPTIONS: > > -s Remove duplicate for single-end reads. By default, the command works for paired-end reads only. > > -S Treat paired-end reads and single-end reads.

しかし，この説明ではオプションを使った場合にどのような挙動をするのかよくわからなかっため，簡単な実験を行いました．

実験

まず，座標が同一な，ペアエンドのリード（ペアの両方または片方のみがマップされているもの）およびシングルエンドのリードから成るテストデータを作成しました．リード名はそれぞれ以下の通り．

test_pe_* はペアエンドでかつ両方のセグメントがマップされているリード（0x1）を表す
test1_peu_* はペアの片方のセグメントのみがマップされているリード（0x1と0x8）を表す
test_se_* はシングルエンドのリードを表す
末尾の0と1は，互いに同じ座標の異なるリードであることを示す

$ samtools view test.bam
test_pe_0     99     chr1     160502     50     101M     =     160543     142     TGGGCTTGCTCTGCTGGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACA     @@@FFFFFHHHGGHHJJIIJJIIJHIIIIJJJGHIJJJJJJJJIJHIJJJJJIJJJIGIJIJJIJIJJJIHHEBB==CBBBBBDDB<BDDDDDD@BDDDD>     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:92A2G5     YT:Z:UU     NH:i:1
test_pe_1     99     chr1     160502     50     101M     =     160543     142     TGGGCTTGCTCTGCTGGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACA     @@@FFFFFHHHGGHHJJIIJJIIJHIIIIJJJGHIJJJJJJJJIJHIJJJJJIJJJIGIJIJJIJIJJJIHHEBB==CBBBBBDDB<BDDDDDD@BDDDD>     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:92A2G5     YT:Z:UU     NH:i:1
test_se_0     0     chr1     160500     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test_se_1     0     chr1     160500     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test1_peu_0     137     chr1     160517     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test1_peu_1     137     chr1     160517     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test_pe_0     147     chr1     160543     50     101M     =     160502     -142     TCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGCCAGTTGATTGACAGCTCAACTGGTTG     >CBBBDDBDDDDBA?CDDDEEDDDDDDDBBBDDBDDDDDDDFFHHFGGHDBJJIIHGIJIIHGGEGIIJJJJIGEIJJIJIJJJIHIJHHGDHDFFFD@C@     AS:i:-12     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:51A2G46     YT:Z:UU     NH:i:1
test_pe_1     147     chr1     160543     50     101M     =     160502     -142     TCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGCCAGTTGATTGACAGCTCAACTGGTTG     >CBBBDDBDDDDBA?CDDDEEDDDDDDDBBBDDBDDDDDDDFFHHFGGHDBJJIIHGIJIIHGGEGIIJJJJIGEIJJIJIJJJIHIJHHGDHDFFFD@C@     AS:i:-12     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:51A2G46     YT:Z:UU     NH:i:1

このテストデータに対してsamtools rmdupを実行します．

$ samtools rmdup test.bam test_rmdup.bam
$ samtools rmdup -s test.bam test_rmdup_s.bam
$ samtools rmdup -S test.bam test_rmdup_S.bam

結果を見てみます．ここではSAM/BAMの１列目のリード名（QNAME）のみを出力しています．

# test.bam
$ samtools view test.bam |awk '{print $1}'
test_pe_0
test_pe_1
test_se_0
test_se_1
test1_peu_0
test1_peu_1
test_pe_0
test_pe_1

# samtools rmdup したもの
$ samtools view test_rmdup.bam |awk '{print $1}'
test_pe_0
test_se_0
test_se_1
test1_peu_0
test1_peu_1
test_pe_0

# samtools rmdup -s したもの
$ samtools view test_rmdup_s.bam |awk '{print $1}'
test_pe_0
test_pe_1
test_se_0
test1_peu_0
test1_peu_1
test_pe_0
test_pe_1

# samtools rmdup -S したもの
$ samtools view test_rmdup_S.bam |awk '{print $1}'
test_pe_0
test_se_0
test1_peu_0
test_pe_0

これらの結果から，

ペアの両方のセグメントがマップされている場合（0x1），デフォルトのときと-Sがついているとき除かれる．
ペアの片方のセグメントのみがマップされている場合（0x1と0x8），-Sがついているとき除かれる
シングルエンドの場合，-sもしくは-Sがついているとき除かれる

ということが分かりました．

なにが挙動を決めているのか

ペアの片方のセグメントのみがマップされている場合に，デフォルトでは除かれないことがわかりました．最初は，FLAG列をみてペアエンドかどうかを判定していると思い込んでいたので驚きました．

テストデータのtest_pe_*とtest1_peu_*を眺めていると，後者では７列目が*となっていました（それまで気付いていませんでした…）． SAM/BAMの7列目のRNEXTはThe SAM Format Specification (PDF)によると，以下のような説明がなされています． > If RNEXT is ‘*’, no assumptions can be made on PNEXT and bit 0x20.

この記述から，rmdupは RNEXTとPNEXTのみを参照しているように思えたのでさらに実験しました．

$ samtools view test2.bam
test_pe_0     99     chr1     160502     50     101M     *     0     0     TGGGCTTGCTCTGCTGGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACA     @@@FFFFFHHHGGHHJJIIJJIIJHIIIIJJJGHIJJJJJJJJIJHIJJJJJIJJJIGIJIJJIJIJJJIHHEBB==CBBBBBDDB<BDDDDDD@BDDDD>     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:92A2G5     YT:Z:UU     NH:i:1
test_pe_1     99     chr1     160502     50     101M     *     0     0     TGGGCTTGCTCTGCTGGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACA     @@@FFFFFHHHGGHHJJIIJJIIJHIIIIJJJGHIJJJJJJJJIJHIJJJJJIJJJIGIJIJJIJIJJJIHHEBB==CBBBBBDDB<BDDDDDD@BDDDD>     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:92A2G5     YT:Z:UU     NH:i:1
test_se_0     0     chr1     160500     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test_se_1     0     chr1     160500     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test1_peu_0     137     chr1     160517     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test1_peu_1     137     chr1     160517     50     101M     *     0     0     GGAATGGAGATTTCGCTCTCAGCTTCTCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGC     @CCFFFFFGFHHHJJJJJJJJIJJJIJIJJJJJIJJJJJJJJJJJJJEHIHIJJJJJJIJJHHHHHHFDDDBDDDDDDD>BDBDDDDDDDDDDDDDDDDDD     AS:i:-10     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:77A2G20     YT:Z:UU     NH:i:1
test_pe_0     147     chr1     160543     50     101M     *     0     0     TCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGCCAGTTGATTGACAGCTCAACTGGTTG     >CBBBDDBDDDDBA?CDDDEEDDDDDDDBBBDDBDDDDDDDFFHHFGGHDBJJIIHGIJIIHGGEGIIJJJJIGEIJJIJIJJJIHIJHHGDHDFFFD@C@     AS:i:-12     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:51A2G46     YT:Z:UU     NH:i:1
test_pe_1     147     chr1     160543     50     101M     *     0     0     TCTGTCTGCCTCCTCTTCTCTGATATCGGTCCGAGTCGTGCGCCGTAAGGGGGTCGAACAGCTCTCCAGCGCAGCCAGTTGATTGACAGCTCAACTGGTTG     >CBBBDDBDDDDBA?CDDDEEDDDDDDDBBBDDBDDDDDDDFFHHFGGHDBJJIIHGIJIIHGGEGIIJJJJIGEIJJIJIJJJIHIJHHGDHDFFFD@C@     AS:i:-12     XN:i:0     XM:i:2     XO:i:0     XG:i:0     NM:i:2     MD:Z:51A2G46     YT:Z:UU     NH:i:1

samtools rmdup test2.bam test_rmdup2.bam
samtools rmdup -s test2.bam test_rmdup_s2.bam
samtools rmdup -S test2.bam test_rmdup_S2.bam

# test2.bam
$ samtools view test2.bam |awk '{print $1}'
test_pe_0
test_pe_1
test_se_0
test_se_1
test1_peu_0
test1_peu_1
test_pe_0
test_pe_1

# samtools rmdup したとき
$ samtools view test_2.bam |awk '{print $1}'
test_pe_0
test_pe_1
test_se_0
test_se_1
test1_peu_0
test1_peu_1
test_pe_0
test_pe_1

# samtools rmdup -s したとき
$ samtools view test__s2.bam |awk '{print $1}'
test_pe_0
test_pe_1
test_se_0
test1_peu_0
test1_peu_1
test_pe_0
test_pe_1

# samtools rmdup -S したとき
$ samtools view test__S2.bam |awk '{print $1}'
test_pe_0
test_se_0
test1_peu_0
test_pe_0

今度は，デフォルトのときでも，ペアの両方のセグメントがマップされているリードが除かれなくなりました．このことから，ペアエンドのリードを除くかどうかの判定には，（FLAGは関係なく）RNEXTおよびPNEXTのみを参照していると考えられました．一方で，-sのときに，ペアの片方のみマップされたリードが除かれないことから，シングルエンドの判定にはFLAGを使っていると考えられます．

まとめ

デフォルト: ペアエンドでペアの両方がマップされるものについて，PCR duplicatesかどうかを判定する
-s: FLAGで0x1 (template having multiple segments in sequencing) がセットされていないリードについて，PCR duplicatesかどうかを判定する
-S: すべてのリードをシングルエンドとして，PCR duplicatesかどうかを判定する

余談ですが，samtoolsの公式には

-S Treat paired-end reads and single-end reads.

と書いてある一方で，例えばこちらのページには

-S Treat paired-end reads as single-end reads.

と書いてあり，この書き方の方が分かりやすい気がしました．

もし勘違いなどあれば，コメントかtwitterでご指摘いただけると幸いです．

Written on July 18, 2013