How to find novel transcripts using GFFcompare?

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.

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.

• Thank you so much for this detailed response, this was life-saving. I had been trying to find a way out for so long on multiple forums! Turns out I just need to look better at the files that the tool had already generated for me. Whats interesting is that using the code you describe, I get about ~80K assembled genes and ~220K transcripts matching annotation! I checked why this could be and found that most entries in the column for assembled gene (query gene ID) are repeated, while there are no repeats in the column for the query ID (used for transcripts matching annotation). I wonder why! Jul 22 '19 at 19:16