Which proteins mark response to melanoma immunotherapy?

The first question is not the scientific one. It is worth knowing what is in a deposition before deciding what to ask of it.

Under twelve seconds, answered from the data profile that was built when the files were registered. The deposition comes back split into two sub-studies that nobody labelled as such: a clinical arm of super-SILAC ratios covering aPD1-1 at roughly 10,011 protein groups over 50 experiments, aPD1-2 at 9,092 and 33, and TIL at 10,152 and 44; and an in-vitro arm of label-free intensities covering four melanoma lines under dichloroacetate and a CRISPR knockout screen of five fatty-acid and lipid genes.

The four hazards it lists at the bottom are the ones that shape the rest of this page. The five matrices are independent searches whose rows do not align, so joins have to run on gene symbols. Reverse and contaminant rows arrive unfiltered. The sample sheet has duplicate header names and lists every raw file twice, once for the specimen and once for the heavy reference. And the ratios and the label-free intensities must never be concatenated. Every one of those becomes a decision further down, and one of them becomes the whole second half of this page.

Then the actual request. It names no method, no software, no gene and no threshold. It asks a question and adds a second question that depends on the first one's answer.

Three minutes and twenty-four seconds later there is a six-step plan across four tracks. Load, filter and normalise; build the patient-level annotation; then a fork into the responder contrast and the Cox and log-rank arm, which rejoin at a cross-cohort consistency step and an integrated shortlist.

The shape of that graph is the answer to the second half of the prompt. Asking whether the same proteins track with survival is a request for a joint claim, and the plan makes the two arms independent and then intersects them, so the place where the two questions meet is a step with a rule in it rather than a sentence at the end.

The plan was approved and executed as run 5a8f439f-38eb-48e7-b8c6-febe171effc2. What comes back is the whole first pass, and the rest of this page unpacks it: the null, the artefact, the survival result that validated nothing, and one candidate protein. The closing half of the same answer sits at the end of section 5, where it does more work than it would here.

Zero at FDR 0.05 and zero at FDR 0.10, in every cohort. A Wilcoxon test on the identical complete cases also returns zero, so the normal model is not what produced the null. A Perseus-style down-shifted imputation was run as a sensitivity analysis and changed no conclusion, which is why the no-imputation fit stays primary. Imputation correlates with the unimputed fold changes at only 0.82 to 0.84, with individual shifts up to 2.79 log2 units, and that distortion lands on the sparse rows where left-censoring bites hardest.

The confound that would have made imputation dangerous is absent, and somebody checked instead of assuming it. Response groups do not differ in depth in any cohort, at Wilcoxon p of 0.84, 0.43 and 0.42 in table order.

A bare null tells you very little, so the calibration is where the information is. All three cohorts show a mild excess of small p-values over the global null, at 1.16, 1.56 and 1.35 times expectation, with null-proportion estimates of 0.939, 0.837 and 0.940. Somewhere between 6% and 16% of the tested proteome is plausibly non-null. Something is in there; these sample sizes cannot say what.

The second anti-PD1 batch gets its own label instead of being folded in. Its patient stratum is marked not estimable with zero samples, and the 8 versus 8 comparison that does exist is drawn entirely from xenografts. Where its patient volcano would have gone there is a placeholder panel, so the absence sits on disk as a file you can open.

Across 4,708 shared symbols the anti-PD1 pair correlates at r = -0.00033 with a 95% interval of -0.029 to 0.028 and p = 0.982. Sign agreement is 0.4868, which a binomial test cannot distinguish from chance at p = 0.073. Meanwhile aPD1-1 against the mechanistically unrelated TIL arm reaches r = 0.183 with sign agreement 0.5576 at p = 2.9e-14.

The explanation is upstream: the second anti-PD1 batch is entirely xenograft. Human stromal and leukocyte content is lost in xenografts and replaced by mouse, so a patient-versus-xenograft split is expected to dominate a batch-versus-batch split. The practical consequence is that these two batches cannot be treated as replicates of one another, and any pooling across them estimates no single underlying effect.

It also spoils the reassuring heterogeneity statistics. Median between-study heterogeneity across the meta-analysis is 0%, which sounds like agreement, but that statistic measures dispersion relative to each study's own standard error and those errors are large here. Two facts drain the 0% of meaning: the genome-wide correlation between the anti-PD1 batches is zero, and the effect-size scales differ nearly twofold, with fold-change standard deviations of 0.402, 0.491 and 0.276.

