A proteogenomic update to Yersinia: enhancing genome annotation

Correcting Annotation Errors

Following the data path outlined in the Methods (Figure 1), ~15 million MS/MS spectra from Yersinia pestis KIM were searched by Inspect and PepNovo against the six-frame translation of the genome. Confident peptide/spectrum matches were mapped onto the genome sequence, and used to infer annotation improvements. In total, we report 30,994 peptides mapping to 1302 proteins when requiring 2 peptides per protein. 1277 proteins (31% of the proteome) contain at least 1 uniquely mapping peptide. Of the 25 proteins lacking a unique peptide, the vast majority are an active transposase (IS285). The finding of an active transposon is intriguing given than transposons have been proposed as a driving force of Y. pestis genome evolution since its recent divergence from Y. pseudotuberculosis[13].

By mapping the proteomic data onto the genome, we are able to objectively determine the quality of the genome annotation. When peptides map outside of predicted proteins, there are two categories of annotation improvement. If the open reading frame lacks a predicted protein, the observed peptides are evidence for a novel gene. If the ORF contains a predicted protein and peptides map upstream, they are evidence for a 5' extension and new start site. Both of these situations necessitate an update to the genome annotation. Working with the RefSeq curators at NCBI, all of the instances discussed below have been updated, and can be accessed seamlessly through any NCBI tool (e.g. blast) or downloaded directly.

We find four ORFs which lack protein annotation, but have at least two uniquely mapping peptides (Figure 2). Three are easily recognizable proteins: major outer membrane lipoprotein between y1943 and y1944, two cold-shock proteins between y1817 and y1818 and y2562 and y2563. We also report an apparent Yersinia specific protein between y2035 and y2036. This protein, now named y5001, lacks homology to any currently described protein domain, and has no significant blast hits outside of the Yersinia genus. Finally, y3734 is predicted to be a pseudogene, an ABC transporter disrupted by insertion elements, indels, and nonsense mutations (Additional File 1, Figure S1). These mutations destroy several functionally relevant motifs, and almost certainly preclude proper biochemical function. However, we find two peptides mapping to the n-terminus of the region, providing evidence for its translation and presence in the cell.

The second major class of proteogenomic correction is start site annotation. Yersiniabactin thioesterase is an enzyme participating in the biosynthesis of a siderophore important for iron acquisition from the host. In our proteomics data we observed several peptides upstream of the current start for this gene, which required a 44 residue extension. Similarly, numerous peptides found upstream of the Y. pestis specific protein y0291 pointed to a 40 amino acid extension. Its lack of homology to any known domain and narrow taxonomic distribution demonstrate the utility of proteomic involvement in gene prediction. A third gene, y2368, highlights a subtle ramification of erroneous start sites, which is that functional elements within the n-terminus are obscured. y2368 is annotated simply as a 'periplasmic protein' with a CDD iron transport domain (cl01377). As expected we found this protein to be highly enriched in periplasmic as compared to cytoplasmic fractions [14], yet the n-terminus lacked localization motifs, e.g. a signal peptide. Six peptides mapped upstream of the annotated start site (Figure 3). Manual inspection of the region immediately upstream from these peptides found a strong signal peptide sequence motif.

Two proteins with erroneous start sites are special cases and widely mispredicted in bacterial genomes. Peptide chain release factor II, prfB, often contains a ribosomal +1 frame shift. As with the KIM genome, erroneous annotations of this protein typically contain only the c-terminal ORF but exclude the true n-terminus (Additional File 1, Figure S2). Our proteogenomics pipeline recognized peptides in both ORFs and the two were stitched together. We also found peptides upstream of infC, protein chain initiation factor IF-3, which utilizes an ultra-rare start codon ATT [15].

Signal Peptides

Proteins exported from the cytoplasm through the Sec-dependant pathway contain a short sequence essential to targeting and export (Figure 4). The ~20 residue motif, or signal peptide, is located at the n-terminus of the full-length protein. The signal peptide helps target the protein to the membrane, where it is temporarily anchored by a patch of hydrophobic residues. A three amino acid motif following the hydrophobic patch is recognized by the signal peptidase enzyme; the protein is cleaved and the c-terminal portion of the protein is released into the periplasm. The signal peptide, still anchored in the membrane is rapidly degraded by the signal peptide peptidase. The Sec-dependent pathway is separate from other export pathways, such as the Type III secretion system. Signal peptides are also present in proteins exported through the Twin Arginine receptor-mediated pathway, but are typically longer and include an additional motif prior to the hydrophobic patch.

