bowtieが終わったあとに表示される以下の様なリード数が表示されるがこれをあとでコマンドから確認する方法
135582741 reads; of these:
135582741 (100.00%) were unpaired; of these:
22383842 (16.51%) aligned 0 times
61518065 (45.37%) aligned exactly 1 time
51680834 (38.12%) aligned >1 times
83.49% overall alignment rate
基本的にフラグが立っているのでそれを利用する。
#reads; of these:これはfastQCした結果のhtmlを見るかbowtieにかける前のfastqをwc -lして4で割れば良い。
#aligned 0 times:マップされなかったリード数の表示は
samtools view -cf 4 filename.bam
#マップされたリード数の表示は
samtools view -cF 4 filename.bam
#aligned >1 times:マルチマップされたリード数の表示は
grep 'XS:' ./filename.sam | wc -l
#aligned exactly 1 time:ユニークにマップされたリード数の表示は
マップされたリード数からマルチマップされたリード数を引くか
sed '/XS:/d' Enh1.sam | wc -lとして、全リード数からマルチマップのリードをsamから取り除いた行カウントを行う。その値からsamのheaderのライン数とマップされなかったリード数を引く。
ちなみに、pairをきちんと組んでいるものを抽出するには
この組み合わせ
samtools view -S -f 99 -S filename.sam > ./filename.99.sam
samtools view -S -f 147 -S filename.sam > ./filename.147.sam
この組み合わせ
samtools view -S -f 83 -S filename.sam > ./filename.83.sam
samtools view -S -f 163 -S filename.sam > ./filename.163.sam
http://ppotato.wordpress.com/2010/08/25/samtool-bitwise-flag-paired-reads/がわかりやすい。
<<<<<追記>>>>>>
XS:による方法は非常に限られた状況でのみ成立するっぽいです.
bowtie2でSEでやった時は少なくともbowtieを走らせたあとの結果と一致しました。
ちなみにPEでやったものを同様に行うと、想定よりも大きい値が出たため、PEにはこの方法は適応できません。
XS:i:<N>
Alignment score for the best-scoring alignment found other than the alignment reported. Can be negative. Can be greater than 0 in --local mode (but not in --end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than AS:i.
ってbowtie2のmanualに記載があるけれどよくわかりません( TДT)
SE/PEどちらの場合も、よくいわれるNH:i:のフラグはbowtie2では立ちませんでした。MAPQの値で行う事も出来ませんでした。
確実にmulti alignを取り除く方法知っている方いたら教えてくださいm(_ _ )m