The agreement that does exist sits in one narrow place. A block-stratified diagnostic, added beyond the original brief because the correlation pattern needed explaining, finds that a co-regulated set of 15 antigen-presentation proteins (the class I and II HLAs, B2M, CD74, TAP1, TAP2, TAPBP, ERAP1 and the immunoproteasome components) carries almost all of it: r = 0.648 within the block against 0.145 in the remaining 4,329 proteins.

That block is displaced the same way in all three cohorts, with Mann-Whitney p values of 7.6e-13, 2.9e-08 and 9.8e-09 against a whole-proteome mean of essentially zero. Pooling the three cohorts turns the run's zero FDR-significant proteins into one, HLA-B, at a random-effects adjusted p of 0.0438. HLA-B belongs to the block.

The block is also displaced, just as coherently, in the xenograft cohort, which has no human immune compartment and cannot host an infiltrate signal at all. The whole second half of this analysis grows out of that observation. It is also why those 15 concordances are worth roughly one degree of freedom between them: the members are co-regulated and co-quantified, so every enrichment p-value is anti-conservative for the one block it most wants to report.

Nobody could read the progression events off a field, so they were derived. Overall survival is always at least progression-free survival, so a strictly shorter progression time means progression was observed, and an equal time means the patient was progression-free when the clock stopped. Since that is inference, it was audited against RECIST before anyone used it: all 26 progressive and stable disease cases come out as events, all 8 complete responses come out censored, and the partial responses split 14 to 14. A correct indicator should behave like that, and this one does.

Overall survival was never fitted. The log states the bounds plainly: the lower bound on events is zero if everyone is alive at last contact, and the upper bound is the number who progressed. There is a sensitivity endpoint that assumes death and progression coincide, labelled as an assumption on every row of every table. It correlates with progression-free survival at 0.89 to 0.96 by construction, so it replicates nothing.

The one analysable cell is TIL progression-free survival, at n = 41 with 32 events. A univariate Cox scan across 4,871 proteins returns zero hits at proteome-wide FDR, with a best adjusted p of 0.2207. Null calibration again shows a real but diffuse excess, 1.41 times expectation with a null proportion of 0.907. Proportional hazards causes no trouble: across 33 leading-protein rows nothing violates, and the smallest test p in the file is 0.134.

Here is the part that earns the section. The 25 proteins carried forward from the differential arm look strongly survival-associated: 22 of 25 reach raw p below 0.05 and 22 of 25 clear a focused FDR over the 25 candidates, with a permutation empirical p of 5e-04 over 2,000 draws. Read at face value, a differential result has just been validated on an independent endpoint.

The endpoint is not independent. Because the event indicator comes out of the survival times, all 21 TIL non-responders are events by construction while only 11 of the 20 responders are. Any protein picked on the responder contrast inherits a hazard ratio for free.

Protein by protein: ACAT1 falls from a hazard ratio of 3.94 at p = 1.1e-04 to 2.375 at p = 0.028 with a focused FDR of 0.692; NUDT4 from p = 3.6e-04 to 0.346; SUPV3L1 from 9.9e-04 to 0.194; HTATIP2 from 4.4e-03 to 0.806; ERLEC1 from 8.0e-03 to 0.963. Every one of them carries a proteome-wide adjusted p between 0.22 and 0.70 regardless of stratification.

The control exists because the discovery and validation endpoints share a variable, and the lesson travels well past this dataset. Select proteins on a response contrast, then validate them on a progression endpoint whose event indicator is itself determined by progression, and you will recover strong-looking survival associations that only restate what you selected on.

5,964 gene symbols enter and none survive. Two counts are reported alongside, because a reader will want them: 1,602 proteins keep the same direction across every cohort in which they were measured, and 2,435 do so across the two patient cohorts. Under independent signs you would expect 25% by chance, and 1,602 out of 4,134 is close to it. Directional consistency on its own is not evidence.

The step also says something most integration steps leave out. Intersecting two FDR-controlled lists does not control the joint error rate in either direction. It neither preserves the false discovery rate nor gives you a valid combined p-value. Had the intersection come back non-empty, the honest label for it would still have been proteins that ranked well on two related criteria.

With no FDR-defined response set, an FDR-defined enrichment is undefined, so enrichment ran on nominal sub-FDR sets, each labelled as such, against that cohort's own quantified proteome instead of the whole genome. Two coherent programmes came out, and over-representation and rank-based methods both find them:

