Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-11T07:33:22.007Z
Answers.Text: Thanks. I work in `Python` and `R`. So python scripts will better help me. So your approach may help others who work in java.
Answers.UserDisplayName: Lot_to_learn
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-10T14:26:20.953Z
Answers.LastActivityDate: 2018-12-10T14:26:20.953Z
Answers.UserDisplayName: Pierre
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-09T09:00:57.797Z
Answers.Text: 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.
Answers.UserDisplayName: Lot_to_learn
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-09T23:51:30.423Z
Answers.Text: This updated answer seems like it is what you are actually after, right?
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-11T07:29:51.110Z
Answers.Text: 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.)?
Answers.UserDisplayName: Lot_to_learn
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-11T07:32:20.987Z
Answers.Text: 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.
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-09T05:35:24.907Z
Answers.LastActivityDate: 2018-12-10T00:02:51.543Z
Answers.LastEditDate: 2018-12-10T00:02:51.543Z
Answers.UserDisplayName: conchoecia
Body: <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>
  
Tags: <phylogenetics><sam><ngs><samtools><read-mapping>
Title: How to count the number of mapped read in 100-bp window from a BAM/SAM file
ThreadText: # 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](https://pysam.readthedocs.io/en/latest/api.html).
  
  ```python
  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.](https://pysam.readthedocs.io/en/latest/faq.html#pysam-coordinates-are-wrong)
  
  ```python
  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](http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html) 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
  
  ```
  
  > **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.  
  
  ---
  
  