Is a gene a fact or a call?

In 2000, a team sequenced Campylobacter jejuni and reported 1,654 genes. In 2007, a team that included two of the original authors went back to the same genome and reported fewer. There was no new strain and no new sequencing. The same 1,641,481 letters were read again, gene by gene, by curators. The DNA did not change. The count did (Figure 3.1).

Somebody printed a microarray with 1,731 spots on it, one for every gene first identified in the genome: real glass, real DNA, an experiment you could run. By the time the genome paper was published, 77 of those spots stood for genes that were no longer in the annotation. The count had been revised down to 1,654, and the array had been built from the longer list. The 77 spots stayed on the glass, and the gene lists in the analysis software flagged each one as a “not annotated gene”. Seven years later the count went to 1,643. From 1,731 to 1,654 to 1,643 is eighty-eight fewer genes in the list.

The same thing still happens every day. Bakta and Prokka are two standard annotation tools, both widely used. Give them the same bacterial genome and, at some places, one says there is a gene here and the other says there is nothing here. Same DNA, the same letters in the same order, read twice. Neither program is broken.

This chapter asks whether that disagreement can be predicted from the sequence alone. It does not ask whether a gene is really there. Nothing in this project can tell you that. The answer being predicted, which machine learning calls the label, was produced by software, so a model trained on it learns the annotator rather than the organism. That is the harder thing about this chapter, and it is the point.

3.1 Where genes come from

Start with what a gene caller is looking for. A protein-coding gene begins with a start codon, runs in groups of three letters and ends at a stop codon. The stretch in between is an open reading frame. Find every open reading frame and you have found every possible gene. That is the easy part, and on its own it is useless.

Three of the sixty-four possible codons mean stop. Campylobacter is AT-rich: most of its letters are A or T, and the three stop codons are written mostly in those two letters, so in its genome they are common. Read a random stretch and you hit a stop every thirteen codons or so. So the genome is full of short open reading frames that mean nothing: millions of letters, six reading frames, stops everywhere. The real genes average 316 codons. They are the long ones. But length alone will not separate them, because a short gene is still a gene and a long accident is still an accident.

So the caller weighs other things. Does the codon usage inside this frame look like the codon usage of known genes in this organism? Is there a ribosome-binding site sitting just upstream? Does the protein it would produce resemble a protein somebody has already described? Each of these is evidence. None of them is proof. And each needs a cut-off: how similar is similar enough, how long is long enough, how close the binding site has to sit. Nobody finds those numbers in the genome. They are chosen by the people who write the software, and two programs can choose them differently. That is the whole of this chapter.

Of the C. jejuni genome, 94.3% codes for protein. The other 5.7% is about 93,564 letters, and it is neither empty nor junk. Some of it is genes that make RNA instead of protein: the ribosomal RNA the cell builds its ribosomes from, and the transfer RNAs that carry amino acids into place. Some of it is control: promoters where transcription starts, ribosome-binding sites, terminators that say stop here. Besides ribosomal and transfer RNA, this genome has four annotated RNAs. Four, in 1.6 million letters. That is a count of how many anybody has written down, not a claim that only four exist.

Which brings me to the word intergenic. It sounds like a description of a place, between the genes. What it actually means is that nothing has been annotated here. It is a residual, defined by absence, and the absence is in the record rather than in the molecule. Move a gene boundary and the intergenic space moves. Withdraw a gene and the intergenic space grows. When two annotation tools disagree about a region, one of them calls it a gene and the other calls it intergenic, and neither is looking at anything the other cannot see.

3.2 Why a random forest

The model in this chapter is a random forest: many decision trees, each a chain of yes-or-no questions learned from a different slice of the data, voting on the answer.

Chapter 2 was built on a single tree, and the best moment in it was reading the thing. I could print the questions it had chosen, and it had found a real binding rule that nobody told it about. A forest gives that up. Three hundred trees will score higher than one, and I will not be able to read any of them. So this chapter has a harder label and a model I cannot read.

