<p>Maybe there's one tool that does everything, but here's one approach that uses a mix of Unix tools and BEDOPS.</p>

<p>At the start, one thing that may be necessary is to fix GTF files that lack an attribute required by specification:</p>

<pre><code>$ awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' annotations.gtf &gt; annotations.fixed.gtf
</code></pre>

<p>This may not be an issue for your GTF file.</p>

<p>With a correctly-formatted GTF, we can filter it by <code>gene_biotype</code> attribute, e.g., <code>snoRNA</code>, and convert those filtered annotations to a sorted BED file with <code>gtf2bed</code>:</p>

<pre><code>$ NEEDLE="snoRNA"
$ grep -v '^#' annotations.fixed.gtf | awk -vneedle=${NEEDLE} '{ match($0, /gene_biotype "([a-zA-Z_]+)"/, a); if ((needle==a[1]) &amp;&amp; ($8=="gene")) { print $0; } }' | gtf2bed - &gt; ${NEEDLE}.bed
</code></pre>

<p>Secondly, we can convert the BAM files to sorted BED files via <code>bam2bed</code>:</p>

<pre><code>$ bam2bed &lt; reads.bam &gt; reads.bed
</code></pre>

<p>Once we have <code>snoRNA.bed</code> and <code>reads.bed</code>, we can use <code>bedmap --count</code> to count the number of reads which overlap each snoRNA-annotated gene from the original annotations:</p>

<pre><code>$ bedmap --count snoRNA.bed reads.bed &gt; answer.txt
</code></pre>

<p>One complication is that the reads should have a chromosome name scheme identical to that of your GTF file, so that mapping can be done between elements on the same chromosome. Mixing Ensembl and UCSC sourced files, for example, could require an extra step to add/remove <code>chr</code> prefix.</p>

<p>In any case, each line of <code>answer.txt</code> is the number of reads that overlap a snoRNA-typed annotation, not the overall sum.</p>

<p>To sum all values, you can add to the <code>bedmap</code> statement an <code>awk</code> command that runs a simple accumulator:</p>

<pre><code>$ bedmap --count snoRNA.bed reads.bed | awk 'START{s=0;}{s+=$0;}END{print s;}' &gt; sum.txt
</code></pre>

<p>The file <code>sum.txt</code> will have the sum of reads overlapping annotations.</p>

<p>This procedure could be repeated for other gene biotype categories, <code>snRNA</code>, <code>rRNA</code>, etc. by making a script that swaps in a category name from a list of such, rerunning the GTF conversion and filtering step, and re-running the <code>bedmap</code> step on the new annotation subset.</p>

CC BY-SA 3.0
2018-02-07T20:04:42.973Z
Hi Alex,
I like the approach and was doing this only but thought there might be a tool already doing this. THanks for your answer. Really appreciate it.I am using featureCounts from subread package to count reads for each gene and then will query that gene to it's biotype from gtf
user3138373
CC BY-SA 3.0
2018-02-07T19:19:10.087Z
2018-02-07T19:19:10.087Z
Alex Reynolds
<p>You may be interested in a Python library I developed called <a href="https://github.com/SACGF/pyreference" rel="nofollow noreferrer">PyReference</a></p>
<pre><code># Install
python3 -m pip install pyreference
# Process GTF into a JSON file
pyreference_gff_to_json.py --gtf gencode.hg38.gtf
# Create a setup file (see website for details)
# Run against your BAM file
pyreference_biotype.py ${BAM_FILE}
</code></pre>
<p>This generates a graph:</p>
<p><a href="https://i.sstatic.net/Tsjr3.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Tsjr3.jpg" alt="enter image description here" /></a></p>
<p>As well as a CSV of read counts per biotype for each length. You can sum up the read lengths for a total.</p>

CC BY-SA 4.0
2022-01-12T03:53:13.553Z
2022-01-13T03:54:04.760Z
2022-01-13T03:54:04.760Z
David Lawrence
<p>I have bam files. I want to find the total read counts associated with all the biotypes, eg snRNA,rRNA,tRNA mRNA,scRNA,snoRNA etc. I can use ht-seq count to get read  counts for the genes, but is there a tool which can directly sum up counts for each of the above category. I have a gencode hg38 gtf file as my annotation file.</p>

<p>Thanks</p>

<sam><reads><gtf>
Counting reads for each biotype
# Counting reads for each biotype

(Asked by: username_0 on 2018-02-07)

I have bam files. I want to find the total read counts associated with all the biotypes, eg snRNA,rRNA,tRNA mRNA,scRNA,snoRNA etc. I can use ht-seq count to get read counts for the genes, but is there a tool which can directly sum up counts for each of the above category. I have a gencode hg38 gtf file as my annotation file.

Thanks

**Tags:** sam, reads, gtf

---

## 2 answers

### Answer by username_1 on 2018-02-07. Score: 1 | Accepted answer

Maybe there's one tool that does everything, but here's one approach that uses a mix of Unix tools and BEDOPS.
At the start, one thing that may be necessary is to fix GTF files that lack an attribute required by specification:
$ awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' annotations.gtf > annotations.fixed.gtf

This may not be an issue for your GTF file.
With a correctly-formatted GTF, we can filter it by gene_biotype attribute, e.g., snoRNA, and convert those filtered annotations to a sorted BED file with gtf2bed:
$ NEEDLE="snoRNA"
$ grep -v '^#' annotations.fixed.gtf | awk -vneedle=${NEEDLE} '{ match($0, /gene_biotype "([a-zA-Z_]+)"/, a); if ((needle==a[1]) && ($8=="gene")) { print $0; } }' | gtf2bed - > ${NEEDLE}.bed

Secondly, we can convert the BAM files to sorted BED files via bam2bed:
$ bam2bed < reads.bam > reads.bed

Once we have snoRNA.bed and reads.bed, we can use bedmap --count to count the number of reads which overlap each snoRNA-annotated gene from the original annotations:
$ bedmap --count snoRNA.bed reads.bed > answer.txt

One complication is that the reads should have a chromosome name scheme identical to that of your GTF file, so that mapping can be done between elements on the same chromosome. Mixing Ensembl and UCSC sourced files, for example, could require an extra step to add/remove chr prefix.
In any case, each line of answer.txt is the number of reads that overlap a snoRNA-typed annotation, not the overall sum.
To sum all values, you can add to the bedmap statement an awk command that runs a simple accumulator:
$ bedmap --count snoRNA.bed reads.bed | awk 'START{s=0;}{s+=$0;}END{print s;}' > sum.txt

The file sum.txt will have the sum of reads overlapping annotations.
This procedure could be repeated for other gene biotype categories, snRNA, rRNA, etc. by making a script that swaps in a category name from a list of such, rerunning the GTF conversion and filtering step, and re-running the bedmap step on the new annotation subset.

> **Comments:**
> **username_0**: Hi Alex, I like the approach and was doing this only but thought there might be a tool already doing this. THanks for your answer. Really appreciate it.I am using featureCounts from subread package to count reads for each gene and then will query that gene to it's biotype from gtf  

---

### Answer by username_2 on 2022-01-12. Score: 0

You may be interested in a Python library I developed called PyReference
# Install
python3 -m pip install pyreference
# Process GTF into a JSON file
pyreference_gff_to_json.py --gtf gencode.hg38.gtf
# Create a setup file (see website for details)
# Run against your BAM file
pyreference_biotype.py ${BAM_FILE}

This generates a graph:

As well as a CSV of read counts per biotype for each length. You can sum up the read lengths for a total.


---

