Transcripts missing from cufflinks

[dnusol]

Hi, I am moving this thread from the RNA Sequencing group, in case anyone here has some ideas. From my search, this is an issue that comes up every now and then with very few comments. I emailed tophat list a while ago but no reply so far.

I am trying to figure out why the following transcript, and some others similarly, is excluded from my comparison analysis. I have an assembled reference fasta file with some 48000 entries that I have indexed with bowtie. I have aligned the reads from both my samples using Tophat and want to find differentially expressed transcripts between both samples (regardless of the power issue).
After running Tophat I get that this transcript I am interested in has been included both in the deletions.bed and the insertions.bed. See below:

#for sample 1
$ grep 'Locus_28169_Transcript_10/10' .
deletions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 817 818 - 1
deletions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 1032 1033 - 1
insertions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 988 988 A 3

#for sample 2
$ grep 'Locus_28169_Transcript_10/10' .
deletions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 140 141 - 1
deletions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 771 772 - 2
insertions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 258 258 T 1
insertions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 453 453 A 5
insertions.bed:Locus_28169_Transcript_10/10_Confidence_0.803_Length_1176 988 988 A 4

I have not found much information as to what columns 2 and 3 indicate but I suspect that they are coordinates for deletions and insertions found after alignment. But what does the last column mean? is it the number of reads supporting the INDEL?

This is the command I run for Tophat for one sample, same for the other one:
$ tophat -r 380 -p 8 --mate-std-dev 18 --solexa1.3-quals --library-type fr-unstranded --min-coverage-intron 2 --min-segment-intron 10 -o ./tophat_s1_blast10-3 /home/Analysis/bowtie-0.12.8/indexes/Match_10e-3 s_1_trim1.fastq s_1_trim2.fastq --max-insertion-length 7 --max-deletion-length 5

After this, I run cufflinks and now I find that this transcript is not included in any further analysis or comparisons, it is missing from the cufflinks output files (genes.fpkm_tracking, isoforms.fpkm_tracking, skipped.gtf and transcripts.gtf), but I am biologically interested in this transcript.

Here my cufflinks command
$ cufflinks -o tophat_s1_blast10-3/cufflinks_s1_blast10-3 tophat_s1_blast10-3/accepted_hits.bam -p 8 --library-type fr-unstranded


Why is this transcript being excluded? How can I rescue it so that I can finally compare it between both samples, knowing that both samples have reads aligning to it? I would like to know the FPKM calculated for both samples and then get a p-value, even if it is not significant (by relaxing the significance flag in cuffdiff).


I am using tophat-2.0.4 and cufflinks-2.1.1 versions

Thanks for any comments

[sdriscoll]

The pipeline you're running was designed and intended for when you've got a complete genome reference and you either 1: need to produce a novel annotation of genes within that genome or 2: want to quantify gene expression for genes that are already annotated. These annotations come in GTF format (or GFF). The reason the output is confusing to you is because cufflinks has produced its own gene annotation relative to your assembled transcripts (that's the 'transcripts.gtf' file) and quantified expression in terms of its own annotation system.

This pipeline isn't appropriate for when you're working with assembled transcripts as a reference. In this case you're aligning reads directly to transcripts which requires a different approach since transcripts tend to share a lot of sequence with each other (if your species has alternative splicing). What that means is the aligner will not necessarily select "correct" alignments between transcripts that share sequence.

There are workflows designed for what you're doing, however. I recommend having a go with RSEM which is specifically designed to work with a assembly type workflow (like output from Trinity) where you've got assembled gene/transcript references and simply need to quantify expression of those for downstream differential expression analysis. RSEM has a pretty good tutorial. You still can't go back to using cufflinks or cuffdiff so you'll have to rely on one of the several R based differential expression software such as DESeq or edgeR. RSEM recommends EBSeq as well which claims to model variance properly for transcript differential expression.