0
$\begingroup$

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!

$\endgroup$
1
  • $\begingroup$ Did you do an RNASeq experiment? It's not clear from what you wrote. $\endgroup$
    – swbarnes2
    Mar 1 at 4:42
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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