Answers.Body: <p><code>Bedtools</code> is not the best tool for this. You can use <code>featureCounts</code>, which gives you the number of read that overlap with features in a genome. In this case, your features are genes.</p>
  <p><code>bedtools coverage</code> 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 <code>featureCounts</code> does.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-02-28T01:45:31.637Z
Answers.LastActivityDate: 2021-02-28T01:51:52.457Z
Answers.LastEditDate: 2021-02-28T01:51:52.457Z
Answers.UserDisplayName: Phoenix Mu
Body: <p>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,</p>
  <pre><code>bedtools coverage -a outCSEQSAligned.out.bam -b GCF_000746645.1_ASM74664v1_genomic.bed
  </code></pre>
  <p>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:</p>
  <pre><code>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
  </code></pre>
  <p>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.</p>
  <p>Thanks for the help!</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2021-03-01T04:42:34.607Z
Comments.Text: Did you do an RNASeq experiment?  It's not clear from what you wrote.
Comments.UserDisplayName: swbarnes2
Tags: <phylogenetics><gene-expression><bedtools>
Title: How can I get bedtools to tell me which genes are being expressed?
ThreadText: # How can I get bedtools to tell me which genes are being expressed?
  
  (Asked by: username_0 on 2021-02-28)
  
  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!
  
  **Tags:** phylogenetics, gene-expression, bedtools
  
  > **Comments:**
  > **username_1**: Did you do an RNASeq experiment?  It's not clear from what you wrote.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2021-02-28. Score: 0 | Accepted answer
  
  `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.
  
  
  ---
  
  