Which Direction Does the Gene Point?

A maddening problem in genome annotation is determining the "sense" strand for a gene, especially when the gene is short and/or the genome has a high GC content (and thus contains few or no stop codons in reverse translation). To convince yourself this is a very real and serious problem, all you have to do is browse a few genomes to see the ridiculously high number of "hypothetical proteins" (over 40% in some genomes), bogus overlaps, genes that score BLAST hits in reverse-complement mode but not frame zero, and other artifacts that are a direct result of the aforesaid problem.

I've presented examples of this problem before, but just so there's no confusion, I want to show you a particularly maddening example so you can see what I'm talking about. (Then I'll suggest a solution.) The following graphic shows a region of E. coli UTI89 in which several genes are shown as overlapping (that is to say, existing on opposite strands of DNA in the same coverage area). Small overlaps sometimes happen between genes, but whole genes rarely, if ever, overlap, and never in clusters. The situation shown below is bogus, but you see it all the time in public genomes. In fact, some of the genes shown below also show up as overlaps in Mycobacterium abscessus M93 (see gene OUW_18941), Citrobacter koseri strain ATCC BAA-895 (gene CKO_00072), and quite a few others. Glimmer has choked on this exact situation many times, in many genomes.  

A region with overlapping genes in the genome of E. coli UTI89.

The big gene on the top strand, middle, is UTI89_C4288 (DNA sequence here). It's annotated as (what else?) a "hypothetical protein." The M. abscessus version of the gene (here) is marked as a "cellobiose phosphorylase," and you can find many BLAST hits (the Rothia version gives an E-value of 6.0×10^-101) for similar "cellobiose phosphorylase" genes in other organisms at UniProt.org and elsewhere. Of course, they're all bogus and represent Glimmer choke points, but the question is how one can determine that, and be sure about it.

E. coli's hypothetical protein (UTI89_C4288) has a wobble-base GC percentage of 64.3%, whereas the gene on the opposite strand (just below it, pointing left), namely UTI89_C4287 (marked as "membrane-bound ATP synthase F1 sector alpha-subunit"), has GC3 = 54.5%. In a much higher-GC organism like Mycobacterium or Pseudomonas, you would find out which gene has the higher GC3 percentage and crown it the winner (and most of the time, you'd be right). In this case, it's not so simple. The gene with the higher GC3 value isn't necessarily the winner.

Of course, in this particular example, you can cheat and look at the identities of the genes in the immediate vicinity of the hypothetical protein, on the bottom strand, and if you do, you'll find that all of the bottom-strand genes are ATPase subunits. Mystery solved, right? Sure, in this particular case. But what about situations where overlapping genes are all shown as "hypothetical protein"? (You can find many such cases in the genome for Burkholderia pseudomallei strain 1710b, for example.) When a hypothetical overlaps a hypothetical in a low-GC genome, then what?

