How to count the number of mapped read in 100-bp window from a BAM/SAM file

Question

Although I know how to get total number of mapped read using samtools flagstat (samtools flagstat file_sorted.bam) but I want to count total number of mapped read in a non-overlapping sliding window of fixed size (let's say 500 bp) only with respect to reference (hg19) (assigning each read only once by its start position).

An example of the desired output would be something like:

Chromosome    Region_Start     Region_Stop    Num_Unique_Reads
chr2          1                500            200
chr2          501              1000           202
chr2          1001             1501           230
chr2          .                 .              .
chr2          .                 .              .
chr2          .                 .              .
chr2          8001              8500           20

Can anyone suggest me how to perform that window operation?

Answer 1

one-liner

Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence.

Note, this one-liner does not double-count reads because of uniq.

samtools view file_sorted.bam | \
    awk '{ if (($3 == "hg19") && ($4 >=1) && ($4 <=500)) { \
  print($1) }}' | sort | uniq | wc -l

python 3

If you need something a little more programmatic to get a known region you should consider trying the pysam package.

import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
region_set = set()
start = 100
stop =  500
counter = 0
for read in samfile.fetch('interesting_contig', 100, 500):
    region_set.add(read.query_name)
    counter += 1
print("{} unique reads in [{}, {}]".format(len(region_set), start, stop))
print("{} alignments in [{}, {}]".format(counter, start, stop))

For my bam file, this program outputs:

1618 unique reads in [100, 500]
2487 alignments in [100, 500]

python 3, non-overlapping windows

It is important to note that the samfile.fetch() function conforms to 1-based indexing, consistent with SAM files, unlike pysam's otherwise default 0-based indexing that matches the BAM-indexing.

import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
print(samfile.references)
print(samfile.lengths)
window = 500
print("ref\tstart\tstop\tuniq_reads")
for i in range(len(samfile.references)):
    refname = samfile.references[i]
    seqlen = samfile.lengths[i]
    for j in range(1, seqlen, window):
        stop = j+window-1 if j+window-1 < samfile.lengths[i] else samfile.lengths[i]
        region_set = set()
        for read in samfile.fetch(refname, j, stop):
            region_set.add(read.query_name)
        print("{}\t{}\t{}\t{}".format(refname, j, stop, len(region_set)))

For the small bam file I looked at the output was:

ref         start  stop  uniq_reads
OdonB       0      500   1668
OdonB       501    1000  2647
OdonB       1001   1500  129400
OdonB       1501   2000  598
OdonB       2001   2500  538
OdonB       2501   3000  872
OdonB       3001   3500  561
OdonB       3501   4000  49120
OdonB       4001   4500  48433
OdonB       4501   4806  28121

Answer 2

using a streaming process with bioalcidaejdk rather than running many random-access jumps.

The script: the.code

final int WINDOW_SIZE=500;

// map position to count number of reads
final Map<Integer,Long> pos2count = new TreeMap<>();
String currChrom=null;

for(;;)
    {
    // get next read in stream
    final SAMRecord rec = iter.hasNext()?iter.next():null;
    //ignore 'bad' reads
    if(rec!=null ) {
        if(rec.getReadUnmappedFlag()) continue;
        if(rec.getDuplicateReadFlag()) continue;
        if(rec.isSecondaryOrSupplementary()) continue;
        }
    // loop over the windows before the current 
    final Iterator<Integer> posIter = pos2count.keySet().iterator();
    while(posIter.hasNext())
        {
        final int start = posIter.next();
        final int end = start + WINDOW_SIZE;
        if(rec!=null && rec.getContig().equals(currChrom) && end >= rec.getStart()) break;
        final long count = pos2count.get(start);
        out.println(currChrom+"\t"+start+"\t"+end+"\t"+count);
        //remove this position from pos2count
        posIter.remove();
        }
    //no more read, exit
    if(rec==null) break;

    //retrieve all the windows for the current read
    for(Integer pos : IntStream.rangeClosed(rec.getStart(), rec.getEnd()).
            mapToObj(POS-> WINDOW_SIZE*(int)(POS/(1.0*WINDOW_SIZE))).
            collect(Collectors.toSet()))
        {
        //increment the count of reads
        final Long count = pos2count.get(pos);
        pos2count.put(pos,1L + (count==null?0L:count));
        }


    currChrom = rec.getContig();
    }

execute:

$java -jar dist/bioalcidaejdk.jar -f the.code  src/test/resources/S1.bam

RF01    0   500 38
RF01    500 1000    61
RF01    1000    1500    69
RF01    1500    2000    80
RF01    2000    2500    73
RF01    2500    3000    56
RF01    3000    3500    20
RF02    0   500 39
RF02    500 1000    76
RF02    1000    1500    73
RF02    1500    2000    72
RF02    2000    2500    59
RF02    2500    3000    13
RF03    0   500 34
RF03    500 1000    74
RF03    1000    1500    85
RF03    1500    2000    67
RF03    2000    2500    51
RF03    2500    3000    7
RF04    0   500 31
RF04    500 1000    68
RF04    1000    1500    90
RF04    1500    2000    75
RF04    2000    2500    29
RF05    0   500 47
RF05    500 1000    100
RF05    1000    1500    53
RF05    1500    2000    5
RF06    0   500 49
RF06    500 1000    89
RF06    1000    1500    29
RF07    0   500 58
RF07    500 1000    65
RF07    1000    1500    7
RF08    0   500 58
RF08    500 1000    69
RF08    1000    1500    3
RF09    0   500 58
RF09    500 1000    69
RF09    1000    1500    2
RF10    0   500 53
RF10    500 1000    39
RF11    0   500 50
RF11    500 1000    32