GFF Utilities: GffRead and GffCompare

Implementation

Both our utilities share a code base built around a C++ class called GffObj that implements many of the common GFF parsing and indexing functions. Because the GFF format has no requirements for grouping and sorting of hierarchically linked genomic features (e.g. a transcript feature can have one of its exons at the beginning of the file and another at the end of the file), the parser has to keep transcript data in memory until the whole file is parsed. Feature identifiers (like transcript IDs) are kept in string hashes for fast identification of hierarchical relationship between features. Reference sequence names and GFF attribute names are also stored in global string hashes with numeric IDs associated, while pointers to the genomic feature objects (GffObj) are stored in dynamic arrays sorted by the genomic location such that a binary search can be used for quick overlap verification. The code shared by these utilities also implements functions to test and classify the structural similarities and overlaps between transcripts in the same location on the genome.

GffRead. We initially implemented the GffRead utility as a fast tool for verification, filtering and conversion of the most popular annotation file formats, GTF and GFF3, and for quick extraction of transcript sequences from the genome sequence. With its many features added over time, GffRead is now a complex and versatile tool that can sort, filter, remap and even cluster transcripts into loci (based on exon overlaps) while optionally discarding "redundant" transcripts from an input GFF data. Different examples for the command lines used to perform all these functions are offered in the Use Cases section below.

GffRead parses the input records given in GTF, GFF3 or BED format, and stores them into an internal collection of GffObj data structures that can be easily sorted and filtered according to different criteria. For instance, GffRead can output only the subset of the input transcripts that are multi-exonic, or do not belong to pseudogenes (see Table 1 for a complete set of filtering options). Besides conversions between different GFF formats, GffRead has many additional output options (see Table 2). Among these is a user-defined tab-delimited format, with a line for each transcript and the columns defined by a custom list of some of the GFF columns and attributes in the input annotation file. If a genome sequence is provided, GffRead can also generate multiple additional sequence data files in FASTA format such as: (1) a file with the transcript sequences produced by extracting and concatenating all of the exon sequences of each transcript; (2) a file with all the protein-coding sequences in each transcript; or (3) a file with the amino-acid translations of the coding sequence of each transcript. If a FASTA index file (such as the one created by the samtools utility²⁶) is not present in the same directory with the genomic sequence, GffRead will first create one in order to accelerate the retrieval of the specific transcript sequences. If the transcripts in the annotation file have coding sequences (represented as CDS features in the file), GffRead can check their validity and add specific annotations to the output file, indicating if either the START or the STOP codons are missing in these transcripts or if there are in-frame STOP codons.

The transcript clustering functions of GffRead can group each set of input transcripts into a locus, where all transcripts in a locus are on the same strand, and any two transcripts in that locus have at least one exonic interval overlap. When clustering is enabled, the GFF output will have a new 'locus' feature for each cluster with attributes listing all the transcript IDs (and gene IDs, if available) that belong to that cluster. Optionally, GffRead can identify transcripts that are structurally "matching" or "equivalent", defined as transcripts that share all their introns, or have more than 80% of their length overlap in the case of single exon transcripts. GffRead can also discard redundant transcripts (either matching or contained within other transcripts) from the output, providing the user with the ability to choose among merging strategies with different levels of stringency when assessing redundancy in such cases.

GffCompare. GffCompare is a generic, standalone tool for merging and tracking transcript structures across multiple samples and comparing them to a reference annotation. Initially written based on the CuffCompare utility program included with the Cufflinks suite²⁷, GffCompare has the following main functions:

The last two purposes require the user to provide a known reference annotation file that GffCompare then uses to classify all the transcripts in the input samples according to the reference transcript that they most closely overlap (Figure 1). To assess the accuracy of transcriptome assemblies, GffCompare reports several accuracy metrics previously employed for gene prediction evaluation²⁸. These metrics include sensitivity and precision as well as the number of novel or missed features, and the metrics are computed at various levels (base, exon, intron chain, transcript, or locus). More details about how to obtain the different reports provided by GffCompare can be found in the Use Cases section.