Proteolytic cleavage of proteins in vivo can be recognized by proteomics by their atypical peptide endpoints. Spectra used in this report were generated from proteins digested with trypsin. Thus, we expect most identified peptides to be fully tryptic. Previously Gupta and colleagues postulated that signal peptides could be discovered simply by identifying non-tryptic peptides [16]. In their analysis, if the first observed peptide in a protein had a non-tryptic n-terminus and was within 17-55 residues of the start site, then it was a considered evidence of signal peptide cleavage. 202 Yersinia proteins fulfill these two requirements. We extend Gupta's criteria to include critical biological motifs within the signal peptide: the hydrophobic patch and cleavage motif [17]. Filtering out proteins lacking these new requirements, we report 82 proteins with observed signal peptide cleavage (Additional File 2, Table S1). These proteins also contain other common signal peptide features: prevalence of LL doublets and early basic residues. Furthermore, we noticed that the hydrophobic patch had a similar placement within the signal peptide for all validated proteins (Figure 4B). This location is consistent with the patch's structural purpose, i.e. membrane anchoring and exposure of the cleavage motif at an appropriate distance from the membrane. Finally, we report that many of the sequences contain not simply early basic residues but contain Met-Lys as the first two residues in the protein sequence. 34 of the 82 proteins start in this manner. An additional 15 have Met-Lys internally which could be the true start, if this trend is seen as a general pattern.

Twenty proteins contain signal peptides longer than 30 residues, which is atypical. When we compared these sequences to close homologs within the gamma-proteobacteria, 13 of them could be better predicted at an alternative start site downstream of the current annotation. This shorter version of the protein would not only have better agreement with homologous genes, but also the characteristics of signal peptides noted above (not just length). An additional four long signal peptides appear to contain the twin arginine motif, which are characteristically longer than those exported through the Sec-dependent pathway.

The set of secreted proteins within a genome is often computationally predicted. To compare our proteomic observations to such predictions, we ran signalp on all proteins in the genome. We plot the score for proteomically validated signal peptide containing proteins against the background of signalp's score distribution (Figure 4C). We also plot proteomically rejected signal peptide containing proteins (see Methods). The proteomically observed and rejected proteins separate very clearly, with the positive set scoring well above the suggested cutoff. Furthermore, proteomics and signalp generally agree on the exact residue of cleavage.

Dubious Genes

In the Yersinia pestis KIM genome, there were over 200 genomic loci with a >50 bp overlap between two protein coding genes. Such a substantial overlap is unusual, especially considering that 10% of the proteome (~400 proteins) falls into this category. We viewed these conflicted loci as unlikely to be correctly annotated. For the 46 loci covered by proteomics, we manually reviewed the evidence supporting the existence of either gene. 38 loci contained a dubious gene (Figure 5). In over half of the instances, the loci contained genes with 100% sequence overlap. The equivocal nature of dubious genes was witnessed by narrow phylogenetic distribution, poor and seemingly random sequence conservation, and weak computational justification noted in the original genome submission (see Methods). The remaining 8 loci were in conflict due of an overly extended 5' on one or both of the genes (Figure 5C). Working with RefSeq curators at NCBI, the dubious genes have been removed. Analysis of the remaining ~150 conflicted loci is being addressed in future work.

Improvement of Related Genomes

As a final step in our analysis, we used our proteomics data to improve the annotation of other Yersinia genomes. Published in 2002, KIM was the second annotated Yersinia and served as a source for subsequent genome projects [18]. Thus errors in KIM's annotation are likely to show up in other genomes. We created one-to-one orthology maps for proteins from 21 Yersinia genomes, and used our peptide mappings from the KIM dataset to evaluate gene models. Transferring proteogenomic improvements is complicated by both legitimate differences between strains/species, but also artifacts from sequencing and assembly. Our approach was to be conservative, and only change gene models where the evidence was clear.

