Answers.Body: <p>First of all, you need gene coordinates. You can retrieve them in GTF format from <a href="https://www.gencodegenes.org/human/" rel="nofollow noreferrer">GENCODE</a> or other sources.</p>
  <p>Next, extract TSS and TES +/- 2kb windows from each gene, taking in mind that genes can be in the forward or reverse strand. Something like this:</p>
  <pre><code>awk -F&quot;\t|\&quot;&quot; -vOFS=&quot;\t&quot; '$3==&quot;gene&quot; {
    TSS=($7==&quot;+&quot; ? $4-2000&quot;\t&quot;$4+2000 : $5-2000&quot;\t&quot;$5+2000)
    TES=($7==&quot;+&quot; ? $5-2000&quot;\t&quot;$5+2000 : $4-2000&quot;\t&quot;$4+2000)
    print $1,TSS,$10&quot;_TSS_2kb&quot;,0,$7
    print $1,TES,$10&quot;_TES_2kb&quot;,0,$7
  }' annotation.gtf &gt; tss_tes.bed
  </code></pre>
  <p>Here I used GNU Awk from Bash, but you can write your own solution using Python, R or whatever you like.</p>
  <p>Finally, count reads falling in these intervals with your tool of choice (you mentioned bedtools, which is fine).</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-08-05T11:05:28.227Z
Answers.LastActivityDate: 2022-08-05T11:05:28.227Z
Answers.UserDisplayName: bepoli
Body: <p>I have BAM files from an RNA-Seq experiment and for all genes (or a subset) I want to get the number of reads in regions around the TSS (e.g. 2kbp) and the TES (e.g. 2 kbp) and calculate the ratio between them. I was thinking of first creating these regions for each gene and store them in a BED file and then use <code>bedtools intersect</code> with the BAM files to get the reads in these regions. My questions is how would I make such a BED file that contains only the specified regions downstream of the TSS and specified regions upstream of the TSS?</p>
  
Tags: <sam><ngs><bedtools><bed>
Title: Count reads at specific gene features
ThreadText: # Count reads at specific gene features
  
  (Asked by: username_0 on 2022-08-02)
  
  I have BAM files from an RNA-Seq experiment and for all genes (or a subset) I want to get the number of reads in regions around the TSS (e.g. 2kbp) and the TES (e.g. 2 kbp) and calculate the ratio between them. I was thinking of first creating these regions for each gene and store them in a BED file and then use `bedtools intersect` with the BAM files to get the reads in these regions. My questions is how would I make such a BED file that contains only the specified regions downstream of the TSS and specified regions upstream of the TSS?
  
  **Tags:** sam, ngs, bedtools, bed
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2022-08-05. Score: 3 | Accepted answer
  
  First of all, you need gene coordinates. You can retrieve them in GTF format from [GENCODE](https://www.gencodegenes.org/human/) or other sources.
  
  Next, extract TSS and TES +/- 2kb windows from each gene, taking in mind that genes can be in the forward or reverse strand. Something like this:
  
  ```
  awk -F"\t|\"" -vOFS="\t" '$3=="gene" {
    TSS=($7=="+" ? $4-2000"\t"$4+2000 : $5-2000"\t"$5+2000)
    TES=($7=="+" ? $5-2000"\t"$5+2000 : $4-2000"\t"$4+2000)
    print $1,TSS,$10"_TSS_2kb",0,$7
    print $1,TES,$10"_TES_2kb",0,$7
  }' annotation.gtf > tss_tes.bed
  
  ```
  
  Here I used GNU Awk from Bash, but you can write your own solution using Python, R or whatever you like.
  
  Finally, count reads falling in these intervals with your tool of choice (you mentioned bedtools, which is fine).
  
  
  ---
  
  