Figure 1. Transcript classification codes based on their relationship to reference transcripts, as generated by GffCompare.

Reference exons and transcripts are shown in black, transcripts to be classified are shown in blue, and hashed regions represent repeated regions in the genome. For example, the transcript in blue on the uppermost left panel is labeled “=” because all of its introns precisely match the annotation in black.

Some pipelines can produce a very large number of transcripts that need to be evaluated; e.g. when merging the transcript assemblies from tens or hundreds of RNA-seq experiments. Because GffCompare always loads the entire transcript data into memory for clustering, running GffCompare on such large GTF/GFF files could be slow and memory intensive. One may be interested only in how these transcripts overlap the reference annotation, and then only wish to further analyze those transcripts that have specific types of overlaps with the reference annotation transcripts. GffCompare also only produces the best match of a transcript to a reference annotation, but for each transcript we might want to know all possible reference matches. In order to address these needs, we built TrMap ("Transcript vs. reference Mapping"), a program that we distribute along with GffCompare and that was designed to avoid using a large amount of memory by streaming the input transcript data. TrMap first loads the reference annotation into an interval tree data structure²⁹, and then for each query transcript it reports all the reference transcripts that overlap it, along with their overlap classification codes. These are the same classification codes described in Figure 1, with the exception of codes p, r, and u which are reserved for transcripts that do not overlap reference transcripts and represent transcripts that are single exon and nearby genes (p), repeats outside of genes (r), and intergenic (u).

Operation

This software can be built on a Linux or MacOS system with no other library dependencies. A GNU C++ compiler (g++) is required for compilation (on Linux at least g++ version 4.5 is required). The release packages on Github include precompiled binaries for Linux and MacOS that can be used directly instead of having to build the programs from source. Linux compatibility goes back as far as RedHat Enterprise Linux 5, while on MacOS the programs can run on systems as old as OS X 10.7 (Lion). We also provide the gffread, gffcompare and trmap executables. These are supposed to be used as command line programs, in a Linux/Unix shell, in a terminal or a script. All programs take GFF3, GTF or BED files as their (main) input files. Both packages require the shared code provided in GCLib (https://github.com/gpertea/gclib³⁰).

Basic usage examples of the GffRead utility

The program GffRead can be used to validate, filter, convert and perform various other operations on GFF files (see Table 1 and Table 2 for the full list of usage options). For instance, GffRead can be used to simply read an annotation file in a GFF format, and print it in either GFF3 (default) or GTF2 format (with the -T option), while optionally discarding any non-essential attributes, and fixing some potential issues with the input file. The command line for such a quick cleanup and a quick visual inspection of a given GFF file would be:

gffread -E annotation.gff -o ann_simple.gff

This will show the minimalist GFF3 re-formatting of the transcript records found in the input file (annotation.gff in this example) which could be given in either GFF3 or GTF2 format. The -E option directs GffRead to "expose" (display warnings about) any potential issues encountered while parsing the input file.

In order to obtain the GTF2 version of the same transcripts, the -T option should be added:

gffread annotation.gff -T -o annotation.gtf

GffRead can be used to generate a FASTA file with the DNA sequences for all transcripts in a GFF file. For this operation a fasta file with the genomic sequences have to be provided as well. This can be accomplished with a command line like this:

gffread -w transcripts.fa -g genome.fa annotation.gff

The file genome.fa in this example would be a multi-fasta file with the chromosome/contig sequences of the target genome. This also requires that every contig or chromosome name found in the 1st column of the input GFF file (annotation.gtf in this example) must have a corresponding sequence entry in the genome.fa file.

Basic usage example of the GffCompare utility

The program GffCompare can be used to compare, merge, annotate and estimate accuracy of one or more GFF files (the “query” files), when compared with a reference annotation (also provided as GFF). A basic command line to compare a list of GTF files to a reference annotation file is:

gffcompare -r annotation.gff transcripts.gtf

The reference annotation is specified in the annotation.gff file and transcripts.gtf represents the query file (more than one query file can be provided). Unless the -o option was provided, the output will be found in multiple files with the prefix “gffcmp.”. A list of the more important options for the GffCompare utility is provided in Table 3.

Transcript accuracy estimation with GffCompare

GffCompare can be used to assess the accuracy of transcriptome assemblies produced by programs like StringTie¹⁹ in respect to a known reference annotation. To this end, GffCompare reports various statistics related to the accuracy of the input transcripts compared to the reference annotation in the <outprefix>.stats file. Among these statistics are sensitivity and precision values computed at various levels (base, exon, intron chain, transcript, locus), which are calculated as:

where TP stands for "true positives", or query features (bases, exons, introns, transcripts, etc.) that agree with the corresponding reference annotation features; FN means "false negatives", i.e. features that are found in the reference annotation but are not present in the input data; FP (“false positives”) are features present in the input data but not confirmed by any reference annotation data. Notice that FP+TP amounts to the whole input set of query features in the input file. If multiple query GTF/GFF files are given as input, these metrics are computed separately for each sample.

Sensitivity and Precision values are estimated at various levels, which are largely an increasingly stringent way of evaluating the accuracy/correctness of a set of predicted transcripts (transfrags), when compared to the reference annotation provided with the -r option. The six different levels that GffCompare uses are described below:

1) Base level. At the base level, TP represents the number of exon bases that are reported at the same coordinate on both the query transcripts and any reference transcript, FN is the number of bases in reference data exons that are not covered at all by any of the query exons, and FP is the number of bases which are covered by predicted transcripts' exons but not covered by any reference transcript exons.