3.3 What the two tools actually do

Bakta and Prokka do not disagree the way I assumed. I assumed two independent gene callers, each trained differently, each finding its own answer. That is not what is happening underneath.

Prokka calls genes with a program called Prodigal. Bakta calls genes with a program called Pyrodigal, and Pyrodigal is a rewrite of Prodigal in a different language: same model, same scoring, faster code, some bugs fixed. The two are close enough to the same caller that most of a genome will come out identical either way.

The disagreement comes from two decisions about what to look for. Bakta adds a search for very short proteins, under thirty amino acids, which Prodigal was never built to find. Bakta also calls pseudogenes: genes broken by a stop codon or a frameshift, still recognisable, still worth recording. Prokka mostly does not. So a real disagreement between them will rarely look like two tools reading the same DNA and reaching different conclusions. It will look like one tool asking a question the other tool never asked.

That changes what the forest should find. If the two callers were independent, disagreement might trace to something subtle in the sequence, a signal one caller catches and the other misses. If they mostly share a caller, disagreement should trace to something blunter. Length, mostly, because short proteins are what one tool looks for and the other does not. If the forest’s top feature, the column of data it leans on most, is not length, I have misunderstood something.

3.4 The rules

The rule for what counts as a disagreement was written into the code before the label was built. Each protein-coding gene Bakta called is one row. If anything Prokka annotated on the same strand overlaps that call, by as little as a single base, the call counts as matched. If nothing does, Bakta made a call there and Prokka did not, and the row is a disagreement.

One base sounds loose, and it is meant to be. The two tools mostly agree on where a gene starts and stops, and when they do not, it is usually over the exact start codon: a difference of a few dozen letters, inside a call both tools clearly made. That is not a disagreement. Both tools found the gene. A rule that demanded a large overlap on both sides would turn some of those start-codon differences into missing genes. So the label asks only whether Prokka claimed the stretch at all. Whether the two calls share a start, and whether they share a stop, are recorded beside the label rather than mixed into it (Figure 3.2).

Strand matters. DNA has two strands, and a gene can sit on either. Two genes on opposite strands can occupy the same letters and still be different genes, so only what Prokka annotated on the same strand is consulted. If the only call Prokka made at that place is on the other strand, the Bakta call counts as one Prokka did not make. The kind of annotation does not matter. If Prokka called a transfer RNA, or anything else that is not a protein-coding gene, on the same strand where Bakta called a protein-coding gene, Prokka did not leave the region empty. It made a different claim, not a missing one, and counting it as silence would put two unlike things in one label. Those calls count as matched, and there are 29 of them.

One more decision was easy to get wrong. When both tools call the same gene but pick different start codons, which coordinates go in the table? Every row carries Bakta’s, matched or not: where Bakta’s call starts and stops, and how long it is. How much of it Prokka’s call covers is kept as a separate count. So the length of every call in the table is the length Bakta chose, and a model that leans on length can learn that tool’s habits instead of the sequence.

The rule did not change once the label existed, and that matters more than it looks. A threshold chosen after seeing which one gives the better story is not a threshold. It is the story wearing a number.

Later in the chapter the question moved from where genes are to what they are called, and that comparison needed a rule of its own. Comparing product names, the names each tool gives to the protein a gene makes, by exact string match is too strict. One tool writes 50S ribosomal protein L7/L12, and the other writes the same thing with a bracket on the end. Those are the same claim, and counting them as different would make the whole exercise meaningless. So names are normalised before comparison. Strip a trailing bracket. Strip a trailing gene symbol, because one tool often appends it and the other does not. Strip a leading hedge word: putative, probable, possible, predicted or conserved. Then ignore capitals and treat punctuation as a space. Under that rule, putative cytosol aminopeptidase and Cytosol aminopeptidase are the same name. The check on that rule was fixed before I ran it.

3.5 Seventy-two rows