We started with the removal of dubious genes, resulting in the deletion of 46 accessions. Remembering that each dubious gene was part of a conflicted locus, we checked to make sure that the true gene was present. At two loci we found that the true gene was missing, excluded from the original annotation by the dubious gene. We created new gene models for each of these. Three additional new genes were created, arising from the transitive annotation of the novel KIM genes. The two cold shock proteins were missing from Y. pseudotuberculosis IP 31758. The major outer membrane lipoprotein, lpp, was missing from Y. pestis Pestoides F.

To correct start sites in other Yersinia genomes, we applied the n-terminal most peptide from a gene in KIM to members of its ortholog set. If an ortholog was too short to include this peptide, we tested whether an extension could be made to accommodate. In 43 instances, this extension was trivial - clear homology up to and including the start codon. For 22, an extension could not be made due to indels or mutations shortening the reading frame. Finally, 10 were difficult cases where we could extend the reading frame to an upstream start codon, but homology was unclear. In these cases, we did not alter the genome annotation. This was particularly problematic in Y. enterocolitica, which is the most divergent genome used in our comparison.

Mass Spectrometry Data

Spectra used in this report have been previously published. Briefly, the Pathogen Functional Genomics Resource Center (PFGRC) at JCVI generated ~15 million spectra from Y. pestis KIM 6+. Cell lysates were fractionated into cytoplasm, membrane and periplasmic components, and then run on 2 D gels, with spectrum acquisition on a Thermo LTQ [14, 22, 23]. All spectral data has been uploaded to NCBI's Peptidome http://www.ncbi.nlm.nih.gov/peptidome/.

Data Processing

Our pipeline can be divided into distinct steps, with the overall goal of mapping peptides onto the genome and then inferring genome annotation (Figure 1). The first step is to prepare protein sequence databases by translating the RefSeq DNA into all six frames: genome NC_004088, and plasmids NC_004836, NC_004837, NC_004838. We do not make any restriction of open reading frames (ORFs). For example, we do not require a minimum length or the presence of a start codon. In addition to these Yersinia sequences, we append common contaminants (trypsin and keratin) to the database. To finish database preparation, we create a decoy database by shuffling each sequence [24].

We use the Inspect software to match spectra to peptide sequences from the protein database [25]. Search parameters: 25 tags/spectrum, cysteine + 57 fixed modification, and 3.0 Da parent mass tolerance. Following Inspect, we rescore peptide/spectrum matches (PSMs) with PepNovo [26]. We assign a p-value to PSMs based on the score distribution of hits to the target and decoy databases and set a 5% local false-discovery rate cutoff [21]. PSMs passing the cutoff are mapped back to their genomic ORF, and subjected to protein-level filters. To leverage the gel-based nature of our experimental set up, we extend the common 2 peptide heuristic to require 2 peptides per protein per gel spot. Manual inspection of the remaining PSMs (now grouped to their ORFs and reported as proteins) revealed a category of questionable results. We noticed PSMs mapping to low sequence-complexity regions, with significant overrepresentation of small amino acids (e.g. glycine or alanine). We believe that such sequences are likely to match noise spectra because the small mass allows for spurious peak matching. We also noticed that none of the PSMs in these objectionable clusters are tryptic. As trypsin was used to generate peptides for the experiment, we reject ORFs that lack any tryptic PSMs.

Proteogenomic Improvement

We use PSMs passing the above filters to improve the RefSeq annotation of NC_004088 in three ways: adding novel gene models, correcting the start site of current models, and deleting dubious protein predictions. The updated KIM genome can be accessed seamlessly through any NCBI tool (e.g. blast) or downloaded directly ftp://ftp.ncbi.nih.gov/genomes/Bacteria/Yersinia_pestis_KIM_10_uid288/. We create a new gene model (novel gene) when at least 2 PSMs uniquely map to an ORF that lacks any protein annotation. We examine start sites when PSMs map upstream of and in frame with a current protein prediction. Through manual inspection, we adopted a simple heuristic, that the peptide supporting an extension must contain at least two amino acids upstream of the currently annotated start. With only a single amino acid upstream, it was possible for the amino acid to have the same mass as common modifications (e.g. glycine and carboxyamidomethyl (CAM) modification have a mass of 57).

