Cufflinks XS tag and -F parameters

[swarbre]

I have RNA-seq data which is strand specific that is the query strand equals the strand of transcription. In the sam output file from Tophat the XS tag is included for spliced transcripts, as I know the strand of the read from the FLAG query strand I went ahead and added the XS tag for unspliced reads. The Cufflinks readme indicates this strand info is utilized in the transcript build therefore i was expecting that the resulting assemblies would all have the strand set in the gtf file. This appears not to be the case and the strand was only set for spliced reads. The XS tags I added appears correct so it seems to be a feature that cufflinks does not output the strand for unspliced assemblies, is this correct?

Also i notice the -F parameter appears to limit two functions --min-isoform-fraction and --pre-mrna-fraction. Ideally I would not want to limit genuine splice variants (i.e. novel splice sites) based on abundance but would like to filter out isoforms with retained introns. Given the F parameter controls both options this doesn't seem possible.

[arthur.yxt]

I also have problem with this tag. I tried to add "XS:A:±" to the end of each line in sam and I tried to replace the original field "XA:i:~" with "XS:A:±", neither worked. cufflink simply prompt these onto the screen:
Processing bundle [ chrY:123123123-234234234] with 5 non-redundant alignments
and nothing written into the 3 result files.

[Cole Trapnell]

Can you send me a gene's worth of SAM alignments with this tag added? This is likely to be a simple fix.

[Cole Trapnell]

Also the "-F" issue mentioned above is a typo in the manual. For pre-mrna-fraction, the short option code is "-j".

[arthur.yxt]

I just found something. I build cufflinks myself and rerun the command like: cufflinks -m 50 -I 1000000 -p 4 -G sample.sam sample.gft [wrongly by mistake]. Error came:
"Error parsing strand (=) from GFF line..."
So I wonder weather GFF or GTF format cufflinks requires?
If is GFF, I tired GFF3 format, didn't work. If is GTF, is it necessary to provide "start_codon" and "stop_codon"? I only have "exon" in the file.

Many thanks!

[Cole Trapnell]

Ah yes, Cufflinks requires GTF (and TopHat is going to switch from GFF3 to GTF soon). You don't need codons in the GTF (they get ignored anyway), just exons.

[swarbre]

sent cole part of the sam file

[marvin.j]

Dear Cole,

First of all: Thank you for cufflinks! I love the idea and am currently trying it out. I also understand the beta-status explains the sparseness of the documentation. But regarding the reference annotations as provided by the -G switch, I had to look at the source-code in order to deduce the details. Are the following observations correct?

* cufflinks expects a two-level hierarchy: exon->mRNA
* mRNA IDs must be unique
* exons can belong to multiple mRNAs
* mRNA entries in the GTF should not have a Parent=... GFF3 attribute
* exons must have this attribute to reference their mRNAs:
e.g. ID=NM_1234.exon1;Parent=NM_1234,NM_11235

After reading the GTF cufflinks seems to group mRNAs into transcription loci, so that a "gene" top-level hierarchy is not necessary or even harmful (as mRNAs should not have Parents).

I would be glad for any clarification on that matter and would volunteer to make an FAQ entry for it and/or a little script to generate cufflinks-compatible GTF from UCSC RefGene table-output (a problem that more people than just me might be confronted with).

Best!
-Marvin

[arthur.yxt]

i guess that exon only belong to one mRNA, but a mRNA may have multiple transcripts

[marvin.j]

Trouble with RefSeq

Well, I tried the stock UCSC->RefSeq Genes->GTF output (mouse ncbi37) after figuring out that cufflinks can also handle the GTF2.2 format (which is what UCSC claims to produce).

Now unfortunately cufflinks crashes with SIGSEGV:

(gdb) run -G ../../../ucsc_mm9_refGene.gtf ../tophat_out/accepted_hits.sam
Starting program: /data/BIO2/mjens/bin/cufflinks -G ../../../ucsc_mm9_refGene.gtf ../tophat_out/accepted_hits.sam
[Thread debugging using libthread_db enabled]
[New Thread 47900392977584 (LWP 28779)]
Counting hits in map
243 duplicate reference transcripts discarded.
        Total map density: 379552491.720572
