<p>using a streaming process with <a href="http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html" rel="nofollow noreferrer">bioalcidaejdk</a> rather than running many random-access jumps. </p>

<p>The script: <strong>the.code</strong></p>

<pre><code>final int WINDOW_SIZE=500;

// map position to count number of reads
final Map&lt;Integer,Long&gt; pos2count = new TreeMap&lt;&gt;();
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&lt;Integer&gt; posIter = pos2count.keySet().iterator();
    while(posIter.hasNext())
        {
        final int start = posIter.next();
        final int end = start + WINDOW_SIZE;
        if(rec!=null &amp;&amp; rec.getContig().equals(currChrom) &amp;&amp; end &gt;= 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-&gt; 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();
    }
</code></pre>

<p>execute:</p>

<pre><code>$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
</code></pre>

CC BY-SA 4.0
2018-12-11T07:33:22.007Z
Thanks. I work in `Python` and `R`. So python scripts will better help me. So your approach may help others who work in java.
Lot_to_learn
CC BY-SA 4.0
2018-12-10T14:26:20.953Z
2018-12-10T14:26:20.953Z
Pierre
<h2>one-liner</h2>

<p>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 <code>1</code> in <code>(<span class="math-container">$4 &gt;=1)</code> and the <code>500</code> in <code>($</span>4 &lt;=500)</code> to set your window. Change <code>"hg19"</code> to your target sequence.</p>

<p>Note, this one-liner does not double-count reads because of <code>uniq</code>.</p>

<pre><code>samtools view file_sorted.bam | \
    awk '{ if ((<span class="math-container">$3 == "hg19") &amp;&amp; ($</span>4 &gt;=1) &amp;&amp; (<span class="math-container">$4 &lt;=500)) { \
  print($</span>1) }}' | sort | uniq | wc -l
</code></pre>

<h2>python 3</h2>

<p>If you need something a little more programmatic to get a known region you <a href="https://pysam.readthedocs.io/en/latest/api.html" rel="noreferrer">should consider trying the <code>pysam</code> package</a>.</p>

<pre class="lang-python prettyprint-override"><code>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))
</code></pre>

<p>For my <code>bam</code> file, this program outputs:</p>

<pre><code>1618 unique reads in [100, 500]
2487 alignments in [100, 500]
</code></pre>

<h2>python 3, non-overlapping windows</h2>

<p><a href="https://pysam.readthedocs.io/en/latest/faq.html#pysam-coordinates-are-wrong" rel="noreferrer">It is important to note that the <code>samfile.fetch()</code> function conforms to 1-based indexing, consistent with SAM files, unlike pysam's otherwise default 0-based indexing that matches the BAM-indexing.</a></p>

<pre class="lang-python prettyprint-override"><code>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)))
</code></pre>

<p>For the small bam file I looked at the output was:</p>

<pre><code>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
</code></pre>

CC BY-SA 4.0
2018-12-09T09:00:57.797Z
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.
Lot_to_learn
CC BY-SA 4.0
2018-12-09T23:51:30.423Z
This updated answer seems like it is what you are actually after, right?
conchoecia
CC BY-SA 4.0
2018-12-11T07:29:51.110Z
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.)?
Lot_to_learn
CC BY-SA 4.0
2018-12-11T07:32:20.987Z
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.
conchoecia
CC BY-SA 4.0
2018-12-09T05:35:24.907Z
2018-12-10T00:02:51.543Z
2018-12-10T00:02:51.543Z
conchoecia
<p>Although I know how to get total number of mapped read using samtools flagstat (<code>samtools flagstat file_sorted.bam</code>) but I want to count total number of mapped read in a <strong>non-overlapping sliding window of fixed size</strong> (let's say 500 bp) only with respect to reference (hg19) (assigning each read only once by its start position). </p>

<p>An example of the desired output would be something like:</p>

<pre><code>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
</code></pre>

<p>Can anyone suggest me how to perform that window operation?</p>

<phylogenetics><sam><ngs><samtools><read-mapping>
How to count the number of mapped read in 100-bp window from a BAM/SAM file
# How to count the number of mapped read in 100-bp window from a BAM/SAM file

(Asked by: username_0 on 2018-12-08)

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?

**Tags:** phylogenetics, sam, ngs, samtools, read-mapping

---

## 2 answers

### Answer by username_2 on 2018-12-09. Score: 10 | Accepted answer

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

> **Comments:**
> **username_0**: 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.  
> **username_2**: This updated answer seems like it is what you are actually after, right?  
> **username_0**: 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.)?  
> **username_2**: 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.  

---

### Answer by username_1 on 2018-12-10. Score: 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 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 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

> **Comments:**
> **username_0**: Thanks. I work in `Python` and `R`. So python scripts will better help me. So your approach may help others who work in java.  

---