We found many peptides which mapped to genomic loci containing multiple protein predictions. A closer inspection revealed that these loci often contained proteins with significant overlap. We marked genomic regions containing two or more genes with a 50 bp overlap as potentially in conflict. In our initial analysis, we divided conflicted loci into three categories defined by the relationship of the overlapping genes: convergent, divergent, unidirectional [27]. This delineation is attractive, because divergent and unidirectional overlaps can potentially be resolved by shortening one or both genes to remove the overlap. However, the extreme length and multiplicity of overlaps (some proteins overlapped two or three others) rendered this classification system awkward. Through visual analysis, we categorized conflicting loci according to the extent of overlap: 100% overlap, > 80% overlap, > 50% overlap, and < 50% overlap. It then became clear that in most instances, one of the overlapping genes was dubious and should simply be deleted.

When both genes had proteomic evidence, we asserted that both genes are real. In these instances, one or both of the genes was always overpredicted, i.e. was longer than homologs from other genomes. After correcting for this, the genes no longer overlapped in their coding regions. When only one gene had proteomic evidence, we tested whether the unobserved protein was a legitimate gene. Equivocal genes were determined using the following criteria. First, their extended sequence overlap (to a gene with proteomic support) is atypical, especially considering that many of them are completely contained within another protein's coding region. Second, they were sporadically predicted in bacterial genomes with poor and uncharacteristic conservation, while the proteomically supported gene had more compelling sequence conservation, was more consistently annotated within the genus/family, and had a generally broader phylogenetic distribution. Third, most of the dubious genes were short (median 95, mean 108). We note that 100 amino acids is the minimum for many genome annotation projects. We finally note that, if present, the computational evidence for asserting the equivocal gene in the original genome submission was weak (e.g. low sequence identity over short regions to phylogenetically distant organisms). Therefore, we term the equivocal gene as dubious, and have deleted it from the annotation.

Mature Protein Annotation

To report a protein as containing a signal peptide, we started with proteins where the first observed peptide was not tryptic on its n-terminus, and was within 15-50 amino acids of the predicted start site [16]. Trypsin cuts after Arg and Lys; if the residue before the peptide is not Arg or Lys, then the peptide has a non-tryptic n-terminus. The region between the initial methionine and the first observed peptide is thus the putative signal peptide. We filtered this set using previously recognized signal peptide characteristics [17]. We required a hydrophobic patch of at least 8 amino acids, using the standard Kyte-Doolittle hydropathy index. We then examined the signal peptide terminus for the expected AxB cleavage motif (where A = [Ile, Val, Leu, Ala, Gly, Ser], B = [Ala, Gly, Ser]). Signal peptides matching these two criteria were considered validated by proteomics.

We ran signalp v3.0 on all proteins in the genome [28]. Proteins listed as containing a computationally predicted signal peptide are those with a D-score above 0.43 as recommended. When comparing signalp against proteomic evidence for signal peptide presence, we simply look gene by gene, whether proteomics or signalp or both list the protein as containing a signal peptide (and at the site of cleavage). Proteins with proteomic validation which have a signalp score below 0.43 were viewed as computational false-negatives. Conversely where peptides localize 5' of the signalp predicted cleave site (at least 4 residues upstream), we viewed these as computational false-positives.

Proteogenomic Improvements in Other Genomes

We created orthology clusters of all genes in the 12 complete Yersinia genomes with the Sybil program [29]. Peptides observed in the KIM dataset were mapped onto the appropriate columns of an aligned cluster, and genome annotation improvements were inferred from this. For example, if a gene model in genome A was shorter than the orthologous member of the KIM genome, and peptide evidence from KIM mapped to the uncalled region, we checked for the ability to move the start site upstream in genome A. When genome improvements were possible (e.g. not precluded by stop codons) they were made and submitted to RefSeq. The same process was performed for the 9 Genbank annotations owned by JCVI. A list of all altered genes is presented in Additional File 3, Table S2. We are actively pursuing the ability to edit other genome annotations in collaboration with the genome owners.