What does a cell show to the immune system?

There is a reason Chapter 1 stopped at logistic regression.

It was already enough. Eighty-four features did no better than one, and adding a more powerful model to a solved problem teaches nothing. So this chapter keeps the task, predicting a category, and changes the problem underneath it. The data is harder, and this time the differences between the models are real. Four models, one to a chapter: decision trees and the rules they learn, random forests and what is lost when a model can no longer be read, gradient boosting and a label that is two numbers together, and support vector machines and the geometry of a decision boundary (Figure 2.1). This chapter is the decision tree.

Right now, every cell in your body is holding up samples of what is inside it. Proteins get chopped into short fragments. Some of those fragments are carried to the cell surface and held there, where passing immune cells inspect them. If a fragment looks foreign, the cell is killed. That is how your immune system finds a virus hiding inside a cell without ever going in. It reads the samples on the outside.

The molecule doing the holding is called Major Histocompatibility Complex class I, or MHC class I. It has a groove. A fragment nine amino acids long sits in that groove like a key in a lock. Two of those nine do most of the gripping, the second and the ninth. Get both right and the fragment is displayed. Get one wrong and it falls out, the immune system never sees it, and a virus with the wrong fragments is invisible.

So here is the chapter. Give the model nine amino acids. Ask whether they bind.

One MHC molecule, HLA-A*02:01, the most studied variant in the world. The first download gave 15,597 laboratory measurements and 9,142 distinct fragments, and 53.7% of those measurements were binders, which is almost exactly a coin flip. In Chapter 1, guessing susceptible every time was right 79.32% of the time before I read a single base of DNA. This time, guessing gets me nothing. Those counts are the download as it arrived, before I checked which unit each measurement was in, and the corrected cohort, smaller and with a lower share of binders, comes later in this chapter.

One more thing, and it should feel familiar from Chapter 1. Binding strength is a number. Somebody drew a line across it and called everything on one side a binder. The database that supplies this data says, in its own documentation, that the line should probably sit somewhere different for every MHC molecule.

2.1 Why a tree

Now look at the shape of the groove’s rule.

Chapter 1 was one position. One change at gyrA 86, and ciprofloxacin stops working. One question, one answer.

This needs two positions to agree at once. And what counts as agreeable at position nine depends on which MHC molecule you have, and people carry different ones. That is a conjunction, not a threshold.

Which is exactly what a decision tree is made of. It asks one question, then asks a different question depending on the answer, then another. Follow the path and you have a chain of conditions. The groove enforces a chain of conditions. The tree is built from them (Figure 2.2).

2.2 Nine letters as a table

GILGFVFTL. Nine letters, positions 58 to 66 of the influenza A matrix protein. If you carry HLA-A*02:01, your cells display this fragment when you have flu, and it is the most studied T-cell target we have. Position two is isoleucine and position nine is leucine. Both sit in their pockets.

Now hand those nine letters to a model. It sees nine characters. No chemistry, no shape, no charge. Only characters.

Machine learning models take numbers, and a peptide sequence is text. Something has to convert one into the other, and that something is a choice.

The obvious choice is to number the alphabet. Alanine becomes 1, cysteine 2, and so on to tyrosine at 20. Nine columns. Done.

It is also wrong. A decision tree splits by asking whether a number falls below some value. Ask whether the amino acid is below ten and you have grouped alanine, aspartate, phenylalanine and lysine. They share nothing except the alphabet.

The honest version is called one-hot encoding, one column for every amino acid at every position. Nine positions and twenty amino acids give 180 columns, each holding a zero or a one. Every column is one yes-or-no question. Is position two a leucine. Is position three a valine. That is the question the groove asks.

The cost is that similarity disappears. Leucine and isoleucine are the same atoms arranged differently. Across 180 columns they are strangers, and the model has to learn each one alone.

The third option is to describe each amino acid instead of naming it. Hydrophobicity, molecular weight and charge become the columns, three numbers for every position and twenty-seven in all. Molecular weight stands in for the size of the residue, and similar residues then get similar numbers. That is fewer columns and more biology, but it is a bet placed before the model has seen anything.

All three go in the repository (Figure 2.3), and I said I would write down which one I expected to win before any of them had a score, because a choice made after seeing the scores is not really a choice. I never wrote it. The pipeline did not run for several days, and when it did, there was nothing on record saying which encoding I had backed. Listing 2.1 shows the three encodings.

def enc_ordinal(sequences):
    return np.array([[ORD.get(aa, 0) for aa in seq] for seq in sequences], dtype=float)

