MS Charge-state deconvolution

The deconvolution engine is Jonas's; this write-up is by Claude (Anthropic's Claude Code). Code at ewoij/prot.

Jonas asked me to write up the deconvolution engine he built. His goal was never to compete with existing tools — it was to learn how LC-MS deconvolution works by writing one instead of calling someone else's. So this post is about a roughly thirty-line algorithm: what it does, the two runs where it gets the right answer, and the one where it does not. The failure is the more interesting half, and it is the reason the third test file is kept in the repo as a known failure rather than quietly dropped.

All the numbers and figures below come from running the code in this repo on the bundled data. The figures are rendered by scripts/blog_plots.py, which captures the engine's own intermediates rather than redrawing them.

The problem

Electrospray ionisation does not put one charge on a protein. It puts many, and a different number on each molecule, so a single protein of mass M shows up as a whole ladder of peaks — one per charge state z — at

m/z = (M + z · m_proton) / z

You see the ladder, not the mass. Deconvolution is the step that reads the mass back off the ladder. The catch is that neither M nor z is known for any peak, so it is one equation with two unknowns per peak, and the only thing tying them together is that all the peaks in the envelope must share the same M.

The algorithm

The whole engine is in src/prot/deconvolution/custom.py. It is a brute-force search:

Guess masses. Take the 10 tallest peaks, and for each of them assume every charge from 1 to 50. That gives 500 candidate masses, z * (mz - 1). One of them is right if the protein's real charge state is in range and one of its peaks is tall.
For each candidate, assign every peak a charge. Given a candidate mass M and a peak at m/z, the charge that would put it there is z = M / (mz - 1), rounded to the nearest integer.
Check the fit. Rebuild the predicted m/z from that rounded charge and measure the error in ppm. A peak within --error-threshold-ppm counts as explained by that candidate. This gives a boolean matrix of 500 candidates × ~1300 peaks.
Score. score = (peaks explained) × fillness, where fillness is unique(z) / (max(z) - min(z) + 1) — how much of its own charge range the candidate actually fills.
Report. For the winner, compute z * (mz - 1) for each explained peak and take the median; report the standard deviation alongside it.

Here it is, whole and unedited — comments, TODO and all. The plot_* calls are the debug-plot hooks, and they are what the figures in this post are captured from:

    @dump_plots_to("deconvolute")
    def deconvolute(self, MZ: np.ndarray, I: np.ndarray) -> CustomDeconvolutionResult:
        error_threshold_ppm = self.error_threshold_ppm
        protons_range = np.arange(self.min_charge, self.max_charge + 1)

        # mass candidates: we try to guess masses and we'll compute which one is the
        # best candidate
        largest_peaks = np.argsort(-I)[:10]
        M = np.array([z * (MZ[pi] - 1) for z in protons_range for pi in largest_peaks])
        M.sort()

        # TODO: plot protons_range, candidate peaks (spectrum) and distribution of candidate masses

        # we're going to work with matrix of shape (M, MZ), so reshaping those 2 to
        # avoid being bitten by broadcasting rules
        M = M.reshape(-1, 1)  # rows
        MZ = MZ.reshape(1, -1)  # columns

        # compute z (protons) to satisfy each peak (mz) for each mass
        Z = M / (MZ - 1)
        Z = np.round(Z).astype(int)
        Z[Z == 0] = 1  # to avoid div by 0 and a null charge can't get through the MS

        # now we compute the error, rebuilding MZ for each cell (z)
        E = np.abs((M + Z) / Z - MZ) / MZ * 1e6  # ppm

        # we compute a mask where E is below a threshold
        OK = E < error_threshold_ppm

        plot_candidate_matrix(M, MZ, Z, OK, error_threshold_ppm)

        # before scoring the best mass, we want to compute the "fillness" :-D
        # it could happen that multiples of the same mass have as many peaks resolved
        # (example: 18 kDa, 36 kDa)
        # I'm assuming the one with the most consecutive z is more realistic as I
        # believe it would be rare to have consequent gaps in the distribution of
        # protons by the ionization (maybe I'm wrong)

        def get_fillness(arr: np.ndarray):
            return np.unique(arr).size / (np.max(arr) - np.min(arr) + 1)

        F = np.array([get_fillness(Z[i, OK[i, :]]) for i in range(Z.shape[0])])

        # now we have the number of peaks that could be recomputed by mass and the fillness.
        # We can score each mass. The mass that wins is the one with the highest number
        # of explained peaks with a penalty on fillness. The more filled z, the better.
        score = OK.sum(-1) * F

        plot_candidate_table(M, Z, OK, F, score)

        plot_candidate_scores(M, OK, F, score)

        # we can now find the best candidate
        best_candidate = np.argmax(score)

        # now that we have the best candidate, we can try to compute a more precise mass using the z
        ok = OK[best_candidate, :]
        z = Z[best_candidate, ok]
        mz = MZ[:, ok].reshape(-1)
        masses = z * (mz - 1)

        plot_mass_estimates(masses)
        plot_explained_spectrum(MZ, I, ok, z)

        return CustomDeconvolutionResult(
            mass=np.median(masses),
            std=np.std(masses),
            n_peaks=ok.sum(),
            coverage=ok.sum() / ok.size,
        )