Warning: large bundle encountered...
Processing bundle [ chr1:9875-24619571 ] with 39148 non-redundant alignments
[New Thread 1075841344 (LWP 31768)]
Warning: nonsense gene merge, skipping

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread 1075841344 (LWP 31768)]
Scaffold::merge (s=@0x401ff870, merged=@0x401ff810, introns_overwrite_matches=<value optimized out>) at scaffolds.cpp:365
365             merged._left = merged._augmented_ops.front().g_left();
(gdb)

Obviously sth. is going very wrong. From some sanity checking on the RefSeq datasets I know that the NM_XYZ transcript IDs are NOT globally unique as required by GTF but can appear multiple times, even referring to loci on different chromosomes. I guess this is because of the cDNA origin of the RefSeq IDs and duplicated genes.

I could fix this by modifying the IDs. But what is the preferred way of letting cufflinks know about exons shared by a number of transcripts (as in isoforms)?
Multiple exon-lines with same coordinates but different transcript_id ?

Please help me out here I guess I just need a hint.

Best!
-Marvin

[bioinfosm]

Hi all,

I am having trouble running cufflinks on tophat output from a lane of RNA-SEQ GA-IIx data. Is there a way I can avoid the tool getting stuck with something that has *so many non-redundant alignments*, maybe just ignore that and move on to complete the analysis?

Any comments are appreciated.

$ ../cufflinks-0.7.0.Linux_x86_64/cufflinks -p 2 berg_42R7WAAXX_300164_41.lane1/accepted_hits.sam
...
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88936024-88936173 ] with 3 non-redundant alignments
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88936699-88936749 ] with 1 non-redundant alignments
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88936940-88936990 ] with 1 non-redundant alignments
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88937389-88942713 ] with 773098 non-redundant alignments
terminate called after throwing an instance of 'std::bad_alloc'
what(): St9bad_alloc
Aborted

[Cole Trapnell]

Hi all,

I am having trouble running cufflinks on tophat output from a lane of RNA-SEQ GA-IIx data. Is there a way I can avoid the tool getting stuck with something that has *so many non-redundant alignments*, maybe just ignore that and move on to complete the analysis?

Any comments are appreciated.

$ ../cufflinks-0.7.0.Linux_x86_64/cufflinks -p 2 berg_42R7WAAXX_300164_41.lane1/accepted_hits.sam
...
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88936024-88936173 ] with 3 non-redundant alignments
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88936699-88936749 ] with 1 non-redundant alignments
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88936940-88936990 ] with 1 non-redundant alignments
Processing bundle [ gi|89161199|ref|NC_000002.10|NC_000002:88937389-88942713 ] with 773098 non-redundant alignments
terminate called after throwing an instance of 'std::bad_alloc'
what(): St9bad_alloc
Aborted

One thing you could try is to increase the value of the collapse-rounds option from it's default of one. Each additional bump should cut the memory use in bundles like this roughly in half (up to a certain point). It carries some risk that Cufflinks will misassemble things if you set it too high, but 2 or 3 should certainly be safe (at least it is in my experience).

[Cole Trapnell]

Well, I tried the stock UCSC->RefSeq Genes->GTF output (mouse ncbi37) after figuring out that cufflinks can also handle the GTF2.2 format (which is what UCSC claims to produce).

Now unfortunately cufflinks crashes with SIGSEGV:

(gdb) run -G ../../../ucsc_mm9_refGene.gtf ../tophat_out/accepted_hits.sam
Starting program: /data/BIO2/mjens/bin/cufflinks -G ../../../ucsc_mm9_refGene.gtf ../tophat_out/accepted_hits.sam
[Thread debugging using libthread_db enabled]
[New Thread 47900392977584 (LWP 28779)]
Counting hits in map
243 duplicate reference transcripts discarded.
        Total map density: 379552491.720572
