Expression: the E of EPIC#

What this page is for. Two of the nine fitted terms are expression, they are the two most often scored at a constant by accident, and the fix is one command. Read the first section and you will know whether your table has the problem.

The one-line version. expr_lvl is what the candidate’s gene is transcribed at in the tumour; expr_norm is the same gene in the tumour’s matched normal tissue. Both are \(\log_2(1 + \mathrm{TPM}/c)\). Both are keyed on a gene symbol, and a row without one takes a mean-imputed constant — so the term measures nothing for that row.

First: does your table carry a gene symbol?#

If it does not, both expression terms are constants and two of the model’s nine coefficients are doing nothing for you. This is the common case, not the edge case: over the neoantigen corpus the symbol is missing on 356,387 of 695,811 rows (51.2 %), and on 5,205 of 5,833 immunogenic candidates (89.2 %). On the VACCIMEL screen that left expr_norm at standard deviation exactly 0.0000 and AUROC exactly 0.5000 — a term present in the model and absent from the answer.

mhcmatch genes recovers it, because a neoantigen is a near-copy of a self peptide: a near-exact proteome search names each parent by its UniProt GN= field.

mhcmatch genes pairs.tsv --species human --out annotated.tsv     # adds a `gene` column
mhcmatch rank pairs annotated.tsv --tumor SKCM --out ranked.tsv   # reads it: no join, no rename

Over that corpus coverage goes to 692,349 of 695,811 rows (99.5 %), 4,511 of the 5,833 positives gain a symbol, and expr_norm’s standard deviation on VACCIMEL goes from 0.0000 to 2.520 (bench/results/epic_gene_repair.md). A tie becomes several rows rather than a refusal, and an unresolved peptide keeps its row with an empty cell — losing the row would be the larger error. See The parent gene, when the deposit does not name one.

The two terms, and why they are two#

term

what it is

expr_lvl

the candidate’s source-gene abundance: the cohort’s own measurement where it has one, else the tumour type’s reference value, else the gene’s matched-normal or cross-tissue level. mhcmatch.rank.expr_level().

expr_norm

the same gene’s median in the tumour’s matched normal tissue, on the same floor, falling back to that gene’s pan-tissue median and never to missing. mhcmatch.rank.expr_norm_level().

They enter free, not as a ratio, and the data says to. A tumour-versus-normal ratio is a difference of logs, which a linear model can express only with equal and opposite coefficients. Entering the two separately lets that ratio be found — and it is not found: both coefficients come back positive. Imposing the ratio would have forced a shape the fit rejects.

The floor c is the tumour type’s own#

\(c\) is the 25th percentile of the tumour type’s non-zero gene medians, so the transform is scaled by the transcriptome the candidate actually comes from rather than by a global constant. It ranges 0.1400 to 0.2400 TPM over 35 cancer types:

from mhcmatch import expression
expression.context_floor(tumor="SKCM")    # 0.1600 TPM
expression.context_floor(tumor="LUAD")    # 0.2000 TPM
expression.context_floor()                # 0.1800 TPM, pooled -- the fallback

The unit does not have to be TPM, because \(c\) is a quantile of the same column and the two cancel — but only while they are the same column. Where a submitted abundance is on some other scale, mhcmatch.expression.batch_scale() estimates the factor by median-of-ratios against the reference and refuses unless the input covers half the context’s expressed genes. A candidate list cannot clear that gate, and should not: a mutation reaches one only where the gene was seen in RNA, so the ratio would measure that conditioning rather than the library.

Picking a context#

--tumor takes a TCGA study code; each is paired to the GTEx normal tissue(s) it is compared against. Nineteen pairings ship, of which 18 currently resolve — CRC is listed in mhcmatch.expression.TUMOR_TISSUE and rejected by mhcmatch.expression.resolve_context(), so use COAD/READ for colorectal:

mhcmatch expression --list-contexts        # every TCGA study code and its matched GTEx tissue(s)
mhcmatch expression PMEL --safety          # where else is this gene expressed?
mhcmatch expression NLVPMVATV --tumor SKCM # has this exact peptide been seen expressed in SKCM?

mhcmatch.expression.resolve_context() also accepts a disease or organ name ("melanoma" and "skin" both reach SKCM) and raises rather than falling back to the pooled reference on anything it cannot place — a silent fallback would return a number from the wrong distribution with no way to tell it had happened.

Two questions, never merged#

The module keys the same table two ways because a ranker asks two different things of it:

key

joined against

answers

gene symbol

a GTEx tissue

Is the source gene transcribed in this lineage? — the ranking read, and what imputes a missing TPM. Also the safety read: a gene expressed everywhere is a toxicity risk, not a target (mhcmatch.expression.safety_profile()).

peptide

a TCGA cancer type

Has this exact neoantigen been seen expressed in this tumour type? Keyed on the peptide deliberately: the TCGA source carries ensp and no ENSP-to-symbol map ships with it, so a gene-level join would be a guess and the peptide-level one is exact.

Missing is encoded, never dropped. mhcmatch.expression.impute() returns the reference value and a flag saying whether it was observed, so a caller carries a missing-indicator column instead of discarding the candidate. That is the standing rule for every partially covered covariate here.

Mouse#

Every function in the module takes species=, defaulting to "human", so no existing call moves. species="mouse" selects different deposits, not a different code path:

from mhcmatch import expression as EX

EX.gene_level("Trp53", tissue="thymus", species="mouse")   # normal rung
EX.gene_level("Tyr", tumor="B16F10", species="mouse")      # tumour rung
EX.context_floor(species="mouse")                          # 0.7174 TPM, pooled
mhcmatch expression Trp53 --species mouse --tissue thymus --safety
mhcmatch expression --species mouse --list-contexts     # tissues and models, listed apart

Four things differ from human, and each is load-bearing:

The normal rung is FANTOM5 CAGE, not GTEx — EBI E-MTAB-3579, 18,830 gene symbols across 35 adult tissues including thymus.

The tumour rung is gene-keyed, where human’s is peptide-keyed. TCGA has per-peptide rows, so human can answer has this exact neoantigen been seen expressed; no mouse deposit has peptide rows anywhere, so the mouse rung answers the gene-level question instead. mhcmatch.rank._expression_for() picks the key from the species rather than growing a branch. lookup(…, tumor=…) covers the six GSE245293 syngeneic models; context_floor and gene_level read the harmonised toil_matrix_mmu.npz — 26,737 genes by 68 contexts on one scale — and reach 33, which is what makes P815, Renca and EMT6 resolvable at all.

Compare floors within a species, never across one. Mouse floors run 0.60 TPM (aorta) to 2.00 (pancreas, stomach, testis) against human’s 0.10–0.40, because CAGE tag density concentrates on fewer genes than RSEM TPM does. A mouse floor read against a human one says nothing.

There is no mouse TUMOR_TISSUE. The tumour-to-matched-normal map is keyed by TCGA study code, so mhcmatch.expression.tissue_floor() with tumor=…, species="mouse" raises, and names context_floor as what to call instead.

Where it is used elsewhere#

Expression is not only a ranking term. mhcmatch.expression.safety_profile() is what the cassette screen consults before a unit is manufactured (Safety, prior evidence, and what goes in the cassette), and mhcmatch.expression.coexpression() is the channel that prices two units firing in the same tissue (Cassette composition).

Full API: mhcmatch.expression in API reference.