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.



2 Answers 2


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:

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

  • $\begingroup$ 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 $\endgroup$ Commented Feb 7, 2018 at 20:04

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:

enter image description here

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.