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

The GGAGG Reflex

A common myth in biology is that genes coding for proteins need to have a Shine Dalgarno sequence upstream of the start codon. Students sometimes spout this as an inarguable fact; a kind of molecular biology catechism. I call it the GGAGG reflex.

In fact, no SD sequence is required. None at all. It's important to be clear on this.

In case you're not a biogeek: In the 1970s,  Australian scientists John Shine and Lynn Dalgarno were the first to notice that the tail end of the 16S bacterial ribosomal RNA contains a short nucleotide sequence whose reverse complement is often found immediately upstream of a protein gene's start codon. The exact sequence varies from organism to organism, but the rRNA trailer sequence is usually pyrimidine-rich. In E. coli, the sequence is CACCTCCTTA. (Here, I am of course talking about the DNA sequence. In RNA it's CUCCUCCUUA.)  If you reverse the sequence, the Watson-Crick complement is TAAGGAGGTG. Some portion of the latter is often found a few nucleotides upstream of a start codon; not 100% of the time, but too often to be by chance.

The key intuition here is that Watson-Crick binding of the tail end of the 16S rRNA to the corresponding antisequence ahead of the start codon helps stabilize the ribosome so that it is more likely to translate the gene. The degree of binding depends, of course, on the fidelity of the SD sequence ahead of the gene. Usually, the purine-rich SD area is not an exact match for the 16S rRNA trailer, and in fact the SD region quite often has no detectable SD signature whatsoever.

How often is "quite often"? In 2002, Ma et al. undertook a survey of 30 organisms representing bacteria from all major taxonomic groups. Somewhat surprisingly, they found that in 17 out of 30 organisms, a Shine Dalgarno sequence was present at fewer than half of all CDS (protein-encoding) genes. Among the bacteria most likely to use SD sequences were Bacillus subtilis and Thermotoga thermophilus, in which 90% of known protein genes have an upstream SD signal. Among those least likely to use SD sequences were low-GC/small-genome organisms (intracellular parasites, Mycoplasmas, and pathogens), with many groups, like the Actinobacteria (47%), falling somewhere in the middle.

Before taking these findings to heart, though, it's worth noting some serious weaknesses in the Ma et al. study. In obtaining the above numbers, Ma et al. used a rather permissive definition of "SD sequence," based on a minimum binding-energy cutoff (∆G) of -4.4 kcal/mol, which means they counted GAGG as a SD sequence (and also GGAG and AGGA). If one were to count only GGAGG (length 5) and longer motifs, the percentages given by Ma et al. would be much lower. (I present some data of my own on this further below.) The reason this is a very serious issue is that the probability of random occurrence of short (length-4) sequences like GAGG is substantial. Ma et al. failed to report the expectation odds for the various "signals" they looked for. Hence, for short motifs, we have no way of knowing, for the various organisms, what the expected rate of occurrence of short signals was. If an organism with genomic GC content of 66% has a putative SD motif of GAGG, AGGA, or GAGG in the 20-bp target region for 20% of its genes, how does that compare with the random occurrence rate for those sequences, given the organism's DNA base composition? We're not told.

Bearing in mind the weaknesses of the study, a number of nonethelesss interesting findings came out of the Ma et al. survey, including:

• A SD sequence is rarely long or canonical; many times it's just GAGG or GGAG or AGGA (putatively) or a corruption of the expected form (e.g., GGTGG instead of GGAGG)
• SD sequences occur more often with highly expressed genes (such as genes for ribosomal proteins and core energy metabolism genes) than with low-expression genes
• In some (not all) organisms, the SD sequence is more likely to occur in conjunction with an ATG start codon and less likely to occur with GTG or TTG
• Vanishingly few SD signals are located further than 14 bases or closer than 4 bases away from a start codon

Anybody with modest JavaScript skills can write scripts that verify some of these findings against public genomes. I took a quick look at the genome for Rothia mucilaginosa DY-18 (a member of the Actinobacteria family and a common inhabitant of the human mouth). First, I determined the most likely SD sequence for Rothia based on the 16S rRNA trailer of CCTCCTTTCT (implying a SD sequence of AGAAAGGAGG), then I had my script scan the genome in both directions, looking for any of the six possible length-5 motifs within the full-length sequence (so, AGAAA, GAAAG, etc.), in the 20 base pairs upstream of every annotated open reading frame start codon. In total, I found 686 putative SD sequences within 4 to 14 bases of an annotated start codon. Since Rothia mucilaginosa has 1905 CDS genes, this means 36.0% of protein genes carry a putative length-5 Shine Dalgarno signal. When I re-ran the check using all possible length-4 SD signal variants (using the relaxed criteria of Ma et al.), I found 1160 positives. Thus, 60.9% of CDS genes in R. mucilaginosa have a length-4 SD signal per Ma et al.

On a probability of abundance basis (given Rothia's actual base composition stats), we would expect to see 203 length-5 SD motifs by pure chance in the genome's 1905 20-bp regions. The actual number (686) is obviously quite a bit higher than expected, tending to validate the notion that these are, indeed, SD motifs we're looking at. For length-6 motifs, the trend is even sharper: The expectation is 40 occurrences by chance; the actual number is 340. So at a length of 6, a motif has high odds (~90% chance) of being real.

By contrast, the statistical expectation for length-4 motifs calculates out at 957, which is only slightly less than the number found (1160). Therefore, in dealing with Shine Dalgarno sequences, at least in Rothia, it's meaningful to deal with length-5 and longer motifs, but probably not meaningful to deal with length-4. When you spot a length-4 motif, odds are very high you're looking at a randomly occurring pattern.

If you enjoyed this post, please share the URL with a friend. Thank you!

Coming to Grips with Pseudogenes

The term pseudogene was coined in 1977, when Jacq et al. discovered a version of the gene coding for 5S rRNA in the African clawed frog (Xenopus laevis) that was truncated yet retained homology with the active gene. Subsequent work has shown that in higher life forms, pseudogenes (genes that have been inactivated through one event or another) are almost as numerous as coding genes, with (for example) the human genome containing 10,000 or more pseudogenes. (A more recent estimate puts the number at 20,000.) Many of these pseudogenes are highly conserved. Looking at pseudogenes in the mouse and human, Svensson et al. found that of a group of 74 such genes that occur in both species, 30 appear to have been conserved since before the evolutionary divergence of mice and humans.

In higher organisms, pseudogenes are sometimes transcribed into RNA, with the RNA filling a regulatory function. For example, Korneev et al. found that simultaneous transcription of neural nitric oxide synthase (nNOS) and the antisense strand of a homologous pseudogene in the same neurons of Lymnaea stagnalis (a snail) leads to the formation of a duplex between the two strands and a reduction in nNOS translation. Further examples can be found in Pink et al. (2007), "Pseudogenes: Pseudo-functional or key regulators in health and disease?"

In bacteria, pseudogenes are somewhat rarer than in eukaryotes, but exist in significant numbers in many pathogens (including many species of Mycobacterium, Shigella, Brucella, Bordetella, and others). A study by Kuo and Ochman (2010) found that pseudogenes are swiftly eliminated from Salmonella. They describe "evidence of a strong deletional bias in Salmonella, such that genes that are not maintained by selection are rapidly inactivated and eliminated by mutational events." In fact, Kuo and Ochman found that pseudogenes are eliminated more rapidly than could be explained by the so-called neutral theory of evolution, indicating that the continued presence of pseudogenes exacts a high cost to the cell.

And yet, many bacteria with slow-evolving genomes (such as Mycobacterium species) retain their pseudogenes with high fidelity across evolutionary timespans. The most celebrated "pseudogene hoarder" of all time, M. leprae (the leprosy bacterium) appears to have acquired its 1000+ pseudogenes 9 to 20 million years ago. Meanwhile, the half-life of pseudogenes in Buchnera aphidicola was measured at 23.9 million years—a staggering number.

So on the one hand, we have work by Kuo and Ochman showing that pseudogenes in bacteria are rapidly eliminated, and on the other hand we have some bacterial lineages in which it seems pseudogenes are not only conserved but actively repaired over periods of tens of millions of years!

In Chapter 5 of Brucella: Molecular Microbiology and Genomics (2012, Caister Academic Press), Garcia-Lobo et al. describe their work with RNA sequence data from the bacterium Brucella abortus:

Twenty-four of the genes selected from the RNAseq data were annotated as pseudogenes in the B. abortus 2308 genome, which was considered a rather unexpected finding. By comparison with other Brucella genomes we can reduce the list of highly expressed pseudogenes to 16 (often, truncated parts of a gene are annotated as different pseudogenes especially in B. abortus 2308). This seems contradictory since high transcription of these genes, which should be not able to translate into functional proteins, will be contrary to biological economy. The high levels of transcription observed for these genes strongly suggest that they could be active genes and their products may perform functions unreported in metabolic reconstructions. High pseudogene expression may also indicate that these are very recently produced pseudogenes that did not turned down transcription yet by accumulation of mutations in their promoter or control regions. It is also possible that these pseudogenes may contain sequencing errors and they are indeed active genes.

It's almost comically obvious from this passage that the authors are troubled by their own finding that some pseudogenes in Brucella are highly transcribed. They try explaining it away by saying it could all be "sequencing errors."

A more parsimonious view is that pseudogenes that haven't been eliminated from a genome are, in fact permanent, legitimate fixtures of the landscape, in microbes just as in higher life forms. And as in higher life forms, pseudogenes in microbes are probably serving perfectly understandable regulatory functions (when they're not actually translated into protein products).

Kuo and Ochman have convincingly shown that useless pseudogenes are quickly eliminated. It follows that any pseudogenes that aren't swiftly eliminated are, in fact, serving a biological purpose, or else they wouldn't be there. This line of reasoning is already well accepted by researchers who study eukaryotic life forms. Those who study bacteria need to take a hint from their up-the-food-chain colleagues.

What could the hundreds of pseudogenes in Bordetella pertussis (or the 1000+ pseudogenes in M. leprae) be doing? First we need to get used to the idea that in bacteria, virtually all genes are transcribed, in both directions. It's been four years since Dornernberg et al. reported finding ~1000 antisense transcripts in E. coli, but no one seems to have gotten the memo.

A section of Rothia mucilaginosa genome (top) and a corresponding portion of Mycobacterium leprae (bottom); click to enlarge. The yellow gene, in each case, is DnaE (error-prone polymerase). Pink bands indicate areas of 65% or more homology between the two organisms. The small-diameter silver genes in the lower panel are M. leprae pseudogenes. "Normal genes" are shown in green. Notice that R. mucilaginosa has open reading frames on both strands of DNA, with many bidirectionally overlapping genes.

A look at the genome of the bacterium Rothia mucilaginosa DY18 shows that a very large proportion of "normal genes" have open reading frames on the opposite strand (see illustration). Bidirectional overlapping genes run throughout the Rothia genome. A massive annotation error? Maybe. Or maybe both strands are transcribed.

If massive wholesale transcription of antisense strands occurs in E. coli, as we know it does, certainly it's no stretch to imagine it occurring in Rothia mucilaginosa. And if it is occurring in Rothia, which is (incidentally) an opportunistic pathogen, how much harder can it be to imagine it occurring in another well-known pathogenic member of the Actinomycetales family, Mycobacterium leprae? We know already that upwards of 40% of M. leprae pseudogenes are transcribed. Antisense transcripts could well be playing a role in silencing certain gene essential genes when attempts are made to grow the organism in defined media. Forward transcripts could be producing nonsense or partial-nonsense/truncated proteins that are excreted as toxins or find their way to the cell wall as surface antigens. Any number of scenarios might be possible.

Some very low-hanging fruit is available to micobiologists who are willing to accept the obvious. Instead of wishing away pseudogenes or imagining them to be useless baggage, we should be looking at them as potential determinants of pathogenicity. We should consider their possible roles in modulating protein expression patterns. We should attempt to learn why they're conserved; what role(s) they're playing in cell physiology. The last thing in the world we should be doing is calling them "junk DNA."

An Example of Antisense Proteogenesis?

The question of how organisms develop entirely new genes is one of the most important open questions in biology. One possibility is that new genes often develop through accidental translation of antisense strands of DNA.

An example of this can be seen with the S1 protein of the 30S bacterial ribosome. If you take the amino-acid sequence for an S1 gene and use it as the query sequence in a blast-p (protein blast), you'll mostly get back hits on other S1 proteins, but you'll also get minor (low-fidelity) hits on polynucleotide phosphorylase. Why? When you do a blast search, the search engine, by default, looks at both DNA strands of target genes (sense and antisense strands) to see if there's a potential sequence match with the query. If there's a match on the antisense strand, it will be reported along with "sense" matches. In the case of the S1 protein, blast-p searches often report weak antisense hits on polynucleotide phosphorylase in addition to strong sense hits on ribosomal S1.

Ribosomal proteins are, of course, among the most highly conserved proteins in nature. It turns out that polynucleotide phosphorylase (PNPase) is very highly conserved as well. It's an enzyme that occurs in every life form (bacteria, fungi, plants, animals), absent only in a scant handful of microbial endosymbionts that have lost the majority of their genes through deletions. While the chemical function of PNPase is well understood (it catalyzes the interconversion of nucleoside diphosphates to RNA), its physiologic purpose is not well understood, although recent research shows that PNPase-knockout mutants of E. coli exhibit lower mutation rates. (Hence, PNPase may actually be involved in generating mutations.)

The bacterium Rothia mucilaginosa, strain DY18, has a (putative) PNPase gene at a genome offset of 1277514. When this gene is used as the query for a blast-p search, the hits that come back include many strong matches for the S1 ribosomal proteins of various organisms. By "strong match," I mean better than 80% sequence identity coupled with an E-value (expectation value) of zero. (Recall that the E-value represents the approximate odds of the match in question happening due to random chance.

If we use the Genome Viewer at genomevolution.org to look at the PNPase gene of Rothia mucilaginosa, we see something extraordinarily peculiar (look carefully at the graphic below). Click to enlarge the following image, or better yet, to see this genome view for yourself, go to this link.

Notice the presence of overlapping sense and antisense open reading frames on a portion of DNA from Rothia mucilaginosa. The top reading frame contains the gene for polynucleotide phosphorylase. The lower (-1 strand) reading frame contains ribosomal S1. To see this in your own browser, go to this link.

Notice that there are overlapping genes. On the top strand is the gene for PNPase; on the bottom strand, in the same location, is a gene for ribosomal S1. These are bidirectionally overlapping open reading frames, something occasionally encountered in virus nucleic acids but rarely seen in bacterial or other genomes.

How do we explain this anomaly? It could be just that: an anomaly, two open reading frames that happen to overlap (but that aren't necessarily translated in vivo). Or it could be that at some point, many millions of years ago, the ribosomal S1 gene of a Rothia ancestor was erroneously translated via the antisense strand, producing a protein with PNPase characteristics. We don't know why PNPase confers survival value (its physiologic purpose is not fully understood), but we do know, with a fair degree of certainty, that PNPase does, in fact, confer survival value—because every organism, at every level of the tree of life, has at least one copy of PNPase. Once Rothia's ancestor, through whatever process, opened up a reading frame on the antisense strand of ribosomal S1, the reading frame stayed open, because it conferred survival value. In this way, the first Rothia PNPase was born. (Arguably.)

At some point in its history, Rothia duplicated its PNPase gene and placed a new copy at genome offset 1650959. Over time, this second copy diverged from the original copy, becoming more like E. coli PNPase (which is also to say, less S1-like). Rothia's second PNPase shows a blast-p similarity of 45% (in terms of AA identities) to E. coli PNPase, with E-value 4.0e-147. It shows a blast-p similarity of 26% (AA identities) with E. coli ribosomal S1 (E-value: 4.0e-17). Neither E. coli PNPase nor Rothia PNPase-2 overlaps an S1 gene. However, both are colocated with the ribosomal S15 protein gene. And you'll find (if you look at lots of bacterial genomes) that PNPase is almost always located immediately next to an S15 ribosomal gene.

Rothia PNPase is an example of an enzyme that may very well have started out as an antisense copy of another protein (the S1 ribosomal protein). Of course, the mere presence of bidirectionally overlapping open reading frames doesn't prove that both frames are actually transcribed and translated in vivo. But the fact that blast-p searches using PNPase as the query almost always turn up faint S1 echoes (in a wide variety of organisms) is highly suggestive of an ancestral relationship between the two proteins.