One of my favorite cheats (but this isn't the final solution!) is to check the gene's AG1 percentage (adenine plus guanine, codon base one). This percentage averages ~60% in something like 90% of protein-coding genes. The problem is, AG1 is often 60% whether you read the gene forward, or backward (off the antisense strand). The reverse complement of a gene usually has high AG1, because the forward AG3 is usually under 50%.

Almost any trick you can dream up will fail under edge cases. GC3 is helpful, but only in high-GC genomes. AG1 is helpful, but only sometimes. Shine Dalgarno signals are not universally used by all organisms, and even in those that do use them, they're usually reserved for highly conserved genes encoding things like ribosomal proteins. Gene context is helpful in some cases but not others.

It turns out, the best clue for positively identifying the correct strand and correct reading frame is codon usage frequency patterns. If you know what the codon frequencies are, genome-wide, for a given organism, you can use this information to good advantage, even if the genome (and therefore the codon table) contains inaccuracies. As long as the codon frequencies are approximately correct, you can use them to verify the reading frame of a protein-coding gene.

The algorithm I came up with is very simple, yet effective. For a given gene, read each triplet of bases sequentially, and score each triplet twice: keep two scores going. First, score it according to its frequency in the codon table for the organism. Then score it according to a second table developed for reverse-complement codons.

The following table shows codon frequencies in Caulobacter crescentus NA1000. If you were to encounter a "hypothetical protein" gene in Caulobacter, and you couldn't decide whether the strand assignment was correct or not, first develop a score for the gene by reading its triplets and adding the frequency value of the corresponding codon to the running total. For example, if you encounter the triplet "CTG," add 6.84 to the score (see table). For every occurrence of CTT, add 0.60, for CTC add 1.70, and so on, using the values in the table.

Codon frequencies for Caulobacter crescentus NA1000.

But you also have to create an anticodon frequency table as follows: For every codon in the original table, apply the same score to the corresponding reverse-complement codon in the "antcodon table." E.g., for CTG, the first table would contain 6.84 (as above), but the value 6.84 would apply to CAG (the reverse complement of CTG) in the second table. I call the first table the "forward" table and the second table the "back" table. One represents the frequencies of codons encountered in protein genes in the forward reading direction. The other represents those same frequencies applied to the reverse-complement of the codons (the same codons read in the reverse direction, off the opposite DNA strand).

When scoring an unknown gene, you tally a "forward table" score, and keep a separate score using the "back" table. When you're done, the gene's "forward table" score should be greater than the "back table" score. If it's not, you're reading the gene off the wrong strand.

When I scored all 3,737 C. crescentus CB15 genes using this technique, I found 136 genes that gave a "back" score higher than the "forward" score. Interestingly, when I checked the identities of those 136 putative "backwards-annotated" genes, 132 of them were listed as "hypothetical proteins." Only four genes with assigned functions gave suspect scores, and one of those (CC_0662) turns out to be a 100%-identity match for the reverse of gene CCNA_00700 in Caulobacter crescenstus NA1000.  The other three are less than 200 bases long and could well be non-coding regions.

For a more challenging test, I turned to the genome of Rothia mucilaginosa DY-18, one of the most disastrous annotation nightmares of all time. In the genome for DY-18 you will find 524 protein-coding genes (out of 1,905 total) that are involved in significant overlaps. (Some overlaps are 2-on-1, some are 1-on-1; but the genome is almost certainly overannotated by at least 260 genes.) I trained my program on the codon usage table of R. mucilaginosa M508 (which contains fewer overlaps than DY-18), then tallied codon and anticodon scores on all of DY-18's CDS genes. In the end, 276 genes gave scores indicative of a reversed reading frame. Of those, 265 were, in fact, involved in overlaps.

Codon scoring is such an effective method, I don't know why programs like Glimmer don't use it. It's quite obvious they're not using it, though, because every genome has reverse-annotated genes (by the hundreds, in some cases) that are easily detected using this simple method.

Here, for the record, are the Caulobacter crescentus CB15 genes that appear to be annotated on the wrong strand:

CC_0023
CC_0048
CC_0073
CC_0099
CC_0149
CC_0354
CC_0480
CC_0546
CC_0564
CC_0605
CC_0662
CC_0666
CC_0676
CC_0677
CC_0680
CC_0681
CC_0687
CC_0728
CC_0739
CC_0775
CC_0782
CC_0786
CC_0825
CC_0850
CC_0853
CC_0913
CC_0987
CC_0996
CC_0997
CC_1020
CC_1022
CC_1031
CC_1032
CC_1050
CC_1069
CC_1073
CC_1084
CC_1094
CC_1123
CC_1127
CC_1161
CC_1174
CC_1212
CC_1222
CC_1238
CC_1245
CC_1274
CC_1312
CC_1322
CC_1340
CC_1349
CC_1392
CC_1393
CC_1394
CC_1395
CC_1414
CC_1416
CC_1513
CC_1561
CC_1648
CC_1789
CC_1793
CC_2000
CC_2086
CC_2116
CC_2163
CC_2184
CC_2193
CC_2240
CC_2256
CC_2308
CC_2334
CC_2338
CC_2351
CC_2376
CC_2413
CC_2424
CC_2442
CC_2445
CC_2450
CC_2452
CC_2471
CC_2475
CC_2499
CC_2519
CC_2525
CC_2571
CC_2574
CC_2597
CC_2602
CC_2621
CC_2624
CC_2665
CC_2698
CC_2705
CC_2718
CC_2719
CC_2720
CC_2731
CC_2732
CC_2738
CC_2739
CC_2756
CC_2769
CC_2800
CC_2850
CC_2865
CC_2875
CC_2878
CC_2907
CC_2916
CC_2949
CC_3050
CC_3055
CC_3251
CC_3302
CC_3318
CC_3342
CC_3360
CC_3429
CC_3437
CC_3438
CC_3451
CC_3453
CC_3463
CC_3479
CC_3517
CC_3519
CC_3547
CC_3548
CC_3553
CC_3554
CC_3608
CC_3665
CC_3671
CC_3700

Biology's Dirty Little Secret

Biology has a dirty little secret. It's a well-known secret among those who deal with sequenced-genome data intensively, but I suspect many non-biologists are unaware of the problem, which is: Much of the existing genome data (for sequenced genomes, ranging from bacteria to human DNA) is either corrupt or misannotated.

"Junk DNA" probably doesn't exist in living cells. But it certainly exists in published genomes.

A substantial portion of published genome data is suspect, at this point, either because of contamination issues, technical problems surrounding DNA sequencing technology, or faulty gene annotation. An example is the Oryza sativa indica (rice) genome, which inexplicably contains at least 10% of the genome of the bacterium Acidovorax citrulli. There's also a Culex (mosquito) genome with a complete copy of Wolbachia embedded. The genome of Rothia mucilaginosa DY-18 contains over 300 genes incorrectly annotated in antisense orientation (as does the genome of Burkholderia pseudomallei strain 1710b, a truly execrable train-wreck of a genome).

Another example of a genome gone wrong (arguably) is that of the bacterium Ktedonobacter racemifer, which is filled with forward and backward copies of transposases. Incredibly, one in 13 Ktedonobacter genes is a transposase, integrase, or resolvase (and that's not counting the many "hypothetical proteins" with "transposase-like" mentioned in the gene ontology notes). Disregarding the 40% of that organism's genes that are marked as hypothetical proteins, one can say that in Ktedonobacter, one in four genes of known function is a transposase, integrase, or resolvase. (Some of the organism's 4000+ "hypothetical proteins" are actually transposases incorrectly annotated in an antisense orientation.) Common sense says something's amiss.

The misannotation problem is getting worse over time. This graph (from Schnoes et al., 2009, PLoS ONE) depicts the number of sequences (left y-axis, bar graph) found to be correctly annotated in green. Sequences found to be misannotated are shown in red. The bars for each year represent only the sequences deposited into the database in that year. The fraction (right y-axis, line plot) of sequences deposited each year into the Genbank NR database that were misannotated is given by the open nodes, connected by the black line to aid in visualizing the overall trend.

The "dark matter" problem in microbial genetics is widespread and openly acknowledged. At least 20% 28.3% (according to the Joint Genome Institute) of bacterial genes are annotated as "hypothetical protein," and most of these are so annotated because they have no sequence similarity match to any known protein. In many cases, there's no match because many of the sequences are in the wrong reading frame, or have an improperly located start codon (or other serious issues). When Ely and Scott (PLoS ONE, 2014) manually reannotated the genome of the bacterium Caulobacter crescentus, they identified 11 new genes, modified the start site of 113 genes, changed the reading frame of 38 genes, and found that 112 "hypothetical proteins" were actually non-coding DNA (not genes at all). A recent transcriptome analysis of the archaeon Sulfolobus solfataricus resulted in correction of 162 gene annotations and the addition of 80 new open reading frames. But these numbers barely hint at the extent of gene misannotation. In examining the Gene Ontology database (GOSeqLite), Jones et al. found:

Annotations made without use of sequence similarity based methods (non-ISS) had an estimated error rate of between 13% and 18%. Annotations made with the use of sequence similarity methodology (ISS) had an estimated error rate of 49%.

Surprisingly, the use of sequence similarity as a guide to function identification is less reliable than non-SS methods. This is no doubt partly a reflection of the fact that gene databases contain  a great deal of aberrant data. Gene-annotation programs like the widely used Glimmer (Gene Locator and Interpolated Markov Modeler) have to be trained, using a training set. If the training set contains faulty data, it's a classic GIGO situation.

The annotation accuracy problem is getting worse by the year (see graph above). Devos and Valencia estimated in 2001 that misannotation levels could be as high as 37%. More recently, Schnoes et al. (2009) concluded that "function prediction error (i.e., misannotation) is a serious problem in all but the manually curated database Swiss-Prot," and yet Artamonova et al. (2005) found that for five types of annotation entries, even the vaunted UniProt/SwissProt database had an error rate between 33% and 43%. So even the best manually curated database is full of errors.

I've spent many hours examining bacterial genomes and it's been my experience that annotation quality is uniformly poor for all but the best-known genes in the best-curated genomes. By "best-known genes," I mean things like ribosomal protein genes, genes for well-known polymerases and chaperones, well-studied metabolic-pathway genes, and so on. The vast majority of genes encode proteins that have not been studied (and may never be studied) experimentally. All of these are suspect and need to be treated with caution, particularly in high-GC genomes where programs like Glimmer frequently can't distinguish between sense and antisense strands.

As an example: Pseudomonas aeruginosa MPAO1 has a gene, O1Q_25367, encoding a phenol hydroxylase, which is an enzyme required for catabolic breakdown of phenols (the kind of thing many Pseudomonas species are good at). If you run a BLAST search of this gene against the UnitProt.org database, you'll find dozens of good-quality hits against phenol hydroxylases of many organisms (including Rhizobium sp. Pop5, Natronolimnobius innermongolicus JCM 12255, Rhodococcus wratislaviensis IFP 2016, and quite a few others). "Good-quality" here means more than 50% amino acid sequence identities, with E-values of ~10^-50. But there's a problem. If you take the DNA sequence for the Pseudomonas phenol hydroxylase and reverse-translate it using the online app at http://web.expasy.org/translate, the resulting protein sequence is a 99% or better match (E=0) for over a hundred Pseudomonas maltose/mannitol transporters. (It's a 100% match in 85 cases.) In other words, the so-called phenol hydroxylase gene is not a phenol hydroxylase gene at all. It's a sugar transporter, backwards. You can do the same trick with the "phenol hydroxylase" of Gordonia terrae C6. Reverse-translate it and it's a much better sugar permease (dozens of E=0 matches) than it is a phenol hydroxylase.

This view of a segment of P. aeruginosa and Streptomyces genomes shows a region of 60% homology (pink band) between the two genomes for two genes. The yellow gene is each case is a "phenol hydroxylase."

How does a program like Glimmer not catch things like this? In this case, Glimmer became confused by an apparent overlapping-gene situation (above). The program found two open reading frames, in the same part of the chromosome, but on opposite strands. Glimmer 3 is supposedly much better at resolving overlaps, but Glimmer 2 (which was used to annotate around half the genomes currently available in public databases) relied on unbelievably crude heuristics to resolve overlap problems. In the above case (perhaps through operator error; these things are configurable, to a degree) Glimmer simply designated each ORF as a legitimate gene, when in fact a BLAST check leaves little doubt the phenol hydroxylase gene is (in reality) a backwards maltose/mannose transporter gene.

My experience has been that Glimmer is very easily confused by high-GC genome data, with organisms like Pseudomonas (typically 63% to 67% GC) showing a greater number of reverse-annotated genes (and no-BLAST-matches "hypothetical proteins") than low-GC organisms like Buchnera, although it should be noted that even E. coli strains tend to contain many suspect annotations. The problem is partly due to the fact that high-GC-content DNA contains relatively few stop codons in the six different reading frames, compared to low-GC DNA (which is rich in stop codons if deciphered the wrong way). Also a problem for Glimmer is the fact that in high-GC genomes, codons tend to contain a purine in the first base position whether they're read forward or in reverse.

I've written before about the fact that codons tend to occur with frequencies roughly equal to that of their reverse-complement twins. This is another confounding factor for Glimmer (codons look the same whether read in the forward direction or the reverse direction), although honestly, I'm beginning to wonder if the codon/anticodon symmetries I've been seeing aren't simply due to widespread reverse annotation of genes (misannotation of a gene's anticoding stand, as with "phenol hydroxylase").

One might naively suppose that it shouldn't be hard to discriminate the "sense" strand accurately, given that so many genes begin with a Shine Dalgarno sequence ahead of the start codon. But in reality, it turns out that not very many organisms make extensive use of SD sequences. Short motifs like GAGG, which occur randomly at a high rate, can be mistaken for a Shine Dalgarno sequence. This only aggravates the false-positive rate.

Erroneous gene annotations are rampant in bacteria, but some authors have suggested that the problem is far worse in eukaryotic genomes. If that's true, we're in trouble.

All of this puts bioinformatics research at a crisis point. Gene discovery algorithms are good (maybe a bit too good) at finding open reading frames but poor at identifying and assigning gene functions, in part because we lack the kind of in-depth understanding of protein folding required for prediction of 3-dimensional structures (and active-site conformations) programmatically, a capability that's sorely needed if we're to progress out of the gene-annotation Stone Age we're now living in. (Protein folding is a Hard Problem requiring supercomputers to untangle.) Maybe in ten years (or twenty?), we'll be able to predict 3D protein structures computationally, and on that basis make better ab initio predictions of protein function. Right now, we have to make do with relatively crude Markov-model pattern recognition software, aided by human intervention, to come up with even a minimally reliable genome annotation. But we have the means, already, of doing much better crosschecks. Some of that can be automated. We just need to have the will to do it.

Are Overapping Genes Real?

Bacteria belonging to the Pseudomonas family are a perennial favorite among bacteriology instructors (and students) because of the curious ability of some of its members to produce pigments that fluoresce under an ultraviolet light. If you're unlucky enough to get an infected cut on the arm while working in the garden, it's possible your cut will fluoresce under a black light. That's enough of a diagnosis to pronounce the infectious agent. 

Fluorescent colonies of Pseudomonas.

Silby and Levy, investigating the adaptation of the bacterium Pseudomonas fluorescens to soil, uncovered the existence of at least ten antisense genes in P. fluorescens. They went on to demonstrate experimentally that one of the genes, cosA, produces not just antisense RNA but an associated protein. Tellingly, Silby and Levy commented:

These findings suggest that current genome annotations provide an incomplete view of the genetic potential of a given organism.

The implication is that additional antitranscriptome genes remain to be found, not only in Pseudomonas but in other organisms.

There's a good reason they haven't been found yet. Overlapping genes are automatically rejected by many of the annotation programs that are commonly used to find, identify, and label genes in genome sequences. (The oft-used freeware Glimmer 2 program allows you to set the overlap-rejection threshold.) Many yet-to-be-discovered antisense genes have been deliberately and systematically obscured in published genomes.

Still, once in a while such genes do surface. For example, in Pseudomonas stutzeri A1501, we find a pair of overlapping genes at an offset of 3035137 on the chromosome (see illustration below).

Overapping genes in Pseudomonas stutzeri.

The top gene is annotated merely as a "hypothetical protein," while the underlying gene on the opposite strand is an aspartyl-tRNA synthetase. One's normal inclination is to dismiss a hypothetical protein as being unimportant, but this may not be wise. Twenty percent or more of bacterial genes are annotated as hypothetical proteins; common sense says they can't all be unimportant. In fact, in "Transcriptome Analysis of Pseudomonas syringae Identifies New Genes, Noncoding RNAs, and Antisense Activity" by Filiatrault et al. (2010), researchers found that 818 out of 1,646 protein genes in P. syringae annotated as "hypothetical proteins" were expressed under iron-limited conditions. Many (probably most) genes annotated as "hypothetical protein" are quite real and should probably be re-annotated as PUF: "protein of unknown function."

In this case, the "hypothetical protein" shown in yellow (above) turns up medium-strength protein-BLAST hits with other "hypothetical proteins" from other organisms, including a hit with an E-value of 3.0×10^-49 in Parasutterella excrementihominis YIT 11859 and a comparable hit on a predicted phosphatase/phosphohexomutase in Rothia mucilaginosa DY-18.

In this particular case, the hypothetical-protein gene lacks a strong upstream Shine Dalgarno sequence (a sequence preceding many genes that helps bind a ribosome to the mRNA). But so too does the gene on the opposite strand. (This is not unusual. The SD sequence is not required for translation and in fact, in about half of bacterial species, a Shine Dalgarno sequence is associated with fewer than 50% of genes.) Hence, the jury's out on whether the antigene is expressed. It could be that no protein is made from the top strand but the gene provides RNA-mediated control of the gene on the bottom strand. We won't know for sure until someone investigates.

In Pseudomonas aeruginosa strain PADK2_CF510, we find another instance of a bidirectional overlapping gene pair (see graphic below). In this case, the gene on the top strand (CF510_06030) encodes the large subunit of an isopropylmalate isomerase. The gene on the bottom strand (CF510_06025, shown in yellow) is annotated as "Flp pilus assembly protein TadG." It could very well be a misannotated non-gene. However, five genes away is FimV (CF510_06060), another pilus-assembly (motility) protein. Moreover, the gene marked TadG has a strong upstream SD sequence containing the canonical GGAGG motif. The gene above it has a weaker GGAAA motif.

P. aeruginosa has an overlap of an isopropylmalate isomerase gene and a gene for a motility protein. The latter is shown in yellow.

In previous posts, I've mentioned (and shown data for) the fact that in the overwhelming majority of protein-encoding genes (across every kind of genome), the first base of a codon tends to be purine-rich. One check of whether a bidi-overlap gene is "real" or not ought to be that the first codon base should be purine rich in both reading directions. This is, in fact, the case for the examples shown above. The aspartyl-tRNA synthetase gene for P. stutzeri has AG1 (1st base, purine) content averaging 59.8%, whereas its bidirectional partner gene ("hypothetical protein") has AG1 = 58.5%. The isopropylmalate isomerase of P. aeruginosa has AG1 = 65.9%, while its antisymmetric partner (TadG) has AG1 = 56.2%.

If you enjoyed this post, please give the URL to your biogeek friends. Thanks!