• Antigen presentation, interferon signalling and allograft rejection. Interferon gamma response reaches 3.0-fold over-representation in the anti-PD1 cohort and a normalised enrichment score of 2.67 on the ranked meta-analytic effect, the strongest set in the analysis.
• Mitochondrial oxidative metabolism. In the TIL cohort, oxidative phosphorylation reaches 3.79-fold at an adjusted p of 3.9e-12, with fatty-acid metabolism, adipogenesis, fatty-acid beta-oxidation and the TCA cycle behind it. The same programme recurs independently on the survival arm.

The run then wrote down which way each one pointed. This is the moment where a re-analysis becomes capable of being confidently wrong.

Those two objections are not buried in a limitations file. They are the closing section of the run's own summary, under the heading "What I'd do next", above the one candidate it was willing to carry.

Five minutes and seventeen seconds went by before a plan appeared, and most of it was not spent planning. After checking which reference panels were available, the agent went back into the first run's own source: it read the load step's summary, grepped its scripts for the ratio column and for anything mentioning polarity or inversion, read the first sixty lines of the loader, then grepped the sample sheet for the label assignment. Then it stopped and said it had found the bug before writing the plan meant to look for it.

The finding is a two-line contradiction inside the first run, and both lines are in the package. At 01_load_filter_normalise.py:49 the loader selects MaxQuant's heavy-over-light normalised ratio columns and applies no reciprocal. At 03_responder_vs_nonresponder_limma.R:11 a comment documents those same values as the log2 of specimen over reference. Heavy over light is heavy over light by definition, and the sample sheet puts the patient specimen in the unlabelled channel and the super-SILAC melanoma mix in the heavy one. The two statements are reciprocals. One of them is wrong, and the one that is wrong is a comment, which is why nothing downstream ever complained.

That is enough to know the answer and not enough to act on. A comment is not evidence, and a convention read out of documentation can be wrong about a particular file. So the documentary reading was converted into a prediction, and the plan built to test it was allowed to contradict it.

The plan narrates itself before asking for approval, and three of its paragraphs are decisions, not descriptions. Polarity gets two panels whose expected signs oppose each other, so no single artefact can satisfy both. The purity scores are kept as separate numbers, computed within cohort, with the all-xenograft cohort as a positive control that is allowed to fail loudly. And the collinearity check is deliberately scheduled before the adjusted models, because CD74 and HLA-DRA belong to both the block being tested and the score doing the testing, and an attenuation you cannot distinguish from arithmetic must not be reported as a refutation.

The last paragraph before the approval prompt is a limitation, stated while the result could still have gone either way: the anti-PD1 cohort is 16 against 5, a covariate costs a degree of freedom it can barely afford, and a null adjusted result there will be weak evidence of absence. Sections 9 and 10 of this page are the collinearity decision and that limitation playing out.

The design is to score marker panels whose direction biology already tells you. If the stored value is specimen over reference, then blood and leukocyte proteins in a blood-perfused, immune-infiltrated tumour biopsy should come out strongly positive against a denominator made of cultured cell lines that contain none of them.

The blood and leukocyte panels are negative everywhere, with medians of -2.3 to -4.7 log2 units. Taken at face value under the documented orientation, that would require the cultured melanoma reference mix to contain 5 to 26 times more haemoglobin, fibrinogen, haptoglobin and CD45 than blood-perfused, immune-infiltrated patient tumour tissue. Nothing marginal about that call. It is biologically impossible.

A further panel was added after two of the tests conflicted, and it was chosen so its expected direction is the mirror of the tissue panels: cell-cycle and replication proteins, far more abundant in exponentially growing cultured lines than in stroma-diluted bulk tumour. Under the correct orientation those should be positive, and they are, in all three cohorts, at medians of +0.88, +0.85 and +0.56 with sign-test p values between 4e-04 and 2.6e-03. TYMS, TK1, RRM1, RRM2, PCNA, MKI67, TOP2A and every MCM subunit point the same way.

A melanocytic marker panel ran too, and its verdict is recorded as non-discriminating by construction, since both channels are melanoma. It shows near-complete detection and no sign bias, which is how an uninformative panel should look. Worth reporting for that reason alone.

