2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ cross posted: biostars.org/p/292247 $\endgroup$ – Pierre Jan 7 '18 at 17:16
  • 1
    $\begingroup$ yes i posted earlier in biostar so i didn't have any response so i posted here hoping to get some solution $\endgroup$ – krushnach Chandra Jan 7 '18 at 17:17
5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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... $\endgroup$ – krushnach Chandra Jan 7 '18 at 16:51
  • $\begingroup$ You can use findOverlaps() with something and itself. $\endgroup$ – Devon Ryan Jan 7 '18 at 16:56
  • $\begingroup$ okay will try your suggestion $\endgroup$ – krushnach Chandra Jan 7 '18 at 16:58
  • $\begingroup$ I get this error " invalid class “GRanges” object: 'mcols(x)' is not parallel to 'x' $\endgroup$ – krushnach Chandra Jan 8 '18 at 19:29
  • $\begingroup$ That should probably be its own question, though my guess is that you have something the wrong length somewhere. $\endgroup$ – Devon Ryan Jan 8 '18 at 19:38
3
$\begingroup$

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.

$\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.