# 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?

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  • Thanks for your reply with proper examples. But, it seems I have missed something in question. By a window, I mean that it is a sliding non-overlapping window of width 500 bp (not just one region). Once again thanks for your answer. Can you extend your answer for sliding non-overlapping windows of given width. Dec 9 '18 at 9:00 • This updated answer seems like it is what you are actually after, right? Dec 9 '18 at 23:51 • Thanks. Can I call number of uniq_reads in a given window as Read Depth of that window (because final results is clearly telling the this many reads are supporting particular region of reference genome.)? Dec 11 '18 at 7:29 • It is pretty similar to "Read Depth", but that makes more sense on a per-position basis. I think it makes sense just to call it what it is -- the number of unique reads mapping to each 500-bp window. Dec 11 '18 at 7:32 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

• Thanks. I work in Python and R. So python scripts will better help me. So your approach may help others who work in java. Dec 11 '18 at 7:33