Question

Counting Distinct Mapped Reads in BAM File of Paired-end Sequencing

0

Entering edit mode

Dario Strbenac ★ 1.5k

@dario-strbenac-5916

Last seen 3 days ago

Australia

Is it possible to determine the distinct number of reads in a BAM file using countBam, if some of the fragments have both read pairs mapped and some have just one of the two reads mapped ? I want each situation to be only counted once towards the grand total.

rsamtools paired-end reads countbam • 1.5k views

ADD COMMENT • link updated 7.6 years ago by Martin Morgan 25k • written 7.6 years ago by Dario Strbenac ★ 1.5k

score 1 · Answer 1 · 2016-09-27

1

Entering edit mode

Martin Morgan 25k

@martin-morgan-1513

Last seen 2 days ago

United States

Perhaps use

param = ScanBamParam(
    flag = scanBamFlag(isFirstMateRead=TRUE))
countBam(..., param=param)

to count only the first (i.e., the only read of a single-end read, or the first of a paired end read) read? Perhaps with additional flags?

ADD COMMENT • link 7.6 years ago Martin Morgan 25k

0

Entering edit mode

It wouldn't be possible to easily count multi-mapping reads only once with this approach. I suppose the most general solution would be to do length(unique(scanBam(bamFile, param = ScanBamParam(what = "qname"))[[1]][[1]])).

ADD REPLY • link 7.6 years ago Dario Strbenac ★ 1.5k