How to find novel transcripts using GFFcompare?

Question

I am trying to find novel transcripts from an RNA-seq database. Based on the advice I got, it seemed that using Stringtie for transcript assembly is a good way to go, and it supports novel transcript discovery even with the reference GTF file provided (which apparently significantly improves the assembly process). So I tried to follow the protocol provided in the Nature Protocols paper that described the use of Stringtie; the paper also suggested using GFFcompare for comparing the assembled transcripts with the reference, to quantify and find the novel transcripts.

I followed the entire protocol, and this was the output of GFFcompare when comparing the reference GTF with that of one of the samples after mapping -

#= Summary for dataset: ./hisat2/ERR188044_chrX.gtf 
#-----------------| Sensitivity | Precision  |
        Base level:    51.7     |    79.5    |
        Exon level:    46.7     |    85.2    |
      Intron level:    47.2     |    93.9    |
Intron chain level:    31.2     |    64.4    |
  Transcript level:    31.6     |    52.3    |
       Locus level:    36.6     |    50.1    |

     Matching intron chains:     580
       Matching transcripts:     664
              Matching loci:     397

          Missed exons:    4395/8804    ( 49.9%)
           Novel exons:     426/4874    (  8.7%)
        Missed introns:    3832/7946    ( 48.2%)
         Novel introns:      83/3992    (  2.1%)
           Missed loci:     610/1086    ( 56.2%)
            Novel loci:     273/795 ( 34.3%)

 Total union super-loci across all input datasets: 749 
1270 out of 1270 consensus transcripts written in 188044_only.annotated.gtf (0 discarded as redundant)

In the paper however, the number of novel genes and transcripts is clearly mentioned -

So, I was wondering how they were able to calculate the number of novel transcripts and genes, from the data that GFFcompare provides (they even mention in the paper that they got it from GFFcompare). I see the numbers related to novel exons and novel introns, but how do I translate it to the number and identities of the novel transcripts/genes?

Cross-posted in Biostars. Let me know if I shouldn't do this.

Answer 1

The output of gffcompare includes several files per run (just like cuffcompare). Example for a run:

$ ls cuffcmp.* | sed 's/\t/\n/'
cuffcmp.combined.gtf
cuffcmp.loci
cuffcmp.output.gtf.refmap
cuffcmp.output.gtf.tmap
cuffcmp.stats
cuffcmp.tracking

From my experience, the easiest file to manipulate is the tmap file:

Tab delimited file lists the most closely matching reference transcript for each query transcript. There is one row per query transcript, and the columns are as follows:

1     Reference gene name    The gene_name attribute of the reference GTF record for this transcript, if present. Otherwise gene_id is used.
2     Reference transcript id     The transcript_id attribute of the reference GTF record for this transcript
3     Class code  The type of relationship between the query transcripts in column 4 and the reference transcript (as described in the Class Codes section below)
4     Query gene id   The query (e.g., Stringtie) internal gene id
5     Query transcript id     The query internal transcript id
6     Number of exons     The number of exons in the query transcript
7     FPKM    The expression of this transcript expressed in FPKM
8     TPM     the estimated TPM for the transcript, if found in the query input file
9     Coverage    The estimated average depth of read coverage across the transcript.
10    Length  The length of the transcript
11    Major isoform ID    The query ID of the gene's major isoform
12    Reference match length  The length of the longest overlap with a reference, ‘-’ if there is no such exonic overlap

The important column to look at for transcript classification is actually the 3rd, "class code". The content is meaningful if you used option -r when running gffcompare.

In the tmap file you can access info on both transcript and gene ID directly, both from the assembler and the reference annotation.

Now what you need is:

1. n. assembled genes
2. n. novel genes
3. n. novel transcripts
4. n. transcripts matching annotation

At this point, it is essential to understand what the authors mean by "novel" (gene/transcript), does it mean "non previously annotated" (= class code u, for "unknown")? Also, what about "matching", does it mean "perfectly matching" (= class code =, for "complete, exact match of intron chain")?

Under these assumptions, you can do some easy stats with some bash help.

1. n. assembled genes:

   $cat cuffcmp.output.gtf.tmap | sed "1d" | cut -f4 | sort | uniq | wc -l

2. n. novel genes:

   $ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f4 | sort | uniq | wc -l

3. n. novel transcripts:

   $ cat cuffcmp.output.gtf.tmap | awk '$3=="u"{print $0}' | cut -f5 | sort | uniq | wc -l

4. n. transcripts matching annotation:

   $ cat cuffcmp.output.gtf.tmap | awk '$3=="="{print $0}' | cut -f5 | sort | uniq | wc -l

I hope this helps. I think it's an easier approach than trying to "translate the numbers related to novel exons and novel introns to the number and identities of the novel transcripts/genes" from the stats file.

Although you have to think a bit what "class codes" to include in each count. Maybe my assumptions are too strict for your use case.