This is a test version of Biostars. For the public version, visit https://www.biostars.org.
subset gtf file based on coordinates

Hi, I have a gtf file based on Ensembl chromosomes and gene IDs. I would like to extract the genes in it based on a second bed file. This Bed file though comes from a different source (lncRNA DB - www.noncode.org).

I would like to extract all the genes from the gtf file which falls within the coordinates of the bed files completely or partially.

I would appreciate any advices on how to do it

thanks

Assa

rna-seq bed gtf subset genomic ranges
• 6,460 views • link
updated by ATpoint ATpoint • written by Assa Yeroslaviz Assa Yeroslaviz

2 answers

Use bedtools:

# Get the GTF entries that overlap with -b (-wa = return entire range of -a interval)
bedtools intersect -wa -a file.gtf -b second_file.bed

Using it without -wa will return only the interval part of -a that partially matches -b.

• 0 views • link
written by ATpoint ATpoint

I'd load the GTF file into R and make a GRangesList out of it (split it by gene name). Then you can subsetByOverlaps() that with the BED file (you'll need to load it too). That should keep all of the genes intact. You can then write the resulting GRangesList back to a GTF file.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks for the suggestion, but I think I am not doing it right. :-)

I load the two files 

# load the gtf file
gtf.gr <- import(con ="~/genomes/Caenorhabditis_elegans/Ensembl/WBcel235.2017/Annotation/Genes/genes.gtf", format = "gtf" )
gtf.gr <- sortSeqlevelsgtf.gr) # make sure it is sorted
gtf.gr <- sortgtf.gr)

# load the bed file
bed.NONCODE <- import("NONCODEv5_ce10.lncAndGene.sorted.ensembl.bed", format="bed")
bed.NONCODE <- sortSeqlevels(bed.NONCODE) # make sure it is sorted
bed.NONCODE <- sort(bed.NONCODE)

and split the gtf based on the geneIDs. 

#split by genes to grangeslist
gtf.grl <- splitgtf.gr, mcolsgtf.gr)$gene_id)

Than I try to subset by overlaps, but i don't think it is done correctly. Do I need to subset it with the GRangesList object of with the GRanges onbject?

# subset by overlap
overlap.subset.gr <- subsetByOverlaps(x=gtf.gr, ranges = bed.NONCODE)
overlap.subset.grl <- subsetByOverlaps(x=gtf.grl, ranges = bed.NONCODE)

Does it make a different?

Exporting the GRanges objects to a gtf file using export

#export the file to a gtf format
rtracklayer::export(object = overlap.subset.gr, con = "overlap.subset.gr.gtf", format = "gtf")
rtracklayer::export(object = unlist(overlap.subset.grl), con = "overlap.subset.grl.gtf", format = "gtf")

But I get two different results here with a different number of rows. Which of the two (if any) is correct than?

thnaks Assa

• 0 views • link
written by Assa Yeroslaviz Assa Yeroslaviz

You'll want to subset the GRangesList object, that way the gene as a whole will get selected. You shouldn't need to sort things (OK, maybe the seqlevels, I always forget about them). Make sure to export the .grl one.

• 0 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this