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?


1 Answer 1


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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