How to get output with read counts and read count percentages from bam file?

Question

[bayramliaziz]

I have a bam file with lots of genes included and I would like get read counts and their percentages as a bam file. I tried samtools idxstats but It doesn't give percentages. Is there any tool that can do what I want? If there is I would like to get you opinion. Also my bam file is derived from pair-end reads.

Thanks in advance.

Answer 1

[WouterDeCoster]

I would suggest using featurecounts to count the number of reads per gene, and calculate the percentages yourself, e.g. using R or Python.

It is not entirely clear whether you aligned to the genome or to the transcriptome, or how you obtained the sequencing data. Often alignment to the genome is the best approach to avoid a bias on (incomplete?) annotation. But since you state that samtools idxstats only gave you the counts it looks like you aligned to a transcriptome, in which case my answer may not be exactly what you need. It is also unclear what you aim to achieve in the long run, so we cannot comment if this is the best strategy.

[bayramliaziz]

I aligned to the genome by using bwa aln after that I used samtools to get bam file. I am new to this thing so I am just trying to explore and experiment with different approaches. I'll check out the featurecounts.