Predicting antibiotic resistance from a genome

Campylobacter has 1.6 million bases of DNA. Change one of them, a single C to a T, and ciprofloxacin stops working. One letter out of 1.6 million, and that is not a figure of speech.

Fluoroquinolones, the family ciprofloxacin belongs to, kill bacteria by jamming the enzymes that untangle DNA during replication. E. coli has two of them, DNA gyrase and topoisomerase IV. Mutate one and the drug still has the other, which is why resistance in E. coli usually needs hits in both. Campylobacter jejuni has no topoisomerase IV. Repeated attempts to find parC have failed, and the genome does not carry it. One target, and no backup.

A single base change at position 257 of gyrA, turning threonine 86 into isoleucine, raises the ciprofloxacin MIC 128-fold. The two numbers name one change: 257 counts bases in the gene, and 86 counts amino acids in the protein it encodes, which forms part of DNA gyrase.

Appearing is one thing. Lasting is another. Most human Campylobacter infection comes from poultry, so what matters is whether the mutation survives inside a chicken, with no drug anywhere near it. It does. Resistant and susceptible strains colonised chickens equally well with no antibiotic present. Nothing pushes the mutation back out, and once it appears, it stays.

So the answer is in the genome, in a known gene at a known position. The clinic still does it the slow way: grow the organism, expose it to the drug, wait. Campylobacter is microaerophilic and slow, so that takes days, and the patient is treated on a guess in the meantime. Sequencing takes a day. Read position 257 and you are done.

Except that you are not. Take two isolates, each a culture of the bacterium grown from one sample, that both carry the change at 257. One dies at a concentration of ciprofloxacin that the other survives. Campylobacter has a pump, CmeABC, that pushes the drug back out of the cell before it ever reaches the gyrase. Disable the pump and the resistance collapses, with the mutation still intact.

Position 257 is not a switch. It is a contribution, and how much it contributes depends on what else the genome is carrying. That is the question this chapter is built on. If resistance is written in the sequence, why can nobody simply read it off? Half of the answer is in the biology. The other half, as it turned out, is in the machine learning.

1.1 Features, labels and a floor

The models in this chapter run off a table, with one row per isolate and one column per resistance gene, holding a 1 where the gene is present and a 0 where it is not. A known resistance mutation, such as threonine to isoleucine at position 86 of gyrA, gets a column of its own and counts among the resistance genes. What a model looks at are its features, and here the features are those columns. What it predicts is its label: for each isolate, susceptible or resistant, as the laboratory measured it. Every isolate arrives with its answer attached, and the model has to learn from those answers how to give the right one for a genome it has not seen. That is supervised learning. Figure 1.1 sets it beside the other possibility, unsupervised learning, where no answers are given and the job is to find the groups. Nothing in this chapter asks for that.

The first model is a rule, with nothing fitted to the data. It asks one question of each genome. For ciprofloxacin, the biology supplies the question: is there a substitution at position 86 of gyrA? It asks nothing else: no pump, no context, no second gene.

The second model is a logistic regression. I gave it all 84 resistance genes and let it weight them freely. For each genome, the model adds the weights of the genes it carries to a starting value shared by every genome, turns the total into a probability of resistance, and calls the isolate resistant when that probability is above one half. A gene with a large weight is one the model leans on, and a gene with a weight near zero is one it barely looks at.

Before either model is believed, it needs something to beat. Most of the isolates in this chapter are susceptible to ciprofloxacin, so a model that reads no DNA at all and guesses susceptible every time is right 79.32% of the time. That is the floor. Find the floor before you believe any number, including mine.

1.2 The rules

Both rules were fixed before any modelling, and each was changed once after I had looked at the data. Each changed version sits beside its original, with the reason for the change, and the sections that follow explain each one.

1.3 One position against eighty-four genes

I took 3,984 Campylobacter genomes from NCBI. Every one had a ciprofloxacin result measured on a plate. The one question about position 86 was right for all but thirteen of them, which is 99.67%. That beats the floor by twenty points. The logistic regression, with all 84 genes to draw on, did not beat the one question. Its score is cross-validated, averaged over repeated divisions of the data and each time measured only on isolates held back from its training. It was right 99.60% of the time, and 99.57% when whole families of related genomes were held back together. The rule has nothing to fit, so it scores 99.67% either way.

The pump is real, and I modelled none of it. A single question was still right 3,971 times. The reason has nothing to do with the bacterium. It is about what the laboratory measured.

