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 |
|---|---|
|
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.
|
|
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.
|
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 ( |
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 |
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.