The verdict is a polarity factor of -1: the stored values are the log2 of reference over specimen, and the comment in the differential-expression script had it backwards. That is the same answer the agent read out of the code before the plan existed, arrived at by a route that never consulted the code. The agreement is what makes it a verdict. A test whose answer was fixed in advance would have been a formality, and this one was not: a third test in the same step came back pointing the other way.

Under the corrected orientation a positive fold change means higher in the responder. The antigen-presentation block now comes out higher in responders in all three cohorts, the canonical direction. Interferon signalling is up in responders. So are oxidative phosphorylation, fatty-acid metabolism and adipogenesis. Higher ACAT1 now means slower progression, not faster.

Enrichment was recomputed from scratch instead of relabelled, which was necessary in principle and dull in practice. On a tie-free continuous statistic, multiplying by -1 is a strict rank reversal, so the leading edge of every one of the 50 hallmark sets comes back identical, no gene enters or leaves any of them, and only the enrichment score's sign moves. Recomputing was still the right call. That identity holds for this transformation and is not a licence to relabel in general.

One row matters more than the others: the comparison with the source publication. The first pass reported that these specimens inverted Harel 2019's mitochondrial model. Corrected, they reproduce it. Oxidative phosphorylation and lipid metabolism are higher in responders in both treatment arms, which is what the original analysis of these same specimens said.

Call that a positive control on the pipeline and stop there. Agreeing with the original analysis of the same data restores internal coherence; it validates nothing externally. The independent literature also makes uncomfortable reading. Patient-derived melanoma lines with deranged oxidative metabolism progress on anti-PD1, checkpoint-resistant models acquire a hypermetabolic phenotype, and a high-oxidative CD8 T-cell state predicts resistance. Steady-state protein abundance in bulk tissue, oxygen consumption rate in derived lines, and intratumoral hypoxia are three different measurements. This analysis has the first one only, so the disagreement stays on the page instead of being talked away.

In the only contrast not confounded with search batch, the area under the curve is 0.5152 with a 95% interval of 0.136 to 0.864 and p = 0.969. The point estimate is on the wrong side of chance and the interval spans it. The three well-powered cross-cohort contrasts, at 33 against 63 specimens, land flat at 0.478 to 0.485.

The diagnosis teaches more than the failure does. Xenograft specimens score high on human immune, human stromal and human blood panels at once, at medians of +2.2, +5.5 and +2.0 against roughly zero for patient tumours. Per protein, the markers most elevated in xenografts are membrane and cytoskeletal proteins with well-conserved mouse orthologues, while the lineage-restricted and MHC class II-dependent members move the other way. A real human leukocyte signal would move them together. What this looks like is a human-only search database absorbing conserved mouse peptides into human protein groups, which is the problem species-aware xenograft proteomics pipelines exist to solve.

Three attempts to rescue the score all fail, and all three are reported. An orientation-free version using only whether a marker was detected returns 0.447. Residualising the immune score on a shared stromal and blood offset returns 0.424, even though that offset separates the xenografts strongly on its own at 0.939. One score in the file does place xenografts low, the global detection rate at 0.136, and all it says is that xenografts yield fewer identifications overall. Nothing immune-specific about it.

The score went forward anyway, flagged not validated, with the flag propagated into every downstream artifact. That was the right call, and it carries a permanent constraint on the language. For the rest of this analysis, saying an effect was attenuated by purity means it shares variance with an unvalidated proteome-internal axis. It never means tumour purity explains it.

The antigen-presentation block correlates with the immune score at 0.824, 0.761 and 0.797, explaining 58 to 70% of its variance, with a maximum variance inflation factor of 3.6. The metabolic block correlates at 0.151, -0.197 and 0.185, with only one of its six intervals against the two purity scores excluding zero, and 79 to 95% of its variance unexplained by purity. The step's own brief had predicted that oxidative phosphorylation would be the complementary side of the purity axis. Its own data say otherwise, and it records that.

Two design choices keep those numbers from being circular. First, the eigengenes use the identical function that built the purity scores, and recomputing the purity scores that way reproduces the published ones at a correlation of 1.0 with a maximum absolute difference below 1e-15. The two quantities being correlated are genuinely commensurable.

Second, CD74 and HLA-DRA sit in both the antigen-presentation block and the immune scoring panel, so shared membership accounts for part of the correlation. Instead of disclaiming it, the step built both stripped variants. Take those two proteins out of both sides and the correlations fall to 0.737, 0.565 and 0.606, with every interval still excluding zero. The relationship costs 0.09 to 0.20 in correlation and survives.