The pipeline ran on 25 bacterial genomes, one from each of 25 species spread across eight phyla. The Campylobacter jejuni genome this chapter opened with is one of them. Across all 25, Bakta made 87,960 gene calls, and 72 of them have no matching call from Prokka. That is a disagreement rate of 0.0008, eight in ten thousand. Of the 87,888 calls both tools made, 87,788 have identical start and stop coordinates, to the base.

The reason was knowable before I downloaded a single genome. With Prokka on Prodigal and Bakta on a rewrite of Prodigal, asking whether Bakta and Prokka disagree about where genes are is mostly asking whether Prodigal agrees with itself. It does. Leave aside the 43 calls Bakta made with its search for very short proteins, and in the other 87,917, there are 38 divergences, which is 0.04%. I expected a small signal. I got 72 rows.

How much overlap the rule demanded hardly mattered. When 87,788 of the 87,888 matched calls are identical to the base, they would match under any threshold, and the 72 with nothing overlapping them would match under none. A stricter threshold could have changed the label for at most the remaining hundred of the 87,960 calls. I could not have known that in advance, and choosing the threshold afterwards, once I could see which one flattered the result, would not have been a choice at all. A rule that turns out not to bind is still worth fixing beforehand.

The reasoning about length held. Asked to pick out the 72 coordinate disagreements, the forest leaned on the length Bakta gave each call more than on anything else. Of those 72, 34 come from one module. After the main gene caller finishes, Bakta runs a second search for short open reading frames, under thirty amino acids. Prokka has nothing equivalent. Bakta made 43 such calls across the whole panel, and 34 of them are disagreements. Nearly every short protein Bakta finds, Prokka misses, because it never looked. Prodigal’s default minimum gene length is 90 base pairs, which is thirty codons. The shortest disagreement here is 45 base pairs.

Now the part that is actually biology. Look at what those short proteins are. They are Type I toxin-antitoxin system Ibs family toxins, and leader peptides from the tryptophan and histidine operons. These are real, well described and functionally important. The first is a toxin whose antitoxin is an RNA that stops it being translated. The second is a short peptide whose ribosome stalls or does not, and in stalling decides whether the operon behind it gets transcribed. This is textbook regulatory biology, and exactly the class a dedicated short open reading frame finder is built to catch.

So the disagreement is real, and it points at something real. But the label is not what I said it was. At the start of this chapter, I described two tools reading the same DNA and reaching different conclusions. What the data holds is one tool running an extra module. That is a fact about software architecture, and the biology is downstream of it.

3.6 Two failures that exited cleanly

One more thing had gone wrong before I could trust any of this. Prokka found zero ribosomal RNAs in all 25 genomes. It reported that, and exited successfully.

Every bacterial genome has ribosomal RNA. It is not optional, and a genome with none is not a genome. The tool did not crash, did not warn and did not return an error code. It printed Found 0 rRNAs and finished. I found it while the annotation was still being set up, before any label had been built.

Here is the chain. Prokka does not find ribosomal RNA itself. It calls a program called barrnap. Barrnap does not do it either: it calls nhmmer. The nhmmer binary in that environment had been compiled against one version of a maths library, and the environment shipped a different version. So nhmmer failed at runtime, every time. Barrnap got nothing back and reported nothing. Prokka reported what barrnap told it. Three programs deep, and nobody raised a hand (Figure 3.3).

Think about what that would have done to the comparison. A disagreement is a region one tool called and the other did not. If Prokka finds no ribosomal RNA anywhere, then every ribosomal RNA Bakta calls looks like one. Not a real one: a missing shared library, which would have gone into the data as biology.

I caught it because the count was zero rather than low. Zero is impossible, and impossible is visible. The first complete Prokka run, 87,901 gene calls and no ribosomal RNA, was thrown away and redone. The rule I took from it is that a program exiting successfully is not evidence that it did its job. The pipeline now records how many of each kind of annotation each tool produced, per genome. I do not expect this bug again. I expect the next one to exit zero as well.

