<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>

CC BY-SA 4.0
2022-08-05T11:05:28.227Z
2022-08-05T11:05:28.227Z
bepoli
<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>

<sam><ngs><bedtools><bed>
Count reads at specific gene features
# 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 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).


---