Warning: large bundle encountered...
Processing bundle [ chr1:9875-24619571 ] with 39148 non-redundant alignments
[New Thread 1075841344 (LWP 31768)]
Warning: nonsense gene merge, skipping

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread 1075841344 (LWP 31768)]
Scaffold::merge (s=@0x401ff870, merged=@0x401ff810, introns_overwrite_matches=<value optimized out>) at scaffolds.cpp:365
365             merged._left = merged._augmented_ops.front().g_left();
(gdb)

Obviously sth. is going very wrong. From some sanity checking on the RefSeq datasets I know that the NM_XYZ transcript IDs are NOT globally unique as required by GTF but can appear multiple times, even referring to loci on different chromosomes. I guess this is because of the cDNA origin of the RefSeq IDs and duplicated genes.

I could fix this by modifying the IDs. But what is the preferred way of letting cufflinks know about exons shared by a number of transcripts (as in isoforms)?
Multiple exon-lines with same coordinates but different transcript_id ?

Please help me out here I guess I just need a hint.

Best!
-Marvin

Can you please send me ([email protected]):

1) The SAM alignments for that bundle
2) The GTF records you're using

I will fix the crash for 0.7.1.

Thanks for the bug report.

[Cole Trapnell]

* cufflinks expects a two-level hierarchy: exon->mRNA

Cufflinks only cares about "exon" records, but they MUST have gene_id and transcript_id attributes

* mRNA IDs must be unique

Correct.

* exons can belong to multiple mRNAs

The Coordinates can, but there must be one record for each one, and the exon records must have distinct transcript_id attributes

* mRNA entries in the GTF should not have a Parent=... GFF3 attribute