What the laboratory measured was a number. The organism goes into a row of wells, each holding twice as much ciprofloxacin as the one before it, and the lowest concentration that stops it growing is the result. That number is continuous. Some isolates stop at a low dose, and some need far more. Then somebody draws a line across it, susceptible below, resistant above.

A good pump buys the organism another tube or two. If the isolate stays on the same side of the line, the pump did something and the label recorded nothing. So I did not solve the problem I opened with, which was how much a change contributes. I answered a newer one: which side of the line an isolate falls on.

That lesson arrives before any algorithm does. The label is a decision somebody made, usually before you arrived, often for convenience. Choose the label and you have fixed the ceiling. Every model after that, however clever, works underneath it.

The weights showed the next problem. I read them from the same kind of model fitted to a narrower table, the 58 genes found in at least three isolates, using the training isolates of a split that kept families of related genomes together. The chapter comes back to both choices. That model gave the isoleucine change at position 86 a weight of 8.845. It gave a real gyrase mutation, the valine change at the same position, 2.364. It also gave a beta-lactamase 2.364: the same number, to three decimals, for a gene that does nothing to fluoroquinolones.

1.4 A gene that cannot touch the drug

A beta-lactamase cuts open a beta-lactam ring. Penicillin has one, and so does amoxicillin. Ciprofloxacin has none, so there is nothing in it for the enzyme to cut.

The gene is blaOXA-493. Five isolates carry it, and all five are resistant to ciprofloxacin, so the model paid it. Now look at those five again. Every one also carries a gyrase mutation, the valine change at position 86, with no exceptions.

Put yourself where the model sits. You see two columns. Wherever the gene is present, the mutation is too, and the isolate is resistant. Which one is responsible? Nothing in the table can tell you. Trained with families kept together, the 58-gene model split the credit, 2.364 each, because it had no other option. That is confounding, the oldest problem in statistics wearing a genomic coat (Figure 1.2).

The part that should worry you comes next. Suppose a sixth isolate turns up carrying blaOXA-493 and no gyrase mutation. The model counts the gene towards resistance and raises the probability it gives that isolate, for a reason that has nothing to do with the drug. There is no such isolate among these 3,984. The mistake is already in the model, and nothing in the data will show it to you.

The model did not learn what the gene does. It learned what the gene travels with. Why the two travel together I cannot tell from five isolates. Shared ancestry, perhaps, or nothing at all. The model cannot tell either, because it only ever saw them arrive as a pair.

This is the limit of the whole method. A model sees co-occurrence, and that is all it ever sees. From inside a table, a cause and its companion look identical. What separates them is knowing what the gene does, and that comes from biology, from somebody who knows that ciprofloxacin has no beta-lactam ring.

1.5 Hiding whole families

If the model is reading company rather than mechanism, there is a way to catch it: hide whole families of isolates from the model, then see whether it still works.

To see why that helps, start with how a model gets tested at all. You hold data back. Train on three quarters of the isolates, then score the model on the quarter it has never seen. If it performs well on that quarter, it has learned something. That only works if the held-back quarter is genuinely new to it.

These isolates are not strangers to one another. NCBI sorts them into families called SNP clusters, a SNP being a difference at a single base. Two genomes in the same cluster differ by a handful of bases out of 1.6 million: the same organism sampled twice, often from two patients in one outbreak. Of my 3,984 isolates, 3,520 fall into 1,128 clusters, and the other 464 have no cluster assignment. Split them at random and 344 of those clusters end up with isolates on both sides. The model learns from one genome and is then tested on its near-identical twin. It scores well and has learned nothing, because it is recognising a genome it has already seen.

So I split again, keeping every cluster whole, with no family on both sides. The 344 became zero (Figure 1.3). Listing 1.1 shows the split.

    # Grouped split. Isolates with no cluster assignment are treated as
    # their own singleton group, which is the conservative reading: we
    # cannot show they are related to anything, so we do not assume it.
    groups = y["snp_cluster"].fillna(
        pd.Series([f"__unclustered_{i}" for i in idx], index=y.index)
    )

    gss = GroupShuffleSplit(n_splits=1, test_size=TEST_FRACTION, random_state=SEED)
    train_g, test_g = next(gss.split(idx, y["resistant"], groups))
    y["split_grouped"] = "train"
    y.loc[test_g, "split_grouped"] = "test"

A single division of the data can be accidentally kind, so I repeated the comparison twenty-five times, dividing the data differently each time. Averaged over those twenty-five, the one question about position 86, asked of the 58-gene table, was right for 0.9965 of the held-back isolates with families allowed on both sides and for 0.9965 with families kept whole, and the 58-gene logistic regression was right for 0.9961 and 0.9960. The rule’s number did not move, and the regression’s barely did.