def enc_onehot(sequences):
    rows = []
    for seq in sequences:
        row = []
        for aa in seq:
            row += [1 if aa == a else 0 for a in AA20]
        rows.append(row)
    return np.array(rows, dtype=float)

def enc_physico(sequences):
    rows = []
    for seq in sequences:
        row = []
        for aa in seq:
            row += [KD.get(aa, 0.0), MW.get(aa, 0.0), CH.get(aa, 0.0)]
        rows.append(row)
    return np.array(rows, dtype=float)

Here is the part that took me a while to see. Chapter 1 never made this choice at all. A tool called AMRFinderPlus handed me a finished table of genes present or absent and substitutions found or not. Somebody else had decided what counted as a feature. I inherited that decision and never noticed I had.

A feature matrix is not a fact about nature. It is a claim about what matters.

2.3 How a tree chooses, and what stops it

Every textbook gives the same rule for HLA-A*02:01. Leucine or methionine at position two. GILGFVFTL has neither. It binds anyway, and it binds so well that your killer T cells go after it before almost anything else in the virus. The rule is not wrong. It is just not the whole story.

Now watch what happens when a machine tries to learn that rule from data.

A decision tree gets all 180 columns, and it may ask only one of them first. It chooses by counting.

Put every measurement in one pile. 15,597 laboratory measurements on this molecule, 8,376 counted as binders and 7,221 not. The pile is a coin flip. Now ask one question and let the pile fall into two. Ask something irrelevant and nothing happens, because both piles come out half and half and you have learned nothing. Ask about leucine at position two and one pile should come out mostly binders. The tree tries all 180 questions, measures how much each one tidies the piles, and keeps the winner. Then it starts the whole search again inside each pile.

The tree cannot look ahead. It takes whichever question pays now. A weaker first question might set up a far stronger second, and the tree never looks. That is called a greedy search. Most of the time it works, because position two really does carry signal on its own, so the tree will find it.

And then the tree’s first question does exactly what the textbook does. It sends GILGFVFTL to the wrong side.

One awkward peptide would be a footnote. Sidney and colleagues counted fourteen A*0201 epitopes in vaccinia virus and reported in 2008 that eight of them did not fit the rule (Figure 2.4).

Chapter 1 never had this problem. One question about one position in Campylobacter was correct for 99.67% of the 3,984 isolates I tested. One question, whole model, done.

Peptides do not hand you that. So the tree asks again. And again. Nothing in the arithmetic tells it when to stop.

Ask enough questions and every peptide ends up alone in its own group. A group of one is never wrong. A near-perfect score on all 15,597 measurements, and worth nothing. A group holding one peptide is not a rule. It is a memory. The tree has recorded that this exact sequence bound. It has learned nothing about a sequence it has never seen.

So what stops the tree? You do. You can cap how many questions deep it may go. You can refuse to split a group below some size. Or you can let it grow wild and cut the branches back afterwards, which is called pruning. Every one of those is a number you picked. Not one of them came out of the data.

The honest way to choose is to hide some peptides, fit the tree on the rest, then ask it about the ones it never saw, and keep the depth that does best on those. Then do it again with a different part held back, five times over, so that every peptide is held back once, and average the five scores. That is cross-validation.

2.4 The rules

Here is what was fixed, and what I changed after looking.

2.5 Where this could go wrong

Cross-validation can still be fooled, and here the data does the fooling.

There are 15,597 measurements and only 9,142 distinct peptides, so 6,455 of those measurements are a repeat reading of a peptide already in the set (Figure 2.5). Hide a peptide from the tree and a copy of it may still be sitting in the data you trained on. The tree is not being tested then. It is being asked to recall. The repair is to collapse the repeats before anything is split, one row for each distinct peptide carrying one label. Choosing that label is a problem of its own, and it comes later in this chapter.

Chapter 1 had the same problem and got away with it. 344 groups of closely related Campylobacter isolates sat on both sides of a random split, and it changed nothing: the one question about position 86 was still right 99.65% of the time. The model had found the one position that decides resistance, and there was nothing left worth memorising. A tree deep enough to memorise has plenty to memorise.

Then comes the harder version. Not copies. Relatives.

Change one letter in GILGFVFTL. Is it still the same peptide? Ask the MHC molecule holding it and the answer may be yes, because the groove grips position two and position nine, and if you change position five, the grip barely notices. Ask a T cell and the answer may be no, because position five points up, out of the groove, straight at the T-cell receptor. Same nine letters, two different answers, depending on who is asking (Figure 2.6).

