text: user2594: 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? 
  
  user776: 
  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.
  
original_text: 2594: <p>I have a file that contains species name, chromosome, start and end locations like the following format
  <code>hg19,chr8,77778733,77779026</code>. I would like to know some information about the regions <code>77778733-77779026</code>. 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? </p>
   776: <blockquote>
    <p><em>How can I quickly find out if it lies in an intron or exon of a gene?</em></p>
  </blockquote>
  
  <p>I have an answer on Biostars that might help: <a href="https://www.biostars.org/p/124515/#124522" rel="nofollow noreferrer">https://www.biostars.org/p/124515/#124522</a> </p>
  
  <blockquote>
    <p><em>Or it is in the upstream or downstream regions of a gene?</em></p>
  </blockquote>
  
  <p>One option is to build a table of annotated bins and then map your regions-of-interest to those annotations. </p>
  
  <p>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.</p>
  
  <p>To start, one could download stranded gene annotations:</p>
  
  <pre><code>$ 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 &gt; ${GENCODE_GENES}
  </code></pre>
  
  <p>Chop up nuclear chromosomes of the genome of interest into 1kb bins (or whatever gradation of bin size you prefer):</p>
  
  <pre><code>$ 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 - &gt; ${GENOME}
  $ bedops --chop ${CHOP_SIZE} ${GENOME} &gt; ${CHOPPED_GENOME}
  </code></pre>
  
  <p>Then assign the midpoints of bins inside a gene, or to the nearest of any particular set of intergenic categories. </p>
  
  <p>Here's a script I wrote which uses the output of <code>closest-features</code> to do that. This could be easily modified for any set of distance categories:</p>
  
  <pre><code>#!/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 &lt; -50000:
      return '-50kb_to_-Inf'
    elif dist &lt; -20000 and dist &gt;= -50000:
      return '-20kb_to_-50kb'
    elif dist &lt; -10000 and dist &gt;= -20000:
      return '-10kb_to_-20kb'
    elif dist &lt; -5000 and dist &gt;= -10000:
      return '-5kb_to_-10kb'
    elif dist &lt; -1000 and dist &gt;= -5000:
      return '-1kb_to_-5kb'
    elif dist &lt; 0 and dist &gt;= -1000:
      return 'TSS_to_-1kb'
    elif dist == 0:
      return 'Gene'
    elif dist &gt; 0 and dist &lt;= 1000:
      return 'TES_to_+1kb'
    elif dist &gt; 1000 and dist &lt;= 5000:
      return '+1kb_to_+5kb'
    elif dist &gt; 5000 and dist &lt;= 10000:
      return '+5kb_to_+10kb'
    elif dist &gt; 10000 and dist &lt;= 20000:
      return '+10kb_to_+20kb'
    elif dist &gt; 20000 and dist &lt;= 50000:
      return '+20kb_to_+50kb'
    elif dist &gt; 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) &lt; 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) &lt; 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) &lt; 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) &lt; 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()))
  </code></pre>
  
  <p>To use it, e.g.:</p>
  
  <pre><code>$ 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 &lt; &lt;(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 - &gt; ${CATEGORIZED_BINS}
  </code></pre>
  
  <p>The output from this file annotates bins with the closest feature to a gene. It looks something like this:</p>
  
  <pre><code>$ 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
  ...
  </code></pre>
  
  <p>For example, the region <code>chr1:33000-34000</code> is closest to the 1kb window downstream of the TES of <code>FAM138A</code>.</p>
  
  <p>Once you have this file, you can quickly map it to your sorted regions-of-interest file:</p>
  
  <pre><code>$ bedmap --echo --echo-map --delim '\t' regions-of-interest.bed categorized_bins.hg38.bed &gt; answer.bed
  </code></pre>
  
  <p>The file <code>answer.bed</code> will contain your regions of interest and the annotated bins which overlap them.</p>
  