That is good news. If the regression had been reading family resemblance, it would have failed on families it had never seen. It did not fail. It predicted resistance in families it had never encountered, and that is what learning a mechanism looks like.

The beta-lactamase is still in there, though. Three of the eight valine isolates carry no blaOXA-493. The random split left valine isolates without the gene on the training side, so the 58-gene model could give the valine change more of the credit: 3.689, against 1.469 for the beta-lactamase. The grouped split, which gave the weights quoted earlier, held all three back for testing, so the model never saw the two apart and paid them 2.364 each. Splitting by family did not remove the gene’s weight, and it could not. Every isolate that carries blaOXA-493 also carries the gyrA mutation, so however the data is divided, the two are still together on whichever side they fall.

So there are two problems here, and they are different problems. Leakage is about how you divide the data, and dividing it better fixes it. Confounding is about what is in the data, and no division fixes it at all.

1.6 The gene was there; the protein was not

A different drug exposed a different weakness. Tetracycline sits down on the ribosome and blocks the seat where each new amino acid has to arrive. Protein synthesis stops, and the cell cannot build anything. Campylobacter’s answer is a gene called tet(O). It makes a protein that binds the ribosome and dislodges the drug.

So the rule is easy. If the genome carries a tet gene, call the isolate resistant. In these genomes, that means tet(O), except in five isolates that carry one of two mosaic genes instead, tet(O/M/O) in four and tet(O/W/32/O) in one, all five of them resistant. Across 3,983 isolates, every one with a tetracycline result measured on a plate, that rule is right 99.32% of the time. A rule that asked about tet(O) alone would miss those five and score 99.20%. But fourteen isolates carry tet(O) and die of tetracycline.

Carrying a gene is one thing. Making a working protein is another. Part of the explanation lay in how the gene had been written down. AMRFinderPlus is the software that reads a genome and reports which resistance genes are in it, and that report is the annotation. It works by comparing the sequence against a reference copy of each known gene, and it does not answer yes or no. When the match covers the whole reference gene, it writes tet(O). When the match covers only part of it, it writes tet(O)=PARTIAL. My rule had been reading both as present.

Of these isolates, 1,890 carry a full-length copy, and 99.6% of them survive tetracycline. Four carry a fragment and nothing else, with no full copy anywhere in the genome. All four died. The protein, Tet(O), has to reach the ribosome and take hold of it. Part of a Tet(O) reaches nothing.

There is a second kind of partial. Fifteen isolates carry tet(O)=PARTIAL_END_OF_CONTIG, and 93% of those are resistant. That flag means something entirely different. Sequencing does not hand you a genome in one piece. It hands you short reads, which software assembles into longer fragments called contigs, and in these isolates the gene ran off the edge of one. The sequence stops because the assembly stops, not because anything is missing from the bacterium. One word, partial, was covering two situations. One concerns the organism, and the other concerns how its genome was pieced together.

Separate them, and 99.32% becomes 99.45%. That is five isolates out of 3,983: the four that carry only a fragment, and one whose only copies are a fragment and a copy broken by an internal stop or a frameshift. I ran this comparison twenty-five times as well. A rule has nothing to fit, so the repeats cannot raise or lower its average. What they show is how far the difference between the two rules moves from one division of the data to the next, and whether a gain of five isolates is larger than that. It is about the same size. The gain is 1.1 times that spread, and anything below about 1 cannot be told apart from noise. The mechanism says a truncated protein cannot work. The data agrees, quietly, and cannot prove it, so the change goes into the repository as a hypothesis. Nine of the fourteen are still unexplained: an intact gene, susceptible on the plate, and no reason I can give.

The lesson underneath does not depend on whether that difference is significant. My table already kept tet(O) and tet(O)=PARTIAL in separate columns. The rule merged them. It counted every tet entry as present, so it read a fragment of tet(O) the same as a whole copy. The annotation had been telling me more than that, and I threw it away before the rule ever ran.

1.7 A threshold chosen in seconds

Back to ciprofloxacin, and to the table. Before any model uses it, you have to decide which columns to keep. Here is why you would drop any. Suppose a gene turns up in two isolates out of 3,984, and both happen to be resistant. In the table, that column looks perfect. It is also worthless. Two isolates tell you nothing about a gene, and the model cannot know that.

So the standard move is to delete rare features. Pick a minimum count and cut everything below it. In an early version of the pipeline, I picked ten. Anything found in fewer than ten isolates was gone before the model ran.