Those 9,142 peptides came out of 2,874 different proteins, just over three peptides per protein. Split them at random and peptides from the same protein land on both sides of the line. The model studies one. Then it gets asked about the neighbour.

The repair is to split by group instead of by row. Every peptide from one protein goes to the same side, training or testing, never both. That is grouped splitting. One line of code, and the accuracy usually drops. That drop is the point. The old number was scoring the model on neighbours it had already studied, and the new one asks about proteins it has never met.

Grouping by source protein is a guess. Two peptides from different proteins can differ by a single letter, and two from the same protein can share almost nothing. The protein is a convenient handle rather than the thing that matters. What matters is how similar two sequences are, and the usual fix is to cluster the peptides by sequence identity and split whole clusters instead. Then you have to choose the identity threshold. Eighty per cent? Seventy? You have swapped one arbitrary line for another.

So you decide, and then you report a number that is partly a report about your decision. There is no neutral split.

2.6 A column with no unit

A 45-degree melting point is not a binding affinity.

IEDB, the database that supplies this data, stores measurements from many assay types in the same column. Radioactive inhibition assays report affinity in nanomolar. Thermal shift assays report melting temperatures in degrees. Stability assays report half-lives in minutes. Structural assays report distances in ångströms. The column is called quantitative measurement. No unit is attached.

A 500-nanomolar threshold on that column does not select binders. It selects numbers below 500.

Of the 15,597 measurements, 3,874 were not nanomolar, and of those, 3,418 fell below 500 and were labelled binders. Among those 3,874 were 2,225 half-lives, 1,087 unitless figures, 333 melting temperatures and 219 interatomic distances. A half-life of 23 minutes is a stable complex, not a tight binder, and it was sitting in the training data as one.

The corrected cohort is 11,723 measurements, 8,346 distinct peptides, and 4,958 that bound below 500 nanomolar. That is 42.29%, not 53.7%. The count of repeat readings falls with it, from 6,455 to 3,377.

2.7 What the tree found

In 1991, somebody stripped the peptides off an MHC molecule and sequenced them all at once. The pattern was unmissable. Two positions out of nine did nearly all the gripping.

The tree found the same two positions. Nobody told it they existed.

The grouped split put 6,971 of the 8,346 peptides on the training side and the other 1,375 on the test side. The tree was given those 6,971 training peptides, nine letters each, and the same 180 columns. That is all. No groove, no pockets, no anchors. Nine letters and a label saying bound or did not.

Its first question was position two, leucine. Then position two again, methionine. Then position nine, valine. Then position two, isoleucine. Then position nine, leucine.

Five questions. Two positions. It had nine to choose from.

Those two are the pockets. An MHC molecule is the clamp on your cell surface, and the second and ninth residues drop into it and hold the fragment still. The tree found the clamp by counting.

It did not discover it. Falk and colleagues published the pattern in Nature in 1991, thirty-five years ago.

A tree insisting that position four decides everything would have told me the pipeline was still broken: contaminated units, duplicated peptides, relatives sitting on both sides of a split. Any of those shows up as a model confidently finding the wrong thing. It found the right thing instead.

Chapter 1 ended the same way. The model walked to the one position in gyrA that microbiologists have known about for decades. Recovering a known answer is a receipt, not a finding.

2.8 What a score is worth

I can build you a model that is right 63% of the time and never finds a single binder. It answers no to everything.

37.2% of the 1,375 peptides in the test set bind. Say no to all of them and you are correct about the other 62.8%. Accuracy 0.6284. Binders found: none.

That 37.2% is a smaller share than the 42.29% above because this one counts peptides, with every repeat reading collapsed into one label.

Accuracy counts how often you were right. It says nothing about what you were right about.

So you need two numbers instead of one. Of the peptides you called binders, how many actually bound? That is precision. Of the peptides that actually bound, how many did you catch? That is recall.

Either one alone is easy to fake. Call a single peptide a binder and be right about it, and precision is perfect while recall is almost nothing. Call everything a binder, and recall is perfect while precision collapses to the share that bind. F1 is one number that only rises when both rise. The model that says no to everything scores 0.0000.

Now the comparison. Every number below comes from the same held-out proteins.

The textbook rule first. Leucine or methionine at position two, predict binder. Accuracy 0.7542. F1 0.6667. One question about one position, and it beats saying no to everything on both counts.

Then the trees, tested on the same proteins they never saw during training. Nine columns of alphabet numbering, F1 0.6983. 180 columns of one-hot, F1 0.7153. 27 columns of hydrophobicity, molecular weight and charge, F1 0.7399.

So yes. The tree beats the rule. One hundred and eighty columns and thirteen levels of questions bought 0.0486 F1 over one rule about one position, and the 27 columns of chemistry bought 0.0732. Chapter 1 went the other way, where eighty-four features bought nothing. Here the model wins, and that is the size of the win.

I expected one-hot to win. It did not. I am saying so after the scores, and without the note I said I would write, so it is not a prediction.

2.9 Why the smallest encoding won

To one of my encodings, leucine and isoleucine were neighbours. To another they were strangers. That difference decided which model won.

I wrote down the cost of one-hot as a caveat before anything ran. It turned out to be the result.

Look at how deep each tree had to go. Chemistry needed seven levels of questions. One-hot needed thirteen.

That is the mechanism. A tree using chemistry asks one question. Is this position hydrophobic. That catches leucine, isoleucine, valine and methionine at once. A tree using one-hot has to ask about each of them separately, four questions to say what the other tree said in one.

The MHC groove does not read letters. It has a pocket with a shape and a chemistry, and any residue that fits, fits. An encoding that describes shape and chemistry sits closer to what the molecule is doing than one that spells names.

One caution before anyone generalises this. Three physicochemical properties per position is a choice I made before seeing any result, and a different three might do better or worse. This is one comparison, on one allele, with one model.

2.10 The depth I reported

Eight questions deep scores a cross-validated F1 of 0.7095. Thirteen questions deep scores 0.7118. Both are the one-hot tree. I reported thirteen. I am not certain that was right.

A decision tree needs a depth limit. Without one, the one-hot tree kept splitting until every group agreed with itself, which took 1,140 groups for 6,971 training peptides, and its accuracy on them was 1.0000. On peptides it had not seen, that same tree scored 0.7789 for accuracy, and the one I kept scored 0.7927.

So I tested every depth from one to twenty, five times each, on peptides held back from different proteins each time, and kept the depth with the best cross-validated F1. Thirteen won (Figure 2.7).

The sweep ran over all 8,346 peptides, and the training and test split was drawn after it, so the 1,375 test peptides sat inside the folds the depth was chosen on. The trees scored earlier in this chapter were fitted on the training peptides alone, and every fold of the sweep held its own peptides back while its tree was fitted, but the depth those trees were given had been chosen over a cohort that included the test peptides. Of everything I report, the depth is the one choice the test set helped make.

Now look at the curve it won on. The score climbs steeply to about depth eight. After that it is almost flat. Depth eight to depth thirteen buys 0.0023. Five extra levels of questions for two thousandths.

There is a standard rule for exactly this situation, and it exists because a flat plateau is mostly noise, and the peak of noise is not a real peak. Under that rule I would have reported depth eight: a tree half as deep, easier to read, almost the same score.

So why did I report thirteen? Because the rule I wrote down beforehand said take the maximum, and I only saw the plateau afterwards. Changing the rule after seeing the answer is how you get results that do not replicate. Even when the new rule is better. Even when you are certain your motives are clean.

Both numbers are in the repository. Depth thirteen is the reported one because it is the pre-specified one.

Notice what this is not. It is not a data problem, or a modelling error, or something a better pipeline would catch. It is a judgement, made by me, that changes the number I publish. There is no neutral depth either.

2.11 Which mistakes it makes

My best model missed 117 real binders. My second best missed 153. There are 511 binders in the test set, so those are both bad numbers. They are bad in different ways, and F1 cannot tell you how.

1,375 peptides. 511 bind. 864 do not.

The chemistry model called 554 of them binders. 394 of those calls were right and 160 were wrong. The one-hot model was more cautious. It called 490 binders, and 358 were right. Being cautious cost it: it missed 153 real binders where chemistry missed 117. So chemistry catches 36 more binders than one-hot and raises 28 more false alarms doing it. That is the trade underneath a gap of 0.0246 in F1.

Precision: chemistry 0.7112, one-hot 0.7306. Recall: chemistry 0.7710, one-hot 0.7006. One-hot is the better guesser. Chemistry is the better finder.

Which one you want depends on what a mistake costs, and nothing in the data can tell you that. A false alarm is one assay you run and throw away. A missed binder is a peptide you never look at again. For screening vaccine candidates, those two prices are nowhere near equal. Twenty-eight more assays would be cheap. Thirty-six peptides you never test could include the one that mattered.

Now the textbook rule, on the same 1,375 peptides. Leucine or methionine at position two catches 338 binders and raises 165 false alarms (Figure 2.8). Chemistry catches 56 more than that, and raises 5 fewer. Better on both counts at once. No trade, no judgement call, just better.