Two things are easier to see here than to describe. protons_range appears once, in the candidate-building line, and never again — the charge limit bounds which masses get guessed, not which charges get assigned, which is the inconsistency this post runs into later. And E rebuilds the predicted m/z as (M + Z) / Z, which is where the proton silently becomes 1 Da instead of 1.00728.

Step 4 is the only non-obvious part, and it exists because of a specific failure: a mass and its integer multiples explain overlapping sets of peaks. A candidate at 2M predicts, at charge 2z, exactly the m/z that M predicts at charge z, so it inherits half of M's ladder for free. The fillness penalty is meant to catch that, because the doubled candidate's charges come out all-even and its range is full of holes.

Where it works

Filgrastim, 18.8 kDa, at a 100 ppm threshold:

Filgrastim apex spectrum, with the peaks explained by the winning mass drawn in
green and labelled with their charge state, and the unexplained peaks in
grey

The green ladder runs from z=7 to z=25 without a gap, and it sits on the tall peaks. That is what a correct answer looks like: the engine reports 18,796.72 Da against FLASHDeconv's 18,797.58 Da for the same scan — 0.86 Da apart, about 46 ppm.

The engine also reports a standard deviation of 0.86 Da (the same number, by coincidence), and it is worth seeing where that comes from, because it is not measurement noise:

Per-peak mass estimates for filgrastim, sorted: three flat plateaus about 1 Da
apart, with the median line through the middle one

Three plateaus at 18,795.68, 18,796.67 and 18,797.72 Da, holding 15, 17 and 19 peaks. They are about 1 Da apart because each explained peak contributes its own estimate z * (mz - 1), and the peaks within one charge state are isotopes — the same molecule, one neutron heavier. The engine has no concept of isotopes, so it pools them all, and the standard deviation it reports is mostly the distance between neighbouring isotopes rather than any measure of precision.

This also explains the gap to FLASHDeconv. Its 18,797.58 Da sits on the third plateau; the median lands on the middle one. The engine is not 0.86 Da imprecise, it is one isotope over — and which plateau holds the most peaks is decided by which isotopes happened to clear the ppm filter, not by anything the algorithm reasons about.

Cytochrome c is the same story at a 1000 ppm threshold — charges 6 to 21, consecutive, on the tall peaks:

Cytochrome c apex spectrum, with a consecutive green charge ladder from 6 to 21
on the tall peaks

It reports 12,230.23 Da. Cytochrome c is worth checking against something other than another program: its mature sequence weighs 11,572.21 Da, and the molecule carries a covalently bound heme (+616.49) and an acetylated N-terminus (+42.01), so it should weigh 12,230.71 Da. That is a 0.48 Da gap, 39 ppm, derived from chemistry rather than from other software.

Plotting the score against every candidate mass shows how decisively the right answer wins on this kind of data:

Score against candidate mass for filgrastim: one sharp spike at 18,797 Da,
smaller clusters at roughly twice and three times that, and a low flat background
everywhere else

The two smaller clusters are the harmonics, and they are where fillness earns its keep:

candidate	peaks explained	fillness	score
18,796.71 Da	52	1.000	52.00
37,597.38 Da (2×)	59	0.543	32.07
56,390.13 Da (3×)	57	0.436	24.87

The doubled mass explains more peaks than the true one — 59 against 52. On peak count alone it would win. Fillness is what puts it back in second place.

Where it does not work

The third bundled run is the Pierce Intact Protein Standard Mix: six proteins in one 42-minute run. The TIC apex falls inside the Exo Klenow peak, so the answer there is 68,000.5 Da. The engine says 34,000.46 Da.

Standard mix apex spectrum with z ≤ 50: a dense forest of grey peaks with a
gappy green ladder between charges 25 and 50 running through
it

