Bioinformatics Asked on July 27, 2021
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!
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.
Correct answer by Phoenix Mu on July 27, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP