GENCODE: The reference human genome annotation for The ENCODE Project

Abstract

Launched in September 2003, the Encyclopedia of DNA Elements (The ENCODE Project Consortium 2011) project brought together an international group of scientists tasked with identifying all functional elements in the human genome sequence. Initially focusing on 1% of the genome (The ENCODE Project Consortium 2007), the pilot project was expanded to the whole genome in 2007. As part of the initiative, the GENCODE collaboration was established whose aim was to annotate all evidence-based gene features on the human genome at high accuracy, again initially focusing on the 1% (Harrow et al. 2006). The process to create this gene annotation involves manual curation, different computational analysis, and targeted experimental approaches. Eight groups in Europe and the United States directly contribute data to this project, with numerous additional sources of evidence also used for the annotation. Figure 1 shows how the different elements of the GENCODE Consortium interact together.

Open in a separate window

Figure 1.

The GENCODE pipeline. This schematic diagram shows the flow of data between the groups of the GENCODE Consortium. Manual annotation is central to the process but relies on specialized prediction pipelines to provide hints to first-pass annotation and quality control (QC) for completed annotation. Automated annotation supplements manual annotation, the two being merged to produce the GENCODE data set and also to apply QC to the completed annotation. A subset of annotated gene models is subject to experimental validation. The Annotrack tracking system contains data from all groups and is used to highlight differences, coordinate QC, and track outcomes.

The ability to sequence genomes has far exceeded the techniques for deciphering the information they encode. Selecting the correct reference gene annotation for a particular project is extremely important for any downstream analysis such as conservation, variation, and assessing functionality of a sequence. The type of gene annotation applied to a particular genome is dependent on its quality; therefore, next-generation sequencing assemblies (Metzker 2010) have had automatic gene annotation applied to them, whereas high-quality finished genomes such as the human (International Human Genome Sequencing Consortium 2004), mouse (Church et al. 2009), and zebrafish (Becker and Rinkwitz 2011) have manual annotation projects associated with them. Publicly available gene sets such as RefSeq (Pruitt et al. 2012), AceView (Thierry-Mieg and Thierry-Mieg 2006), and GENCODE are generated by a combination of manual and automatic annotation and have developed different methods to optimize their annotation criteria. For example, RefSeq annotates cDNAs rather than genomic sequence to optimize full-length gene annotation and is thus able to ignore sequencing errors in the genome. This publication will describe the generation of the GENCODE gene set and its strengths over other publicly available human reference annotation and the reasons it has been adopted by the ENCODE Consortium (The ENCODE Project Consortium 2011), The 1000 Genomes Project Consortium (2010), and The International Cancer Genome Consortium (2010) as their reference gene annotation.

