Counting reads for each biotype

Question

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

Answer 1

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.

Answer 2

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.