0
$\begingroup$

I have a file that contains species name, chromosome, start and end locations like the following format hg19,chr8,77778733,77779026. I would like to know some information about the regions 77778733-77779026. How can I quickly find out if it lies in an intron or exon of a gene? Or it is in the upstream or downstream regions of a gene? Is there a tool or web service from where I can get this information?

$\endgroup$
1
  • $\begingroup$ you could probably use bedtools intersect or bedops bedmap $\endgroup$ – user3479780 Apr 28 '20 at 10:18
1
$\begingroup$

How can I quickly find out if it lies in an intron or exon of a gene?

I have an answer on Biostars that might help: https://www.biostars.org/p/124515/#124522

Or it is in the upstream or downstream regions of a gene?

One option is to build a table of annotated bins and then map your regions-of-interest to those annotations.

These annotations would consist of whether the bin is inside a gene, or, after picking the nearest gene by strand, locate the bin to an upstream region of the gene TSS or to an downstream region of the gene TES.

To start, one could download stranded gene annotations:

$ GENCODE_URL=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.basic.annotation.gff3.gz
$ GENCODE_GENES=${PWD}/../data/gencode.v28.genes.bed
$ wget -qO- "${GENCODE_URL}" | gunzip -c | awk '($3=="gene")' | gff2bed - | awk -v FS="\t" -v OFS="\t" '{ n=split($10,a,";"); p=a[4]; n=split(p,a,"="); $4=a[2]; print $0; }' | cut -f1-6 > ${GENCODE_GENES}

Chop up nuclear chromosomes of the genome of interest into 1kb bins (or whatever gradation of bin size you prefer):

$ ASSEMBLY=hg38
$ GENOME=${PWD}/../data/${ASSEMBLY}.bed
$ CHOP_SIZE=1000
$ CHOPPED_GENOME=${PWD}/../data/chopped.${ASSEMBLY}.bed
$ fetchChromSizes ${ASSEMBLY} | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | grep -ve "\.*_\.*" | grep -v "chrM" | sort-bed - > ${GENOME}
$ bedops --chop ${CHOP_SIZE} ${GENOME} > ${CHOPPED_GENOME}

Then assign the midpoints of bins inside a gene, or to the nearest of any particular set of intergenic categories.

Here's a script I wrote which uses the output of closest-features to do that. This could be easily modified for any set of distance categories:

#!/usr/bin/env python

import sys
import numpy as np

'''
Input should be bin midpoints run through closest-elements, against some set of BED6+ (stranded) gene annotations
'''

def assign_to_category(dist):
  if dist < -50000:
    return '-50kb_to_-Inf'
  elif dist < -20000 and dist >= -50000:
    return '-20kb_to_-50kb'
  elif dist < -10000 and dist >= -20000:
    return '-10kb_to_-20kb'
  elif dist < -5000 and dist >= -10000:
    return '-5kb_to_-10kb'
  elif dist < -1000 and dist >= -5000:
    return '-1kb_to_-5kb'
  elif dist < 0 and dist >= -1000:
    return 'TSS_to_-1kb'
  elif dist == 0:
    return 'Gene'
  elif dist > 0 and dist <= 1000:
    return 'TES_to_+1kb'
  elif dist > 1000 and dist <= 5000:
    return '+1kb_to_+5kb'
  elif dist > 5000 and dist <= 10000:
    return '+5kb_to_+10kb'
  elif dist > 10000 and dist <= 20000:
    return '+10kb_to_+20kb'
  elif dist > 20000 and dist <= 50000:
    return '+20kb_to_+50kb'
  elif dist > 50000:
    return '+50kb_to_+Inf'

