Answers.Body: <p>You almost had the correct python code already, you just need to filter out secondary alignments:</p>
  
  <pre><code>def get_reads_hist(bam): 
      bam = pysam.AlignmentFile(bam, 'rb')
      counts = Counter()
      for query in bam.fetch():
          if query.is_secondary:
              continue
          nh_count = Counter(dict(query.get_tags())['NH'])
          counts += nh_count            
      return counts
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-30T09:08:11.583Z
Answers.LastActivityDate: 2017-05-30T09:08:11.583Z
Answers.UserDisplayName: Devon Ryan
Body: <p>I have been using <a href="https://github.com/alexdobin/STAR" rel="noreferrer">STAR</a> for our RNA-Seq samples. The <code>final.out</code> log file reports percentage of uniquely mapped reads along with percentage of reads that map to <code>multiple loci</code> (less than or equal to 10) and percentage of reads mapping to <code>too many loci</code> (greater than 10). However, I want to break down the <code>multiple loci</code> part to individual counts: Reads mapping to 2 locations, 3 locations, 4 locations .. 10 locations.</p>
  
  <p>The <code>NH</code> tag seems to be used by <code>STAR</code>. However a naive read counting approach results in it reporting more number of reads than total reads.</p>
  
  <p>For example, my <code>final.out</code> looks like this:</p>
  
  <pre><code>   Mapping speed, Million of reads per hour |       1403.36
  
                        Number of input reads |       53015978
                    Average input read length |       26
                                  UNIQUE READS:
                 Uniquely mapped reads number |       368916
                      Uniquely mapped reads % |       0.70%
                        Average mapped length |       26.45
                     Number of splices: Total |       1057
          Number of splices: Annotated (sjdb) |       0
                     Number of splices: GT/AG |       802
                     Number of splices: GC/AG |       1
                     Number of splices: AT/AC |       0
             Number of splices: Non-canonical |       254
                    Mismatch rate per base, % |       0.31%
                       Deletion rate per base |       0.00%
                      Deletion average length |       1.45
                      Insertion rate per base |       0.00%
                     Insertion average length |       1.00
                           MULTI-MAPPING READS:
      Number of reads mapped to multiple loci |       45766732
           % of reads mapped to multiple loci |       86.33%
      Number of reads mapped to too many loci |       3757890
           % of reads mapped to too many loci |       7.09%
                                UNMAPPED READS:
     % of reads unmapped: too many mismatches |       0.00%
               % of reads unmapped: too short |       5.89%
                   % of reads unmapped: other |       0.00%
  </code></pre>
  
  <p>Counting histogram of number of positions a read maps to using <code>pysam</code>:</p>
  
  <pre><code>def get_reads_hist(bam): 
      bam = pysam.AlignmentFile(bam, 'rb')
      counts = Counter()
      for query in bam.fetch():
          nh_count = Counter(dict(query.get_tags())['NH'])
          counts += nh_count            
      return counts
  </code></pre>
  
  <p>results in </p>
  
  <pre><code>Counter({1: 330606,
       2: 86772164,
       3: 329,
       4: 38083,
       5: 31,
       6: 1094,
       7: 129,
       8: 50,
       10: 50})
  </code></pre>
  
  <p>The count <code>1</code> reads are fine even though they do not match the counts in <code>final.out</code> file since I am counting a certain category of reads (say those mapping to <code>tRNA</code> only), but the reads mapping to 2 locations are highly overestimated. Why is that?</p>
  
Tags: <rna-seq><sam><samtools><star><pysam>
Title: Quantifying reads mapping to multiple loci
ThreadText: # Quantifying reads mapping to multiple loci
  
  (Asked by: username_0 on 2017-05-30)
  
  I have been using [STAR](https://github.com/alexdobin/STAR) for our RNA-Seq samples. The `final.out` log file reports percentage of uniquely mapped reads along with percentage of reads that map to `multiple loci` (less than or equal to 10) and percentage of reads mapping to `too many loci` (greater than 10). However, I want to break down the `multiple loci` part to individual counts: Reads mapping to 2 locations, 3 locations, 4 locations .. 10 locations.
  
  The `NH` tag seems to be used by `STAR`. However a naive read counting approach results in it reporting more number of reads than total reads.
  
  For example, my `final.out` looks like this:
  
  ```
     Mapping speed, Million of reads per hour |       1403.36
  
                        Number of input reads |       53015978
                    Average input read length |       26
                                  UNIQUE READS:
                 Uniquely mapped reads number |       368916
                      Uniquely mapped reads % |       0.70%
                        Average mapped length |       26.45
                     Number of splices: Total |       1057
          Number of splices: Annotated (sjdb) |       0
                     Number of splices: GT/AG |       802
                     Number of splices: GC/AG |       1
                     Number of splices: AT/AC |       0
             Number of splices: Non-canonical |       254
                    Mismatch rate per base, % |       0.31%
                       Deletion rate per base |       0.00%
                      Deletion average length |       1.45
                      Insertion rate per base |       0.00%
                     Insertion average length |       1.00
                           MULTI-MAPPING READS:
      Number of reads mapped to multiple loci |       45766732
           % of reads mapped to multiple loci |       86.33%
      Number of reads mapped to too many loci |       3757890
           % of reads mapped to too many loci |       7.09%
                                UNMAPPED READS:
     % of reads unmapped: too many mismatches |       0.00%
               % of reads unmapped: too short |       5.89%
                   % of reads unmapped: other |       0.00%
  
  ```
  
  Counting histogram of number of positions a read maps to using `pysam`:
  
  ```
  def get_reads_hist(bam): 
      bam = pysam.AlignmentFile(bam, 'rb')
      counts = Counter()
      for query in bam.fetch():
          nh_count = Counter(dict(query.get_tags())['NH'])
          counts += nh_count            
      return counts
  
  ```
  
  results in
  
  ```
  Counter({1: 330606,
       2: 86772164,
       3: 329,
       4: 38083,
       5: 31,
       6: 1094,
       7: 129,
       8: 50,
       10: 50})
  
  ```
  
  The count `1` reads are fine even though they do not match the counts in `final.out` file since I am counting a certain category of reads (say those mapping to `tRNA` only), but the reads mapping to 2 locations are highly overestimated. Why is that?
  
  **Tags:** rna-seq, sam, samtools, star, pysam
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2017-05-30. Score: 9 | Accepted answer
  
  You almost had the correct python code already, you just need to filter out secondary alignments:
  
  ```
  def get_reads_hist(bam): 
      bam = pysam.AlignmentFile(bam, 'rb')
      counts = Counter()
      for query in bam.fetch():
          if query.is_secondary:
              continue
          nh_count = Counter(dict(query.get_tags())['NH'])
          counts += nh_count            
      return counts
  
  ```
  
  
  ---
  
  