text: user161: I have been using 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?
  
  user77: 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
  
  
original_text: 161: <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>
   77: <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>
  