for in_line in sys.stdin:
  (bin_str, upstream_gene_str, upstream_gene_dist, downstream_gene_str, downstream_gene_dist) = in_line.rstrip().split('|')
  try:
    upstream_gene_dist = int(upstream_gene_dist)
  except ValueError:
    upstream_gene_dist = np.nan
  try:
    downstream_gene_dist = int(downstream_gene_dist)
  except ValueError:
    downstream_gene_dist = np.nan
  upstream_gene_elems = upstream_gene_str.split('\t')
  downstream_gene_elems = downstream_gene_str.split('\t')

  '''
  We have six possible scenarios, depending on the minimum distance and strandedness of closest elements
  '''

  '''
  1) If either distance value is NaN, the bin is assigned to the non-NaN distance element
  '''
  if upstream_gene_dist is np.nan and downstream_gene_dist is np.nan:
    raise SystemError("Bad NaN input on:\n[{}]\n".format(in_line.rstrip()))
  elif upstream_gene_dist is np.nan:
    gene_id = downstream_gene_elems[3]
    if downstream_gene_elems[5] == '+':
      # bin is closer to downstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      sys.stdout.write('{}\t{}\t\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue
  elif downstream_gene_dist is np.nan:
    gene_id = upstream_gene_elems[3]
    if upstream_gene_elems[5] == '+':
      # bin is closer to upstream TES
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to upstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    continue

  '''
  2) If either distance value is zero, the bin is assigned to "Gene" to state overlap with a gene
  '''
  if upstream_gene_dist == 0:
    gene_id = upstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  if downstream_gene_dist == 0:
    gene_id = downstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  '''
  3) If the upstream strand is positive, if the downstream strand is positive, then: 

   -- the upstream distance is a signed (-) distance to the upstream gene's TES 
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  4) If the upstream strand is negative, if the downstream strand is positive, then: 

   -- the upstream distance is a signed (-) distance to the upstream gene's TSS 
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  5) If the upstream strand is positive, if the downstream strand is negative, then: 

   -- the upstream distance is a signed (-) distance to the upstream gene's TES 
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  6) If the upstream strand is negative, if the downstream strand is negative, then: 

   -- the upstream distance is a signed (-) distance to the upstream gene's TES 
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  Sanity check: If we get here, something went wrong
  '''
  raise SystemError("Bad logic on:\n[{}]\n".format(in_line.rstrip()))

To use it, e.g.:

$ MP_LEFT=500
$ MP_RIGHT=-499 
$ MP_LEFT_UNDO=-500
$ MP_RIGHT_UNDO=499
$ CATEGORIZED_BINS=${PWD}/../data/categorized_bins.${ASSEMBLY}.bed
$ ./assign_chops.py < <(bedops --range ${MP_LEFT}:${MP_RIGHT} --everything ${CHOPPED_GENOME} | closest-features --dist - ${GENCODE_GENES}) | sort-bed - | bedops --range ${MP_LEFT_UNDO}:${MP_RIGHT_UNDO} --everything - > ${CATEGORIZED_BINS}

The output from this file annotates bins with the closest feature to a gene. It looks something like this:

$ more categorized_bins.hg38.bed
...
chr1  32000 33000 TES_to_+1kb RP11-34P13.3
chr1  33000 34000 TES_to_+1kb FAM138A
chr1  34000 35000 Gene  FAM138A
chr1  35000 36000 Gene  FAM138A
chr1  36000 37000 Gene  FAM138A
chr1  37000 38000 TSS_to_-1kb FAM138A
chr1  38000 39000 -1kb_to_-5kb  FAM138A
...

For example, the region chr1:33000-34000 is closest to the 1kb window downstream of the TES of FAM138A.

Once you have this file, you can quickly map it to your sorted regions-of-interest file:

$ bedmap --echo --echo-map --delim '\t' regions-of-interest.bed categorized_bins.hg38.bed > answer.bed

The file answer.bed will contain your regions of interest and the annotated bins which overlap them.

$\endgroup$
0
$\begingroup$

How can I quickly find out if it lies in an intron or exon of a gene? Or it is in the upstream or downstream regions of a gene? Is there a tool or web service from where I can get this information?

Since you emphasise speed and suggest a web service, I'm guessing you just want to check a few individual loci. In that case, just use a genome browser like those from UCSC or Ensembl. If you look up those coordinates in UCSC you can see the sequence falls within the last exon of the ZFHX4 protein coding gene and other information.

enter image description here

$\endgroup$
0
$\begingroup$

There is already some very useful answers out there but in my experience Homer does a very efficient job at annotating bed files. This gives you a lot of details about nearest transcription start sites or exonic regions and more.

Be sure to convert your data to bed format first

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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