2) Exon level. We define the TP, FN, and FP values at the exon level similar to the base level, but now the unit of comparison is the exon interval on the genome, i.e. if an exon of the predicted transcript overlaps and matches the boundaries of a reference transcript exon, then it is counted as a TP.

3) Intron Level. Intron intervals are the units that are matched at the intron level, therefore each intron of the predicted transcript is checked against any introns of the reference transcripts in the same region and if there is one with the same exact start-end coordinates, it is counted as a TP.

4) Intron chain level. At this level we count as a TP any query transcript for which all of its introns can be found, with the same exact intron coordinates as in a reference transcript that has the same number of introns. Matching all the introns at this level implies that all the internal exons also match, but this might not be true for the external boundaries of the terminal exons.

5) Transcript level. Note that intron chain level values are calculated only by looking at multi-exon transcripts, so it completely ignores the single-exon transcripts, which can be quite numerous in a RNA-Seq experiment (possibly due to a lot of transcriptional and alignment noise). The transcript level considers single-exons as well. A TP at this level is defined as a full exon chain match between the predicted transcript and a reference transcript, where all internal exons match and the outer boundaries of the terminal query exons can only slightly differ from the reference exons (with at most 100 bases by default). Also GffCompare considers single-exon transcripts as matching an overlapping single-exon reference transcript if there is a significant overlap between the two (more than 80% of the longer transcript by default).

6) Locus level. At this level GffCompare considers that an observed locus, defined as a cluster of exon-overlapping transcripts, matches a similarly built reference locus if at least one predicted transcript has a transcript level match with a reference transcript in the corresponding reference locus.

Other statistics reported by GffCompare are the number of missed or novel exons, missed or novel introns and missed or novel loci. Note that in order to properly evaluate precision and sensitivity when comparing two sets of transcripts, special care must be taken for duplicated (or redundant) entries within each set. GffCompare uses different levels of stringency of what to consider duplicated transcripts, depending on the option given in its input (see options -D, -S, -C, -A, -X in Table 3).

Merging structurally equivalent transcripts with GffCompare

When multiple input GTF/GFF files are provided, GffCompare reports a GTF file named <outprefix>.combined.gtf containing the union of all transfrags in each sample. If a transfrag with the same exact intron chain is present in both samples, it is thus reported only once in the output file.

The "super-locus" concept

