How can I get bedtools to tell me which genes are being expressed?

Question

I'm trying to align the Acinetobacter baumanii genome to a genome. I've already done the alignment, and I want to use bedtools to see which genes are being expressed exactly. When I try running the following command,

bedtools coverage -a outCSEQSAligned.out.bam -b GCF_000746645.1_ASM74664v1_genomic.bed

unfortunately, it doesn't tell me which genes are being expressed at each place but instead tells me that the sequence aligns to the single Acinetobacter chromosome which isn't too helpful:

NZ_CP009257.1   282134  282254  GWNJ-1012:340:GW200702637th:1:1105:15519:10019/1        255     +       282134  282254  0,0,0   1       120,    0,4       120     120     1.0000000

I'm assuming that 282134-282254 are the positions on the chromosome that correspond to that alignment. I suppose I could extract these positions and figure out which genes they correspond to, but I'm wondering if there's an easier way to use bedtools to figure out which genes are being expressed. For reference, outCSEQSAligned.out.bam is my bamfile which I aligned the Acinetobacter genome to, and GCF_000746645.1_ASM74664v1_genomic.bed is a bedfile corresponding to the Acinetobacter genome.

Thanks for the help!

Answer 1

Bedtools is not the best tool for this. You can use featureCounts, which gives you the number of read that overlap with features in a genome. In this case, your features are genes.

bedtools coverage is used to calculate the overlap between each read in BAM and each line in BED. It is not very straightforward to quantify gene expression. Plus bedtools does not account for exon in each gene and summarize them to gene level, like featureCounts does.