Now look at what that removed. The change this chapter opened with turns threonine at position 86 of gyrA into isoleucine. There is a second change at the same position, threonine to valine. Eight isolates in this cohort carry it, and every one of them is resistant. Those eight sit in seven different SNP clusters, counting an isolate with no cluster as a cluster of its own, as the split does. So they are seven independent observations rather than one outbreak counted repeatedly. Eight is fewer than ten.

With the threshold at ten, the rule asking about position 86 could only find the isoleucine version, and it was right 99.55% of the time. Lower the threshold to five, so that the valine version survives the cut, and the same rule is right 99.75% of the time. The gap is eight isolates, exactly the eight carrying the change I had deleted (Figure 1.4). Both thresholds also cut the two rarest changes at that position, which come later in the chapter, so neither figure is the 99.67% that the rule scores when it counts every change.

Nothing warned me. There was no error, no missing file and no collapse in accuracy. Both numbers sit above 99%, and both look like a model that works. Ten was a number I chose in a few seconds because it sounded reasonable, and I did not test it. It removed the best-evidenced resistance mutation in the dataset after the common one.

I have tested it since. The pipeline now scores the rule and the logistic regression with the threshold at 1, 3, 5, 10 and 20, and at no setting does the regression beat the rule. The default is three, the smallest count that removes a two-isolate column like the one above. That is why the valine change, in eight isolates, and blaOXA-493, in five, both carry weights in the model I read them from. Three is a choice too, and it still cuts the rarest change at position 86, which a single isolate carries.

So the model I trained was the algorithm plus every decision I made before it started: which columns to keep, what counts as one sample, what the label is allowed to mean. None of those show up in an accuracy figure. All of them decide what that figure is measuring.

1.8 Four changes at one position

Position 86 is the amino acid in DNA gyrase that sits exactly where ciprofloxacin has to bind. Change it and the drug loses its grip. Four different changes turn up there in these 3,984 genomes. Threonine normally sits at 86. In 810 isolates it has become isoleucine, in eight valine and in one lysine. That is 819 isolates, and 817 of them survive ciprofloxacin.

The fourth change is alanine. Four isolates carry it, and none of the four is resistant.

The explanation came from a second drug, nalidixic acid, the original quinolone, in clinical use since the 1960s. Ciprofloxacin came later. It is a fluoroquinolone, a far more potent version of the same idea, aimed at the same enzyme. The NCBI release that supplied these genomes also holds nalidixic acid results for 3,816 sequenced isolates, the four alanine isolates among them. All four are resistant to nalidixic acid. Not one of them is resistant to ciprofloxacin.

The other three substitutions defeat both drugs. Against nalidixic acid, counting only isolates with a result for it, the count for isoleucine is 772 resistant out of 775, for valine eight out of eight and for lysine one out of one. Alanine defeats only the old drug (Figure 1.5).

Look at what each change does to the enzyme. Isoleucine and valine are bulky. Lysine is long and carries a charge. Alanine is the smallest substitution available, barely a change at all. It is enough to shake off nalidixic acid and not enough to shake off ciprofloxacin.

None of this is new. Jesse and colleagues reported exactly this in 2006, in isolates that resisted high concentrations of nalidixic acid and were killed by low concentrations of ciprofloxacin. My four isolates reproduced a twenty-year-old result without being asked to.

Now the part that concerns the model. My feature was any substitution at position 86 of gyrA. I chose it from the mechanism, and the opening of this chapter argued for it. That feature covers four substitutions, and one of them points the other way. Drop alanine from the rule and it goes from 99.67% correct to 99.77%.

That is four isolates. Three of the four sit in the same SNP cluster, so it is closer to two independent observations than four. The improvement is a fraction of the variation between cross-validation folds, the repeated divisions of the data into training and test sets. The biology is published. The statistics here are silent. Both of those are true at once.

The discussion of the label, earlier in the chapter, was about the resolution of the answer, where susceptible or resistant is a line drawn across a continuous measurement. This is about the resolution of the question. How finely you define a feature is a choice, and taking it from the mechanism does not protect you. The mechanism at position 86 turned out to be four mechanisms.

1.9 Changing the rule after looking

I found the four alanine isolates by reading through the data, and then changed the rule to account for them. That sequence is how results go wrong.

The number went up because I moved the goalposts after watching where the ball landed. Do that often enough and you can improve any model to any level you like. Every dataset contains patterns that mean nothing, and a rule shaped around them will not survive contact with the next dataset.

