How to count the number of mapped read in 100-bp window from a BAM/SAM file 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?
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 &gt;=1) and the 500 in ($4 &lt;=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") &amp;&amp; ($4 &gt;=1) &amp;&amp; ($4 &lt;=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 &lt; 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