A super-locus is a region of the genome where predicted transcripts and reference transcripts get clustered together by exon overlaps. When multiple GFF files with are provided as input to GffCompare, this clustering is performed across all the input files. Due to the transitive nature of this clustering, these super-loci can occasionally get very large, sometimes merging a few distinct reference gene regions together, especially if there is a lot of transcription or alignment noise around the individual gene regions. For each super-locus, GffCompare assigns a unique identifier with the XLOC_ prefix.

Annotating transcripts with GffCompare

One can run GffCompare on a single GTF/GFF input file using with the -r option (which provides a reference annotation), and without any specific options to remove redundant transfrags (such as the -D, -S, -C, -A, -X options) to produce an GTF file called <outprefix>.annotated.gtf that contains all the input transcripts annotated with several additional attributes: xloc, tss_id, cmp_ref, and class_code. The xloc attribute specifies the super-locus a specific transcript belongs to. The tss_id attribute uniquely identifies the transcription start for that transcipt, and using this value the user can quickly see which transcripts use the same transcription start, or how many different transcription starts are present in a locus. The cmp_ref gives the closest reference transcript (where applicable), while the relationship to this reference transcripts is given by the class_code attribute. The possible values for the class_code attribute are listed in Table 4.

Tracking transcripts with GffCompare

GffCompare can also be used to track all transcripts that are structurally equivalent among the different input files. GffCompare considers transcripts matching (or structurally equivalent) if all their introns are identical. Note that matching transcripts are allowed to differ on the length of the first and last exons, since these lengths can usually vary across samples for the same biological transcript. A list of all matching transcripts is reported in a file called <outprefix>.tracking in which each row represents a transcript. The first column in this file represents a unique id assigned to that transcripts. The second file represents the super-locus that contains that transcript. If GffCompare was run with the -r option, the 3^rd and 4^th columns contain the reference annotation transcript that was found to be closest to the transcript and the classification code (as specified by Table 4) that specifies the relationship between these two transcripts, respectively. The rest of the columns show the corresponding transcript from each input file in order. An example and a brief description for each column are given in Table 5.

In order to quickly see which reference transcripts match which transcripts from a sample file, two other files, called <outprefix>.<input_file>.refmap and <outprefix>.<input_file>.tmap are also created for each query <input_file>. The <outprefix>.<input_file>.refmap file is a tab-delimited file that has a row for each reference transcript that either fully or partially matches a transcript from the given input file. Its columns are described in Table 6. Conversely, the <outprefix>.<input_file>.tmap file has a row for each input transcript, while the columns in this file (as detailed in Table 7) describe the most closely matching reference transcript for that transcript.

Overlap classification for a large set of transcripts with TrMap

The utility TrMap was designed for large scale overlap analysis of streaming transcript prediction data (millions of transcripts) with a reference annotation data set. Particularly, TrMap performs detection and classification of all the overlaps found between the streamed transcripts and the reference annotation transcripts.

The program trmap is distributed with GffCompare and a basic usage for it is shown below:

trmap [-S] [-o ] <ref_gff> <query_gff>
Positional arguments:
  <ref_gff>    reference annotation file name (GFF/BED format)
  <query_gff>  query file name (GFF/BED format) or "-" for stdin
Options:
  -o <outfile> write output to <outfile> instead of stdout
  -S           report only simple reference overlap percentages, without
               classification (one line per query)

The default output is a pseudo-FASTA format showing a record for each query transcript that had at least one reference overlap. The query transcript is shown in the header of the record, with space delimited fields showing the genomic location and strand. Each reference overlap follows, as a line with tab delimited fields, starting with the "classification code" for the overlap and then providing the genomic location of the transcript (chromosome, strand, transcript-start, transcript-end, reference_transcriptID, exons).

The exons for both query and reference transcripts are shown as comma delimited lists of intervals. These are all 1-based coordinates like in the GTF/GFF format (even when input is BED).

Underlying data

All data underlying the results are available as part of the article and no additional source data are required.