That is the clearest thing this chapter has produced, and it took 27 columns and seven levels of questions.

2.12 Peptides that disagree with themselves

127 peptides in this dataset disagree with themselves. They were measured more than once, and one result says binder while another says not. I expected that number to be far worse.

8,346 distinct nine-mers sit in the corrected cohort, and 2,660 of them were measured more than once, in different labs, in different years, with different reagents. Nearly a third of the dataset. For most of them the repeats agree. Two readings of 15 and 40 nanomolar both say binder. Two at 6,000 and 12,000 both say not. 127 do not agree, because one reading falls below 500 nanomolar and another falls above it. Same nine amino acids. Same molecule. Two answers. That is 4.8% of the repeat-measured peptides.

So what do you do with 127 peptides that contradict themselves? You cannot keep both readings, because the model would see the same nine letters labelled two ways and learn nothing from either. You cannot drop them either, because that deletes the peptides sitting nearest the line, and those are the ones a binding model most needs to get right.

I took the median of every peptide’s measurements and labelled that. Three readings of 200, 400 and 900 become 400, and 400 is a binder (Figure 2.9). It is a defensible choice. It is still a choice, and a different one would give a different training set.

Now the part worth sitting with. A binding affinity comes out of an assay, run by people, with reagents that vary. Run it twice and you get two numbers. Draw a line anywhere through the middle of them and some pairs will land on opposite sides. Those 127 peptides are what measurement looks like when you stop rounding it off. Move the threshold and a different set of peptides straddles it. The disagreement does not go away. It relocates.

2.13 The allele filter

My human dataset contained 263 measurements from a rhesus macaque. Also horse. Also cattle. Also chicken.

None of that was a download error. I wrote one line of code that let them in.

The allele I am studying is HLA-A*02:01. To pull its measurements out of a nine-gigabyte file, I asked for every row whose allele name contains 02:01. That is the natural thing to write, and it is wrong.

Mamu-A1*002:01 is a macaque molecule. Eqca-1*002:01 is horse. BoLA-3*002:01 is cattle. Gaga-BF2*002:01 is chicken. All of them contain 02:01, because the naming is done per species. The number after the asterisk identifies an allele within its own species, not across them. Two different species can both have an 02:01 and share nothing but the label.

Then there were eleven more. HLA-A*02:01 K66A and HLA-A*02:01 A150P are point mutants, engineered in a lab, one amino acid changed inside the groove. Those are human, and they are the same allele, and they are still the wrong molecule. The whole point of mutating position 66 is to change what binds there.

Together they added 302 measurements and 262 peptides that do not belong. Matching the allele name exactly removed all of them and reproduced my feasibility scan to the row. The cohort of 11,723 is what the exact match leaves, so none of those rows are in it. Listing 2.2 shows the filter.

ALLELE     = "HLA-A*02:01"  # Exact string, deliberately not a substring match.
                            # Matching on "02:01" also picks up Mamu-A1*002:01
                            # (macaque), Eqca-1*002:01 (horse), BoLA-3*002:01,
                            # Gaga-BF2*002:01 and eleven point mutants of
                            # HLA-A*02:01 itself (K66A, A150P and so on). Those
                            # are different grooves on different molecules.

Unit contamination announced itself. Melting points and half-lives sit in the wrong column, and a check catches them. This did not announce itself. Every one of those rows is a real peptide, measured properly, against a real MHC molecule. Nothing is corrupt. They are simply not the molecule I said I was studying, and a model would have learned from them happily. The nomenclature was doing exactly what it was designed to do, and my filter was the thing making the claim.

A substring is not an identifier.

That is where this chapter ends. Not on the model that worked, but on the filter that almost let a macaque define a human molecule. Everything is in the repository: the corrected cohort, the three encodings, every metrics file and every number I got wrong on the way. Chapter 3 asks what happens when two annotation tools read the same genome and disagree about what is in it.

This chapter is rewritten from the posts for Days 12 to 25 of the series.

References

Falk, K., Rötzschke, O., Stevanović, S., Jung, G. and Rammensee, H.-G. (1991). Allele-specific motifs revealed by sequencing of self-peptides eluted from MHC molecules. Nature, 351(6324), 290-296. https://doi.org/10.1038/351290a0

Sidney, J., Assarsson, E., Moore, C., Ngo, S., Pinilla, C., Sette, A. and Peters, B. (2008). Quantitative peptide binding motifs for 19 human and mouse MHC class I molecules derived using positional scanning combinatorial peptide libraries. Immunome Research, 4, 2. https://doi.org/10.1186/1745-7580-4-2