# Parse out exon coordinates from bed file for each gene

This is my bed file for all the exon coordinates I want to take out all the exons of a given gene. Let say I have gene in chr 1 which starts from chr1 11868 12227 so I want to parse out all the exons that comes in between 11868 12227.

This is my small subset:

cat exon.bed | head -10
chr1    11868   12227   +   exon
chr1    11871   12227   +   exon
chr1    11873   12227   +   exon
chr1    12009   12057   +   exon
chr1    12178   12227   +   exon
chr1    12594   12721   +   exon
chr1    12612   12697   +   exon
chr1    12612   12721   +   exon
chr1    12612   12721   +   exon
chr1    12974   13052   +   exon


How do I parse out?

I use mostly R and bit of shell script but I'm not sure if I can use R, may a few lines of perl or shell script can help me solve my problem.

• cross posted: biostars.org/p/292247 Jan 7 '18 at 17:16
• yes i posted earlier in biostar so i didn't have any response so i posted here hoping to get some solution
– kcm
Jan 7 '18 at 17:17

## 2 Answers

There are multiple ways to go about it. On the command line you can make a 1 line BED file:

chr1   11868   12227


And then bedtools intersect with it.

In R, you could load your original BED file and use GenomicRanges:

library("GenomicRanges")
bed = read.delim("foo.bed", header=F)  # Rename this
# N.B., BED files use 0-based coordinates, I've switched to 1-based here
exons = GRanges(seqnames=levels(bed$V1), ranges=IRanges(start=bed$V2+1, end=bed$V3), strand=bed$V4)
# Again, I switched to 1-based coordinates...you don't have to do that
region = GRanges(seqnames="chr1",
ranges=IRanges(start=11869, end=12227),
strand="+")
theExonsYouWant = subsetByOverlaps(exons, region, type="within")


You'll need to either specify a strand in your region or use ignore.strand=T in subsetByOverlaps(). I hope you don't have any genes that overlap on the same strand, since they'll prove problematic.

• thanks for the reply so what I want is I want to parse out the exon coordinates then 11868 12227 + exon chr1 11871 12227 + exon chr1 11873 12227 + exon chr1 12009 12057 + exon chr1 12178 12227 + exon as i you can see 11868 to 12227 there are like 5 exons so i want only the first and the last exon to be taken and then i want to use it for divergent primer..so the way you suggested so i have to kind of define the exon i want to parse every time, is there a better way, because , i want to do the same for all the genes with their exons...
– kcm
Jan 7 '18 at 16:51
• You can use findOverlaps() with something and itself. Jan 7 '18 at 16:56
• okay will try your suggestion
– kcm
Jan 7 '18 at 16:58
• I get this error " invalid class “GRanges” object: 'mcols(x)' is not parallel to 'x'
– kcm
Jan 8 '18 at 19:29
• That should probably be its own question, though my guess is that you have something the wrong length somewhere. Jan 8 '18 at 19:38

Via BEDOPS bedops -n and Unix I/O streams:

$echo -e "chr1\t11868\t12227" | bedops -n 1 exon.bed - > answer.bed  Or, if you have your genes in a BED file called genes.bed: $ bedops -n 1 exon.bed genes.bed > answer.bed


If you have your genes in some other format, like GFF or GTF, you can use gff2bed or gtf2bed, e.g.:

$bedops -n 1 exon.bed <(gff2bed < genes.gff) > answer.bed  Or: $ bedops -n 1 exon.bed <(gtf2bed < genes.gtf) > answer.bed


The file answer.bed will contain exons that do not overlap a gene annotation.