Cufflinks support GFF3, so no Parent, ID attributes. Cufflinks ignores mRNA, gene, CDS records, etc. Only the exon records are used. Any code that you see examining these records is due to the fact that the GTF parser came from an earlier GFF3 library, and such code may be removed in the future.
[QUOTE

* exons must have this attribute to reference their mRNAs:
e.g. ID=NM_1234.exon1;Parent=NM_1234,NM_11235

As noted above, exons need to be attached to their parents transcripts, but through the transcript_id attribute, not the ID/Parent tree.

[swarbre]

Just to bring this back to my original query, modifying the sam file so that the XS tag was included for non-spliced reads did not result in the resulting cufflinks gtf file containing any strand info.

For example the following gtf file is produced from the sam file below.

Chr1 Cufflinks transcript 3655 3836 1000 . . gene_id "CUFF.1"; transcript_id "CUFF.1.0"; RPKM "5384615.3846153850"; frac "1.000000"; conf_lo "0.000000"; conf_hi "0.000000"; cov "9.692308";
Chr1 Cufflinks exon 3655 3836 1000 . . gene_id "CUFF.1"; transcript_id "CUFF.1.0"; exon_number "1"; RPKM "5384615.3846153850"; frac "1.000000"; conf_lo "0.000000"; conf_hi "0.000000"; cov "9.692308";

***********************************************


SRR013427.10015398 16 Chr1 1197 255 36M * 0 0 TGAAAAAAGTTGTAATTATTAATGATAGTTCTGTGA 5IIIIIIII0EI?II=II5:IIIII6IIE>IIIIFI NM:i:0 XS:A:-
SRR013429.2412261 0 Chr1 3655 255 36M * 0 0 AGAAAACAAATACATAATCGGAGAAATACAGATTAC IIIIIIIIIIIIII5DI9@(99@)6-.9,0*+(,.% NM:i:0 XS:A:+
SRR013425.9887919 0 Chr1 3658 255 36M * 0 0 AAACAAATACATAATCGGAGAAATACAGATTACAGA IIIIIIIIIIIIIIIIIIII:IIIIIID1IA66=72 NM:i:0 XS:A:+
SRR013412.15349639 0 Chr1 3659 255 36M * 0 0 GCCAAATACATAATCGGAGAAATACAGATTACAGAG 4IIIIIIIIIIIIIIII/I/50A@;1.%+:=/&*%. NM:i:2 XS:A:+
SRR013413.11338591 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGAAATACAGATTACAGAG IIIIIIIIIIIIIIDII?I2IIEBC<I5IBC:F>4E NM:i:0 XS:A:+
SRR013413.5298488 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGAATTACAGATTACAGAG IIIIII6IIIIIII=$IIC%IIII&II@+I'B/54" NM:i:1 XS:A:+
SRR013414.11501860 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGAAATACAGATTACAGAG IIIIIIIIIIIIIIIIIIII-IIIB@I@FII3I</+ NM:i:0 XS:A:+
SRR013414.11712788 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGAAATACAGATTACAGAG IIIIIIIIIIIIIIIIIII>FIIIB>@IEII0<G6F NM:i:0 XS:A:+
SRR013415.7231575 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGAAATACAGATTACAGAG IIIIIIIIIIIIIIIIIIIII?IIIIIIIIIIIC?B NM:i:0 XS:A:+
SRR013415.9176134 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGTAATACAGATTACAGAG IIIIIIIIIIIIIII,I+/+EIDI(II;:4>)1/I? NM:i:1 XS:A:+
SRR013416.3752313 0 Chr1 3659 255 36M * 0 0 AACAAATACATAATCGGAGGAATACAGATTACAGAG IIBIIIIIIIIIIIIII2I%9I/I.FIICIH2/I/I NM:i:1 XS:A:+
SRR013411.14229422 0 Chr1 3661 255 36M * 0 0 CAAATACATAATCGGAGAAATACAGATTACCGAGAG II@-ICI.<FII48=00$*(1,2/+&/#'&%#&!#" NM:i:1 XS:A:+
SRR013413.10519592 0 Chr1 3661 255 36M * 0 0 CAAATACATAATCGGAGTAATACAGATTAAAGAGAG IIIIIIIIIIIII5I-04,IEI5I7I&0I-I)H2I1 NM:i:2 XS:A:+
SRR013413.10886083 0 Chr1 3661 255 36M * 0 0 CAAATACATAATCGGAGAAATACAGATTACAGAGAG IIIIII-IIIIIHIIII@<IIII1I&D?I/I?;IFH NM:i:0 XS:A:+
SRR013415.6280320 0 Chr1 3661 255 36M * 0 0 CAAAAACATAATCGGAGAAATACAGATTACAGAGAG ):;5$I-B9II6III0<"I9:0'./*'&5$$4++(, NM:i:1 XS:A:+
SRR013416.394478 0 Chr1 3661 255 36M * 0 0 CAAATACATAATCGGAGCAATACAGATTACAGAGAG I;IIAI4III2%I<I>A549I13;I3%I)*,++&-' NM:i:1 XS:A:+
SRR013423.58140 0 Chr1 3661 255 36M * 0 0 CAAATACATAATCGGAGAAATATAGATTACAGAGAG IIIIIIIIIIIIIIIIIE=II=IIIIIIF4C.5938 NM:i:1 XS:A:+
SRR013419.12428410 0 Chr1 3662 255 36M * 0 0 AAATACATAATCGGAGAAATACAGATTACAGAGAGC [email protected]?(71=,;($ NM:i:0 XS:A:+
SRR013424.9910945 0 Chr1 3662 255 36M * 0 0 AAATACATAATCGGAGAAATACAGATTACAGAGAGC IIIIIIIIIIIIIIII/I5FI-0I&6A3(0726)8% NM:i:0 XS:A:+
SRR013424.9913205 0 Chr1 3662 255 36M * 0 0 AAATACATAATCGGAGAAATACAGATTACAGAGAGC 1IIIIIIIIIIIIIII.I64I-0I&6A3(0726)8) NM:i:0 XS:A:+
SRR013425.7971507 0 Chr1 3662 255 36M * 0 0 AAATACATAATCGGAGAAATACAGATTACAGAGAGC IIIIIIIIIIIIIIIIIIIIIIIIIIIIII53I8I' NM:i:0 XS:A:+
SRR013419.2519606 0 Chr1 3663 255 36M * 0 0 AATACATAATCGGAGAAATACAGATTACAGATAGCG >@9IIIIIII=IH9I2%%>H/C//C?<$2'+%/'"* NM:i:1 XS:A:+
SRR013420.10715881 0 Chr1 3663 255 36M * 0 0 AATACATAATCGGAGAAATACAGATTACAGAGAGTG IIIIIIIIIIIII>I0F<II4FII&II>5;2=+0)I NM:i:1 XS:A:+
SRR013428.5275064 0 Chr1 3663 255 36M * 0 0 AATACATAATCGGAGAAATACAGATTACAGAGAGCG IIIIIIIIIIIIIII<IIIIIII=CI?04>8=863. NM:i:0 XS:A:+
SRR013423.4644761 0 Chr1 3670 255 36M * 0 0 AATCGGAGAAATACAGATTACAGAGAGCGAGAGAGA IIIIIIIIIIIIIIIIIIIIII8BA+=,9$7+.-A' NM:i:0 XS:A:+
SRR013424.3035743 0 Chr1 3670 255 36M * 0 0 AATCGGAGAAATACAGATTACAGAGAGCGAGAGAGA IIIIIIII:IIIIIIIDIA8B455F/6"7*1%-&1* NM:i:0 XS:A:+
SRR013428.4636073 0 Chr1 3683 255 36M * 0 0 CAGATTACAGAGAGCGAGGGAGATCGACGGCGAAGG ;?I)IHI.ID5I*I"3#2!6#:,)",#')($(!!'" NM:i:2 XS:A:+
SRR013429.4979690 0 Chr1 3710 255 36M * 0 0 CGGCGAAGCTCTTTACCCGGAAACCATTGAAATCGG IIIIII?IIIIIIIII>FEI5/8+5.@86)3,;$*) NM:i:0 XS:A:+
SRR013429.5258574 0 Chr1 3730 255 36M * 0 0 AAACCATTGAAATCGGACGGTTTAGTGAAAATGGAG IIIBIIIIIIIIIII6<8I;.183%***2-,'*(%% NM:i:0 XS:A:+
SRR013425.6452022 0 Chr1 3731 255 36M * 0 0 AACCATTGAAATCGGAAGGTTTAGTGAAAATGGAGG IIIIIIIIIIIIIIIC%II7FI:54G91.6+;+%I+ NM:i:1 XS:A:+
SRR013425.8655977 0 Chr1 3732 255 36M * 0 0 ACCATTGAAATCGGACGGTTTAGTGAAAATGGAGGA IIIIIIIIIIIIIIIIII1AI;92<&68?:+*&1,! NM:i:0 XS:A:+
SRR013428.332860 0 Chr1 3735 255 36M * 0 0 ATTGAAATCGGACGGTTTAGTGAAAATGGAGGATCG AIIIIIIIIII>IIIIIIHEH73DI9A/A+HB+/$/ NM:i:1 XS:A:+
SRR013412.15231224 0 Chr1 3737 255 36M * 0 0 TGAAATCGGACGGTTTAGTGAAAATGGAGGGTCAAG IIIIIIIII:IIIIIIIII68-61F7:&>*#%**+( NM:i:1 XS:A:+
SRR013413.1398893 0 Chr1 3737 255 36M * 0 0 TGAAATCGGACGGTTTAGTGAAAATGGAGGAGCAAG IAIIIIIICIAI:$I;:&+9,6>5*%-+'/%"!!$( NM:i:1 XS:A:+
SRR013413.5384349 0 Chr1 3737 255 36M * 0 0 TGAAATCGGACGGTTTAGTGAAAATGGAGGATCAAG IIII7IIIIIHII08II5->6+I*&,..%"!!+" NM:i:0 XS:A:+
SRR013415.1168240 0 Chr1 3737 255 36M * 0 0 TGAAATCGGACGGTTTAGTGAAAATGGAGGGTCAAG IIIIIIIIIIIIIIIII2IEI<I9GD3):5#4&-)+ NM:i:1 XS:A:+
SRR013416.6865013 0 Chr1 3737 255 36M * 0 0 TGAAATCGGACGGTTTAGTGAAAATGGAGGATCAAG IIIIIIIIIIIII2IIE@CD86H9?C6.:9%+"$.2 NM:i:0 XS:A:+
SRR013418.1277542 0 Chr1 3737 255 36M * 0 0 CGAAATCGGACGGTTTAGTGAAAATGGGGGATCAAG &IIB:III=6III-??393++4304/+"+*"+"&"( NM:i:2 XS:A:+
SRR013425.3965135 0 Chr1 3737 255 36M * 0 0 TGAAATCGGACGGTTTAGTGAAAATGGAGGATCAAG IIIEIIIIIIFII6II=7<A?;??:/5-2+$#!!*5 NM:i:0 XS:A:+
SRR013425.3858880 0 Chr1 3739 255 36M * 0 0 AAATCGGACGGTTTAGTGAAAATGGAGGATTAAGTT IIIIIIIIIIIIIIIIIIIIIIIII+II"(#.+<$- NM:i:1 XS:A:+
SRR013428.8141273 0 Chr1 3740 255 36M * 0 0 AATCGGACGGTTTAGTGAAAATGGAGGATAAAGTTG ?@4IIH3BII%II9-%I32'0368#,'&3"*)4.(& NM:i:1 XS:A:+
SRR013424.5843681 0 Chr1 3741 255 36M * 0 0 ATCGGACGGTTTAGTGAAAATGGAGGATCAAGTTGG IIIIIIIII/IIICD=+9E1I3;,CH%$%-,/)4-, NM:i:0 XS:A:+
SRR013419.12061437 0 Chr1 3754 255 36M * 0 0 GTGAAAATGGAGGATCAAGTTGGGTTTGTGTTCCTT I6I0@>9IID*A<&%-1+;./+-0%'++"%#+"&"$ NM:i:2 XS:A:+
SRR013418.7028110 0 Chr1 3765 255 36M * 0 0 GAATCAAGTTGGGTTTGGGTTCCGTCCGAACGACGA ,I+IBB5GFIII:/9C)2&+.&%&)&'+($$,&%+& NM:i:1 XS:A:+
SRR013419.5339155 0 Chr1 3765 255 36M * 0 0 GGATCAAGTTGGGTTTGGGTTCCGTCCGAACGACGA IIEIIIIIIIIIIIIIDI49I+.>/&+6','-)%:" NM:i:0 XS:A:+
SRR013417.1566108 0 Chr1 3799 255 36M * 0 0 GAGGAGCTCGTTGGTCACTATCACCGTAACAAAATA IH>I5B7337*1+5+;6+++',!%5)"6&.8%.$!" NM:i:2 XS:A:+
SRR013419.11051550 0 Chr1 3800 255 36M * 0 0 AGGAGCTCGTTGGTCACTATCTCCGTAACAAAATCG DII>IIFIC:CG<2865I/809//+)+'+)++(%%( NM:i:0 XS:A:+
SRR013424.8602330 0 Chr1 3800 255 36M * 0 0 AGGAGCTCGTTGGTCACTATCTCCGTAACAAAATCG IIIIIIIIIIIIIIIIII=IIIIC;-35:6>)C.)1 NM:i:0 XS:A:+
SRR013411.15124438 0 Chr1 3801 255 36M * 0 0 GGAGCTCGTTGGTCACTATCTCCGTAACAAAATCGA II@IIEIIIIII2!=4I3D88.4,(/5+8/1.'7,+ NM:i:0 XS:A:+
SRR013412.1388062 0 Chr1 3801 255 36M * 0 0 GGAGCTCGTTGGTCCCTATCTCCGTAACAAAATCGA CI/IBIII<=?54E3DC*@(?;,(/&&,%&&*$)"" NM:i:1 XS:A:+