tophat / cufflinks and annotation

[jwhite]

Hello,

Both tophat and cufflinks accept annotation files as arguments (in the form of GTF files). I am processing 30 samples from mouse RNA-seq data, but find NO annotation in my output files. I have checked the following:

1) that the BAM files and GTF files use the same chromosome nomenclature, ie chrN
2) that the reference genome files use the same nomenclature, ie chrN
3) that the GTF file is usable by tophat, ie gtf_juncs <gtf_name> yields appropriate results.

I expected to see gene symbols or refseq IDs in the output, but found neither. Below is and example of one of the transcripts.gtf output files, and the first few lines of the GTF file. Shouldn't the annotation be added to the transcripts.gtf file if the GTF option is provided? Or am I missing something here?

Is there a tool in the tophat or cufflinks packages that will add the annotation? Or will I have to write that myself?

Thanks for any help you can supply.
Cheers,
Joe


transcripts.gtf

chr1 Cufflinks transcript 4481088 4483630 1000 - . gene_id "CUFF.1"; transcript_id "CUFF.1.1"; FPKM "0.1707050811"; frac "1.000000"; conf_lo "0.098790"; conf_hi "0.240801"; cov "2.589206";
chr1 Cufflinks exon 4481088 4482749 1000 - . gene_id "CUFF.1"; transcript_id "CUFF.1.1"; exon_number "1"; FPKM "0.1707050811"; frac "1.000000"; conf_lo "0.098790"; conf_hi "0.240801"; cov "2.589206";
chr1 Cufflinks exon 4483181 4483630 1000 - . gene_id "CUFF.1"; transcript_id "CUFF.1.1"; exon_number "2"; FPKM "0.1707050811"; frac "1.000000"; conf_lo "0.098790"; conf_hi "0.240801"; cov "2.589206";
chr1 Cufflinks transcript 4766459 4775776 1000 - . gene_id "CUFF.2"; transcript_id "CUFF.2.1"; FPKM "1.0690828163"; frac "0.790339"; conf_lo "0.749463"; conf_hi "1.370969"; cov "15.665388";
chr1 Cufflinks exon 4766459 4766882 1000 - . gene_id "CUFF.2"; transcript_id "CUFF.2.1"; exon_number "1"; FPKM "1.0690828163"; frac "0.790339"; conf_lo "0.749463"; conf_hi "1.370969"; cov "15.665388";
chr1 Cufflinks exon 4767606 4767729 1000 - . gene_id "CUFF.2"; transcript_id "CUFF.2.1"; exon_number "2"; FPKM "1.0690828163"; frac "0.790339"; conf_lo "0.749463"; conf_hi "1.370969"; cov "15.665388";
[...]

mm9/genes.gtf

chr1 unknown exon 3204563 3207049 . - . gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown stop_codon 3206103 3206105 . - . gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown CDS 3206106 3207049 . - 2 gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown CDS 3411783 3411982 . - 1 gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown exon 3411783 3411982 . - . gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown CDS 3660633 3661429 . - 0 gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown exon 3660633 3661579 . - . gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown start_codon 3661427 3661429 . - . gene_id "Xkr4"; gene_name "Xkr4"; p_id "P1298"; transcript_id "NM_001011874"; tss_id "TSS1978";
chr1 unknown exon 4280927 4283093 . - . gene_id "Rp1"; gene_name "Rp1"; p_id "P11404"; transcript_id "NM_001195662"; tss_id "TSS19117";
chr1 unknown stop_codon 4283062 4283064 . - . gene_id "Rp1"; gene_name "Rp1"; p_id "P11404"; transcript_id "NM_001195662"; tss_id "TSS19117";
[...]

[jphekman]

I am pretty new to Cufflinks as well, but since no one else has jumped in yet, I will take a shot at this.

I do see gene_id values in my transcripts.gtf files:

13 Cufflinks transcript 37401708 37410554 1000 -. gene_id "EPPK1"; transcript_id "XM_532347.2"; FPKM "1.6826560109"; frac "1.000000"; conf_lo "1.436585"; conf_hi "1.928727"; cov "2.236443"; full_read_support "yes";

My command line was:

cufflinks -p 20 -g Annotation/Genes/genes.gtf -o cufflinks-out tophat-juncs/accepted_hits.bam

Did you use -g or -G? When I initially used -G I also found that my output had no gene IDs in it.

Hope this helps,
Jessica

[SrCardgage]

think of the annotation file as a guide for what to quantify

I don't think the point of supplying the annotation files is to copy the annotation data into the output file(s). The reason I use them is to stop cufflinks from attempting to build new transcripts. Instead, I want it to stick to the known genes and quantify all their transcripts. Supplying this annotation is how to tell tophat and cufflinks what are the known genes are.

Running it this way (with GTF flag and file) should also decrease the time it takes both programs to execute. However, that should not be the reason you choose to use these flags, unless you are just prototyping the programs on a new system.

[jphekman]

Interesting. I'm curious -- how would you annotate cufflinks' output?

Jessica

[SrCardgage]

programmer's overview of cufflinks annotation

Just answering off the cuff; personally, I would study the file structure and write a script to cross-reference the gene-name or genomic location (which ever is in the Cuff* files) to RefSeq database (or gene database of your choice).

If you're not proficient with scripting, there are likely tools to do this. Look into cummeRbund, briefly described in this paper on the Tuxedo pipeline:

http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html

If you plan on doing a lot of analysis with Cufflinks et. al., I highly suggest becoming proficient with unix and scripting (python is a good language).

[jwhite]

I don't think the point of supplying the annotation files is to copy the annotation data into the output file(s). The reason I use them is to stop cufflinks from attempting to build new transcripts. Instead, I want it to stick to the known genes and quantify all their transcripts. Supplying this annotation is how to tell tophat and cufflinks what are the known genes are.

Running it this way (with GTF flag and file) should also decrease the time it takes both programs to execute. However, that should not be the reason you choose to use these flags, unless you are just prototyping the programs on a new system.

Thank you for pointing out that cufflinks may attempt to build new transcripts. However, I was looking for annotation to be added to my cufflinks output. Its clear that this can be done; I just needed to learn how. The -g option seems to be the answer.

jwhite

[SrCardgage]

I have to disagree. From the documentation:

"-g/--GTF-guide <reference_annotation.(gtf/gff)> Tells Cufflinks to use the supplied reference annotation (GFF) to guide RABT assembly. Reference transcripts will be tiled with faux-reads to provide additional information in assembly. Output will include all reference transcripts as well as any novel genes and isoforms that are assembled. "

RABT assembly -> Reference Annotation Based Transcript assembly

Which is described in this link.

http://cufflinks.cbcb.umd.edu/howitworks.html#hrga

I'm sorry to sound arrogant, but you really need to sit down and pore over the documentation. Most of what you need to know is in there. The cufflinks step is more like bookkeeping for later stages. Annotation is not going to happen until the very end, when running the cummeRbund package in R.

If you are proficient in scripting and databases, you can perform your own annotation by cross-referencing the unique IDs in the cuff* output. But, why do this when the cummeRbund packages has already been written.

This is a very complex program. It's unlikely a person can just breeze through doc and get things to come out correctly. Give yourself a week to just READ. It will pay off in a huge way down the road.

Best of luck!