Question

Is It Possible To Map Transcripts To Genes Given That We Have A .Gtf File?

0

Entering edit mode

11.7 years ago

Jordan ★ 1.3k

Hi,

I'm doing RNA-seq analysis and I have a .gtf file produced by cufflinks which includes transcripts with their locations. I wanted to know if there is a way to map these transcript positions to gene id's?

Thanks!

cufflinks gtf gene rna-seq • 5.0k views

ADD COMMENT • link updated 11.7 years ago by Alex Reynolds 36k • written 11.7 years ago by Jordan ★ 1.3k

score 2 · Answer 1 · 2013-03-28

First, to demonstrate the annotation process I will describe, we need some genes. Let's assume you're working with human data and that you visit the GENCODE site to grab GENCODE v15 records in GTF format. You can use BEDOPS gtf2bed and BEDOPS sort-bed to convert this to a sorted UCSC BED-formatted file containing GENCODE genes:

$ wget ftp://ftp.sanger.ac.uk/pub/gencode/release_15/gencode.v15.annotation.gtf.gz
$ gunzip --stdout gencode.v15.annotation.gtf.gz \
    | gtf2bed - \
    | sort-bed - \
    | grep "\tgene\t" \
    > gencode.v15.annotation.gtf.genes.bed

You can convert your Cufflinks-sourced GTF file to BED with the same process: Run it through BEDOPS gtf2bed, sorting it with BEDOPS sort-bed.

Now that you have both datasets in BED format, you can now map the Cuffinks elements against the GENCODE-sourced BED file with BEDOPS bedmap, which returns Cufflinks transcripts that overlap GENCODE genes by one or more bases, e.g.:

$ gtf2bed < cufflinks_transcripts.gtf \
    | sort-bed - \
    | bedmap --echo --echo-map-id - gencode.v15.annotation.gtf.genes.bed \
    > answer.bed

(To make this analysis faster and neater, I merged the Cufflinks-conversion and -mapping steps.)

The file answer.bed contains the annotation results. This file is a sorted BED file that contains lines of the format:

[ transcript1 ] | [ gene1-overlapping-transcript1 ] ; ... ; [ geneA-overlapping-transcript1 ]
[ transcript2 ] | [ gene1-overlapping-transcript2 ] ; ... ; [ geneB-overlapping-transcript2 ]
...
[ transcriptN ] | [ gene1-overlapping-transcriptN ] ; ... ; [ geneZ-overlapping-transcriptN ]

The pipe character (|) delimits BED-formatted transcript elements from the IDs or names of overlapping genes. The semi-colon (;) delimits names of overlapping genes for a given transcript.

You can replace these delimiters if you need to with the --delim and --multidelim operators, but basically this lets you easily post-process the results with Perl, Python, awk — whatever scripting or other language you prefer.

If you want the entire gene record and not just the gene name, use --echo-map in place of --echo-map-id. Other --echo-map-* operators are available to recover different attributes of the mapped element. For example, if you want a non-redundant list of gene name annotations, use --echo-map-id-uniq to remove duplicates.

The default criteria for overlap between transcript and mapped gene is one or more bases. If you want this to be more stringent, review the overlap criteria section of the BEDOPS bedmap documentation.

This demo assumes that we are working with human data, but the process for annotation is identical regardless of organism: Get the reference and mapping inputs (here, transcripts and genes, resp.) into sorted BED files, and then process them with BEDOPS bedmap.