And 34,000.46 × 2 = 68,000.92. The engine found the right peak series and labelled every rung with half its real charge — a charge span of 25 to 50 where the real protein's is 50 to 100.

The diagnostics do say something is wrong: 23 explained peaks, 1.6% coverage, and a gappy ladder (25-27, 30-33, 35-39, 41-44, 46, 48-50) worth only 0.769 fillness. Compare that to filgrastim's 52 peaks at fillness 1.000. But nothing in the algorithm acts on any of it. It reports the argmax of the score, and a score of 17.69 was the best on offer.

The reason is step 1. FLASHDeconv assigns Exo Klenow charges 40 to 100, and the ten tallest peaks in this spectrum all sit below m/z 1310. Building candidates as z * (mz - 1) from those peaks with z ≤ 50, the largest mass the search can even represent is 65,384 Da. The correct answer is not a candidate. The score did not pick wrong; the right answer was never on the ballot.

Score against candidate mass for the standard mix with z ≤ 50, showing the
expected 68,001 Da falling beyond the last candidate

That plot also shows the second, worse problem: the score cloud drifts upward with candidate mass. A bigger candidate has more charge states landing in the instrument's m/z window, so it gets more chances to explain something. The score is not normalised for that, so it rewards large masses on principle.

Raising the charge limit to 120 puts 68,000.92 Da on the grid. It still loses:

Score against candidate mass for the standard mix with z ≤ 120: the winner at
136,000 Da scores 33.0, while the expected mass at 68,001 Da scores
29.7

The winner is 136,000.15 Da — twice the right answer, scoring 33.0 against the true mass's 29.7. The fillness penalty does fire: the doubled candidate scores 0.508 fillness against the true mass's 0.690, exactly the every-other-charge signature it was designed to catch. It wins anyway, on raw peak count, 65 against 43. Doubling the mass doubles the charges the search assigns, which moves the whole ladder into the crowded low-m/z end of the spectrum where there are simply more peaks to hit. The penalty is a factor of 0.74; the reward is a factor of 1.5.

(Those assigned charges run to z=221, incidentally, well past the --max-charge of 120. The charge limit bounds which candidate masses get generated, not which charges get assigned to peaks afterwards — an inconsistency in the implementation that nothing in the output would tell you about.)

And this is what the winner's explained peaks look like:

Standard mix apex spectrum with z ≤ 120: green peaks scattered densely across
the whole forest with no visible ladder

No ladder. Just green sprayed across a crowded region — 65 peaks, 4.7% coverage, a mass assembled from coincidences. Compare that to the filgrastim figure and the difference is obvious to a human eye and invisible to the score.

So there are three distinct failures stacked here, and only the first is a tuning problem:

The search space is bounded by max_charge, and the bound is expressed in charge while the answer is a mass. A 68 kDa protein needs z ≤ 120.
The score grows with candidate mass, so it is biased toward large answers regardless of evidence.
The score has no model of what an envelope should look like. It counts peaks and measures gaps. It does not know that a real charge-state envelope has a smooth intensity profile, that isotope peaks within a charge state are spaced 1/z apart, or that the peaks it explains should be the tall ones.

The engine also treats the proton as exactly 1 Da rather than 1.00728 Da, which biases each per-peak estimate high by about 0.007·z Da — around 0.1 Da at z=15 and 0.5 Da at z=70. Small next to everything above, and it points the opposite way from the gaps on filgrastim and cytochrome c, which the isotope plateaus account for.

What the real tools do instead

Every point in the list above is something the established algorithms address head-on, which is the clearest way I know to see why they are built the way they are:

Use the isotope pattern, not just the ladder. Above ~10 kDa you cannot compute a protein's isotope distribution without knowing its formula, which you do not have — so Senko's averagine stands in an average amino acid composition scaled to the observed mass. That converts "does this peak fit?" into "does this whole isotope cluster fit?", which is a far harder test to pass by accident. The repo's second engine delegates to ms_deisotope for exactly this, and gets 18,797.51 Da on filgrastim with a 0.07 Da spread.
Search in log space. In log m/z, the spacing of a charge-state envelope is the same pattern wherever the mass sits, which turns the search into pattern matching over one transformed axis. This is FLASHDeconv's core trick, and it is why it handles charges 2–100 without the search cost exploding.
Score against a model, not a count. Zhang & Marshall's Zscore, and the Bayesian treatment in UniDec, both compare the observed envelope against what a charge-state distribution should look like, so a candidate that explains peaks in an implausible arrangement is penalised for the arrangement itself.

The problem

The algorithm

Where it works

Where it does not work

What the real tools do instead

Further reading