The GENCODE reference gene set is a combination of manual gene annotation from the Human and Vertebrate Analysis and Annotation (HAVANA) group (http://www.sanger.ac.uk/research/projects/vertebrategenome/havana/) and automatic gene annotation from Ensembl (Flicek et al. 2011). It is updated with every Ensembl release (approximately every 3 mo). Since manual annotation of the whole human genome is estimated to take until the end of 2012, the GENCODE releases are a combination of manual annotation from HAVANA and automatic annotation from Ensembl to ensure whole-genome coverage.

Protein-coding genes were annotated automatically using the Ensembl gene annotation pipeline (Flicek et al. 2012). Protein sequences from UniProt (Apweiler et al. 2012) (only “protein existence” levels 1 and 2) were included as input, along with RefSeq sequences. Untranslated regions (UTRs) were added using cDNA sequences from the EMBL Nucleotide Archive (ENA) (Cochrane et al. 2011). Long intergenic noncoding RNA (lincRNA) genes were annotated using a combination of cDNA sequences and regulatory data from the Ensembl project. Short noncoding RNAs were annotated using the Ensembl ncRNA pipelines, using data from mirBase (Griffiths-Jones 2010) and Rfam (Gardner et al. 2011) as input.

This process of combining the HAVANA and Ensembl annotation is complex. During the merge process, all HAVANA and Ensembl transcript models are compared, first by clustering together transcripts on the same strand which have any overlapping coding exons, and then by pairwise comparisons of each exon in a cluster of transcripts. The merge process is summarized in the Supplemental Figures and Tables, including the rules involved in each step. Ensembl have developed a new module, HavanaAdder, to produce this GENCODE merged gene set. Prior to running the HavanaAdder code, the HAVANA gene models are passed through the Ensembl health-checking system, which aims to identify any inconsistencies within the manually annotated gene set. Annotation highlighted by this system is passed back to HAVANA for further inspection. In addition, the HAVANA transcript models are queried against external data sets such as the consensus coding sequence (CCDS) (Pruitt et al. 2009) gene set and Ensembl's cDNA alignments of all human cDNAs. If annotation described in these external data sets is missing from the manual set, then this is stored in the AnnoTrack system (see below) (Kokocinski et al. 2010) so that a record is kept for the annotators to inspect these loci.

The genes in the GENCODE reference gene set are classified into three levels according to their type of annotation. Level 1 highlights transcripts that have been manually annotated and experimentally validated by RT-PCR-seq (Howald et al. 2012), as well as pseudogenes that have been validated by three-way consensus, namely, that have been independently validated by three different strategies. Level 2 indicates transcripts that have been manually annotated. Some Level 2 transcripts have been merged with models produced by the Ensembl automatic pipeline, while other Level 2 transcripts are annotated by HAVANA only. Level 3 indicates transcripts and pseudogene predictions arising from Ensembl's automated annotation pipeline. GENCODE 7 consists of 9019 transcripts at Level 1, 118,657 transcripts at Level 2, and 33,699 transcripts at Level 3. Many of the protein-coding genes in Level 3 are contributed by Ensembl's genome-wide annotation in regions where HAVANA has not yet provided manual annotation.

Manually annotated GENCODE gene features are subdivided into categories on the basis of their functional potential and the source of the evidence supporting their annotation. Annotated gene models are predominantly supported by transcriptional and/or protein evidence. Once the structure of a model has been established, it is classified into one of three broad locus level biotypes: protein-coding gene, long noncoding RNA (lncRNA) gene, or pseudogene. In addition, more detailed biotypes are associated with transcripts to attempt to assign a functionality, for example, protein-coding or subject to nonsense mediated decay (NMD) (see landscape Supplemental Tables).

To provide a more complete description of the gene model, a “status” is assigned at both the locus and transcript level. Loci can be assigned the status “known,” “novel,” or “putative” depending on their presence in other major databases and the evidence used to build their component transcripts. In brief, loci have the status “known” if they are represented in the HUGO Gene Nomenclature Committee (HGNC) database (Seal et al. 2011) and RefSeq (Pruitt et al. 2012); loci with the status “novel” are not currently represented in those databases but are well supported by either locus-specific transcript evidence or evidence from a paralogous or orthologous locus. Finally loci with status “putative” are supported by shorter, more sparse transcript evidence. A similar status categorization is employed at the transcript level (see Supplemental Figures and Tables).

In addition to the information captured by biotype and status, controlled vocabulary attributes are attached to both transcripts and loci. They are used to describe other features relevant to the structure or functional annotation of a transcript. Attributes may be subdivided into three main categories: those that explain features related to splicing, those related to the translation of the transcript, and those related to the transcriptional evidence used to build the transcript model. For a comprehensive list of all attributes along with the definitions used in the GENCODE annotation, see the landscape Supplemental Tables. Where further explanation of annotation is required, free text remarks are added. New controlled vocabulary is developed wherever possible so that annotation text strings can be searched computationally.

Over the last decade, evidence from numerous high-throughput array experiments has indicated that evolution of the developmental processes regulating complex organisms can be attributed to the noncoding regions and not only to the protein-coding regions of the genome (Bertone et al. 2004; Mattick 2004; Kapranov et al. 2007; Clark et al. 2011). The GENCODE gene set has always attempted to catalog this noncoding transcription utilizing a combination of computational analysis, human and mammalian cDNAs/ESTs alignments, and extensive manual curation to validate their noncoding potential. GENCODE 7 contains 9640 lncRNA loci, representing 15,512 transcripts, which is the largest manually curated catalog of human lncRNAs currently publicly available. All the lncRNA loci in the catalog originate from the manual annotation pipeline and are initially classified as noncoding due to the lack of homology with any protein, no reasonable-sized open reading frame (ORF; not subject to NMD), and no high conservation, confirmed by PhyloCSF (see later section), through the majority of exons. The transcripts are not required to be polyadenylated but 16.8% are, and chromatin marks have been identified for 13.9% (Derrien et al. 2012). These lncRNAs can be further reclassified into the following locus biotypes based on their location with respect to protein-coding genes:

1. Antisense RNAs: Locus that has at least one transcript that intersects any exon of a protein-coding locus on the opposite strand, or published evidence of antisense regulation of a coding gene.
2. LincRNA: Locus is intergenic noncoding RNA.
3. Sense overlapping: Locus contains a coding gene within an intron on the same strand.
4. Sense intronic: Locus resides within intron of a coding gene but does not intersect any exons on the same strand.
5. Processed transcript: Locus where none of its transcripts contain an ORF and cannot be placed in any of the other categories because of complexity in their structure.

The GENCODE lncRNA data set is larger than other available lncRNA data sets, and it shows limited intersection with them. Forty-two percent (44 out of 96) of the lncRNA database lncRNAdb (Amaral et al. 2011) are represented in GENCODE lncRNAs. We checked the same strand overlap against recent lncRNA catalogs: GENCODE v7 lncRNAs contain 30% of Jia et al. (2010) lncRNAs, 39% of Cabili et al. lincRNAs (Cabili et al. 2011), and 12% of vlincs (Kapranov et al. 2007) (for more details, see Derrien et al. 2012). While this level of overlap between data sets shows how lncRNA annotation is improving, it also shows that substantial additional work is still required. There are likely to be a number of reasons for the limited overlap between the published lincRNAs and GENCODE, not least that a substantial fraction of transcript annotations are currently incomplete (see below). Another reason is that some of the published transcripts are single exons, which up to now have not been annotated in GENCODE unless there is additional support, for example, polyA features, conservation, submitted sequence, or publications. We are addressing this weakness and re-examining single exons lincRNAs based on annotation from Jia et al. (2010) in collaboration with the Lipovich group, and the data will be incorporated into GENCODE 10.

Although the current definition of lncRNAs requires the transcript to be >200 bp (Wang and Chang 2011), the GENCODE ncRNA set contains 136 spliced transcripts <200 bp (all of them single transcript loci) to highlight that there is evidence of expression at that position in the genome. We currently group the transcripts into loci, which is different compared with other lncRNA analysis groups, for example, the Fantom Consortium (Katayama et al. 2005). Multiple lncRNA transcripts appear to start from the same transcription start site (TSS), for example, the DLX6-AS1 locus shown in Supplemental Figure 2. To estimate the completeness of the lncRNA transcripts, we took advantage of CAGE tags from 12 different cell lines and manually annotated polyA features to assess the TSS and 3′ end of transcripts (Djebali et al. 2012). The beginning and end of 15% and 16.8% of lncRNA are supported, respectively, indicating that the majority of transcripts are incomplete. Interestingly lncRNA transcripts have an unusual exon structure compared with protein-coding transcripts, with their distribution peaking at two and five exons, respectively (see Fig. 2). This lower number does not appear to be an artifact or the product of incomplete annotation but most probably is a bona fide characteristic of the lncRNAs, as it is also observed in potential lncRNA models identified using the Illumina Human Body Map 2.0 (HBM) RNA-seq data generated in 2010 on HiSeq 2000 instruments (described below), which are not built from partial ESTs.

Open in a separate window

Figure 2.

Analysis of exon number of protein-coding and noncoding RNA transcripts. The numbers of exons for each individual transcript annotated at protein-coding and lncRNA loci are plotted for GENCODE 3c (red lines) and GENCODE 7 (blue lines). For each release, darker lines indicate protein-coding transcripts, and lighter lines indicate lncRNA transcripts. The 5′ and 3′ UTR exons of protein-coding transcripts are included.

Small noncoding RNAs are automatically annotated from the Ensembl pipeline and included within the GENCODE gene set. The number has remained relatively stable at 8801 since release 4. Protein-coding and noncoding transcripts that contain a small ncRNA within at least one intron or exon will be annotated with the attribute ncrna_host. Thirty-three percent of small ncRNAs map within the boundaries of a GENCODE gene, the majority of which reside in introns. The GENCODE 7 release contains 1679 protein coding and 301 lncRNA genes with ncrna_host attributes, and there is a sixfold enrichment of small nuclear RNAs (snoRNAs) within exons of lncRNAs (Derrien et al. 2012).

In summary, the lncRNAs data set in GENCODE 7 consists of 5058 lincRNA loci, 3214 antisense loci, 378 sense intronic loci, and 930 processed transcripts loci. Manually evaluating the RNA-seq models generated from HBM data and ENCODE data could potentially double this number in later releases of GENCODE and produce a uniform data set.

Within most gene catalogs, pseudogenes have been annotated as a byproduct of protein-coding gene annotation, since a transcript has been identified with a frameshift or deletion, rather than an important entity in its own right. However recent analysis of retrotransposed pseudogenes such at PTENP1 (Poliseno et al. 2010) and DHFRL1 (McEntee et al. 2011) have found some retransposed pseudogenes to be expressed and functional and to have major impacts on human biology. The GENCODE catalog is unique in its annotation of the comprehensive pseudogene landscape of the human genome using a combination of automated, manual, and, more recently, experimental methods.

The assignment process for pseudogenes is described in detail by Pei et al. (2012). Briefly, in silico identification of pseudogenes is obtained from routine implementation of Yale's Pseudopipe (Zhang et al. 2006) with every new release of Ensembl. Pseudopipe identified 18,046 pseudogenes based on the human genome release in Ensembl 61. These pseudogenes were compared to a recent run of UCSC's RetroFinder, which included 13,644 pseudogenes, and HAVANA's latest annotations of 11,224 pseudogenes based on GENCODE 7, level 2. A three-way Yale, UCSC, and HAVANA pseudogene consensus set was obtained by using an overlap criteria of 50 bp and was developed for the annotation of 1% ENCODE Regions (Zheng et al. 2007). This resulted in a consensus set of 7183 pseudogenes, which are tagged level 1. The functional paralog of a pseudogene is often referred to as the “parent” gene. Currently, we have successfully identified parents for 9369 pseudogenes of the manually annotated pseudogenes, whereas the parents for the remaining 1847 pseudogenes are still ambiguous and may require further investigation. It is important to note, however, that it is not always possible to identify the true parent of a pseudogene with certainty, for example, when a pseudogene is highly degraded and is derived from a parent gene with highly similar paralogs or when the parent contains a common functional domain. We have added this information to the pseudogene annotation if known (based on protein alignments), and it is also available from the pseuodgene decorated resource (psiDR described in Pei et al. 2012) http://www.pseudogene.org/psidr/psiDR.v0.txt.

A pseudogene ontology was created to associate a variety of biological properties—such as sequence features, evolution, and potential biological functions—to pseudogenes and is incorporated into the GENCODE annotation file. The hierarchy of these properties is shown in Supplemental Figure 3. The ontology allows not only comprehensive annotation of pseudogenes but also automatic queries against the pseudogene knowledge database (Holford et al. 2010). The breakdown of the different biotypes within the GENCODE data set can be seen in Supplemental Table 4. A schematic to describe the different manually annotated pseudogene biotypes is presented in Figure 3. For example, unitary pseudogenes (i.e., genes that are active in mouse but pseudogenic in the human lineage) were all manually checked for false positives due to genomic sequencing errors or incorrect automated gene predictions in the mouse (Zhang et al. 2010).

Open in a separate window

Figure 3.

A schematic showing the structural annotation of different pseudogene biotypes. The schematic diagram illustrates the categorization of GENCODE pseudogenes on the basis of their origin. Processed pseudogenes are derived by a retrotransposition event and unprocessed pseudogenes by a gene duplication event in both cases, followed by the gain of a disabling mutation. Both processed and unprocessed pseudogenes can retain or gain transcriptional activity, which is reflected in the transcribed_processed and transcribed_unprocessed_pseudogene classification. Polymorphic pseudogenes contain a disabling mutation in the reference genome but are known to be coding in other individuals, while unitary pseudogenes have functional protein-coding orthologs in other species (we have used mouse as a reference) but contain a fixed disabling mutation in human.

Computational approaches followed by experimental validation were implemented to examine how many GENCODE pseudogenes appeared to be transcribed (Pei et al. 2012). Briefly, transcribed pseudogenes were identified manually and tagged by the HAVANA team examining locus-specific transcription evidence (by aligning of mRNAs or ESTs). This identified 171 transcribed processed and 309 unprocessed pseudogenes. The locus-specific transcriptional evidence must indicate a best-in-genome alignment and clear differences compared with the parent locus. Interestingly, there was over one-third more unprocessed pseudogenes annotated as transcribed compared with processed pseudogenes, even though there are approximately four times as many processed pseudogenes present in the genome than unprocessed pseudogenes (see Supplemental Table 4). In addition, automated pipeline analysis of RNA-seq data from the total RNA of ENCODE cell line GM12878 and K562 plus HBM RNA-seq resource (Pei et al. 2012) generated an additional 110 and 344 transcribed processed and unprocessed pseudogenes, respectively. Specific primers could be designed for 162 potentially transcribed pseudogenes and have been subjected to experimental validation of transcription by the RT-PCR-seq pipeline within the GENCODE Consortium (Howald et al. 2012). After the validation experiments, 63 pseudogenes were found to be transcribed within at least one of eight tissues.

To summarize, we have manually annotated 11,224 pseudogene loci that are in GENCODE 7 level 2 category. These have been compared with automatic models to produce a consensus set of 7183 loci that have been elevated to level 1. Of these, a total of 480 pseudogenes have been manually tagged as “transcribed” based on publicly available EST or cDNA evidence, and an additional 454 pseudogenes have been identified as transcribed due to alignment of HBM data (Pei et al. 2012).

The GENCODE gene set has developed substantially between releases 3c and 7 (see Fig. 4). Release 3c was the first complete merge set containing all the CCDS transcripts and used by the 1000 Genomes Consortium as its reference annotation. GENCODE release 7 is the reference for the analysis of ENCODE project data carried out in 2011. First-pass manual annotation has been done on 18 chromosomes (chr), and HAVANA still has chr14–19 to complete before the whole genome has been fully manually annotated. Supplemental Figure 4 demonstrates how the number of lncRNAs has increased dramatically with the full manual annotation of the chromosomes. The number of protein-coding loci has decreased significantly between GENCODE releases 3c and 7; however, this is almost entirely due to the removal of poorly supported automatic annotation models, particularly between releases 3c and 3d, where 1004 models were removed from the automatic annotation set. All GENCODE small noncoding RNAs are Level 3 and, as such, show a different pattern to other locus biotypes with their numbers dropping between 3d and 4 as incorrect automated gene models are removed and remaining stable thereafter.

Open in a separate window

Figure 4.

Analysis of GENCODE annotation in 3c through 7. (A) The content of the GENCODE 3c to 7 at the locus level for four broad biotypes: protein-coding, pseudogene, long noncoding RNA (lncRNA), and small RNA (sRNA). The yellow section of each column indicates the proportion of loci classified as Level 1 (validated), the blue part as Level 2 (manually annotated), and the red part as Level 3 (automatically annotated). (B) The analysis of the content of the GENCODE 3c to 7 at the level of the individual transcript. Again, the yellow section of each column indicates the proportion of transcripts classified as Level 1, the blue part as Level 2, and the red part as Level 3.

The patterns of change in the GENCODE loci across releases 3c to 7 are reproduced at the transcript level. It is clear from the data that the vast majority (>75%) of transcripts are associated with protein-coding loci. While the total number of protein-coding loci is decreasing, the number of coding locus transcripts is increasing with each release. The lncRNA transcript numbers show less stability than protein-coding loci and pseudogenes because of the novel status of the whole area of lncRNAs and the method of identification has changed, for example, based on chromatin signatures or position relative to a protein coding gene. A key change in lncRNA transcripts between releases 3c and 7 is the introduction of a more refined set of biotypes for Level 1 and 2 transcripts (see Supplemental Table 4), specifically the number of transcripts with the biotype processed_transcript reduced significantly and the number of antisense, lincRNA, noncoding and sense_intronic biotypes correspondingly increases.

To assess whether a locus or transcript is full-length, it is necessary to identify the TSS and transcription termination site (TTS). TSSs may be tested by determining overlap with CAGE tags (Takahashi et al. 2012) and the TTS by the presence of polyadenylation feature (polyA signals and polyA sites) (Ara et al. 2006). The number of protein-coding genes with at least one polyadenylation feature annotated increased by ∼10% between 3c and 7, and in the later release, the majority (∼62%) of protein-coding loci have one or more polyA sites annotated. However, the mean number of individual transcripts with a polyA site decreased between versions 3c and 7 as the overall number of annotated transcripts increased (Supplemental Table 6). Transcripts at both coding and lncRNA loci show a reduction in the percentage with polyA site of ∼2.5%. Figure 5 shows the number of polyA features increases in every annotated chromosome between release 3c and 7 for protein-coding gene loci.

Open in a separate window

Figure 5.

Comparison of polyA features annotated across all chromosomes. The mean number of polyA features (sites plus signals) for all protein-coding loci are plotted for every chromosome for GENCODE 3c (red columns) and 7 (blue columns).

While around 62% of GENCODE coding loci have support for identification of the TTS of the locus, the picture is less clear at the TSS. Around 65% of GENCODE annotated TSSs do not overlap with CAGE tag clusters generated by the ENCODE transcript subgroup (Djebali et al. 2012). One likely reason for the disparity between the validation of the two ends of the gene is the limited tissue coverage of the TSS set, which is derived from 24 whole-cell line polyA+ CAGE experiments (Djebali et al. 2012).

To assess the consistency in the structure of annotated gene models between releases, the number of exons per transcript was plotted for all splice variants at protein-coding and lncRNA loci in releases 3c and 7 (see Fig. 2). It is clear that, although their numbers have increased in release 7, the distribution of the numbers of exons per transcript suggests that the models themselves are very consistent in structure. Transcripts annotated at protein-coding loci demonstrate a peak at four exons per transcript, while lncRNAs show a very similar pattern given the large increase in their numbers between 3c and 7, with a distinct peak at two exons. This analysis confirms a high degree of homogeneity in the structure of transcripts annotated between releases 3c and 7. While the structure of annotated transcripts is invariant, there is a difference between the annotation of UTRs in those models in releases 3c and 7 (Fig. 6). Both the mean 5′ UTR and 3′ UTR length increase with each release between 3c and 7, with the mean 5′ UTR more than 41 bases longer in release 7 and the mean 3′ UTR 180 bases longer. Both increases are likely to be due to an increase in the amount of transcript data (ESTs and mRNAs) available to support the extension of transcripts.

Open in a separate window

Figure 6.

Examining the length of 5′ and 3′ UTRs between GENCODE 3c and 7. (A) The length of 5′ UTR sequence (in 50-bp bins) for each protein-coding transcript. 5′ UTR annotation from GENCODE 3c (red) and 7 (blue). (*) A cutoff was made at 949 bases; longer 5′ UTRs do exist. (B) The length of 3′ UTR sequence (in 250-bp bins) for each protein-coding transcript. 3′ UTR annotation from GENCODE 3c (red) and 7 (blue).

We compared the composition of annotation across the five major gene sets publicly available in UCSC, GENCODE, CCDS, RefSeq, and AceView. Both the number of protein-coding loci and transcripts at those loci were investigated. The CCDS set has the lowest number of protein-coding loci and alternatively spliced transcripts since it is a high-quality conservative gene set derived from RefSeq and Ensembl/HAVANA gene merge (Pruitt et al. 2009). In CCDS, every splice site of every transcript must agree in both the RefSeq and Ensembl/Havana gene set and all transcripts must be full-length. While the number of protein-coding loci in RefSeq, GENCODE, and UCSC is comparable, AceView has ∼20,000 more coding loci. One likely source of inflation is the predisposition for AceView to add a CDS to transcript model and hence create novel loci from lncRNAs and pseudogenes (e.g., PTENP1). AceView predicts 31,057 single exon loci compared with 1724 in GENCODE, 3234 in RefSeq, and 4731 in UCSC genes. Excluding single exon loci predicted by AceView from this analysis, the number of AceView gene loci is much closer to the number in other gene sets (Fig. 7A).

Open in a separate window

Figure 7.

(A) Comparing different publicly available gene sets. The protein-coding content of five major publicly available gene sets— GENCODE, AceView, consensus coding sequence (CCDS), RefSeq, and UCSC—were compared at the level of total gene number, total transcript number, and mean transcripts per locus. (Blue) GENCODE data; (orange) AceView; (yellow) CCDS; (green) RefSeq; (red) UCSC. The lncRNA content of three of these gene sets—GENCODE, RefSeq, and UCSC—were also compared at the level of total gene number, total transcript number, and mean transcripts per locus. Again, GENCODE data are shown in blue, RefSeq in green, and UCSC in red. (B) Overlap between GENCODE, RefSeq, and UCSC at the transcript and CDS levels. Both protein-coding and lncRNA transcripts of all data sets were compared at the transcript level. Two transcripts were considered to match if all their exon junction coordinates were identical in the case of multi-exonic transcripts, or if their transcript coordinates were the same for mono-exonic transcripts. Similarly, the CDSs of two protein-coding transcripts matched when the CDS boundaries and the encompassed exon junctions were identical. Numbers in the intersections involving GENCODE are specific to this data set, otherwise they correspond to any of the other data sets.

The GENCODE gene set contains 140,066 annotated alternative transcripts at coding loci compared with 66,612 in UCSC genes and 38,157 RefSeq. However it must be noted that not all GENCODE transcripts are full length, and if an annotated transcript is partial, it is tagged with start_not_found or end_not_found to highlight this to the user. The GENCODE gene set has 9640 lncRNA loci compared with 6056 in UCSC genes and 4888 in RefSeq. The three transcript data sets (UCSC, RefSeq, and GENCODE) were compared computationally to see how many transcripts were contained in all data sets and how many were unique to each data set (Fig. 7B). As expected, the majority (89%) of CDSs from RefSeq matched in all data sets exactly since Ensembl and UCSC genes use RefSeq cDNA in their automatic pipelines. However, GENCODE has 33,977 unique coding sequences outside RefSeq compared with 18,712 in UCSC genes. Of these unique transcripts, there are only 9319 exact matches in both these sets, indicating the different methods of annotation and the way they interpret EST data.

We analyzed the annotated CDS in GENCODE 7 using the data in the APPRIS database (http://appris.bioinfo.cnio.es/). APPRIS defines the principal variant by combining protein structural, functional, and conservation information from related species in order to determine the proportion of transcripts that would generate functional isoforms with changes to their protein-coding features relative to the constitutional variant. Of the 84,408 transcripts annotated as translated in the GENCODE 7 release, 30,148 (35.7% of all transcripts or 47.3% of alternative transcripts) would generate protein isoforms either with fewer Pfam functional domains (Finn et al. 2010) or with damaged Pfam domains with respect to the constitutional variant for the same gene. Twenty-six thousand nine hundred fifty-five isoforms (31.9% of all isoforms or 42.3% of alternative isoforms) would have lost or damaged structural domains, based on alignments with Protein Data Bank (PDB) structures, and 16,540 isoforms (19.6% of all isoforms or 26% of alternative isoforms) would lose functionally important residues.

In total, 44.9% of the translated isoforms (59.9% of the alternative isoforms) would lose either functional or structural domains or functional residues relative to the constitutive isoform. For the “putative” and “novel” biotypes and for those isoforms predicted to be NMD targets, these figures were much higher: 71.5% of all putative transcripts, 62.7% of all novel transcripts, and 79.5% of all predicted NMD transcripts would give rise to isoforms with loss of protein functional or structural information.

The location of potential alternative isoforms would also be affected for some isoforms. A total of 1287 genes would generate isoforms with differing numbers of cellular or mitochondrial signal sequences, and 2697 genes would generate isoforms with differing numbers of trans-membrane helices. In addition, 8842 genes contained at least one transcript that appeared to be subject to non-neutral evolution. We were able to select a constitutive variant for 16,553 protein-coding genes (80%) based on the output of the methods in APPRIS (Tress et al. 2008).

It is important that users understand how to assess transcript annotations that they see in GENCODE. Some transcript models have a high level of support through the entirety of their structure, while other transcripts are poorly supported with one or two ESTs. A new method is being developed to differentiate this level of support. This method relies on mRNA and EST alignments supplied by UCSC and Ensembl. The mRNA and EST alignments are compared to the GENCODE transcripts, and the transcripts are scored according to how well the alignment matches over its full length (for more details, see the Supplemental Methods). Figure 8 shows support level statistics on a set broken down by annotation method for GENCODE 7. The annotations are partitioned into those produced only by the automated process, and those only from the manual method and the merged annotations, where both processes result in the same annotation. For the analyzed subset of GENCODE, 38% of the transcripts have good full-length support and an additional 39% have a lower level. The different distribution of support levels between the annotations is consistent with the underlying approaches. The merged subset is heavily weighted toward well-supported transcripts. The higher fraction of manual-only annotations supported only by ESTs is due to Ensembl not directly using ESTs in their pipeline.

Open in a separate window

Figure 8.

Quality of evidence used to support automatic, manually, and merged annotated transcripts. The level of supporting evidence for automatic only (A), manual only (B), and merged (C) annotated transcripts is shown for each chromosome. (Yellow) The proportion of models with good support; (dark blue) those supported by suspect mRNAs from libraries with known problems with quality; (light green) those with multiple EST support; (orange) those with support from a single EST; (red) those supported by ESTs from suspect libraries; (pale blue) those lacking good support. The number of transcripts across all chromosomes represented in A is 23,855; B, 89,669; and C, 22,535.

We have focused on removing nonconsensus introns, those not matching the known splicing patterns of GT..AG, GC..AG, and AT..AC. We found ∼2200 such introns in GENCODE 7, which is a reduction from 3c, where there were ∼3300 nonconsensus introns. For version 7, 13% of these show good cDNA support. A small number of nonconsensus splice sites are believed to be other donor/acceptor combinations recognized by U12 spliceosomes. Others are suspected polymorphisms, some with SNP support. Only 2% could be converted to canonical splice sites by known SNPs from dbSNP version 132 in some members of the population, but this proportion has doubled since 3c, reflecting the recent increase in SNP discovery from genome-wide projects. There may be more of these splice sites that are polymorphic but the SNPs are as yet unknown. There are a number of nonconsensus splice sites that are only one base different from canonical, suggesting that they could have been formed due to mutation, but some may represent low frequency polymorphisms. As the 1000 Genomes Project (http://www.1000genomes.org) and other human sequencing projects progress, more SNPs will be discovered, and therefore, the number of known polymorphic splice sites should increase.

The aim of the experimental validation group of GENCODE is to identify gene models that have limited or lower confidence transcribed evidence and to systematically experimentally validate them. Predicted exon–exon junctions were evaluated by RT-PCR amplification in eight different tissues followed by highly multiplexed sequencing readout, a method referred to as RT-PCR-seq (Howald et al. 2012). Eighty-two percent of all assessed junctions (n = 5871) are confirmed by this evaluation procedure, demonstrating the high quality of the annotation reached by the GENCODE gene set. RT-PCR-seq was also efficient at screening gene models predicted using the HBM RNA-seq data. We validated 73% of these predictions, thus confirming 1168 novel genes, mostly noncoding, which will further complement the GENCODE annotation (Howald et al. 2012). Our RT-PCR-seq–targeted approach can also be exploited to identify novel exons. We discovered unannotated exons in ∼10% of assessed introns. We thus estimate that at least 18% of loci contain as yet unannotated exons.

This novel experimental validation pipeline is significantly more effective than unbiased transcriptome profiling through RNA sequencing, which is becoming the norm. Exon–exon junctions specific to a single GENCODE annotated transcript are five times more likely to be corroborated with our targeted approach than with extensive large human transcriptome profiling such as the HBM and ENCODE RNA-seq (validation rates of 82% and 16%, respectively), as random sampling of RNA molecules leads to poor assessment of low expressed transcripts. It should be noted that the RNA-seq samples were deeply sequenced with the resulting data sets and contained a total of 4.99 and 5.55 billion sequence reads, respectively.

Our work demonstrates that the cataloging of all the genic elements encoded in the human genome will necessitate a coordinated effort between unbiased and targeted approaches, like RNA-seq and RT-PCR-seq, respectively (Howald et al. 2012).

To identify novel coding and noncoding genes represented in RNA-seq data, we studied transcript models reconstructed using Exonerate (Howald et al. 2012) and Scripture (Guttman et al. 2010), based on the high-depth HBM transcriptomic data from 16 tissues made publicly available from Illumina (ArrayExpress accession: E-MTAB-513; ENA archive: ERP000546).

We analyzed the resulting transcripts that did not overlap any GENCODE loci for coding potential using PhyloCSF (Lin et al. 2011), which examines evolutionary signatures within UCSC vertebrate alignments, including 33 placental mammals. There were 136 Ensembl HBM models with positive PhyloCSF scores out of a total of 3689 loci, although only five of these had sufficient support for manual reannotation as coding genes (see Supplemental Table 8). The remaining 131 transcripts showed varying quality and evidence; ∼50% overlap novel processed transcripts and could be a result of misalignment of reads or actual expressed pseudogenes. Two hundred Scripture transcript predictions that were outside GENCODE but had high PhyloCSF scores were also manually examined. Of these, 15 were added as novel loci, and only nine were annotated as coding genes (see Supplemental Table 9) and will be added to the next release of GENCODE). Considering the depth of reads of the HBM data (averaging over a billion read depth) from the 16 different tissues, we have not identified many missing coding genes based on PhyloCSF. Indeed, since 3127 HBM Ensembl genes consist of only two exons, it is highly likely these constitute new lncRNAs we have not yet annotated and will be merged into a later release of GENCODE.

A pipeline has been set up to verify the annotation of gene models with mass spectroscopy data from human proteomics experiments (M Tress, P Maietta, I Ezkurdia, A Valencia, J-J Wesselink, G Lopez, A Pietrelli, and JM Rodriguez, in prep.). The data from tandem mass spectrometry experiments are stored in two huge proteomics data repositories, the GPM (Craig et al. 2004) and Peptide Atlas (Desiere et al. 2006). Peptides are detected by mapping spectra from individual proteomics experiments to the gene products from the GENCODE annotation using the search engine X!Tandem (Craig and Beavis 2004). A single peptide may be detected in many different experiments, though only once per experiment. We generate P-values for all detected peptides by combining the X!Tandem P-values for each individual experimental peptide-spectrum match, and a target-decoy approach is used to determine false-discovery rates.

The pipeline was able to map peptides to almost 40% of the protein-coding genes in the GENCODE 7 release. The 83,054 tryptic peptides detected at a false-discovery rate of 1%, mapped unambiguously to 8098 of the 20,700 annotated protein-coding genes. We were able to detect the translation of multiple splice isoforms for 194 genes, and within this set of genes, we validated the expression of eight isoforms that were tagged as candidates for NMD degradation. We found peptide support for 33 transcripts annotated as “putative” and another 50 transcripts annotated as “novel.” With the mass spectroscopy data we generated for the GENCODE 3c release, we also detected the expression of peptides that mapped to a pseudogene (MST1P9) (Pei et al. 2012).

The AnnoTrack software system (Kokocinski et al. 2010) was developed as part of the GENCODE collaboration to facilitate the processing and tracking of the HAVANA annotation and the heterogeneous sources of information used in the genome annotation. It integrates data from distributed sources via DAS (Distributed Annotation System) servers set up by the GENCODE partners following a defined format (http://www.gencodegenes.org/gencodedasformats.html), a direct database connection to Ensembl, and flat file adaptors. The system highlights conflicts and facilitates the quick identification, prioritization, and resolution of problems during the process of genome annotation. Using controlled terms for the solutions chosen by the manual annotators when resolving conflicts allows a retrospective assessment in order to improve the analysis methods or external data sources. More than 4000 issues have been resolved this way in the last year. AnnoTrack also helps to track the progress of the overall HAVANA annotation and to inform external scientists about current issues about the genes they might be working with. The interface can be accessed at http://annotrack.sanger.ac.uk/human.

As the GENCODE gene set is now built partly through the Ensembl pipeline, the GENCODE data release cycle is coupled to the trimonthly Ensembl releases. Dates and release notes, as well as more details of the data sets and formats, are listed at http://www.gencodegenes.org. The GENCODE releases contain updated gene sets where either new data from the manual annotation has been integrated as described above or additionally the automated gene set was rebuilt or refined. Users can view GENCODE data in the UCSC browser (Fig. 9), and also it is the default gene set shown in Ensembl. All genes and pseudogenes within the release have stable Ensembl (ENS) identifications and the manual annotated genes have additional Vega (OTT) IDs (Wilming et al. 2008). All OTT identifications are also versioned so the user can identify when a transcript was last manually updated.

Open in a separate window

Figure 9.

Accessing the GENCODE gene set through UCSC and Ensembl. (A) The composite of screenshots from the UCSC browser shows GENCODE gene annotation displayed in the basic and comprehensive display mode, along with the GENCODE pseudogenes, CCDS models, and a subset of histone modification tracks, DNaseI hypersensitivity clusters, and transcription factor binding site tracks. (B) The configuration display where the user can filter on biotype, annotation method, and transcript type (C). (D) The transcript page in UCSC where the different identifications and version of the transcript can be seen, as well as the evidence used to build the transcript. From the page, the user can click on the Ensembl identification and immediately jump to the Ensembl gene view page (E) and see an overview of the different transcripts in the locus as well as which is a CCDS.

Figure 9 shows the UCSC browser with the “basic” view, which only displays full-length transcripts for every loci unless a partial transcript is the only representative of that locus. This is in demand from users that want a more compact view when browsing the genome. However, the “comprehensive” view shows everything annotated and is designed for the user that is interested in all annotated transcripts in a particular region of the genome. The data set can also be filtered by biotype and annotation method to display the transcript, and submitted nucleotide evidence that each particular transcript was built on is also available on the transcript page. GENCODE data are released and accessible in various formats, which are described in Supplemental Table 8. The gene sets and supporting data are submitted as GTF files to the ENCODE Data Coordination Center (DCC) for integration in the UCSC genome browser and for redistribution.

Discussion

As sequencing technology steadily improves and becomes cheaper a thousand dollar genome will soon become a reality. However, utilizing this genomic data efficiently is still dependent on which reference annotation is chosen to highlight its gene landscape and variation between other individuals. The GENCODE human reference set is a merge of automatic and manual annotation, thus producing the most comprehensive gene set publicly available. Its regular release cycle of every 3 mo ensures that models are continually refined and assessed based on new experimental data deposited in the public databases. In addition, users can send updates or queries concerning individual genes using the Annotrack system or via the Gencodegenes.org website to trigger investigation and updates of specific issues alongside those flagged by global quality-control (QC) checks.

Interestingly, there is still uncertainty about the number of protein-coding and long noncoding loci in the genome. Assessing how many protein-coding loci are missing from the catalog is difficult, but our analysis of coding potential using conservation indicates that the number is likely to be small, namely, around 100 protein-coding genes. A similar figure was suggest by Lindblad-Toh et al. (2011), who recently reported the sequencing and comparative analysis of 29 eutherian mammals. However the recent publication by Ingolia et al. (2011) suggests that there is a new class of small “polycistronic” ribosome-associated coding RNAs encoding small proteins that can now be detected using ribosome profiling. They highlight that the majority of predicted lncRNAs in the mouse from Guttman et al. (2009) actually show comparable translatability to that of protein-coding genes. In addition, Cabili et al. (2011) have found 2798 lincRNAs not in GENCODE 4 using a combination of HBM RNA-seq and additional RNA-seq from eight additional cell lines and tissues totaling 4 billion reads. This indicates that there are still many thousands of lncRNA loci to add to the GENCODE catalog, and completeness will be dependent on the depth and variety of tissues and cell lines sequenced.

The GENCODE catalog is highly specific in its subcategorization of protein-coding and noncoding transcripts highlighting transcripts subject to phenomena such as NMD and nonstop decay (Mazzoni and Falcone 2011). Interestingly mass spectrometry analysis within the project has identified four transcripts, annotated to be subject to NMD, that produce peptides. A recent publication by Bruno et al. (2011) has identified a brain-specific microRNA (miR-128) that represses the NMD pathway by binding to the RNA helicase UPF1. Thus suggesting certain miRNAs can induce cell-specific transcription/translation during development by inhibition of the NMD pathway.

Analysis of transcriptomic and proteomic data is also revealing pseudogenes that are potentially expressed. The GENCODE reference set aids such analysis since it is the only gene set to contain comprehensively manually annotated pseudogenes to the same level as protein and noncoding genes. We currently predict a total of around 10,000 pseudogenes within the human genome. Recent publications highlight the implications of pseudogenes as regulators of gene expression (Han et al. 2011) and specifically a role in tumor biology (Poliseno et al. 2010), and thus we will have to rethink the classification of pseudogenes as nonfunctional entities on the genome.

The GENCODE gene set gives the most comprehensive overview of alternative splicing than of any gene set available. However improved transcriptome sequencing is being facilitated by high-depth RNA sequencing (Mercer et al. 2011), which reveals that the full extent at which different transcripts are expressed in different cells and tissues and at different developmental stages is still to be fully characterized. The RT-PCR-seq methodology developed within the GENCODE project has identified novel exons within 10% of annotated introns being targeted (Howald et al. 2012). Therefore, even when the first pass of whole-genome manual annotation is finished in 2012, reannotation updates will be required on a large percentage of loci to correctly classify the new alternatively spliced variants being revealed by next-generation transcriptomics. However the major drawback with using next-generation short reads to assemble transcripts de novo is that the correct structure of the transcript is hard to predict, as was investigated in the RNA-seq Genome Annotation Assessment Project (RGASP) (J Harrow, T Steijger, F Kokocinski, JF Abril, C Howald, A Reymond, A Mortazavi, B Wold, T Gingeras, R Guigó, et al., in prep.). We therefore look forward to “third-generation” sequencing methods such as Pacific Biosciences (Schadt et al. 2010) that show promise of generating longer reads of 1–2k, which will enable improvements in transcriptome annotation, facilitating the investigation of expression of each transcript structure within a cell during its development.

Methods

Data access

Dates and release notes are listed on the website http://www.gencodegenes.org, as well as more details of the data access, which are listed in Supplemental Table 11.

Acknowledgments

Footnotes

References