With the coordinates settled, the same regions could be asked a different question. At 87,859 places, both tools called a protein-coding gene, almost always with the same start and stop. Do they agree on what those genes are? The second failure belongs to that question, and this time the code was mine.

My code decided that GTPase Era and GTPase Era were different names, 23 times (Figure 3.4). Same string: same letters, same spaces, same capitals, scored as a disagreement between two tools that had written exactly the same thing. The way it happened is common.

The gene-symbol strip is where it broke. To strip a trailing gene symbol, you need to know what the gene symbols are, and my first version read each record’s own symbol field and stripped that. Bakta fills that field on 19,216 of these regions. Prokka fills it on 48,201. So for thousands of pairs the strip fired on one side and not the other. Two identical strings went in, one came out shortened, and the comparison said different.

The fix is one line of thinking rather than one line of code. Build the symbol list once per pair, from both tools, and apply it to both sides. Listing 3.1 is the comparison as it now stands.

def compare(product_bakta, gene_bakta, product_prokka, gene_prokka, level="primary"):
    """-> (agree, normalised_bakta, normalised_prokka)."""
    symbols = (frozenset() if level == "strict"
               else symbol_candidates(product_bakta, gene_bakta,
                                      product_prokka, gene_prokka))
    a = normalise(product_bakta, symbols, level)
    b = normalise(product_prokka, symbols, level)
    return a == b, a, b

Then I added an assertion, because a fix you cannot verify is a hope. If two raw strings are identical, they can never be scored as a disagreement. The pipeline checks this on every run and refuses to write the file if it fails. Across all 87,859 regions both tools called, 29,808 pairs have identical raw strings. There are zero violations.

Note the direction of the bug. It inflated the disagreement rate. It made my result more interesting.

3.7 Half the names

Bakta and Prokka read the same 50,210 genes. They give different names to 25,977 of them. That is 51.7%, about half. These are the same regions, with both tools confident enough to write a name, and the names do not match.

The obvious objection is that this is a formatting artefact, and the check fixed in advance answers it. Compare the strings as written, ignoring only capitals, spacing and a full stop at the end, and 56.0% disagree. Compare them under my normalisation rule and 51.7% disagree. Compare them under the loosest rule I was willing to defend, dropping generic words and ignoring word order, and 50.6% disagree. Five points across the entire range from strictest to loosest. The disagreement is not punctuation.

A second objection is worth answering. Bakta has a naming style that guarantees a mismatch. When it can only match a structural domain, it writes a name ending in domain-containing protein. Prokka does not do that, so those are near-certain disagreements. There are 4,210 such regions, and 284 more where Bakta calls the protein uncharacterised and adds only the name of its open reading frame, which Prokka never matches either. Together that is 4,494 regions, 9% of the set, disagreeing almost every time. Take them out and the remaining 45,716 regions disagree 47.0% of the time.

The shorthand disagrees too. Where both tools assign a gene symbol, the symbols differ 41.1% of the time. Where both assign an enzyme number, the numbers differ 16.7% of the time.

Pangenome analyses group genes by name. Functional enrichment counts categories built from names. Resistance reports quote a gene product. Every one of those inherits whichever tool was run. None of this is a claim that either tool is wrong. I have no way to know which name is right, and neither does anybody working from annotation output alone.

3.8 A model that reads nothing

Can the naming disagreement be predicted? The first answer was a warning. A model that reads nothing won.

It uses no sequence, no database matches and no information at all. It always guesses the same answer, and it got the highest score in the whole comparison. Here is how a model that knows nothing beats models that know something. Half the regions in this data are disagreements between the two tools, so a model that always says yes, this is a disagreement, is right about half the time by luck alone. The score it won on was F1, which Chapter 2 also used. It rewards catching every true case, and it gives no credit for a region correctly left alone. Say yes to everything and you catch every true case. So the model that never looks at the data scored 0.664 on F1, higher than any model that actually learned something.