So the safe answer would be to throw the improvement away and keep the original number. That answer is wrong too. Refusing to learn anything from data you have already seen looks like rigour. It is superstition. The real test is where the justification for a change comes from.

I have two such reasons, both set out above. The nalidixic acid result came from a different experiment on the same four bacteria, not from the accuracy I was trying to raise. In the 2006 paper, Jesse and colleagues found the same thing in bacteria I have never seen. Neither of those could have come out of tuning.

So both numbers go into the repository with their labels attached, one rule fixed before any modelling and the other changed after looking (Figure 1.6). Nobody has to guess which is which.

Looking at your data is not the sin. Pretending you did not is.

1.10 The list behind the verdict

Two genes in this dataset look identical from the numbers. The first appears in five isolates, all five resistant to ciprofloxacin, in four separate SNP clusters, and never once without a mutation at position 86 of gyrA, the change that actually causes the resistance. The second appears in six isolates, all six resistant, in five separate clusters, and never once without a mutation at that position either. One of them has nothing to do with ciprofloxacin at all. The other pumps ciprofloxacin out of the cell. Nothing in the table tells you which is which.

The first gene is blaOXA-493, the beta-lactamase from earlier in the chapter.

So I wrote code to find the other passenger genes, the ones that travel with the cause and do nothing themselves. The code starts from the rule this chapter has been testing, a change at position 86 of gyrA. It then looks for genes that appear only in isolates that rule already calls resistant, and that have no known way of acting on the drug. The first half of that is arithmetic. The second half is a list of gene names I typed into the script: gyrA, parC, tet, 23S, erm. Five entries, covering the mechanisms I had in mind that afternoon (Listing 1.2).

# Prefixes with an accepted mechanistic link to some drug in this chapter.
# Anything outside this set that tracks the rule perfectly is a passenger.
MECHANISTIC = ("gyrA_", "parC_", "tet(", "23S_", "erm(")

The second gene is cmeB. It is part of CmeABC, the pump Campylobacter uses to throw fluoroquinolones back out of the cell before they reach their target. The annotation reported cmeB in only six of the 3,984 genomes, the six counted above. cmeB was not on my list. So the code did exactly what it was told, and accused the pump (Figure 1.7).

Look at the two genes again. Statistically they are the same gene. What separates them is knowing what the protein does, and here that knowledge was in a list. Mine had five entries and was missing one.

Being more careful next time would not fix this. Every method that sorts real causes from things that merely travel beside them uses knowledge from outside the dataset: a pathway database, a curated catalogue, a prior, somebody’s judgement. That knowledge is part of the method, and its gaps do not announce themselves.

1.11 What survives

I set out to ask whether a genome can tell you what a laboratory plate would. For ciprofloxacin it can. One feature is right 99.67% of the time across 3,984 isolates, and eighty-four features do no better, at 99.60% under cross-validation. That part took an afternoon. The rest of the work went on finding out what the number does not mean (Figure 1.8).

The label made the problem easy, because the efflux pump moves an isolate a tube or two along the measurement without carrying it across the line, and so never reaches the answer. The regression on the 58 genes found in at least three isolates paid a gene that cannot touch the drug, because the gene never appears without a real gyrase mutation. Hiding whole families did not fix it. With families kept together, the regression paid the two the same. A rarity threshold I chose in seconds, in an early version of the pipeline, deleted a resistance mutation, and no accuracy figure noticed. The feature I built from the mechanism turned out to be four mechanisms. The code I wrote to find passenger genes named the efflux pump, because the list it checks against was mine, and incomplete.

What survives all of that is narrow, and worth stating carefully. In Campylobacter, ciprofloxacin resistance is explained well enough by one position in one gene that nothing I built improved on it. That is a claim about this organism and this drug, and not about resistance prediction in general.

Everything is in the repository: thirteen scripts, one pinned NCBI release, the metrics file behind every result in this chapter, and a limitations section naming what the data could not establish. Nine tetracycline isolates in there carry an intact resistance gene and die anyway. I do not know why. Chapter 2 needs a harder problem.

This chapter is rewritten from the posts for Days 2 to 11 of the series.

References

Jesse, T. W., Englen, M. D., Pittenger-Alley, L. G. and Fedorka-Cray, P. J. (2006). Two distinct mutations in gyrA lead to ciprofloxacin and nalidixic acid resistance in Campylobacter coli and Campylobacter jejuni isolated from chickens and beef cattle. Journal of Applied Microbiology, 100(4), 682-688. https://doi.org/10.1111/j.1365-2672.2005.02796.x