text: user2040: 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
  
  user776: 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.
  
original_text: 2040: <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>
   776: <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>
  