This had happened before, in Chapters 1 and 2. A score means nothing until you know what a model that knows nothing would score.

So the honest comparison uses a different measure, MCC, with every model and baseline scored on the same regions, all from five genomes held out of training. The model that always says yes scores 0.000, as any model that gives the same answer every time must. A model that looks only at how long the protein is scores 0.028. One that looks only at whether Prokka found a similar protein in UniProtKB, a large public protein database, scores 0.141. It stands in for the baseline I had planned, which asked whether either tool had left the protein unnamed: every region here was named by both tools, so that question has the same answer everywhere. A single decision tree scores 0.261. The random forest scores 0.303. Gradient boosting, a more advanced technique, scores 0.288. The most elaborate model on that list did not win; the random forest did (Figure 3.5).

What does the random forest get right and wrong? It was tested on 12,283 regions it had never seen, and it flagged 6,105 of them as disagreements. 3,965 of those flags were correct. 2,140 were not. So for every three regions it flags, one is a false alarm, against one in two if it flagged at random. That is real signal. It is nowhere near reliable enough to trust on its own.

3.9 What the forest was reading

I ran the same random forest three times, with the same folds, the same random seed and the same way of choosing its settings. Only the columns it was allowed to see changed (Figure 3.6). The seventeen columns that describe the raw DNA include the GC content of the gene and of the hundred letters on either side, how repetitive the sequence is, how often stop codons fall in its emptiest reading frame, and three properties of the whole genome. Thirty-five come from the gene caller: the boundaries it chose, the reading frame and strand they imply, and everything computed from the protein those boundaries translate into, such as its length and the share of each amino acid. The other nine describe what each tool found in its databases.

With all 61 features, the forest scores 0.303 on MCC on the five held-out genomes. Take away everything that comes from the gene caller itself, leaving 26 features, and it scores 0.259. It barely moved. Take away the database features too, leaving the 17 that describe the raw DNA, and it scores 0.142.

That last number is generous. A test on five genomes is small, and its score moves with which five genomes happen to be in it. The steadier measure comes from the twenty training genomes: split them into five groups, hold each group out in turn while the forest trains on the other four, and average the five scores. On that measure the DNA-only forest scores 0.070, which is close to nothing. That is the honest answer to the question this chapter actually asked. Does the DNA sequence predict disagreement? Barely.

So which features had been carrying the score? There is a way to check this honestly. Shuffle one column, so that it no longer lines up with the right answers, and score again to see how much was lost. Do it on genomes the model never saw, so the answer is about generalisation. Then add up the damage by category. Shuffling the nine features that describe database cross-references cost 0.165. Shuffling the thirty-five features that describe the gene call itself cost 0.019. Shuffling the seventeen features that describe raw DNA cost negative 0.001. Negative means that shuffling those columns made the model slightly better. That is what noise looks like when you measure it honestly.

The top three features were how many cross-references Bakta attached, whether Prokka assigned an enzyme number, and whether Prokka matched a protein motif.

The obvious reading is that the names disagree where the evidence runs thin: neither database knows the protein well, each tool guesses, and the guesses differ. I tested that rather than assuming it, and it holds only in part. It fits the 4,494 regions where Bakta fell back on a domain name or an open reading frame name, which disagree almost every time. Among the other 45,716, it holds on Prokka’s side. Where Prokka assigned no enzyme number, the names disagree 54.0% of the time, against 42.1% where it did.

On Bakta’s side it runs the other way. Where Bakta attached two cross-references or fewer, as it does about nine times in ten, the names disagree 44.2% of the time. Where it attached three, they disagree 79.8% of the time, and at four or more, 62.8%. A third cross-reference is usually an enzyme number or a match in a more specialised database, and it comes with a more specific name from Bakta, which is then more likely to differ from Prokka’s more general one. That is why the number of cross-references Bakta attached is the forest’s strongest feature: on Bakta’s side, more evidence goes with more disagreement. Either way, what the model learned is a fact about two reference databases and how each tool names from them. It has nothing to do with the DNA sequence itself.

