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
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.)
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
--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-* 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
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