A background comparison makes the same point sharper. Against all quantified genes that belong to neither the block nor the scoring panel, the antigen-presentation block's median member sits at the 99.5th to 99.9th percentile of the background distribution of correlation with the immune score, and in one cohort all 15 members individually clear the background 95th percentile. The metabolic block sits between the 44th and 73rd percentile. Unremarkable.

What that means for the adjustment cuts both ways. No variance inflation factor anywhere reaches even the conventional threshold of 5, so the adjusted model is identifiable and a surviving coefficient means something. But the antigen-presentation block's standard error inflates 1.84 to 1.90 times at n = 21 and n = 41, while 58 to 70% of its variance goes away by construction, so you cannot tell its attenuation apart from lost power. Cannot separate is the honest phrase here; purity explains it is not. For the metabolic block, adjustment costs 1.06 to 1.16 times, close enough to free that any change there is real confounding.

In TIL the metabolic family keeps 5 of 5 sets at a median enrichment cost of 1.84%, with oxidative phosphorylation still at an adjusted q of 1.2e-18. The interferon and antigen-presentation family keeps 0 of 5. In the anti-PD1 cohort the same asymmetry appears more weakly: metabolic keeps 1 of 5, interferon and antigen presentation keeps 0 of 5.

The adjustment is a scalpel, not a hammer, which is what makes the asymmetry readable. It barely moves the average protein, with an adjusted-on-unadjusted slope of 0.90 and 0.97, while stripping a median 48.6% and 47.7% off the antigen-presentation block against a background attenuation of -3.7% and -0.3%. The block's immune coefficient runs six to seven times the proteome median.

Inside that block, the two proteins that also belong to the scoring panel are the two most attenuated, at 88% and 76%. Circularity doing precisely what you would predict. It is flagged on each row; dropping those proteins would have made the problem invisible instead of fixing it.

Two named proteins go opposite ways, and both are instructive. GBP1 ranks first in the anti-PD1 responder contrast at a corrected fold change of +2.98, then splits by cohort under adjustment: 24.7% attenuation in anti-PD1, 79.4% in TIL, where its adjusted p reaches 0.76. It never approaches FDR significance in either. The literature puts GBP1 mostly in macrophages, so an infiltrate read-out that fails to replicate across arms and collapses under a composition adjustment is doing what you would expect of it.

ACAT1 does the opposite. Its response contrast moves by 0.47%. Its immune coefficient is +0.0086 per standard deviation at p = 0.926. On the survival side its hazard ratio goes from 0.2541 to 0.2540 under adjustment, a change of 0.045% in the log hazard, and proteome-wide the adjusted and unadjusted log hazard ratios correlate at 0.988. The purity covariate carries no survival signal of its own, with a hazard ratio of 0.918 and p = 0.61.

Then comes the control that does bite. Stratify the same Cox model on the response label and the prior 25 candidates collapse from 22 to 1 at raw p below 0.05, and to 0 at every FDR threshold, while the proteome-wide correlation between stratified and unstratified log hazard ratios drops to 0.361. ACAT1 is the lone survivor at p = 0.026, with a focused FDR of 0.64 and a proteome-wide FDR of 0.89. Purity turns out to be nearly orthogonal to the protein effects. The response label is anything but.

Proportional hazards was checked across 92 diagnostic rows, and 10 violate at p below 0.05. All 10 are in the output. ACAT1's own protein term satisfies the assumption in every model, and the violation in its adjusted fit sits on the purity covariate. B2M and HLA-DPA1 violate on the protein term itself, so their single-number hazard ratios are not valid constant-hazard summaries and nobody should read them as such. Both are non-significant anyway.

Three habits from these runs carry over to any dataset.

Test the sign; do not read it off a comment. An orientation you took on trust is a coin flip that propagates into every fold change, every hazard ratio and every enrichment direction downstream, and it is the one error that yields a perfectly coherent wrong answer. Two panels with opposite expected signs cost almost nothing, and no centring artefact can satisfy both at once.

Build the control that could falsify your covariate, and publish the result when it fails. The purity score here is beautiful on the inside and did not do what it claimed. Saying so made the analysis more useful than a clean adjustment nobody could check.

Work out what a correction can touch and what it cannot before you apply it. The ledger on this page has outlived the flips it records, because it tells a reader who doubts the verdict which rows would revert and which would stay put.