2
$\begingroup$

I run a cuffcompare:

cuffcompare -r transcripts.gtf TransStandard_1.0.gtf

I got the following result:

#= Summary for dataset: /home/tedwon/Sources/QA/data/trans/TransStandard_1.0.gtf :
#     Query mRNAs :     165 in      78 loci  (157 multi-exon transcripts)
#            (69 multi-transcript loci, ~2.1 transcripts per locus)
# Reference mRNAs :     165 in      78 loci  (157 multi-exon)
# Super-loci w/ reference transcripts:       77
#--------------------|   Sn   |  Sp   |  fSn |  fSp  
    Base level:     100.0   100.0     -       - 
    Exon level:     110.6   110.6   100.0   100.0
  Intron level:      99.6    99.6    99.6    99.6
Intron chain level:     100.0   100.0   100.0   100.0
  Transcript level:     100.0   100.0    99.4    99.4
       Locus level:     100.0   100.0   100.0   100.0

 Matching intron chains:     157
          Matching loci:      78

      Missed exons:       0/872 (  0.0%)
       Novel exons:       0/872 (  0.0%)
    Missed introns:       3/756 (  0.4%)
     Novel introns:       0/756 (  0.0%)
       Missed loci:       0/78  (  0.0%)
        Novel loci:       0/78  (  0.0%)

 Total union super-loci across all input datasets: 78 

What I don't understand is that why the sensitivity and specificity for the exon level are over 100? The values in this example are both 110.6. Shouldn't it be from 0% to 100%? The website doesn't address it.

$\endgroup$
  • $\begingroup$ Can you please format the output properly (adjust the tabs so that the lines are aligned) $\endgroup$ – WYSIWYG Sep 8 '15 at 5:07
1
$\begingroup$

I don't know if this is still needed but I went through the source code of cuffcompare to see if this phenomenon is indicative of the false positive rate or an issue with the tool (and to provide a more thorough answer). It seems like it may be a bug in the code for cuffcompare (which is a massive amount of code). If you look at line ~1518 in the source code at this URL:

https://github.com/cole-trapnell-lab/cufflinks/blob/master/src/cuffcompare.cpp

There are several statements that suggest Sn values > 100% would be possible and should be reassigned to 100%. These results may occur from keeping "redundant" transcripts if they start with the same 5' intron. How do your "fuzzy" Sn and Sp compare to this result? are they are similar?

The developers are no longer maintaining this software in favor their new faster version, Hisat/Stringtie. However I find the output to be more complete from Tophat and Cufflinks. To get a better understanding of the output, you could compare the number of potentially novel features in your cuffcmp.combined.gtf to the total number of features in your new cuffcmp.combined.gtf using overlapping and contained features as a proxy for false positive rate. Comparing these values with the total number of lines in your reference annotation may provide you with another Sn estimate without performing qPCR.

# to count the number of potentially novel
grep -c "\"j\"" cuffcmp.combined.gtf 

# to count the overlapping and contained features
grep -c "\"[c,o]\"" cuffcmp.combined.gtf 

# to count the number of exact matches
grep -c "\"=\"" cuffcmp.combined.gtf

# To count how many features are identified for each gene
cut -f4  -d ";" cuffcmp.combined.gtf |   uniq -c

# To count how many known genes gene were identified 
cut -f4  -d ";" cuffcmp.combined.gtf |   uniq -c | wc -l

# To count number of features in the reference gtf
wc -l reference_annotation.gtf

or if you prefer to use Bioconductor packages like rtracklayer in an R environment, you could read the cuffcmp.combined.gtf file and analyze the data further

library(rtracklayer) # if you have the package installed
rtracklayer::readGFF("~/full/path/to/cuffcmp.combined.gtf")
$\endgroup$
  • 2
    $\begingroup$ is that science direct an open access article? if not, some users may not have access. Also, you may want to check out the how to best answer a question in the help center biology.stackexchange.com/help/how-to-answer it would definitely improve your answer if you added some background and links to other users to follow along - thanks for your efforts! $\endgroup$ – Vance L Albaugh Jun 22 '16 at 16:59
  • $\begingroup$ @Drew The Evaluation of Gene Structure Prediction Programs paper does give the definitions but nowhere it define how sensitivity can be over 100%. $\endgroup$ – SmallChess Jun 23 '16 at 0:31
  • $\begingroup$ @Drew I believe this is a bug as sensitivity over 100% is totally unacceptable. This is like probability over 1. $\endgroup$ – SmallChess Jun 23 '16 at 0:31
  • $\begingroup$ I thought that this might be the result of including novel features in the cuffcmp.stats output which are added to the numerator but not to the denominator of the Sn and Sp calculations. However, given that your output lacks novel features and produces greater than 1 Sn and Sp calculations, that assumption is incorrect. It may be a small enough task to modify the source code on your local workstation. $\endgroup$ – Drew J-H Aug 8 '16 at 9:03
  • $\begingroup$ github.com/cole-trapnell-lab/cufflinks/blob/master/src/… $\endgroup$ – Drew J-H Aug 8 '16 at 9:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.