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

CC BY-SA 3.0
2017-05-30T09:08:11.583Z
2017-05-30T09:08:11.583Z
Devon Ryan
<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>

<rna-seq><sam><samtools><star><pysam>
Quantifying reads mapping to multiple loci
# Quantifying reads mapping to multiple loci

(Asked by: username_0 on 2017-05-30)

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?

**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


---