One decision made this visible. Every feature was given flags at creation saying where it came from, and the check reads those flags automatically rather than from a list I typed by hand. Listing 3.2 shows how the database features were declared.

        "bakta_has_uniref": (
            0, 1, "Bakta's Dbxref contains a UniRef cluster"),
        "bakta_has_pfam": (
            0, 1, "Bakta's Dbxref contains a Pfam accession"),
        "bakta_has_ec": (
            0, 1, "Bakta's Dbxref contains an EC number"),
        "bakta_has_is": (
            0, 1, "Bakta's Dbxref contains an ISFinder insertion-sequence match"),
        "bakta_n_dbxref": (
            0, 1, "number of cross-references Bakta attached"),
        "prokka_has_similarity_hit": (
            0, 1, "Prokka's inference cites a UniProtKB protein similarity"),
        "prokka_has_protein_motif": (
            0, 1, "Prokka's inference cites an HMM protein motif"),
        "prokka_has_cog": (
            0, 1, "Prokka assigned a COG"),
        "prokka_has_ec": (
            0, 1, "Prokka assigned an EC number"),
        "mobile_neighbours_5kb": (
            0, 1, f"Bakta CDS within {MOBILE_WINDOW} bp carrying an ISFinder "
                  "match; the only mobile-element context derivable from what "
                  "is already annotated"),

Without those flags, the 26 features left once the gene-caller columns were gone would have passed for features of the DNA, and the score of 0.259 would have looked like real biology. It was not.

3.10 What this chapter got wrong

I built an experiment around a disagreement that turned out to be 72 rows. The reason was sitting in the tools’ own documentation, and I had not read it first.

Bakta ran on db-light, the smaller of its two reference databases, because the full one would not fit on my machine. I predicted that the smaller database would leave Bakta naming fewer genes than Prokka. I had it backwards. Of the 87,859 regions both tools called, Bakta left 12.5% with a placeholder such as hypothetical protein instead of a name. Prokka left 40.2%.

My own code counted identical text as a disagreement, 23 separate times, and that mistake made my result look better than it was.

One feature in the first experiment, on where genes are, never changed value across its 87,960 rows, so it could tell the model nothing. Nobody caught it until the second experiment, on names.

3.11 Fact or call

The chapter asked whether a gene is a fact or a call, and the answer split in two.

Where does a gene start? Bakta called 87,960 genes, and the two tools disagreed about only 72 of them. That looks remarkable until you see why. Both tools find genes with the same underlying program, so this was never two independent opinions. It was one opinion, copied. Where a bacterial gene starts looks like a fact, but it is one decision, made once and inherited by everyone who has used the same software since. A fact that spread is still a copy.

What is that gene, once you have found it? Here the story turns over. For the same 50,210 genes, both tools were confident enough to give a name, and the names disagree 51.7% of the time. Without one tool’s habit of naming genes after a generic protein domain, or after an open reading frame alone, it is still 47.0%. A model can predict which genes will disagree, but not from the DNA. Give the model only the features that describe the DNA, and the score falls to 0.142 on the five held-out genomes, and to 0.070 in cross-validation on the twenty training genomes. What carries the prediction is what each tool found in its reference library. Where Prokka found less, the names disagree more often. Where Bakta found more, it gives a more specific name, which is more likely to differ from Prokka’s. The DNA never changed. The books each tool checked did.

So a gene is both. Where it starts is a call so old and so widely copied that it behaves like a fact. What it does is a call made fresh, by whoever has the newer library, and it disagrees over half the time. Every number from the two experiments, right and wrong, is in the repository. Chapter 4 asks a harder kind of question, not whether something happens but when, and some of its patients are still alive when the study ends.

This chapter is rewritten from the posts for Days 26 to 39 of the series.