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?
$\begingroup$
$\endgroup$
1
3 Answers
$\begingroup$
$\endgroup$
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.
$\begingroup$
$\endgroup$
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.
$\begingroup$
$\endgroup$
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
bedtools intersectorbedops bedmap$\endgroup$