tcren package#

Structure I/O#

tcren.structure.model module#

Lightweight structure data model.

These dataclasses wrap the parsed contents of a PDB/mmCIF file in the shape the rest of the pipeline needs: per-chain residue lists carrying both a 0-based sequential index (matching the legacy mir residue.index) and the original author numbering, plus heavy-atom coordinates for contact computation.

class tcren.structure.model.Atom(name, element, coord)[source]#

Bases: object

A single (heavy) atom.

Parameters:

name (str)
element (str)
coord (ndarray)

name: str#

element: str#

coord: ndarray#

class tcren.structure.model.Residue(seq_index, pdb_index, insertion_code, aa, resname, atoms)[source]#

Bases: object

A polymer residue.

Variables:

seq_index (int) – 0-based sequential index over the chain’s resolved polymer residues (the legacy residue.index); independent of author numbering gaps.
pdb_index (int) – Author residue number (residue.index.pdb).
insertion_code (str) – Author insertion code ('' when absent).
aa (str) – One-letter amino-acid code ('X' for unknown).
resname (str) – Three-letter residue name (HIS, MSE …).
atoms (tuple[tcren.structure.model.Atom, ...]) – Heavy atoms of the residue.

Parameters:

seq_index (int)
pdb_index (int)
insertion_code (str)
aa (str)
resname (str)
atoms (tuple[Atom, ...])

seq_index: int#

pdb_index: int#

insertion_code: str#

aa: str#

resname: str#

atoms: tuple[Atom, ...]#

property ca: ndarray | None#: Cα coordinate, or None if the residue has no Cα atom.

property cb: ndarray | None#: Cβ coordinate, or None if the residue has no Cβ atom (e.g. glycine).

property cb_or_ca: ndarray | None#: Cβ coordinate, falling back to Cα (glycine / missing Cβ); None if neither.

class tcren.structure.model.Chain(chain_id, residues, chain_type=None, chain_supertype=None, allele_info=None, regions=<factory>)[source]#

Bases: object

A polymer chain and its annotations.

Parameters:

chain_id (str)
residues (list[Residue])
chain_type (str | None)
chain_supertype (str | None)
allele_info (str | None)
regions (list[RegionMarkup])

chain_id: str#

residues: list[Residue]#

chain_type: str | None#

chain_supertype: str | None#

allele_info: str | None#

regions: list[RegionMarkup]#

sequence()[source]#

One-letter sequence in residue order.

Return type:: str

by_seq_index(seq_index)[source]#

Return the residue at a given sequential index, or None.

Parameters:: seq_index (int)
Return type:: Residue | None

class tcren.structure.model.RegionMarkup(region_type, start_seq_index, end_seq_index, sequence, residues)[source]#

Bases: object

An annotated region (CDR/FR for TCR, groove regions for MHC).

Parameters:

region_type (str)
start_seq_index (int)
end_seq_index (int)
sequence (str)
residues (list[Residue])

region_type: str#

start_seq_index: int#

end_seq_index: int#

sequence: str#

residues: list[Residue]#

class tcren.structure.model.Structure(pdb_id, chains, complex_species=None, cell_type=None)[source]#

Bases: object

A parsed complex: a set of annotated chains.

Parameters:

pdb_id (str)
chains (list[Chain])
complex_species (str | None)
cell_type (str | None)

pdb_id: str#

chains: list[Chain]#

complex_species: str | None#

cell_type: str | None#

chain(chain_id)[source]#

Return the chain with the given id (raises KeyError if absent).

Parameters:: chain_id (str)
Return type:: Chain

by_type(*types)[source]#

Return chains whose chain_type is in types.

Parameters:: types (str)
Return type:: list[Chain]

tcren.structure.io module#

Parse PDB / mmCIF files into the tcren.structure.model data model.

Accepts plain .pdb/.ent/.cif/.mmcif files, their gzip-compressed forms (.pdb.gz/.cif.gz …), and — for batches — directories or .tar/.tar.gz archives of any of those (see iter_structures()). Structure identifiers are resolved from the file name by structure_id_from_path().

tcren.structure.io.is_structure_file(name)[source]#

True if name is a (optionally gzipped) PDB/mmCIF structure file.

Parameters:: name (str | Path)
Return type:: bool

tcren.structure.io.structure_stem(path)[source]#

The file name with a trailing .gz and the structure extension removed.

Parameters:: path (str | Path)
Return type:: str

tcren.structure.io.structure_id_from_path(path)[source]#

Resolve a structure identifier from a file name.

Strips a trailing .gz and the structure extension, then takes the part before the first _ (so 4x5w_renumbered.cif and 1ao7.pdb.gz and 6uk4_TCRpMHCmodels.pdb all resolve to their PDB id).

Warning

Lossy by design, and not unique for cohorts whose file names encode metadata after an underscore – VDJdb_Model_603_min.pdb and VDJdb_Model_604_min.pdb both give VDJdb. Prefer iter_structures(), which detects that case for a whole set and falls back to the full stem rather than silently collapsing rows.

Parameters:: path (str | Path)
Return type:: str

tcren.structure.io.resolve_structure_ids(paths)[source]#

Map each path to an id: the PDB-id prefix when that is unique over the set, else the stem.

The prefix rule is what makes 4x5w_renumbered.cif come back as 4x5w, and it is right for RCSB-derived files. It is wrong – silently, and in a way that destroys rows downstream – for any set whose names carry metadata after the first underscore. Deciding per SET rather than per file keeps the convenience where it is unambiguous and refuses it where it is not.

Return type:: dict[str, str]

tcren.structure.io.parse_structure(path, pdb_id=None, model=0, keep_hydrogens=True)[source]#

Parse a structure file into a Structure.

Residues are taken in author order; only amino-acid residues (standard or modified, via the extended three→one table) are kept — waters, ions and ligands are dropped. Each kept residue receives a 0-based sequential seq_index per chain, matching the legacy mir residue.index.

Parameters:

path (str | Path) – Path to a .pdb/.ent or .cif/.mmcif file.
pdb_id (str | None) – Structure identifier; defaults to the file stem.
model (int) – Model index to read (default 0 — the first model).
keep_hydrogens (bool) – Keep hydrogen atoms (default True — the legacy mir contact definition counts hydrogens when a structure provides them).

Returns:

The parsed Structure.

Return type:

Structure

tcren.structure.io.import_structure(path, pdb_id=None, model=0, keep_hydrogens=True, trim_c_gene=True, keep_c_gene=False, min_constant_score=80.0)[source]#

Parse a structure and prepare it for interface analysis.

Wraps parse_structure(), records the αβ/γδ cell type from the TCR constant region, and — by default — trims that constant region so downstream analysis works on the variable domains and the interface.

Parameters:

path (str | Path) – as in parse_structure().
pdb_id (str | None) – as in parse_structure().
model (int) – as in parse_structure().
keep_hydrogens (bool) – as in parse_structure().
trim_c_gene (bool) – Trim the TCR constant domain (default True).
keep_c_gene (bool) – Retain the constant domain even if trim_c_gene is set. Use this for molecular-dynamics / FlexPepDock and any workflow that needs the full chain — those depend on the presence of the C-gene.
min_constant_score (float) – Minimum constant-region alignment score to trim on.

Returns:

The imported Structure with cell_type set.

Return type:

Structure

tcren.structure.io.structure_paths(src)[source]#

List structure files for src, sorted.

src may be a single structure file, a directory (scanned for structure files), a glob pattern (models/*.pdb.gz), or a manifest: a .txt/.list/.lst file holding one path per line, # comments and blank lines ignored, relative paths resolved against the manifest’s own directory. Recognises plain and gzipped PDB/mmCIF (.pdb, .cif.gz, …). For archives or streaming, use iter_structures().

Parameters:: src (str | Path)
Return type:: list[Path]

tcren.structure.io.resolve_sources(sources)[source]#

Split a CLI -s spec (or several) into individual sources for iter_structures().

Accepts one spec or an iterable of them, and splits each on commas, so -s a.pdb.gz,b.pdb.gz, a repeated -s, and a shell glob all mean the same thing. Directories, tar archives, globs and manifests are passed through untouched — they are expanded downstream, where the tar streaming lives.

Parameters:: sources (str | Path | Iterable[str | Path])
Return type:: list[str]

tcren.structure.io.iter_structures(src, importer=<function import_structure>, on_error='raise', **kwargs)[source]#

Yield (pdb_id, Structure) for a file, directory, or .tar/.tar.gz archive.

Handles plain and gzipped PDB/mmCIF (.pdb/.cif/.pdb.gz/.cif.gz …); a directory is scanned for those; a tar archive is streamed member-by-member (each member materialised to a temp file so the path-based importer works unchanged). The identifier is resolved per file by structure_id_from_path().

Parameters:

src (str | Path) – structure file, directory, or tar archive.
importer (Callable[[...], Structure]) – per-file parser — import_structure() (default, trims the C-gene) or parse_structure() (parity-pure). Extra kwargs are forwarded to it.
on_error (str) – "raise" (default) or "skip" to ignore files that fail to parse.

Return type:

Iterator[tuple[str, Structure]]

tcren.structure.io.pdb_lines(structure, transform=None, keep_hydrogens=True)[source]#

ATOM/TER/END record lines for structure (optionally coordinate-transformed).

One conformer per atom name per residue (drops duplicate altlocs). transform is an optional coord -> coord callable (e.g. for an oriented frame); identity if None. Author residue numbers + insertion codes are preserved.

Parameters:

structure (Structure)
keep_hydrogens (bool)

Return type:

list[str]

tcren.structure.io.cif_lines(structure, transform=None, keep_hydrogens=True)[source]#

Minimal mmCIF atom_site loop for structure (optionally transformed).

Same atom selection as pdb_lines() (one conformer per atom name per residue). Only the _atom_site category is written — enough to round-trip coordinates + chain/residue identity through the Biopython MMCIF parser, which is all tcren consumes.

Parameters:

structure (Structure)
keep_hydrogens (bool)

Return type:

list[str]

tcren.structure.io.write_pdb(structure, path, transform=None, keep_hydrogens=True)[source]#

Write structure to a PDB file; return the path.

A .gz suffix (foo.pdb.gz) transparently gzip-compresses the output.

Parameters:

structure (Structure)
path (str | Path)
keep_hydrogens (bool)

Return type:

Path

tcren.structure.io.structure_output_path(directory, pdb_id, mmcif=False, compress=False)[source]#

Build an output path <directory>/<pdb_id>.<ext> from format flags.

.pdb by default, .cif if mmcif, with a trailing .gz if compress.

Parameters:

directory (str | Path)
pdb_id (str)
mmcif (bool)
compress (bool)

Return type:

Path

tcren.structure.io.write_structure(structure, path, transform=None, keep_hydrogens=True)[source]#

Format-dispatch writer: PDB or mmCIF, optionally gzipped (by the path suffix).

Parameters:

structure (Structure)
path (str | Path)
keep_hydrogens (bool)

Return type:

Path

Annotation#

tcren.annotation.arda_adapter module#

TCR chain annotation via the arda library.

Extracts each chain’s amino-acid sequence, runs arda’s AIRR annotation, and projects the returned region coordinates (1-based, end-inclusive into the input sequence) back onto structure residues as RegionMarkup. Region names are mapped to the legacy mir vocabulary (FR1/CDR1/…/FR4).

tcren.annotation.arda_adapter.annotate_chain(chain, organism)[source]#

Annotate one chain with arda; return the AIRR record if it is a TCR chain.

Mutates chain in place when arda recognises it as TRA/TRB: sets chain_type, chain_supertype ("TRAB"), allele_info and regions. Returns the arda record (for any locus) or None if arda produced no locus.

Parameters:: organism (str)
Return type:: dict | None

tcren.annotation.arda_adapter.apply_records(chains, by_id)[source]#

Project a cached {chain_id: record} map onto chains in place (no arda call).

Parameters:: by_id (dict[str, dict])
Return type:: None

tcren.annotation.arda_adapter.score_records(chains, by_id)[source]#

(receptor_ids, summed mmseqs2_score) from already-computed records.

Parameters:: by_id (dict[str, dict])
Return type:: tuple[list[str], float]

tcren.annotation.arda_adapter.annotate_chains(chains, organism)[source]#

Annotate a batch of chains in a single arda call; apply records in place.

One mmseqs invocation for all chains (the per-call process/DB overhead dominates the actual alignment, so batching is ~hundreds× faster than per-chain calls). Returns a {chain_id: record} map for chains that had a sequence.

Parameters:: organism (str)
Return type:: dict[str, dict]

tcren.annotation.arda_adapter.annotate_tcr_chains(structure, organism='human')[source]#

Annotate all chains; return ids recognised as antigen-receptor (TCR/BCR) chains.

Parameters:: organism (str)
Return type:: list[str]

tcren.annotation.arda_adapter.annotate_tcr_chains_scored(structure, organism='human')[source]#

Annotate all chains; return (receptor_ids, summed mmseqs2_score).

The summed mmseqs2 alignment score over the receptor chains measures how well the structure’s TCR/BCR chains match this organism’s germline reference — the signal used to pick the correct species when annotating against human vs mouse.

Parameters:: organism (str)
Return type:: tuple[list[str], float]

tcren.annotation.chains module#

Chain classification: receptor chains (via arda), PEPTIDE, and (provisional) MHC.

classify_chains tags MHC chains with the generic type "MHC" so the TCR↔peptide scoring path is complete. Precise MHC sub-typing (MHCa/MHCb/B2M, class I/II, allele) is a separate step in tcren.mhc (annotate_mhc / mhc.mapper).

tcren.annotation.chains.classify_chains(structure, organism='human', peptide_max_len=30, autodetect_species=True, precomputed_records=None)[source]#

Classify every chain of structure in place.

Receptor chains are assigned from arda’s locus call (TCR TRA/TRB and, for BCR inputs, IGH/IGK/IGL); the shortest remaining chains (length ≤ peptide_max_len) become PEPTIDE; longer remaining chains are tagged "MHC".

Parameters:

structure (Structure) – Structure to annotate (mutated in place).
organism (str) – arda organism ("human"/"mouse").
peptide_max_len (int) – Maximum residue count for a chain to be called PEPTIDE.
autodetect_species (bool) – Annotate against both supported species (human and mouse) and keep whichever gives the higher total mmseqs alignment score over the receptor chains. TCR/BCR germlines are organism-specific, so the wrong species scores measurably lower (e.g. mouse BM3.3 scores ~435 vs ~197 under human); this avoids mis-typing a chain under the wrong reference. Ties keep the requested organism. Disable to force organism verbatim.
precomputed_records (dict[str, dict[str, dict]] | None) – Optional {organism: {chain_id: record}} of arda records for this structure’s chains, to reuse instead of calling arda (the batch path in annotate_structure_set() annotates the whole dataset in one mmseqs call per organism and injects the per-structure slices).

Return type:

None

tcren.annotation.cgene module#

αβ vs γδ T-cell classification from the TCR constant region (C-gene).

arda annotates the variable V(D)J region but not the constant domain. When a structure includes an ordered constant domain, aligning each chain to the TCR constant references (TRAC/TRBC1/TRBC2 → αβ; TRGC/TRDC → γδ) identifies the chain (α/β/γ/δ) unambiguously and therefore the cell type. This is authoritative for αβ-vs-γδ and independent of the (occasionally ambiguous, e.g. TRAV/DV) V-gene call. Variable-domain-only chains carry no constant region and yield no call (cell type "unknown").

tcren.annotation.cgene.MIN_CONSTANT_SCORE = 80.0#: Minimum local-alignment score to accept a constant-domain match. The V domain alone scores ~30-45 against any constant; a real constant domain scores in the hundreds.

class tcren.annotation.cgene.ConstantCall(chain_id, gene, chain_class, cell_type, score)[source]#

Bases: object

A TCR constant-region identification for one chain.

Parameters:

chain_id (str)
gene (str)
chain_class (str)
cell_type (str)
score (float)

chain_id: str#

gene: str#

chain_class: str#

cell_type: str#

score: float#

tcren.annotation.cgene.classify_chain_constant(sequence, min_score=80.0)[source]#

Identify the constant region of a single chain sequence, if one is present.

Parameters:

sequence (str)
min_score (float)

Return type:

ConstantCall | None

tcren.annotation.cgene.constant_span(sequence, min_score=80.0)[source]#

Return the (start, end) query span aligning to the best TCR constant region.

start/end are 0-based half-open indices into sequence. Returns None if no constant domain is present (score below min_score). The constant region is C-terminal, so callers trim residues with seq_index >= start.

Parameters:

sequence (str)
min_score (float)

Return type:

tuple[int, int] | None

tcren.annotation.cgene.classify_constants(structure, min_score=80.0)[source]#

Identify the constant region of every chain that has one.

Parameters:: min_score (float)
Return type:: list[ConstantCall]

tcren.annotation.cgene.cell_type(structure, min_score=80.0)[source]#

Return "ab", "gd" or "unknown" from the constant regions present.

γδ wins if any γ/δ constant is found; otherwise αβ if any α/β constant is found; otherwise "unknown" (no ordered constant domain — e.g. variable-only chains).

Parameters:: min_score (float)
Return type:: str

MHC mapping#

tcren.mhc.imgt module#

Download and parse MHC allele references (IMGT/HLA + UniProt mouse H-2 + B2M).

Produces MhcAllele records labelled with species, MHC class and chain role (MHCa = class-I heavy or class-II alpha; MHCb = class-II beta; B2M). Human alleles come from IMGT/HLA (hla_prot.fasta); mouse H-2 and beta-2-microglobulin come from reviewed UniProt entries.

class tcren.mhc.imgt.MhcAllele(allele, locus, mhc_class, chain_role, species, sequence)[source]#

Bases: object

A reference MHC allele sequence with its functional labels.

Parameters:

allele (str)
locus (str)
mhc_class (str)
chain_role (str)
species (str)
sequence (str)

allele: str#

locus: str#

mhc_class: str#

chain_role: str#

species: str#

sequence: str#

tcren.mhc.imgt.download_human(cache_dir, force=False)[source]#

Download IMGT/HLA hla_prot.fasta into cache_dir.

Parameters:

cache_dir (Path)
force (bool)

Return type:

Path

tcren.mhc.imgt.download_mouse(cache_dir, force=False)[source]#

Download mouse H-2 / B2m and human B2M into cache_dir.

Parameters:

cache_dir (Path)
force (bool)

Return type:

tuple[Path, Path]

tcren.mhc.imgt.parse_human(path)[source]#

Parse IMGT/HLA, keeping classical loci collapsed to two-field resolution.

Parameters:: path (Path)
Return type:: list[MhcAllele]

tcren.mhc.imgt.parse_mouse(mouse_path, human_b2m_path)[source]#

Parse reviewed mouse H-2 / B2m and human B2M from UniProt FASTA headers.

Parameters:

mouse_path (Path)
human_b2m_path (Path)

Return type:

list[MhcAllele]

tcren.mhc.reference module#

Build and load the curated MHC reference shipped under database/mhc/.

The committed reference is a single FASTA (alleles.aa.fasta) whose headers encode the metadata (allele|locus|mhc_class|chain_role|species) plus a metadata.tsv mirror. The mmseqs search index is built on demand into a gitignored cache (mirroring arda’s commit-FASTA / build-index-on-demand split).

tcren.mhc.reference.build(species=('human', 'mouse'), cache_dir=PosixPath('/home/runner/work/tcren/tcren/data/mhc_cache'), out_dir=PosixPath('/home/runner/work/tcren/tcren/database/mhc'), force_download=False)[source]#

Download, curate and write the committed MHC reference.

Parameters:

species (tuple[str, ...]) – Which species to include.
cache_dir (Path) – Where raw downloads are cached (gitignored).
out_dir (Path) – Where the curated alleles.aa.fasta + metadata.tsv are written.
force_download (bool) – Re-download even if cached files exist.

Returns:

Path to the written alleles.aa.fasta.

Return type:

Path

tcren.mhc.reference.reference_fasta(out_dir=PosixPath('/home/runner/work/tcren/tcren/database/mhc'))[source]#

Path to the committed reference FASTA (raise if the reference is not built).

Parameters:: out_dir (Path)
Return type:: Path

tcren.mhc.reference.reference_db(cache_dir=PosixPath('/home/runner/work/tcren/tcren/data/mhc_cache'))[source]#

Path to a compiled, pre-indexed mmseqs DB of the allele reference (built once, cached).

mmseqs easy-search otherwise rebuilds the target DB and its k-mer prefilter index from the ~28k-allele FASTA on every call. Caching createdb saves little; the dominant cost is the prefilter index, so we also run createindex once. Reusing this DB cuts a single-structure MHC search from ~4.5 s to ~0.9 s. Built into the gitignored data/mhc_cache when missing or older than the FASTA.

The build is serialized through arda._locking.build_lock(): tcren is routinely run concurrently against the same cache (one process per SLURM-array task / Nextflow sample), and an unguarded createdb + createindex into the shared path would let every other process search a half-written index. The createindex marker is written last, so its freshness gates completeness.

Parameters:: cache_dir (Path)
Return type:: Path

tcren.mhc.reference.parse_header(header)[source]#

Parse a reference FASTA header back into its metadata fields.

Parameters:: header (str)
Return type:: dict[str, str]

tcren.mhc.mapper module#

Map a structure’s MHC chains to allele / class / role via mmseqs.

Searches each not-yet-typed chain against the curated MHC reference and assigns the best hit’s class (MHCI/MHCII), chain role (MHCa/MHCb/B2M), locus and allele. Class is reconciled across the complex (B2M ⇒ class I; a class-II beta chain ⇒ class II).

class tcren.mhc.mapper.MhcCall(chain_id, chain_role, mhc_class, allele, locus, species, identity, bits, qstart, qend, tstart, tend, cigar)[source]#

Bases: object

Result of mapping one chain to the MHC reference.

Parameters:

chain_id (str)
chain_role (str)
mhc_class (str)
allele (str)
locus (str)
species (str)
identity (float)
bits (float)
qstart (int)
qend (int)
tstart (int)
tend (int)
cigar (str)

chain_id: str#

chain_role: str#

mhc_class: str#

allele: str#

locus: str#

species: str#

identity: float#

bits: float#

qstart: int#

qend: int#

tstart: int#

tend: int#

cigar: str#

tcren.mhc.mapper.map_mhc(structure, sensitivity=5.7)[source]#

Map the structure’s MHC chains against the curated reference.

Parameters:

structure (Structure) – A structure whose TCR/peptide chains are already typed.
sensitivity (float) – mmseqs search sensitivity.

Returns:

One MhcCall per chain that produced a reference hit.

Return type:

list[MhcCall]

tcren.mhc.mapper.calls_from_hits(candidates, best, key=None)[source]#

Build reconciled MhcCall`s for ``candidates` from precomputed mmseqs hits.

key(chain) -> str maps a candidate chain to its key in best (default the chain id; a batched search uses "<struct_idx>|<chain_id>"). Lets one mmseqs search over many structures’ chains be sliced back per structure — no per-structure mmseqs call.

Parameters:: best (dict[str, dict])
Return type:: list[MhcCall]

tcren.mhc.mapper.apply_mhc_calls(structure, calls)[source]#

Write MHC calls onto the structure’s chains in place.

Parameters:

structure (Structure)
calls (list[MhcCall])

Return type:

None

tcren.mhc.domains module#

Canonical MHC groove region definitions.

Loads the bundled mhc_canonical.json: for each "<class>|<role>" key it holds a canonical mature chain sequence and the 0-based positions of each groove region (HELIX_A1/HELIX_A2 for class I, HELIX_A1/HELIX_B1 for class II, and GROOVE_FLOOR). Region boundaries follow established mature-numbering ranges for the α1/α2 (class I) and α1/β1 (class II) groove domains. Regions are projected onto query chains in tcren.mhc.regions.

tcren.mhc.domains.canonical_groove()[source]#

Return the bundled canonical groove definitions.

Return type:: dict

tcren.mhc.domains.groove_for(mhc_class, chain_role)[source]#

Canonical groove definition for a (class, role), or None if none exists.

Parameters:

mhc_class (str)
chain_role (str)

Return type:

dict | None

tcren.mhc.regions module#

Project canonical MHC groove regions onto a structure’s MHC chains.

Each MHC chain is aligned (global, BLOSUM62) to the canonical chain for its class/role; the canonical region positions are then mapped through the alignment onto the chain’s residues, producing RegionMarkup entries (HELIX_A1, HELIX_A2/HELIX_B1, GROOVE_FLOOR) in the same schema as the TCR region markup.

tcren.mhc.regions.partition_chain(chain, mhc_class, chain_role)[source]#

Return groove RegionMarkups for one MHC chain (empty for B2M / unknown roles).

Parameters:

chain (Chain)
mhc_class (str)
chain_role (str)

Return type:

list[RegionMarkup]

tcren.mhc.regions.partition_mhc(structure, calls)[source]#

Assign groove regions to every mapped MHC chain in the structure (in place).

Parameters:

structure (Structure)
calls (list[MhcCall])

Return type:

None

tcren.mhc.regions.annotate_mhc(structure)[source]#

Map and partition the MHC chains of an (already chain-typed) structure.

Returns the MhcCall list and, in place, sets each MHC chain’s type (MHCa/MHCb/B2M), class supertype, allele and groove regions.

Parameters:: structure (Structure)
Return type:: list[MhcCall]

tcren.mhc.regions.annotate_mhc_batch(structures, sensitivity=5.7, threads=1)[source]#

MHC-annotate many (chain-typed) structures with a SINGLE mmseqs search.

Gathers every candidate MHC chain across all structures, runs one easy_search (mmseqs parallelises internally — no Python threads, no per-structure call), then slices the hits back and applies the calls + groove partitioning to each structure in place. This is the batched equivalent of calling annotate_mhc() per structure, for dataset-scale work.

Parameters:

structures (list[Structure])
sensitivity (float)
threads (int)

Return type:

None

tcren.mhc.linker module#

Detect and split covalently linked (single-chain) peptides via MHC alignment.

Engineered single-chain pMHC constructs fuse the peptide to an MHC chain through a flexible (usually Gly/Ser-rich) linker, so the peptide is not a separate chain. Aligning each chain to the MHC reference reveals this: the MHC domain aligns, leaving an unaligned terminal segment that — after stripping the linker — is the peptide. This module provides the alignment check and a splitter that lifts such peptides into their own chain.

No covalently linked peptides occur in the bundled TCR3D / PDB datasets (all conventional, separate-chain complexes); this is robustness for engineered and predicted structures.

class tcren.mhc.linker.MhcAlignmentCheck(chain_id, best_ref, score, query_start, query_end, n_term_extra, c_term_extra)[source]#

Bases: object

Result of aligning a chain onto the MHC reference (a chain-identity check).

Parameters:

chain_id (str)
best_ref (str)
score (float)
query_start (int)
query_end (int)
n_term_extra (int)
c_term_extra (int)

chain_id: str#

best_ref: str#

score: float#

query_start: int#

query_end: int#

n_term_extra: int#

c_term_extra: int#

property is_mhc: bool#

tcren.mhc.linker.check_against_mhc(chain)[source]#

Align a chain onto the MHC reference and report coverage / terminal extensions.

Parameters:: chain (Chain)
Return type:: MhcAlignmentCheck

tcren.mhc.linker.detect_linked_peptide(chain, min_len=7, max_len=25)[source]#

Return the residues of a peptide fused to an MHC chain, or None.

Looks for a peptide-length segment (after stripping an adjacent Gly/Ser linker) at the N- or C-terminus of a chain whose remainder aligns to the MHC reference.

Parameters:

chain (Chain)
min_len (int)
max_len (int)

Return type:

list | None

tcren.mhc.linker.split_linked_peptides(structure, peptide_chain_id='p')[source]#

Split covalently linked peptides off their MHC chains, in place.

For each chain carrying a fused peptide, the peptide residues are removed and added as a new PEPTIDE chain. Returns the list of chain ids that were split (empty if none).

Parameters:

structure (Structure)
peptide_chain_id (str)

Return type:

list[str]

tcren.mhc.pseudo module#

MHC pseudosequence (MPS) annotation.

NetMHCpan defines a 34-residue “pseudosequence” per allele — the polymorphic groove positions that contact the peptide (class I: α1/α2 of MHCa; class II: α1 of MHCa + β1 of MHCb). The committed mhci_pseudo.fa / mhcii_pseudo.fa (see scripts/build_pseudo_fasta.py) hold the unique pseudosequences.

annotate_pseudo() adds an MPS region to a chain-typed + MHC-annotated structure, on demand. The 34 pseudo positions are scattered along the chain (not a contiguous motif), so an mmseqs/local search can’t find them — there is no shared k-mer to seed on. Instead we thread each candidate 34-mer through the chain with a fitting alignment (gaps in the chain are free, the pseudosequence may not gap), which recovers the positions because NetMHCpan lists them N→C. The best-scoring pseudosequence is chosen (one hit), and its identically-matched residues are marked — across MHCa only for class I, split across MHCa+MHCb for class II, never β2m. Scoring all ~5.4k pseudosequences this way is ~0.1 s, so no prebuilt index is needed.

tcren.mhc.pseudo.annotate_pseudo(structure)[source]#

Add an MPS region to each groove chain from the best-matching pseudosequence.

structure must already be chain-typed + MHC-annotated. Returns the chosen pseudosequence id (or None if there is no MHC). The best hit is selected once over the class groove sequence (MHCa for class I; MHCa+MHCb for class II) and its residues marked per chain.

Parameters:: structure (Structure)
Return type:: str | None

Contacts#

tcren.contacts.geometry module#

Atom-level contact and Cα-distance computation.

Ports the legacy mir compute-pdb-contacts / compute-pdb-geom steps using a scipy.spatial.cKDTree for the all-atom neighbour search.

tcren.contacts.geometry.all_atom_contacts(structure, cutoff=5.0, count_atoms=False)[source]#

Closest inter-chain atom contact for each residue pair within cutoff Å.

For every pair of residues on different chains that have at least one heavy-atom pair within cutoff (inclusive, matching the legacy dist <= 5), the row with the minimum atom–atom distance is kept.

Parameters:

structure (Structure) – The (parsed) structure.
cutoff (float) – Contact distance threshold (Å, inclusive).
count_atoms (bool) – When True add an extra n_atom_contacts column = the count of heavy-atom pairs within cutoff for that residue pair (always >= 1 for any kept row, and >= the single closest-atom row this function keeps). Default False keeps the schema and every value byte-identical to the legacy output.

Returns:

chain.id.from, residue.index.from, chain.id.to, residue.index.to, residue.aa.from, residue.aa.to, atom.from, atom.to, dist (plus n_atom_contacts when count_atoms is set). Each unordered residue pair appears once, in (chain.id, residue.index) lexicographic order.

Return type:

Columns

tcren.contacts.geometry.representative_atom_contacts(structure, kind='ca', cutoff=12.0)[source]#

Inter-chain residue contacts by a single representative atom per residue.

kind="ca" uses Cα (default cutoff 12 Å); kind="cb" uses Cβ with a glycine/ missing-Cβ fallback to Cα (default cutoff 8 Å). Mirrors all_atom_contacts()’ residue-pair schema (atom.from/atom.to carry the representative atom kind).

Parameters:

structure (Structure)
kind (str)
cutoff (float)

Return type:

DataFrame

tcren.contacts.geometry.ca_distance_matrix(structure)[source]#

Pairwise Cα–Cα distance matrix over all residues with a Cα atom.

Returns:: (matrix, keys) where matrix[a, b] is the Cα distance and keys[a] is the (chain_id, seq_index) of row/column a.
Parameters:: structure (Structure)
Return type:: tuple[ndarray, list[tuple[str, int]]]

tcren.contacts.definitions module#

Flexible multi-threshold contact definition.

Beyond the legacy single 5 Å all-atom contact (the TCRen parity default, d1), this adds two coarser residue-level layers: d2 over Cβ atoms (Cα for glycine) and d3 over Cα atoms. The layers nest from tight side-chain proximity to backbone neighbourhood, giving the 2D maps and scoring a tunable contact model without changing the 5 Å default.

class tcren.contacts.definitions.ContactDefinition(d1=5.0, d2=8.0, d3=12.0)[source]#

Bases: object

Three nested contact thresholds (Å).

Variables:

d1 (float) – closest heavy-atom distance (all-atom contact).
d2 (float) – closest Cβ distance (Cα for glycine / missing Cβ).
d3 (float) – closest Cα distance.

Parameters:

d1 (float)
d2 (float)
d3 (float)

d1: float#

d2: float#

d3: float#

tcren.contacts.definitions.multi_contacts(structure, definition=ContactDefinition(d1=5.0, d2=8.0, d3=12.0))[source]#

Stacked inter-chain residue contacts across the three layers.

Returns the union of the d1/d2/d3 residue-pair tables with a layer column ("d1"/"d2"/"d3") and the layer’s distance. A residue pair can appear in several layers; callers filter by layer as needed.

Parameters:

structure (Structure)
definition (ContactDefinition)

Return type:

DataFrame

tcren.contacts.table module#

Annotate and symmetrise the residue contact table.

Joins per-residue annotations (chain type, region type, region start, amino acid) onto the raw contacts and mirrors the R rbind(contacts, swapped) symmetrisation, yielding the fully annotated, bidirectional contact table the contact map is built from.

tcren.contacts.table.residue_annotation(structure)[source]#

Per-residue annotation table for joining onto contacts.

Columns: chain.id, residue.index, chain.type, chain.supertype, region.type, region.start, residue.aa. region.type/region.start are null for residues without a region annotation.

Parameters:: structure (Structure)
Return type:: DataFrame

tcren.contacts.table.symmetrize(contacts)[source]#

Return contacts plus their from/to-swapped mirror (R rbind semantics).

Parameters:: contacts (DataFrame)
Return type:: DataFrame

tcren.contacts.table.tidy_contacts(structure, cutoff=5.0, count_atoms=False)[source]#

Symmetrised, fully annotated contact table for a structure.

Each inter-chain residue contact appears in both directions, with chain type, region type, region start and amino acid attached on both the from and to sides — the input to tcren.contactmap.ContactMap.

When count_atoms is set, the n_atom_contacts per-residue-pair heavy-atom count is carried through (it is symmetric, so it survives the from/to swap unchanged). Default False keeps the table byte-identical to the legacy output.

Parameters:

structure (Structure)
cutoff (float)
count_atoms (bool)

Return type:

DataFrame

tcren.contactmap module#

Residue-level contact map and interface partitioning.

A ContactMap wraps the annotated, symmetrised contact table and exposes the three biological interfaces (TCR↔peptide, TCR↔MHC, peptide↔MHC). The TCR↔peptide interface is the central object for scoring and reproduces the schema of data/contact_maps_PDB.csv once chains and regions are annotated.

tcren.contactmap.TCR_REGIONS: dict[str, set[str] | None] = {'all': None, 'cdr': {'CDR1', 'CDR2', 'CDR3'}, 'cdr+fr': {'CDR1', 'CDR2', 'CDR3', 'FR1', 'FR2', 'FR3'}}#: TCR region sets selectable on the from (TCR) side of an interface. "all" (no filter) is the default and reproduces the legacy behaviour byte-for-byte; "cdr" keeps only the three CDRs; "cdr+fr" adds the FR1–FR3 framework regions (FR4 excluded).

class tcren.contactmap.ContactMap(pdb_id, contacts, peptide_length=None)[source]#

Bases: object

Annotated, symmetrised residue contacts for one structure.

Parameters:

pdb_id (str)
contacts (DataFrame)
peptide_length (int | None)

pdb_id: str#

contacts: DataFrame#

peptide_length: int | None#

classmethod from_structure(structure, cutoff=5.0, count_atoms=False)[source]#

Build a contact map from an (annotated) structure.

When count_atoms is set, the annotated table carries an n_atom_contacts per-residue-pair heavy-atom count column (needed for atomic-weighted scoring). Default False keeps the contacts table byte-identical to the legacy output.

Parameters:

structure (Structure)
cutoff (float)
count_atoms (bool)

Return type:

ContactMap

interface(which, tcr_regions='all')[source]#

Return the contacts of one interface with within-region positions.

Parameters:

which (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – "tcr_peptide", "tcr_mhc" or "peptide_mhc".
tcr_regions (str) – which TCR regions to keep on the from (TCR) side — "all" (default, no filter; legacy behaviour), "cdr" (CDR1–CDR3 only), or "cdr+fr" (CDR1–CDR3 plus FR1–FR3). Has no effect on "peptide_mhc" (no TCR side).

Returns:

Filtered contacts with added pos.from/pos.to columns.

Return type:

DataFrame

tcr_peptide()[source]#

Convenience accessor for the TCR↔peptide interface.

Return type:: DataFrame

to_csv(path)[source]#

Write the full annotated contact table to CSV.

Parameters:: path (str | Path)
Return type:: None

class tcren.contactmap.ModeCentroid(apex_x, y_alpha, y_beta, sigma_sum, footprint_width_alpha, footprint_width_beta, n_contacts_alpha, n_contacts_beta)[source]#

Bases: object

Length-invariant binding-mode centroid in graphon coordinates (see binding_mode()).

Parameters:

apex_x (float)
y_alpha (float)
y_beta (float)
sigma_sum (float)
footprint_width_alpha (float)
footprint_width_beta (float)
n_contacts_alpha (int)
n_contacts_beta (int)

apex_x: float#: contact-weighted mean loop position, pooled over both loops (~0.47)

y_alpha: float#: contact-weighted mean peptide position read by CDR3alpha (~0.41)

y_beta: float#: contact-weighted mean peptide position read by CDR3beta (~0.62)

sigma_sum: float#: y_alpha + y_beta — the sigma involution reads ~1.0

footprint_width_alpha: float#

footprint_width_beta: float#

n_contacts_alpha: int#

n_contacts_beta: int#

tcren.contactmap.registered_map(structure, *, grid=8, target='peptide', metric='distance', cutoff=5.0)[source]#

CDR3×``target`` Cα map resampled onto a fixed grid``×``grid graphon grid.

Bilinearly resamples the variable-length CDR3α/β × peptide Cα distance map onto normalised coordinates (i/(L+1), j/(M+1)) so loops of different length are directly comparable, averageable and linear-model-ready. Featurisation, not a binder score — the discriminative signal is epitope-identity provenance and collapses to chance under epitope matching (algem monograph E13).

Parameters:

structure – a chain-typed structure (classify_chains run).
grid (int) – output grid size G per axis.
target (str) – "peptide" (only target currently supported).
metric (str) – "distance" (Cα distances) or "contact" (binary at cutoff Å).
cutoff (float) – contact threshold in Å when metric="contact".

Returns:

(2, G, G) stacking the α then β blocks when both loops resolve; (G, G) when only one does; None when neither loop or the peptide is available.

tcren.contactmap.binding_mode(structure, *, contact=5.0)[source]#

Length-invariant binding-mode centroid of the CDR3α/β × peptide contact map.

Reduces each loop’s contact map to graphon-coordinate centroids: apex_x = i/(L+1) (where on the loop the contacts sit, ~0.47 = the apex) and y = j/(M+1) (where on the peptide each loop reads). The σ involution shows up as y_alpha + y_beta ≈ 1 (α reads the N-terminal half, β the C-terminal). The two loops are never pooled for y — the split is the point — but apex_x is contact-weighted over both.

Characterises the mode, does not discriminate specificity — the gross mode is universal across HLA-A*02:01 9-mers, so do not use it as a same-epitope classifier (algem monograph E19).

Parameters:

structure – a chain-typed structure.
contact (float) – closest-Cα contact threshold in Å (the reference used an 8 Å Cα proxy; 5 Å is tighter).

Returns:

A ModeCentroid, or None if neither loop makes >= 3 contacts with the peptide.

Return type:

ModeCentroid | None

tcren.contact_types module#

Chemical typing of interface contacts — a DSSP-style annotation layer for TCR:pMHC contact maps.

Classifies each atom–atom contact in a ContactMap interface into one chemical type from heavy-atom geometry alone (no hydrogens, no external DSSP binary — models and many crystals lack explicit H), by priority:

salt_bridge — cationic (Lys/Arg/His N) ↔ anionic (Asp/Glu O), ≤ 4.0 Å
hydrogen_bond— two polar N/O atoms, ≤ 3.5 Å
aromatic — two aromatic-ring atoms of aromatic residues (Phe/Tyr/Trp/His), ≤ 5.0 Å
hydrophobic — two carbon atoms of apolar residues, ≤ 4.5 Å
other — anything else within the contact-map cutoff

This replaces ad-hoc one-off H-bond counters: contact_type_counts() gives per-type contact and residue-pair counts (the documented, reproducible source of an n_hbond-style feature).

tcren.contact_types.classify_contacts(interface_df)[source]#

Return interface_df with an added contact.type column (one of TYPES per row).

Parameters:: interface_df (DataFrame) – an interface frame from ContactMap.interface(), carrying residue.aa.from/to, atom.from/to and dist.
Returns:: The same frame with a contact.type string column.
Return type:: DataFrame

tcren.contact_types.contact_type_counts(cm, interface='tcr_peptide', tcr_regions='all')[source]#

Per-type contact counts + distinct residue-pair counts for one interface.

Parameters:

cm – a ContactMap.
interface (str) – interface name ("tcr_peptide", "tcr_mhc", "peptide_mhc").
tcr_regions (str) – passed through to ContactMap.interface().

Returns:

Mapping with n_<type> (atom-pair contacts of each type) and pairs_<type> (distinct residue–residue pairs having ≥1 contact of that type), e.g. pairs_hydrogen_bond is the documented n_hbond feature.

Return type:

dict[str, int]

tcren.geometry module#

Closed-form CDR3 Ω-loop reachability — the feasibility filter.

The CDR3 loop is a chain of n = L + 1 virtual Cα–Cα bonds of length BOND_CA, pinned at its neck separation R (the cdr3{a,b}_ext frame descriptor). Pure geometry then bounds how far the loop’s apex can reach from the neck axis, independent of sequence or energy:

reach_max() — the maximum apex reach for a loop of length L pinned at R.
reachability_floor() — the shortest loop that can span a target at distance d.
span_saturation() — the fraction of maximal reach a structure actually uses, reach / reach_max; > 1 is geometrically impossible and real repertoires sit well below 1 (~0.3–0.5 at CDR3 length 10–20), i.e. with large conformational slack.

These are feasibility / mode descriptors, not binder classifiers — backbone Cα geometry does not discriminate binders (every apparent signal was a provenance / epitope confound). Bill them as structure -> geometry featurisation only.

Ported from model/src/Feasibility.jl of the 2026-tcren-algem monograph.

tcren.geometry.BOND_CA = 3.8#: Virtual Cα–Cα bond length (Å) used for the loop chain.

tcren.geometry.reach_max(length, neck, bond=3.8)[source]#

Maximum apex reach sqrt(((L+1)b/2)^2 - (R/2)^2) of a length-length loop pinned at neck.

Parameters:

length (int) – CDR3 loop length L (number of residues).
neck (float) – neck separation R in Å (the ext frame descriptor).
bond (float) – virtual Cα–Cα bond length (default BOND_CA).

Returns:

Maximum reach in Å, or 0.0 if the loop is too short to span the neck.

Return type:

float

tcren.geometry.reachability_floor(distance, neck, bond=3.8)[source]#

Shortest loop length that can reach a target at distance from a neck of separation neck.

Inverts reach_max(): the minimum L with reach_max(L, neck) >= distance.

Parameters:

distance (float)
neck (float)
bond (float)

Return type:

int

tcren.geometry.span_saturation(structure)[source]#

Per-loop reach / reach_max for cdr3a and cdr3b of a chain-typed structure.

Composes the existing CDR3 frame descriptors: reach (loop-centroid distance from the groove origin) and ext (neck separation R), with L the resolved CDR3 length. > 1 flags an infeasible pose (the loop cannot reach that far); real structures sit ~0.3–0.5. Returns NaN for a loop whose frame or CDR3 span is undefined.

Note

The structure must be chain-typed (classify_chains) so its CDR3 regions are populated.

Return type:: dict[str, float]

tcren.geometry.cdr3_internal_coords(structure, loop='cdr3b')[source]#

Virtual-bond internal coordinates (θ, τ, neck R) of a CDR3 loop’s Cα trace.

A thin structure -> coordinates extractor: the loop’s model (torsion sampling, closure, ensemble entropy) lives in the algem monograph’s Julia code and does not migrate — tcren only hands off the intrinsic internal coordinates of the observed loop. Reuses the CDR3 region markup (classify_chains) for the span and the shared _dihedral primitive for torsions.

Parameters:

structure – a chain-typed structure.
loop (str) – "cdr3a" (TRA) or "cdr3b" (TRB).

Returns:

A LoopInternalCoords, or None if the loop has fewer than 4 resolved Cα (a torsion needs 4 points).

Return type:

LoopInternalCoords | None

class tcren.geometry.LoopInternalCoords(bonds, angles, torsions, neck, n_ca)[source]#

Bases: object

Virtual-bond internal coordinates of a CDR3 Cα trace (the input to the Julia loop model).

A Cα chain of N points has N-1 bonds, N-2 bond angles (interior vertices) and N-3 pseudo-torsions — the inverse of a discrete-Frenet reconstruction. neck is the end-to-end span |p_last - p_first| (the ext frame descriptor). Angles and torsions are in radians.

Parameters:

bonds (object)
angles (object)
torsions (object)
neck (float)
n_ca (int)

bonds: object#: (N-1,) virtual Cα–Cα bond lengths, Å

angles: object#: (N-2,) Cα–Cα–Cα pseudo bond angles θ, radians

torsions: object#: (N-3,) Cα pseudo-dihedrals τ, radians in (−π, π]

neck: float#: end-to-end span R = |p_last − p_first|, Å

n_ca: int#

tcren.cohort module#

Cohort-relative recognition scores — the recommended, fit-free screening layer.

Prefer these over the fitted tcren.binder.binder_score() (p_bind) and tcren.recognition.forced_pose_score() (p_forced). Those carry trained coefficients; the functions here carry none — no logistic, no fit, no training set — so they cannot leak, cannot go stale, and there is nothing to re-derive. The benchmark repo settled the trade-off empirically (ledger C24/C25/C26):

q_score() matches or beats the fitted p_bind and, unlike it, generalises across cohorts — a logistic trained on one cohort learns that cohort’s epitope composition and does not transfer, whereas Q has nothing to transfer. With ipTM it reproduces the headline synergy fit-free: z(ipTM) + z(Q) reaches macro ROC 0.83 on TCRvdb against ipTM’s 0.79.
strain_z() grades pose forcedness (crystal < AF-real < AF-decoy) reproducibly, unlike FORCED_POSE_MODEL whose training rows are lost.

They are cohort-relative by default: each standardizes a feature over the set being ranked. For a candidate set, score the whole batch together (tcren recognize over a directory). For a single structure, or a small/heterogeneous user set where the batch is not a fair reference, pass reference=native_reference() (with features=Q_FEATURES_GEOM): the descriptors are then standardized against the shipped Native2026 crystal manifold, so Q is defined for one structure and transfers across inputs. The descriptors are counts and bounded ratios (mildly non-normal), so method="rank" gives a robust, assumption-free percentile standardization; on the benchmarks it agrees with the default z to ρ≈0.98. The division of labour is scores in tcren, evaluation (ROC/PR/CI) downstream.

All functions take the table tcren recognize --full emits (a mapping of column name to sequence, a polars/pandas frame, or a dict of arrays) and return one value per row.

Sign convention: every term is oriented so that higher = more binder-like for q_score(), and higher = more forced/strained for strain_z().

Note

phi_bind() is deprecated — every term it adds to Q lowers ranking accuracy (benchmark ledger C19b), and its z(-pitch) term is both below chance on its own and derived from an AlphaFold-contaminated angle. Use q_score().

tcren.cohort.zscore(x, reference=None, method='z')[source]#

NaN-aware standardization. reference calibrates against another cohort.

Passing reference is what makes strain_z() crystal-calibrated: the mean and sd come from the crystallographic ensemble, so the score reads ~0 on crystals by construction and grows as a pose departs from the natural manifold. Without it, a cohort of uniformly forced poses would standardize to zero mean and the shift would be invisible.

Parameters:

x – values to standardize.
reference – cohort defining the location/scale; defaults to x itself (cohort-relative).
method – "z" (mean/sd, the default) or "rank" — the percentile of each x against the reference, mapped to [-1, 1] (2·percentile − 1). "rank" is scale-free and makes no normality assumption, so it is the robust choice for the bounded/count descriptors of Q (chain balance, H-bond and contact counts are not normal). On the benchmarks z and rank agree to Spearman ρ≈0.98 and differ by <0.005 AUROC, so z is kept as the default; use rank when a heavy-tailed user descriptor could distort the mean/sd.

Return type:

ndarray

tcren.cohort.q_score(table, reference=None, features=('burial', 'n_pep_contacted', 'chain_balance', 'n_hbond'), method='z', decorrelate=True)[source]#

Interface-quality score Q — fit-free, single-structure-capable; the default binder score.

The default is the directional, decorrelated one-class score over k = len(features) terms

\[Q(x) \;=\; z(x)^{\top}\, C^{-1}\, \mathbf{1}, \qquad z(x)_k = \frac{d_k(x)-\mu_k}{\sigma_k},\]

where each descriptor is standardized against the native crystal reference (reference, default native_reference() — its $\mu_k,\sigma_k$), $C$ is the native descriptor correlation matrix, and $\mathbf 1$ is the biophysical every-descriptor-higher-is-better direction. Whitening by $C^{-1}$ stops correlated descriptors from double-counting. The score carries no fitted coefficient (only the native covariance is estimated; $\mathbf 1$ is fixed), needs no negative set, is calibrated on natives so it transfers across inputs, and is defined for a single structure. It reduces to the equal-weight mean when the descriptors are uncorrelated ($C=I$). See the manuscript Methods (S Scores).

decorrelate=False recovers the legacy equal-weight mean $Q=\frac1k\sum_k z(d_k)$.

Parameters:

table – the tcren recognize --full table (dict / pandas / polars).
reference – cohort defining $\mu,\sigma,C$. None uses native_reference() when decorrelate (so the covariance is defined for any input, incl. one structure); with decorrelate=False it means cohort-relative (the table itself).
features – the k descriptors. Default Q_FEATURES_GEOM (k=4, geometry only) — the validated default. Adding the pp_combo energy term (k=5, Q_FEATURES) degrades ranking on generated poses, because that term inverts on forced ones (ledger C42).
method – per-descriptor standardization, "z" (default) or "rank" — see zscore().
decorrelate – whiten by the native covariance and project onto $\mathbf 1$ (default); else the equal-weight mean.

Return type:

ndarray

tcren.cohort.q_iptm(table, iptm, reference=None, features=('burial', 'n_pep_contacted', 'chain_balance', 'n_hbond'), decorrelate=True)[source]#

Fit-free synergy score z(ipTM) + z(Q) — the interface-quality score composed with the generator’s own confidence.

Q (interface geometry) and the AlphaFold/TCRmodel2 ipTM are near-orthogonal (they fail in different pose regimes, benchmark ledger C26/C35), so their standardized sum out-ranks either alone: macro ROC 0.83 / PR 0.83 on TCRvdb vs ipTM 0.79, and on well-modelled epitopes it beats raw-AF ipTM on both metrics (ledger C42). Both terms are standardized over the same candidate set, so pass an iptm vector aligned row-for-row with table. Use features=Q_FEATURES_GEOM for the geometry-only Q_geom variant that is robust to the forced-pose energy inversion.

Parameters:

table – the tcren recognize --full table (dict / pandas / polars).
iptm – per-structure ipTM, aligned to table rows. Structures whose ipTM is missing (NaN) fall back to z(Q) alone, so the score always ranks; an all-missing iptm returns plain z(Q) — i.e. rank by the model geometry when no generator confidence is available.
reference – optional cohort to standardize against (see zscore()).
features – descriptors for Q; defaults to the five Q_FEATURES.

Return type:

ndarray

tcren.cohort.f_score(table, reference=None, terms=('F_tcr_pep', 'F_tcr_mhc'))[source]#

Binder-oriented TCRen contact energy F = z(-(F_tcr_pep + F_tcr_mhc)) — the chemistry channel.

The standardized, sign-flipped sum of the F_TERMS contact energies, so higher = more binder-like and it is on the same z-scale as q_score(). Unlike Q (interface geometry), F reads the actual contact chemistry — and unlike Q it is pose-conditional: it works on well-modelled poses and inverts on forced ones (benchmark ledger C27/C42). Do not use it unconditioned on pose quality; see F_TERMS and q_f().

Cohort-relative (standardized over the ranked set); pass reference to standardize against another cohort (see zscore()).

Return type:: ndarray

tcren.cohort.q_f(table, reference=None, sign=1.0, features=('burial', 'n_pep_contacted', 'chain_balance', 'n_hbond'), terms=('F_tcr_pep', 'F_tcr_mhc'), decorrelate=True)[source]#

Pure-tcren combiner z(Q_geom) + sign * z(F) — geometry plus contact energy, no deep learning.

With sign=+1 this is z(Q)+z(F); with sign=-1 it is z(Q)-z(F). On clean (template-covered) poses z(Q)+z(F) beats raw-AF ipTM on both ROC and PR with no DL term (benchmark ledger C42: macro 0.759 ROC / 0.725 PR vs ipTM 0.692 / 0.693). On forced poses the energy inverts, so z(Q)-z(F) is the one that ranks (C27: on the forced GLCTLVAML pose z(Q)-z(F)``=0.71 vs ``z(Q)+z(F)``=0.52). Pick the sign from pose quality — grade it with :func:`strain_z` — or prefer :func:`q_iptm` (``z(ipTM)+z(Q)), the geometry channel that is robust to the inversion without needing the energy at all.

Parameters:

table – the tcren recognize --full --scores table (dict / pandas / polars).
reference – optional cohort to standardize against (see zscore()).
sign – +1 for z(Q)+z(F) (clean poses), -1 for z(Q)-z(F) (forced poses).
features – Q descriptors; defaults to the geometry-only Q_FEATURES_GEOM.
terms – energy terms for F; defaults to F_TERMS.
decorrelate – passed to q_score(). False recovers the legacy equal-weight Q.

Return type:

ndarray

tcren.cohort.q_f_iptm(table, iptm, threshold=0.5, reference=None, features=('burial', 'n_pep_contacted', 'chain_balance', 'n_hbond'), terms=('F_tcr_pep', 'F_tcr_mhc'))[source]#

Pose-adaptive z(Q) + s·z(F) where the F sign s is chosen per structure from ipTM.

Automates the forced-pose inversion: a confident pose (ipTM >= threshold) keeps +z(F) because the contact energy is trustworthy there; a forced pose (ipTM < threshold) flips to -z(F) because the energy inverts on forced poses (benchmark ledger C27/C42). A structure with no ipTM (NaN) keeps +z(F) — nothing marks it as forced. See f_invert_by_iptm() for the boolean flag alone.

ipTM is a pose-confidence proxy, not a calibrated forced-pose detector — grading forced-ness with strain_z() is the principled alternative (C27), and q_iptm() (z(ipTM)+z(Q)) sidesteps the energy entirely. Provided because it is the single-call pose-adaptive combiner.

Parameters:

table – the tcren recognize --full --scores table (dict / pandas / polars).
iptm – per-structure ipTM aligned to table rows.
threshold – ipTM below which a pose is treated as forced and F is inverted (default 0.5).
terms (reference / features /) – as in q_f().

Return type:

ndarray

tcren.cohort.f_invert_by_iptm(iptm, threshold=0.5)[source]#

Boolean per-structure flag: invert F where ipTM < threshold (a forced pose). NaN ipTM is not inverted. This is the F_invert column q_f_iptm() acts on.

Return type:: ndarray

tcren.cohort.phi_bind(table, reference=None)[source]#

Deprecated screening score Phi_bind = Q + 0.5 * [z(-pitch) + z(-F_tcr_mhc)].

Deprecated since version Use: q_score(). Both terms this adds to Q lower ranking accuracy — on TCRvdb macro ROC falls from Q’s 0.795 to 0.653, and z(-pitch) alone is below chance (0.43) (benchmark ledger C19b). The pitch axis also carries AlphaFold-confidence leakage (ledger C19). It is retained only to reproduce older figures; do not use it for new work.

Return type:: ndarray

tcren.cohort.q_coupled(q, energy)[source]#

Parameter-free binder score: interface geometry and coupling-weighted contact energy.

\[S(x) \;=\; \tfrac14\Big[1+\operatorname{erf}\tfrac{z(Q(x))}{\sqrt2}\Big] \Big[1+\operatorname{erf}\tfrac{r\,z(\Delta\Phi(x))}{\sqrt2}\Big], \qquad r=\operatorname{corr}(Q,\Delta\Phi)\]

Each bracket is $2\times$ a Gaussian tail probability, written with erf rather than the normal-CDF symbol so nothing collides with the potential $\Phi$. Three biophysical statements, no free parameter — $z$ is standardization and $r$ is measured, not chosen.

Binding needs both. A complex forms only if there is an interface and the residues in it are favourable. Each factor is the one-class probability that the candidate is native-like on that channel, and the product is the conjunction of two pieces of evidence — the smooth AND, with no threshold and no softness constant.
The energy is admitted in proportion to its coupling. Under joint normality $\mathbb E[z(Q)\mid z(\Delta\Phi)] = r\,z(\Delta\Phi)$, so $r\,z(\Delta\Phi)$ is exactly the part of the energy that is evidence about interface nativeness. Nothing is discarded and nothing is over-trusted.
A forced pose disarms itself. $r<0$ on a fabricated cohort (coupling()), which flips the energy’s sign automatically; $r\approx0$ shrinks the factor to $\Phi(0)=\tfrac12$, a constant, leaving the geometry alone. The failure mode that makes raw $\Phi$ inverting and dangerous (ledger C27) is handled by the same $r$ that measures it.

On TCRvdb this reaches macro ROC 0.799 / PR 0.817 / precision-at-10 %-recall 0.949, ahead of every TCRmodel2 confidence (best 0.795 / 0.800 / 0.916) with no generative term, and it is balanced across the two epitopes rather than trading one for the other.

Parameters:

q – interface-quality scores for the cohort, e.g. q_score() output.
energy – binder-oriented referenced contact energy for the same rows — for receptor ranking use the TCR-referenced $\Delta_{\mathrm{TCR}}\Phi$ (the peptide is fixed there, so the peptide reference carries no signal); for peptide ranking use $\Delta_{\mathrm{pep}}\Phi$ (tcren.reference_delta()).

Returns:

Scores in $(0,1)$; higher is more binder-like. Cohort-relative — rank within the set you scored, do not compare across cohorts.

Return type:

ndarray

tcren.cohort.coupling(q, energy)[source]#

Interface–energy coupling $r=\mathrm{corr}(Q,\,\Delta\Phi)$ over a cohort — the label-free forced-pose diagnostic.

In a genuine complex the two channels are physically tied: a larger, better-packed interface holds more contacts, so favourable contact energy and good interface geometry rise together and $r>0$. A structure generator that manufactures a pose optimises contacts without the interface, breaking the tie — the two channels decouple or run opposite, and $r<0$.

So the sign and size of $r$ say how far the energy of this cohort can be trusted, using no labels and no reference set. It is the weight q_coupled() admits the energy with.

Parameters:

q – interface-quality scores (e.g. q_score() output) for the cohort.
energy – binder-oriented referenced contact energy for the same rows.

Returns:

Pearson $r$ over the rows where both are finite; 0.0 if fewer than three remain (an uninformative cohort contributes no energy evidence rather than a spurious weight).

Return type:

float

tcren.cohort.strain_z(table, reference=None)[source]#

Crystal-calibrated interface strain; higher = more forced. The recommended forced-pose score.

Directional mean-z of STRAIN_TERMS with fixed physical signs. Pass the crystal cohort as reference to reproduce the provenance gradient (crystal +0.02 < generated-real +0.40 < generated-decoy +0.81); without it the score is only relative within the input set.

Prefer this over tcren.recognition.forced_pose_score() (p_forced): it is unfitted — no logistic, no coefficients, just signed standardization — so it carries no training set and is fully reproducible, whereas FORCED_POSE_MODEL’s coefficients were frozen from a training set that no longer exists (benchmark ledger C23). It also grades forced-ness continuously, which is what pairs with q_score() to catch the forced poses where the contact energy inverts (ledger C27).

Return type:: ndarray

tcren.cohort.native_reference()[source]#

The interface-geometry descriptors over the 374 Native2026 crystal complexes, shipped so a single user structure (or any small cohort) can be standardized against the natural interface manifold instead of against itself — the deployment path for generic input:

from tcren import cohort
q = cohort.q_score(user_table, reference=cohort.native_reference(),
                   features=cohort.Q_FEATURES_GEOM)

Use Q_FEATURES_GEOM (the four geometry terms) for one structure: the fifth term pp_combo is a within-cohort z-contrast and is undefined for a single row. Returns a dict of column arrays (burial, n_pep_contacted, chain_balance, n_hbond, F_cdr12, F_cdr3a) usable as the reference argument. Provenance: tcren recognize --full over the Native2026 set.

Return type:: dict

tcren.cohort.Q_FEATURES = ('burial', 'n_pep_contacted', 'chain_balance', 'n_hbond', 'pp_combo')#: The five interface-quality descriptors, equal-weighted in q_score(). Each is oriented positive-is-better as given. pp_combo is the CDR1/2-vs-CDR3alpha TCRen contrast — the one energy term robust to the forced-pose inversion (benchmark ledger C27), since it is a contrast rather than an absolute contact energy. Per-term macro AUROC on TCRvdb: burial 0.73, n_hbond 0.69, pp_combo 0.66, n_pep_contacted 0.62, chain_balance 0.61; the terms are near-independent (mean absolute Spearman 0.20).

tcren.cohort.Q_FEATURES_CORE = ('burial', 'chain_balance', 'n_hbond', 'pp_combo')#

it is the weakest term and removing it raises macro AUROC 0.795 -> 0.801 on TCRvdb (benchmark ledger, energy memo). Pass features=Q_FEATURES_CORE to q_score() for the simpler, marginally better score.

Type:: The four load-bearing descriptors. n_pep_contacted is dropped

tcren.cohort.Q_FEATURES_GEOM = ('burial', 'n_pep_contacted', 'chain_balance', 'n_hbond')#: The four geometry-only descriptors — Q_FEATURES without the pp_combo energy contrast. This is Q_geom, the AF-orthogonal channel that survives the forced-pose regime where the contact energy inverts (benchmark ledger C27/C42): z(ipTM) + z(q_score(..., features=Q_FEATURES_GEOM)) beats raw-AF ipTM on well-modelled (“template-covered”) epitopes on both ROC and PR, while the energy term is used only conditioned on pose quality. Pass to q_score() / q_iptm().

tcren.cohort.F_TERMS = ('F_tcr_pep', 'F_tcr_mhc')#: The TCRen contact-energy terms summed into the binder-oriented f_score(). F_tcr_pep is the TCR:peptide TCRen energy, F_tcr_mhc the TCR:MHC energy; both are emitted by tcren recognize. They are raw energies (lower = tighter), so f_score() negates the sum to make higher = more binder-like. This term is pose-conditional — it reads real binding chemistry on well-modelled (crystal-templated) poses and inverts on forced ones (benchmark ledger C27/C42): on the forced GLCTLVAML TCRvdb pose -F_tcr_pep ranks binders at AUROC 0.36 (backwards), on the clean YLQPRTFLL pose at 0.59. Use it only conditioned on pose quality — gate with strain_z(), or read z(Q)-z(F) on forced poses and z(Q)+z(F) on clean ones.

tcren.cohort.STRAIN_TERMS = (('cdr3b_topep', 1.0), ('cdr3b_reach', 1.0), ('extent_per_ct', 1.0), ('chain_balance', -1.0))#: Crystal-calibrated interface-strain terms with their physical signs. A forced pose reaches further from the peptide with a thinner, less balanced interface.

Potentials#

tcren.potential.model module#

Pairwise residue-level statistical potentials.

A Potential is a long-form table of pairwise amino-acid energies keyed on (residue.aa.from, residue.aa.to). The “from” side is conventionally the TCR residue and the “to” side the antigen (peptide) residue, matching the orientation of the legacy R pipeline. Potentials can be loaded from the two CSV layouts shipped with the project (wide and long) and exported to a dense matrix for fast scoring.

tcren.potential.model.AA20: tuple[str, ...] = ('L', 'F', 'I', 'M', 'V', 'W', 'Y', 'C', 'H', 'A', 'G', 'P', 'T', 'S', 'Q', 'N', 'D', 'E', 'R', 'K')#: 20 standard amino acids (one-letter), TCRen ordering used in the paper.

tcren.potential.model.AA21: tuple[str, ...] = ('A', 'I', 'L', 'V', 'R', 'H', 'K', 'C', 'M', 'S', 'T', 'D', 'E', 'N', 'Q', 'G', 'P', 'Y', 'F', 'W', '-')#

21 amino acids plus the gap symbol.

Type:: Alphabet of the alignment-matrix variant

class tcren.potential.model.Potential(name, matrix, alphabet)[source]#

Bases: object

A pairwise amino-acid potential in long form.

Variables:

name (str) – Identifier of the potential (e.g. "TCRen", "MJ", "Keskin").
matrix (polars.dataframe.frame.DataFrame) – Long-form table with columns residue.aa.from, residue.aa.to, value.
alphabet (tuple[str, ...]) – Amino-acid symbols present on each axis.

Parameters:

name (str)
matrix (DataFrame)
alphabet (tuple[str, ...])

name: str#

matrix: DataFrame#

alphabet: tuple[str, ...]#

value(aa_from, aa_to)[source]#

Return the energy for an ordered residue pair.

Parameters:

aa_from (str) – One-letter code of the “from” (TCR) residue.
aa_to (str) – One-letter code of the “to” (antigen) residue.

Returns:

The pairwise energy.

Raises:

KeyError – If the pair is absent from the potential.

Return type:

float

as_matrix()[source]#

Return a dense (n, n) matrix and an amino-acid → index map.

Rows are indexed by residue.aa.from, columns by residue.aa.to. Missing pairs are filled with nan. The dense form is cached (the table is immutable), so repeated scoring/energy calls over one potential rebuild it once. Callers treat the returned array as read-only.

Return type:: tuple[ndarray, dict[str, int]]

to_csv(path)[source]#

Write the potential to a long-form CSV (from, to, value).

Parameters:: path (str | Path)
Return type:: None

classmethod from_csv(path, name=None, value_col=None)[source]#

Load a potential from a CSV, auto-detecting wide vs long layout.

Two layouts are supported:

wide — residue.aa.from, residue.aa.to, <name> (e.g. TCRen_potential.csv with a TCRen value column).
long — residue.aa.from, residue.aa.to, potential, value (e.g. MJ_Keskin_potentials.csv); load a single named potential from it.

Parameters:

path (str | Path) – Path to the CSV file.
name (str | None) – Which potential to select (long layout) or the name to assign (wide layout). Defaults to the value-column name (wide) and is required when a long file holds more than one potential.
value_col (str | None) – Override the value column name for the wide layout.

Returns:

The loaded Potential.

Return type:

Potential

tcren.potential.model.tcren()[source]#

Load the bundled classic TCRen potential (cached; treat as read-only).

Return type:: Potential

tcren.potential.model.mj()[source]#

Load the bundled Miyazawa–Jernigan potential (cached; treat as read-only).

Return type:: Potential

tcren.potential.model.keskin()[source]#

Load the bundled Keskin contact potential (cached; treat as read-only).

Return type:: Potential

tcren.potential.derive module#

Derivation of the TCRen statistical potential from observed contact maps.

This is a direct port of the R derivations in code_paper/2_TCRen_derivation.Rmd (variant="classic") and tcren_am/tcren_am.Rmd (variant="am"). The classic variant reproduces TCRen_potential.csv; the alignment-matrix variant reproduces tcren_am/tcren.txt.

tcren.potential.derive.symmetrize_counts(counts)[source]#

Fold a directed aa-pair count table onto its transpose: N + Nᵀ.

TCRen counts are directed — from is a TCR residue and to a peptide residue — so N[a,b] and N[b,a] are different observations and the derived matrix is asymmetric. Adding the transpose treats each contact as an unordered pair, which is the convention Miyazawa–Jernigan uses. Diagonal cells double, as they must: a C–C contact is one unordered pair observed from both sides.

Symmetrising here — on the raw counts, before the log-odds — is not the same as averaging the finished potential. The marginals (total.from / total.to) are recomputed from the folded counts, so the expected term of the log-odds changes too; averaging the energies afterwards leaves the asymmetric background in place. On the Native2026 derivation set the two disagree by 0.29 on average (max 0.82), so the distinction is not cosmetic.

Cysteine. The classic directed derivation drops from == "C" because free Cys is essentially absent from CDR loops — on Native2026 only 4 of 8062 contacts (0.05 %) have a TCR-side Cys, against 32 (0.40 %) on the peptide side. Folding grafts those peptide-side observations onto the Cys row instead of discarding the column, so the symmetric matrix keeps a full 20×20 alphabet at no cost: the row that would have been dropped for having no data inherits the data the other axis did have.

Parameters:: counts (DataFrame) – Long table with residue.aa.from, residue.aa.to and count.
Returns:: The folded table, with one row per unordered pair-cell (still stored in both orientations, so it is a full symmetric matrix).
Return type:: DataFrame

Example

>>> import polars as pl
>>> c = pl.DataFrame({"residue.aa.from": ["A"], "residue.aa.to": ["W"], "count": [3.0]})
>>> out = symmetrize_counts(c)
>>> sorted((r["residue.aa.from"], r["residue.aa.to"], r["count"]) for r in out.iter_rows(named=True))
[('A', 'W', 3.0), ('W', 'A', 3.0)]

tcren.potential.derive.derive_tcren(contacts, include=None, exclude=None, pseudocount=1, variant='classic', beta=44.0, drop_cys=None, weights=None, symmetric=False)[source]#

Derive a TCRen potential from a table of residue contacts.

Parameters:

contacts (DataFrame) – Long table of TCR↔peptide contacts with at least residue.aa.from, residue.aa.to and (for filtering) pdb.id.
include (list[str] | None) – If given, keep only contacts whose pdb.id is in this list.
exclude (list[str] | None) – If given, drop contacts whose pdb.id is in this list.
pseudocount (int) – Added to every amino-acid pair count (default 1).
variant (str) – "classic" (natural-log log-odds over 20 aa, Cys dropped from the “from” axis) or "am" (log2/beta over 21 symbols including a gap, Cys retained).
beta (float) – Temperature divisor used by the "am" variant.
drop_cys (bool | None) – Override the per-variant default for dropping from == "C" rows. Forced to False when symmetric is set (dropping one axis would un-symmetrise the result).
weights (dict[str, float] | None) – Optional per-structure weights {pdb.id: weight}. When given, each structure’s contributions to the aa-pair counts are multiplied by its weight (rows whose pdb.id is absent from the map default to weight 1.0); this down-weights redundancy while keeping all data (see tcren.potential.redundancy.cluster_weights()). None (default) is unweighted and byte-identical to the legacy derivation.
symmetric (bool) – Fold the raw counts onto their transpose (symmetrize_counts()) before the log-odds, yielding a symmetric value[a,b] == value[b,a] potential over an unordered amino-acid pair — the same convention as the bundled Miyazawa–Jernigan matrix, and therefore directly comparable to it. Default False keeps the directed TCR→peptide potential, which is the shipped TCRen_potential.csv.

Returns:

The derived Potential. For "am" the long matrix additionally carries a count column.

Return type:

Potential

tcren.potential.derive.derive_tcren_loo(contacts, pdb_ids, **kwargs)[source]#

Leave-one-out TCRen: derive once per structure, excluding it each time.

Parameters:

contacts (DataFrame) – Contact table (see derive_tcren()).
pdb_ids (list[str]) – Structures to leave out one at a time (also the inclusion set).
**kwargs – Forwarded to derive_tcren().

Returns:

Long table residue.aa.from, residue.aa.to, TCRen.LOO, pdb.id stacking the per-structure potentials.

Return type:

DataFrame

tcren.potential.redundancy module#

Non-redundancy clustering for TCRen derivation inputs.

The TCRen potential is derived from a non-redundant set of αβ TCR–pMHC structures: near-duplicate complexes (same/similar CDR3α + CDR3β + peptide) are collapsed to a single representative so the contact statistics are not dominated by repeated entries.

This is a verbatim lift of the nonredundant/alphabeta helpers from notebooks/natcompsci2022/01_nonred_and_derivation.ipynb: a complete-linkage hierarchical clustering on the summed Damerau–Levenshtein distance of the CDR3α+CDR3β+peptide strings, cut at distance t (default 6.0), keeping one representative per cluster.

tcren.potential.redundancy.nonredundant_ids(markup, t=6.0, fields=('cdr3a', 'cdr3b', 'peptide'), linkage_method='complete')[source]#

Non-redundant representative pdb.id for each cluster of similar complexes.

Clusters structures by the summed Damerau–Levenshtein distance over fields (CDR3α + CDR3β + peptide) using hierarchical clustering, then returns one representative pdb.id per cluster (the first in row order).

Parameters:

markup (DataFrame) – Per-structure table with a pdb.id column and the fields columns.
t (float | None) – Distance cutoff for fcluster (criterion="distance"). None turns redundancy filtering off and returns every pdb.id unchanged.
fields (Sequence[str]) – Sequence columns whose per-pair distances are summed (default cdr3a, cdr3b, peptide).
linkage_method (str) – Linkage method for scipy.cluster.hierarchy.linkage (default "complete").

Returns:

Representative pdb.id values, one per cluster, in first-seen order.

Return type:

list[str]

tcren.potential.redundancy.cluster_weights(markup, t=6.0, fields=('cdr3a', 'cdr3b', 'peptide'), linkage_method='complete')[source]#

Inverse-cluster-size (Henikoff-style) weight for each pdb.id.

Clusters structures exactly as nonredundant_ids() (same distance, same linkage), then assigns every structure the weight 1 / cluster_size. A unique structure gets weight 1.0; each member of a size-k redundancy cluster gets 1/k, so the cluster contributes a total weight of 1 to the derivation. Feed the result to tcren.potential.derive.derive_tcren()’s weights argument to down-weight redundancy while keeping every structure’s data.

Parameters:

markup (DataFrame) – Per-structure table with a pdb.id column and the fields columns.
t (float) – Distance cutoff for fcluster (criterion="distance").
fields (Sequence[str]) – Sequence columns whose per-pair distances are summed.
linkage_method (str) – Linkage method for scipy.cluster.hierarchy.linkage.

Returns:

Mapping {pdb.id: 1 / cluster_size}.

Return type:

dict[str, float]

tcren.potential.redundancy.alphabeta_ids(contacts)[source]#

pdb.id of complexes whose receptor contacts are exclusively TRA/TRB (αβ TCRs).

Parameters:: contacts (DataFrame) – Contact table with pdb.id and chain.type.from columns.
Returns:: The pdb.id values whose set of chain.type.from is a subset of {"TRA", "TRB"} (i.e. no γδ chains).
Return type:: list[str]

Scoring#

tcren.scoring module#

Candidate-peptide scoring by amino-acid substitution.

Ports the second half of run_TCRen.R: for each candidate peptide, substitute its amino acids at the contacted peptide positions of a structure’s contact map and sum the pairwise potential over all contacts. Lower scores indicate more favourable interactions.

tcren.scoring.score_peptides(contact_map, candidates, potential, interface='tcr_peptide', require_same_length=True, substituted_side=None, tcr_regions='all', contact_weight='residue')[source]#

Score candidate peptides against a structure’s contact map.

Parameters:

contact_map (ContactMap) – The structure’s contact map.
candidates (Iterable[str]) – Candidate peptide sequences (one-letter).
potential (Potential) – Pairwise potential to score with.
interface (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – Which interface to score over (default "tcr_peptide").
require_same_length (bool) – Only score candidates whose length matches the structure’s peptide length (mirrors the legacy length join). Ignored when the contact map has no recorded peptide length.
substituted_side (str | None) – "to" or "from" — which contact side the candidate is threaded onto. Defaults to the peptide side of interface.
tcr_regions (str) – which TCR regions to keep on the TCR side ("all" default = no filter = legacy behaviour; "cdr" or "cdr+fr" to restrict).
contact_weight (str) – "residue" (default, legacy) gives every contacting residue pair unit weight; "atomic" weights each residue pair by its n_atom_contacts heavy-atom-pair count, so the energy tracks the LJ+Coulomb atom-pair sum more closely. "atomic" requires the contact map to have been built with count_atoms=True.

Returns:

Columns complex.id, peptide, potential, score sorted by complex.id then ascending score.

Return type:

DataFrame

tcren.scoring.score_structures(contact_maps, candidates, potential, **kwargs)[source]#

Score candidates against several structures and stack the results.

Parameters:

contact_maps (Iterable[ContactMap])
candidates (Iterable[str])
potential (Potential)

Return type:

DataFrame

class tcren.scoring.RecognitionMatrix(positions, aa, energy, side, interface)[source]#

Bases: object

Per-position × amino-acid substitution-energy landscape (see recognition_matrix()).

energy[i, a] is the summed pairwise potential over position positions[i]’s contacts when amino acid aa[a] sits there and the other side is held fixed. Lower = more favourable, so a per-position preference is -energy (higher = preferred). Positions with no contact are omitted. Entries are NaN for amino acids the potential leaves undefined (e.g. cysteine pairs in TCRen), exactly as score_peptides() drops those contacts — so reduce columns with np.nan* ops.

Parameters:

positions (list)
aa (tuple)
energy (object)
side (str)
interface (str)

positions: list#: one (chain_type, region, pos, native_aa) tuple per row of energy

aa: tuple#: the 20 amino-acid column order

energy: object#: (n_positions, 20) float ndarray of substitution energies

side: str#: which side was scanned ("from" = TCR, "to" = peptide)

interface: str#

tcren.scoring.recognition_matrix(contact_map, potential, *, interface='tcr_peptide', side=None, tcr_regions='all')[source]#

The per-position × 20-AA substitution-energy matrix for one interface side.

For interface="tcr_peptide", side="from" scans the TCR/CDR3 (the motif-matrix analog) and side="to" scans the peptide (the CPL-matrix analog). Each entry is the F energy summed over that position’s contacts with the given amino acid substituted in, the opposite side fixed — the same virtual-substitution path as score_peptides(), resolved per position rather than summed over the whole sequence.

Parameters:

contact_map (ContactMap) – the structure’s contact map.
potential (Potential) – pairwise potential (TCRen for TCR:peptide).
interface (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – which interface to score over.
side (str | None) – "from" or "to"; defaults to the non-peptide side for tcr_peptide (i.e. the TCR), and to the peptide side for the presentation interfaces.
tcr_regions (str) – TCR-region filter ("all"/"cdr"/"cdr+fr") — use "cdr" to restrict a TCR-side scan to the CDRs.

Returns:

A RecognitionMatrix. Rows are the contacted positions in (chain, region, pos) order; columns are RecognitionMatrix.aa.

Return type:

RecognitionMatrix

tcren.scoring_rank module#

Percentile rank of a native peptide’s energy against a random pMHC background.

For a given structure, scores the native peptide together with a background of random (or epitope-sampled) peptides of the same length and reports where the native score falls in that distribution. Lower TCRen energy = better binder, so a small rank_pct means the native peptide scores at least as well as only a small fraction of the background.

tcren.scoring_rank.background_peptides(length, n=1000, seed=0, source=None, anchors=None, pool=None)[source]#

Build a background set of n peptides of the given length.

Parameters:

length (int) – Required peptide length.
n (int) – Number of background peptides to return.
seed (int) – Random seed for reproducibility.
source (str | None) – Optional FASTA/text file of epitopes. When given, peptides of the requested length are sampled from it (with replacement); each sampled sequence is also position-permuted so the background is a shuffled-epitope distribution rather than a verbatim copy. When None, peptides are drawn uniformly at random over the 20 amino acids.
anchors (dict[int, str] | None) – 1-based position -> fixed residue, applied after sampling. Pass the cognate peptide’s anchor residues (typically {2: p[1], length: p[-1]}) to hold MHC presentation constant so the background varies only where the TCR reads. Without this, a “random peptide” background is dominated by peptides that would not be presented at all, which inflates the apparent discrimination. Required to reproduce the published cognate-vs-decoy benchmark.
pool (list[str] | None) – Explicit candidate peptides to draw from (with replacement) instead of generating them – e.g. the native epitopes of the other structures in a cohort, which is the decoy set used by the crystal peptide-swap benchmark.

Returns:

A list of n upper-case peptide strings, each of length length.

Return type:

list[str]

tcren.scoring_rank.percentile_rank(contact_map, peptide, potential, *, interface='tcr_peptide', n_background=1000, seed=0, tcr_regions='all', background=None, contact_weight='residue')[source]#

Percentile rank of a peptide’s energy against a random pMHC background.

Scores peptide together with a background set (supplied or generated) and returns the fraction of background peptides whose score is <= the native score. Because lower TCRen energy means a better binder, a smaller rank_pct indicates the native peptide is among the strongest binders.

Parameters:

contact_map (ContactMap) – The structure’s contact map.
peptide (str) – The native peptide to rank.
potential (Potential) – Pairwise potential to score with.
interface (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – Which interface to score over (default "tcr_peptide").
n_background (int) – Size of the generated background (ignored if background is given).
seed (int) – Random seed for background generation.
tcr_regions (str) – Which TCR regions to keep on the TCR side ("all" default, "cdr" or "cdr+fr"); passed through to score_peptides.
background (Iterable[str] | None) – Explicit background peptides; when None a uniform-random background of length len(peptide) is generated.
contact_weight (str) – "residue" (default) or "atomic"; passed through to score_peptides ("atomic" needs an n_atom_contacts contact map).

Returns:

Mapping with keys peptide, score (native energy), rank_pct (fraction of background with score <= native), and n_background.

Return type:

dict

tcren.ddg module#

Fast ΔΔG of peptide point mutations (virtual-matrix path).

Implements the paper’s fast ΔΔG: no atoms move and no re-docking is performed. A mutation’s effect is read straight off the substitution potential by re-scoring the mutant sequence on the same contact map. ddg = E(native) - E(mutant), and lower energy is a more favourable interface throughout tcren, so a positive value flags a stabilising mutation: the mutant scores below the native and binds better. A negative value is the destabilising one.

tcren.ddg.ddg(contact_map, native, mutant, potential, *, interface='tcr_peptide', tcr_regions='all', contact_weight='residue')[source]#

ΔΔG of a peptide mutation as E(native) - E(mutant).

Parameters:

contact_map (ContactMap) – The structure’s contact map.
native (str) – Native peptide sequence.
mutant (str) – Mutant peptide sequence (same length as native).
potential (Potential) – Pairwise potential to score with.
interface (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – Which interface to score over (default "tcr_peptide").
tcr_regions (str) – Which TCR regions to keep on the TCR side (passed through to score_peptides).
contact_weight (str) – "residue" (default) or "atomic"; passed through to score_peptides.

Returns:

E(native) - E(mutant); positive means the mutant has the LOWER energy, i.e. the mutation is stabilising. Negative means destabilising. Always 0.0 for interfaces that do not contain the peptide (e.g. "tcr_mhc"), since a peptide mutation cannot affect them.

Return type:

float

tcren.ddg.alanine_scan(contact_map, native, potential, *, interface='tcr_peptide', tcr_regions='all', contact_weight='residue')[source]#

Alanine scan of the native peptide.

Mutates each position of native to alanine in turn and reports the ΔΔG of that single substitution. One row per peptide position.

Parameters:

contact_map (ContactMap) – The structure’s contact map.
native (str) – Native peptide sequence.
potential (Potential) – Pairwise potential to score with.
interface (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – Which interface to score over (default "tcr_peptide").
tcr_regions (str) – Which TCR regions to keep on the TCR side.
contact_weight (str) – "residue" (default) or "atomic"; passed through to score_peptides.

Returns:

Columns pos (0-based), wt_aa (native residue at that position) and ddG (E(native) - E(Ala@pos)). Positions without TCR contacts yield ddG == 0.0. For interfaces that do not contain the peptide (e.g. "tcr_mhc") every position yields ddG == 0.0.

Return type:

DataFrame

tcren.ddg.neoantigen_ddg(contact_map, native, mutants, potential, **kw)[source]#

ΔΔG of candidate neoantigen mutants relative to a native peptide.

Parameters:

contact_map (ContactMap) – The structure’s contact map.
native (str) – Native peptide sequence.
mutants (Iterable[str]) – Candidate mutant peptides (each the same length as native).
potential (Potential) – Pairwise potential to score with.
**kw – Forwarded to ddg() (interface, tcr_regions).

Returns:

Columns native, mutant and ddG (E(native) - E(mutant); positive means the mutant has the lower energy, i.e. the substitution is stabilising), one row per mutant – so ranking candidates by descending ddG puts the best first.

Return type:

DataFrame

tcren.ddg.reference_delta(contact_map, peptide, potential, *, interface='tcr_peptide', reference_aa='A', tcr_regions='all', contact_weight='residue')[source]#

Poly-alanine reference difference ΔΦ = Φ(peptide) − Φ(reference) on this contact map.

ΔΦ is the full-peptide alanine-scan difference — the sum of the per-position native→Ala ΔΔGs of alanine_scan(). It subtracts the interface’s identity-independent baseline Φ(reference), i.e. what the pose geometry scores when every peptide residue is reference_aa, leaving the sequence-specific part.

On a fixed contact map this is Φ(peptide) minus a constant, so it does not change the ranking of candidates threaded onto one structure. It differs from raw Φ only across candidates that each have their own structure (e.g. AlphaFold peptide-swap models), where it normalises out the per-pose interface geometry. That normalisation rescues forced / wrong-register poses whose geometry corrupts the raw contact energy (the CPL ila1 case: TCR-ranking ROC 0.35 → 0.83), at a small cost on clones where the generated geometry is itself informative — so it is a scoring mode for generated poses, not a default. It is not an affinity ΔΔG: a dimensionless contact-preference difference, not a free energy (see tcren.refine.register for the geometry defect it corrects). Empirically both raw Φ and ΔΦ are within-receptor ranking scores, not binding constants — on the ATLAS SPR set they correlate with ΔG/Kd/koff/kon only at ρ ≤ 0.3 in magnitude (off-rate comes from tcren.mechanics).

Parameters:

contact_map (ContactMap) – The candidate’s own contact map.
peptide (str) – The candidate peptide sequence.
potential (Potential) – Pairwise potential to score with.
interface (Literal['tcr_peptide', 'tcr_mhc', 'peptide_mhc']) – Which interface to score over (default "tcr_peptide").
reference_aa (str) – The amino acid the reference peptide is made of (default alanine).
tcr_regions (str) – Which TCR regions to keep on the TCR side.
contact_weight (str) – "residue" (default) or "atomic".

Returns:

ΔΦ = Φ(peptide) − Φ(reference); more negative = the sequence adds more favourable contacts than the reference baseline. 0.0 for interfaces without the peptide (e.g. "tcr_mhc").

Return type:

float

tcren.cpl module#

Predict a combinatorial-peptide-library response matrix from one template TCR:pMHC structure.

A positional-scanning combinatorial peptide library (CPL) measures a T-cell clone’s peptide preference one position at a time. For a peptide of length L, each of the L x 20 sublibraries fixes position i to amino acid a and leaves every other position an equimolar 1/20 mixture, so the measured cell is an ensemble mean,

R[i, a] = E[ response | x_i = a ] .

This module predicts that matrix from a single deposited or modelled complex. Position i has a fixed set of contact partners in the template, so threading each of the twenty residues through the same contact map and re-reading the potential costs one batched call per interface. Nothing is re-docked, nothing moves, and nothing is fitted to any assay.

WHAT A CELL IS SCORED WITH. The assay reads activation, which needs the peptide presented and the receptor engaged: a substitution that abolishes MHC binding abolishes the response whatever the receptor would have done. Every cell therefore carries the sum of both peptide-bearing interfaces,

Phi = Phi(TCR:peptide) + Phi(peptide:MHC),

TCRen over the first and Miyazawa–Jernigan over the second by default. The two channels are statistically uncorrelated over these cells, so they add rather than duplicate. A position the receptor never touches is an anchor: its TCR term is identically zero, so scoring by the sum degrades gracefully to presentation alone, and the older “TCRen at receptor-facing positions, MJ at anchors” partition is this rule’s special case rather than a separate mode.

TWO REFERENCE STATES, BOTH USEFUL. A raw Phi carries a large per-position offset that says only how many contacts the position makes, so a cell is meaningful only relative to the other residues that could sit there. Two references are offered and both are exact differences of whole-peptide energies, so the “rest of the peptide” cancels:

"wild_type": Phi(x_{i->wt}) - Phi(x_{i->a}) – a mutation scan off the residue the template carries. This is the neoantigen / epitope-design question: is this substitution better than what is there? It is tcren.ddg.ddg() resolved per position.
"equimolar" (default): mean_b Phi(x_{i->b}) - Phi(x_{i->a}) – referenced to the 1/20 mixture, which is the assay’s own null: the CPL background at position i is exactly that mixture. Use this to compare against measured CPL cells. It is also the only one of the two under which the template’s own residue is an ordinary measurement rather than a forced zero.

The two differ by a per-position constant, and that constant is not noise: it is how far the template’s residue sits above its column’s mean. Referencing to the wild type folds that between-position quantity into every cell of the column, which is the right thing for a mutation scan and the wrong thing for a comparison against an assay whose background is the mixture.

SIGN. Lower energy is a better binder throughout tcren, and both references are written as reference - candidate, so positive means favourable: a positive wild_type value says the substitution improves on the template residue, and a positive equimolar value says the residue is better than the average residue at that position.

Example

>>> from tcren import ContactMap, parse_structure, response_matrix, position_scan
>>> from tcren.annotation import classify_chains
>>> from tcren.mhc import annotate_mhc
>>> s = parse_structure("3HG1.pdb", pdb_id="3HG1")
>>> classify_chains(s, organism="human"); annotate_mhc(s)
>>> rm = response_matrix(ContactMap.from_structure(s, cutoff=5.0))
>>> rm.peptide
'ELAGIGILTV'
>>> position_scan(rm, 5).head(3)          # every residue at position 5

tcren.cpl.AA20: tuple[str, ...] = ('A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y')#: Column order of every matrix this module returns.

tcren.cpl.INTERFACES: tuple[str, ...] = ('tcr_peptide', 'peptide_mhc')#: The two peptide-bearing interfaces, and the potential family each is scored with by default. tcr_mhc is absent on purpose: a peptide substitution cannot change it, so its contribution to every cell of the matrix is identically zero.

tcren.cpl.REFERENCES: tuple[str, ...] = ('equimolar', 'wild_type')#: Accepted values of the reference argument.

class tcren.cpl.ResponseMatrix(peptide, positions, interface_class, aa, phi, phi_tcr, phi_mhc, structure_id)[source]#

Bases: object

A predicted CPL response matrix: peptide positions x twenty amino acids.

phi[i, a] is the total interface energy of the template with amino acid aa[a] threaded at positions[i], everything else held at the template sequence and geometry. It is a whole-complex energy, not a per-position contribution, so differences within a row are exact and differences between rows are meaningless – which is why every accessor references a row against something in that same row.

Use referenced() (or mutation_effect() / position_scan() / equimolar_effect()) rather than reading phi directly.

Parameters:

peptide (str)
positions (tuple[int, ...])
interface_class (tuple[str, ...])
aa (tuple[str, ...])
phi (object)
phi_tcr (object)
phi_mhc (object)
structure_id (str)

peptide: str#: the template’s own peptide sequence

positions: tuple[int, ...]#: 1-based peptide positions that contact either interface

interface_class: tuple[str, ...]#: "receptor" or "anchor", one per row

aa: tuple[str, ...]#: the twenty column labels, AA20

phi: object#

MHC

Type:: (n_positions, 20) total energy, TCR
Type:: peptide + peptide

phi_tcr: object#

peptide component alone

Type:: the TCR

phi_mhc: object#

MHC component alone

Type:: the peptide

structure_id: str#: the template’s pdb_id

row_of(position)[source]#

Row index of a 1-based peptide position, or raise if it contacts nothing.

Raising is deliberate. A position that touches neither interface has an entirely flat row, and returning that silently would report “this position tolerates everything” for what is really “this template says nothing about this position”.

Parameters:: position (int)
Return type:: int

column_of(aa)[source]#

Column index of a one-letter amino acid code.

Parameters:: aa (str)
Return type:: int

wild_type_at(position)[source]#

The residue the template carries at a 1-based peptide position.

Parameters:: position (int)
Return type:: str

referenced(reference='equimolar')[source]#

The matrix as effects, (n_positions, 20), positive = favourable.

Parameters:: reference (str) – "equimolar" references each row to the mean over its twenty residues – the 1/20 mixture the assay holds the other positions at. "wild_type" references it to the residue the template carries, making the matrix a mutation scan whose wild-type column is identically zero.
Returns:: reference_energy - phi, so a positive entry is a favourable residue.

to_frame(reference=None)[source]#

Long form: one row per (position, amino acid) cell.

Parameters:: reference (str | None) – emit only this reference’s effect column. None (default) emits both as effect_equimolar and effect_wild_type, which is what a side-by-side comparison against a measured matrix wants.
Returns:: Columns structure_id, pos, wt_aa, aa, is_wt, interface_class, phi, phi_tcr, phi_mhc, and the requested effect column(s).
Return type:: DataFrame

tcren.cpl.response_matrix(contact_map, peptide=None, *, tcr_potential=None, mhc_potential=None, tcr_regions='all', contact_weight='residue')[source]#

Predict the CPL response matrix of a template TCR:pMHC complex.

Every peptide position that contacts either peptide-bearing interface gets a row; each of the twenty residues is threaded there in turn on the template’s own contact map and scored with both potentials, and the two interface energies are summed. See the module docstring for why both interfaces enter and what the two reference states mean.

Parameters:

contact_map (ContactMap) – the template’s contact map, e.g. ContactMap.from_structure(s, cutoff=5.0). The structure must have been chain-typed (tcren.annotation.classify_chains) and MHC-annotated (tcren.mhc.annotate_mhc) first, or the peptide:MHC interface comes out empty.
peptide (str | None) – sequence to thread substitutions off. Defaults to the template’s own peptide, taken from the contact map. Must match the template’s peptide length.
tcr_potential (Potential | None) – potential for TCR:peptide. Default: the bundled TCRen.
mhc_potential (Potential | None) – potential for peptide:MHC. Default: Miyazawa–Jernigan, which is what the presentation interface is scored with throughout tcren – TCRen is a directed TCR-to-peptide recognition potential and does not describe the groove.
tcr_regions (str) – which TCR regions contribute on the TCR side ("all"/"cdr"/"cdr+fr").
contact_weight (str) – "residue" (default) or "atomic"; passed through to tcren.scoring.score_peptides().

Returns:

A ResponseMatrix.

Raises:

ValueError – if the complex has no peptide contacts at all, or if peptide has the wrong length. Both are silent-empty failures otherwise: an empty matrix is indistinguishable from a genuinely flat one.

Return type:

ResponseMatrix

tcren.cpl.mutation_effect(rm, position, aa, *, reference='equimolar')[source]#

Effect of ONE substitution at ONE position, positive = favourable.

Parameters:

rm (ResponseMatrix) – a predicted ResponseMatrix.
position (int) – 1-based peptide position.
aa (str) – the one-letter residue to put there.
reference (str) – "equimolar" (default) scores it against the 1/20 mixture at that position; "wild_type" scores it against the residue the template carries, i.e. the classical mutation-scan ddG.

Returns:

A single number. Under "wild_type" the template’s own residue returns exactly 0.0.

Return type:

float

tcren.cpl.position_scan(rm, position, *, reference='equimolar')[source]#

Effect of EVERY substitution at one position – one column of the response matrix.

Parameters:

rm (ResponseMatrix) – a predicted ResponseMatrix.
position (int) – 1-based peptide position.
reference (str) – as mutation_effect().

Returns:

Twenty rows – pos, aa, wt_aa, is_wt, interface_class, phi, effect – sorted best residue first, so the head of the frame is what the receptor prefers at that position.

Return type:

DataFrame

tcren.cpl.equimolar_effect(rm, position, aa=None)[source]#

Effect of replacing a residue with a random 1/20 mixture at that position.

This is the sublibrary construction read backwards. The CPL background at position i is the equimolar mixture, so the cost of giving up a defined residue for that mixture is

mean_b Phi(x_{i->b}) - Phi(x_{i->a})

– exactly the "equimolar" reference. A positive result means aa is better than the average residue at that position, so scrambling it to the mixture loses that much; a negative result means the mixture is an improvement, i.e. the template residue is disfavoured there.

Parameters:

rm (ResponseMatrix) – a predicted ResponseMatrix.
position (int) – 1-based peptide position.
aa (str | None) – the residue being given up. Defaults to the one the template carries, which is the usual question: how much does this position’s identity matter?

Returns:

A single number in the potential’s units.

Return type:

float

tcren.binder module#

TCR binder/non-binder classification from AF-orthogonal interface geometry.

The shipped model (binder_score(), BINDER_MODEL) scores a TCR-pMHC complex from native interface descriptors (tcren._geom) plus the CDR1/2-vs-CDR3α TCRen potential term — signal that beats AlphaFold/TCRmodel2 confidence for ranking candidate TCRs against a fixed pMHC (denoised AUC 0.928 vs AF 0.872). Feature extraction (binder_features()) is added once its native potential term is validated; the frozen classifier is available now.

tcren.binder.binder_score(feats)[source]#

P(binder) from the 5 native descriptors (keys = FEATURES).

Parameters:: feats (dict[str, float])
Return type:: float

tcren.binder.is_real_interface(n_contacts, scanning_angle, pitch_angle)[source]#

False if the TCR:pMHC interface is assay noise / a failed dock.

NaN (or None) in any input => False (undocked). Thresholds are the ~p01-p99 range of real (binder) interfaces from the TCRvdb AF set.

Parameters:

n_contacts (float | None) – number of TCR:peptide residue-pair contacts at the interface.
scanning_angle (float | None) – TCR crossing (scanning) angle in degrees (DockingAngles.crossing_angle).
pitch_angle (float | None) – TCR incident (pitch) angle in degrees (DockingAngles.incident_angle).

Returns:

True only if n_contacts >= N_CONTACTS_MIN and scanning_angle lies in SCANNING_RANGE and pitch_angle lies in PITCH_RANGE; False otherwise, including when any input is missing.

Return type:

bool

Example

>>> is_real_interface(25, 45.0, 5.0)
True
>>> is_real_interface(0, 45.0, 5.0)
False

tcren.shuffle module#

Shuffled decoys: wrong-TCR-on-real-pMHC negatives for TCR-recognition models.

Given a set of co-framed (oriented) TCR-pMHC complexes, keep each complex’s pMHC (peptide + MHC) intact and graft on a different complex’s TCR — a within-MHC-class derangement, so no complex keeps its own TCR. The result is a physically-implausible recognition complex: a real, correctly-presented pMHC with the wrong TCR docked over it. Both the TCR:peptide and TCR:MHC interfaces are therefore mismatched, while the peptide:MHC interface is untouched (a clean internal control: peptide:MHC energy is invariant under the graft).

These decoys are the true negatives that a one-class “plausible complex” density lacks. A classifier trained on real (label 1) vs shuffled (label 0) learns TCR-recognition compatibility from structures alone, with no binding-assay labels — useful as a general, label-free recognition prior and as a supplementary benchmark.

Inputs MUST be canonically oriented (all superposed into the common MHC frame): run tcren.orient.run_folder() (tcren orient) or tcren.orient.superimpose() (tcren superimpose) first. The graft is then a direct chain replacement with no per-pair alignment — deliberately unlike tcren.orient.graft.substitute_tcr(), which superposes the donor MHC onto the host MHC pairwise. Because every complex already sits in the one canonical frame, dropping in the donor TCR as-is lets it keep its own native docking angle relative to the canonical MHC, so the decoy set spans the real MHC–TCR docking-angle variance across the whole database rather than forcing every TCR onto one host’s pose. Chains are typed by tcren.annotation.classify_chains() + tcren.mhc.annotate_mhc(), so the graft is by chain type.

CLI: tcren shuffle -s oriented/ -o shuffled/ --n 10.

tcren.shuffle.mhc_class(s)[source]#

MHC class ("MHCI"/"MHCII") from the annotated MHC-chain supertype, or None if unknown.

Parameters:: s (Structure)
Return type:: str | None

tcren.shuffle.graft_tcr(pmhc_source, tcr_source, pdb_id=None)[source]#

Build a decoy: the pMHC chains of pmhc_source + the TCR chains of tcr_source.

Both structures must already share a coordinate frame (be oriented into the canonical MHC frame). No coordinate transform is applied — chains are copied as-is, so the grafted TCR keeps its own native docking angle relative to the canonical MHC (see the module docstring for why this beats pairwise alignment). TCR chain ids that collide with a pMHC chain id are reassigned.

Parameters:

pmhc_source (Structure) – the complex whose peptide + MHC (and their coordinates) are kept.
tcr_source (Structure) – the complex whose TCR is grafted on.
pdb_id (str | None) – id for the decoy (default "<pmhc>__tcr_<tcr>").

Returns:

The decoy Structure.

Raises:

ValueError – if pmhc_source has no pMHC or tcr_source has no TCR.

Return type:

Structure

tcren.shuffle.make_decoys(structures, n_per=10, within_class=True, seed=0)[source]#

Yield decoy structures: each input pMHC paired with n_per distinct other TCRs.

Within each MHC class (if within_class) the TCR sources are a random selection excluding the pMHC’s own complex, so no decoy reproduces a real pairing. Reproducible for a given seed.

Parameters:

structures (Iterable[Structure]) – co-framed, chain-typed + MHC-annotated TCR-pMHC structures.
n_per (int) – decoys generated per input pMHC (capped at the pool size).
within_class (bool) – only graft TCRs from complexes of the same MHC class.
seed (int) – RNG seed.

Yields:

Decoy Structure objects.

Return type:

Iterator[Structure]

tcren.shuffle.run_shuffle(struct_dir, out, n=10, seed=0, within_class=True, organism='human', compress=False)[source]#

Load a folder of oriented complexes, generate n decoys per pMHC, write them, return the count.

Parameters:

struct_dir (str | Path)
out (str | Path)
n (int)
seed (int)
within_class (bool)
organism (str)
compress (bool)

Return type:

int

tcren.recognition module#

Gaussian Bayesian-network classifier: real vs shuffled TCR-pMHC complexes.

A conditional-linear-Gaussian Bayes net. A DAG is learned (BIC hill-climbing) over the standardized interface features on the within-class-centred data, so the edges capture genuine feature-feature dependence rather than the class shift. The binary class y (real = 1 / shuffled = 0) and the MHC class are then added as discrete parents of every feature node, shifting its conditional mean. Classification is the Gaussian log-likelihood ratio log P(x | y=1) - log P(x | y=0) (plus the class-prior log-odds if not balanced).

Pure numpy (dep-light). Trained parameters serialise to gzipped JSON (GaussianBNClassifier.save() / load()); to_dot() renders the network with graphviz. Trained on the Shuffled2026 decoys from tcren.shuffle.

This module also provides BayesianLogisticRecognizer — a frozen distribution-aware Bayesian logistic regression (fit externally with PyMC): each feature enters via its family’s canonical form (encode_features() — circular angles as cos/sin, bounded ratios as logit, counts/continuous linearly), so unlike the Gaussian BN it does not mis-specify the count and angle features.

class tcren.recognition.GaussianBNClassifier(feature_names, max_parents=3)[source]#

Bases: object

Conditional-linear-Gaussian BN classifier (see the module docstring).

Parameters:

feature_names (list[str])
max_parents (int)

fit(X, y, mhc_class=None)[source]#

Return type:: GaussianBNClassifier

decision_function(X, mhc_class=None)[source]#

Log-likelihood ratio log P(x|y=1) - log P(x|y=0) (balanced; add prior log-odds separately).

Return type:: ndarray

predict_proba(X, mhc_class=None, balanced=True)[source]#

Parameters:: balanced (bool)
Return type:: ndarray

marginal_decision(X, keep, mhc_class=None)[source]#

LLR log P(x_G|y=1) - log P(x_G|y=0) after marginalizing out every feature not in keep.

keep is a list of feature names (e.g. the geometry features, energy marginalised out). Because the covariance is shared across classes the marginal LLR is linear in the kept features.

Return type:: ndarray

marginal_proba(X, keep, mhc_class=None)[source]#

Return type:: ndarray

to_dict()[source]#

Return type:: dict

classmethod from_dict(d)[source]#

Parameters:: d (dict)
Return type:: GaussianBNClassifier

save(path)[source]#

Parameters:: path (str | Path)
Return type:: Path

classmethod load(path)[source]#

Parameters:: path (str | Path)
Return type:: GaussianBNClassifier

to_dot(coef_threshold=0.15)[source]#

Graphviz DAG: feature-feature edges (partial slopes) + class/MHC covariate edges above threshold.

Parameters:: coef_threshold (float)
Return type:: str

tcren.recognition.encode_features(X, feature_names)[source]#

Distribution-aware design matrix (pre-standardization).

dock_torsion (circular, wraps) -> its von Mises sufficient statistics (cos, sin); chain_balance ([0, 0.5] Beta) -> logit(2x); n_hbond dropped (exact duplicate of ct_tp_hydrogen_bond); everything else (counts + continuous + unit-vector cos/sin components) enters linearly.

Parameters:

X – (n, len(feature_names)) raw feature array.
feature_names – column names of X.

Returns:

(Z, encoded_names) — the encoded matrix and its column names.

Return type:

tuple[ndarray, list[str]]

class tcren.recognition.BayesianLogisticRecognizer(feature_names, encoded_names, mean, sd, alpha, beta, prior='normal')[source]#

Bases: object

Frozen distribution-aware Bayesian logistic (posterior-mean coefficients) — dep-light numpy predictor.

Applies encode_features(), standardizes with the stored training statistics (nan -> train mean), and returns sigmoid(alpha + Z @ beta). Fit externally by PyMC (logistic_stan/build.py in the technical appendix, which lives with the manuscript, not in this repo) and frozen here; serialises to gzipped JSON.

Parameters:: prior (str)

decision_function(X)[source]#

Return type:: ndarray

predict_proba(X)[source]#

Return type:: ndarray

to_dict()[source]#

Return type:: dict

classmethod from_dict(d)[source]#

Parameters:: d (dict)
Return type:: BayesianLogisticRecognizer

save(path)[source]#

Parameters:: path (str | Path)
Return type:: Path

classmethod load(path)[source]#

Parameters:: path (str | Path)
Return type:: BayesianLogisticRecognizer

tcren.recognition.RECOGNITION_FEATURES = ('extent', 'chain_balance', 'pitch', 'crossing', 'crossing_signed', 'dock_d', 'dock_torsion', 'dock_tcr_uy', 'dock_tcr_uz', 'dock_mhc_uy', 'dock_mhc_uz', 'F_cdr12', 'F_cdr3a', 'F_cdr3b', 'F_tcr_pep', 'F_tcr_mhc', 'F_pep_mhc', 'dF_tcr_pep', 'dF_pep_mhc', 'n_contacts_tp', 'n_pep_contacted', 'n_contacts_tm', 'ct_tp_salt_bridge', 'ct_tm_salt_bridge', 'ct_tm_hydrogen_bond', 'ct_tp_aromatic', 'ct_tm_aromatic', 'ct_tp_hydrophobic', 'ct_tm_hydrophobic', 'ct_tp_other', 'ct_tm_other', 'n_hbond', 'burial', 'mhc_class_bin')#

The core descriptor block recognize emits, and the vector the frozen recognizers consume.

Every statistical-potential energy is named F_* — there is one potential per interface (TCRen on TCR:peptide, MJ on the two presentation interfaces), so the potential’s name does not belong in the column’s. Two exact duplicates were dropped in the 2026-07-28 audit: e_tcr_mhc (the same number as F_tcr_mhc) and ct_tp_hydrogen_bond (the same number as n_hbond, which is the name Eq. Q uses). The frozen models still ask for the old names and get the same values through _FROZEN_ALIASES, so their predictions are unchanged.

tcren.recognition.INTERFACE_SYMMETRY_FEATURES = ('cdr3_dominance', 'cdr3_ab_imbalance', 'chain_cdr_imbalance')#

peptide contact counts (not energies), emitted as extra recognize output columns — not part of RECOGNITION_FEATURES (the frozen models’ 35-vector is fixed). cdr3_dominance = CDR3(α+β) share of CDR contacts (higher = CDR3-dominated, oriented positive); cdr3_ab_imbalance = |CDR3α−CDR3β| normalised (absolute); chain_cdr_imbalance = |α−β| whole-CDR normalised (absolute). See _interface_symmetry().

Type:: Interface-symmetry descriptors from per-loop TCR

tcren.recognition.FULL_FEATURES = ('extent', 'chain_balance', 'pitch', 'crossing', 'crossing_signed', 'dock_d', 'dock_torsion', 'dock_tcr_uy', 'dock_tcr_uz', 'dock_mhc_uy', 'dock_mhc_uz', 'F_cdr12', 'F_cdr3a', 'F_cdr3b', 'F_tcr_pep', 'F_tcr_mhc', 'F_pep_mhc', 'dF_tcr_pep', 'dF_pep_mhc', 'n_contacts_tp', 'n_pep_contacted', 'n_contacts_tm', 'ct_tp_salt_bridge', 'ct_tm_salt_bridge', 'ct_tm_hydrogen_bond', 'ct_tp_aromatic', 'ct_tm_aromatic', 'ct_tp_hydrophobic', 'ct_tm_hydrophobic', 'ct_tp_other', 'ct_tm_other', 'n_hbond', 'burial', 'mhc_class_bin', 'cdr3a_reach', 'cdr3a_ou', 'cdr3a_ow', 'cdr3a_on', 'cdr3a_au', 'cdr3a_aw', 'cdr3a_an', 'cdr3a_topep', 'cdr3a_ext', 'cdr3b_reach', 'cdr3b_ou', 'cdr3b_ow', 'cdr3b_on', 'cdr3b_au', 'cdr3b_aw', 'cdr3b_an', 'cdr3b_topep', 'cdr3b_ext')#

the core recognition descriptors + the 18 CDR3-frame descriptors.

The 12 “matrix-swap” columns (tcren_{g}/mj_{g}/d_{g} for g ∈ {tp, cdr12, cdr3a, cdr3b}) were removed in the 2026-07-28 audit. The tcren_* four were exact duplicates of the F_* energies; the mj_* four scored TCR:peptide contacts under the generic MJ potential, which is not the potential this method uses on that interface, and the d_* four were their difference. Nothing consumed them.

Type:: The full feature vector

tcren.recognition.DESCRIPTORS: dict[str, tuple[str, bool]] = {'F_cdr12': ('physics', True), 'F_cdr3a': ('physics', True), 'F_cdr3b': ('physics', True), 'F_pep_mhc': ('physics', False), 'F_tcr_mhc': ('physics', True), 'F_tcr_pep': ('physics', True), 'K_shear': ('kinetics', True), 'K_tens': ('kinetics', True), 'S_tot': ('kinetics', True), 'aniso': ('kinetics', True), 'burial': ('geometry', True), 'cdr3_ab_imbalance': ('geometry', True), 'cdr3_dominance': ('geometry', True), 'cdr3a_an': ('geometry', True), 'cdr3a_au': ('geometry', True), 'cdr3a_aw': ('geometry', True), 'cdr3a_ext': ('geometry', True), 'cdr3a_on': ('geometry', True), 'cdr3a_ou': ('geometry', True), 'cdr3a_ow': ('geometry', True), 'cdr3a_reach': ('geometry', True), 'cdr3a_topep': ('geometry', True), 'cdr3b_an': ('geometry', True), 'cdr3b_au': ('geometry', True), 'cdr3b_aw': ('geometry', True), 'cdr3b_ext': ('geometry', True), 'cdr3b_on': ('geometry', True), 'cdr3b_ou': ('geometry', True), 'cdr3b_ow': ('geometry', True), 'cdr3b_reach': ('geometry', True), 'cdr3b_topep': ('geometry', True), 'chain_balance': ('geometry', True), 'chain_cdr_imbalance': ('geometry', True), 'clash_score': ('geometry', True), 'couple_mhc': ('kinetics', True), 'couple_pep': ('kinetics', True), 'couple_tcr': ('kinetics', True), 'couple_total': ('kinetics', True), 'crossing': ('geometry', True), 'crossing_signed': ('geometry', True), 'ct_tm_aromatic': ('geometry', True), 'ct_tm_hydrogen_bond': ('geometry', True), 'ct_tm_hydrophobic': ('geometry', True), 'ct_tm_other': ('geometry', True), 'ct_tm_salt_bridge': ('geometry', True), 'ct_tp_aromatic': ('geometry', True), 'ct_tp_hydrophobic': ('geometry', True), 'ct_tp_other': ('geometry', True), 'ct_tp_salt_bridge': ('geometry', True), 'dF_pep_mhc': ('physics', False), 'dF_tcr_pep': ('physics', True), 'dock_d': ('geometry', True), 'dock_mhc_uy': ('geometry', True), 'dock_mhc_uz': ('geometry', True), 'dock_tcr_uy': ('geometry', True), 'dock_tcr_uz': ('geometry', True), 'dock_torsion': ('geometry', True), 'exp_lost': ('kinetics', True), 'extent': ('geometry', True), 'frac_robust': ('kinetics', True), 'lam_max': ('kinetics', True), 'lam_min': ('kinetics', True), 'mean_margin': ('kinetics', True), 'mhc_class_bin': ('geometry', False), 'n_clashes': ('geometry', True), 'n_contacts_tm': ('geometry', True), 'n_contacts_tp': ('geometry', True), 'n_hbond': ('geometry', True), 'n_interface': ('kinetics', True), 'n_pep_contacted': ('geometry', True), 'n_spring': ('kinetics', True), 'p_bind': ('score', True), 'p_forced': ('score', True), 'p_real': ('score', True), 'p_real_bn': ('score', True), 'pitch': ('geometry', True), 'q_bind': ('score', True), 'rupture_force': ('kinetics', True), 'rupture_work': ('kinetics', True), 's_strain': ('score', True)}#

Family of each descriptor, and whether the TCR enters its definition.

Three families, matching the three physical channels the method reports:

geometry — coordinates, docking angles, and the contact topology and chemistry read off them. This is the kind of quantity Eq. Q is built from.
physics — statistical-potential interface energies F and their poly-alanine references dF. Lower is more favourable.
kinetics — the interface as a network of breakable springs: stiffness, anisotropy, strain, rupture, and the residues that couple the pre-formed scaffold to the interface.

involves_tcr is False for a quantity computed from the peptide and the MHC alone. Such a column is a property of the cohort, not of the receptor: two structures of the same epitope on the same allele share its value whatever their TCR. A model handed one can reach a cohort-level label through epitope or allele identity instead of through interface physics, so any analysis whose question is about receptors must select with tcr_only=True (descriptors()).

Fitted and cohort-relative composites (p_real, p_real_bn, p_forced, p_bind, q_bind, s_strain) are listed under score. They are outputs built from the descriptors above and must never be fed back in as inputs; descriptors() excludes them by default.

tcren.recognition.descriptors(family=None, *, tcr_only=False, with_scores=False)[source]#

Descriptor names from DESCRIPTORS, filtered by family and receptor involvement.

Parameters:

family (str | None) – keep one of FAMILIES ("geometry", "physics", "kinetics"), or all of them if None.
tcr_only (bool) – keep only descriptors the receptor enters. Set this whenever the question being asked is about receptors — a peptide- or MHC-only column carries cohort identity.
with_scores (bool) – also return the fitted/cohort-relative composites of the score family. Off by default: they are model outputs, not inputs.

Returns:

The matching names, in catalogue order.

Return type:

tuple[str, …]

Example

>>> descriptors("physics", tcr_only=True)
('F_tcr_pep', 'F_tcr_mhc', 'F_cdr12', 'F_cdr3a', 'F_cdr3b', 'dF_tcr_pep')

tcren.recognition.FORCED_POSE_MODEL = {'coef': (-0.46517433874162056, 0.14437146872011086, -0.31411562068257676, -2.114810136001524, 1.198769596894963, -0.6237422800760706), 'cv_auc': 0.762, 'features': ('dock_d', 'cdr3b_reach', 'cdr3b_topep', 'cdr3a_ext', 'extent_per_ct', 'chain_balance'), 'intercept': 26.11747560652168}#

P(this pose is an AF-forced interface rather than a crystal-natural one). A raw-feature logistic (no standardization) over interface strain — stretched CDR3 loops and thin contacts. Trained ONLY on provenance (Canonical2026 crystals = 0 vs AF/TCRmodel2 models = 1; n=2681, 268 crystal / 2413 forced), so it is independent of any binder label; 5-fold CV AUC 0.762. High p_forced marks a “too-good-to-be-true” pose; the score grades crystal < AF-real < AF-decoy.

Note

For new work prefer the fit-free tcren.cohort.strain_z() (S_strain). It grades the same crystal < AF-real < AF-decoy provenance gradient by signed standardization of the strain terms, with no training set — so it is fully reproducible, unlike the coefficients below.

Warning

These coefficients are frozen and not re-derivable – the n=2681 training set no longer exists. models/fit_frozen.py::forced_pose in the benchmark repo recovers the procedure (unstandardized L2 logistic, C=0.1, which reproduces the 0.762 CV above to within 0.001) but not the coefficients. Refitting on the surviving 1168-row fixture gives a better in-sample ROC (0.769 vs 0.745), which is how we know these were fit on different rows rather than overfit to what survives. Do not replace them with a refit without re-basing the benchmarks.

Type:: Frozen “forced-pose” classifier

tcren.recognition.recognition_features(source, *, organism='human', potential=None, full=False, annotate=True)[source]#

Extract the core recognition vector from a TCR–pMHC structure (path or parsed).

Returns a dict keyed by RECOGNITION_FEATURES (degenerate/undefined terms are NaN): docking geometry, per-interface energies (raw F and poly-alanine ΔF), contact-type tallies, interface ΔSASA burial, and the mhc_class_bin indicator. The structure is chain-typed and MHC-annotated in place. DESCRIPTORS gives each column’s family and whether the receptor enters its definition.

Feed the result to frozen_recognizers() (or BayesianLogisticRecognizer) for P(real) — the probability the complex looks like a genuine TCR–pMHC recognition interface.

With full=True the row is extended with the 18 CDR3-frame descriptors (CDR3_FRAME_FEATURES) — the complete FULL_FEATURES vector.

Parameters:

organism (str)
full (bool)
annotate (bool)

Return type:

dict[str, float]

tcren.recognition.recognition_table(items, *, organism='human', full=False, scores=False, with_p_real=True, threads=1, chunk=64, autodetect_species=True, mechanics=False, _mmseqs_threads=0, _cohort_scores=True)[source]#

Batched feature (+score) extraction for a whole set of TCR–pMHC structures.

items is an iterable of (id, structure-or-path). The set is annotated with a single arda call per organism (tcren.paper.helpers._batch_annotate()) and a single mmseqs MHC search (tcren.mhc.annotate_mhc_batch()) — the dataset-scale path that avoids the per-structure annotation cost — then recognition_features() (full=) is extracted for each. With with_p_real the p_real / p_real_bn recognizer columns are added; with scores the fit-free cohort scores q_bind / s_strain (recommended, see tcren.cohort) plus the fitted p_forced / p_bind (retained for reproducibility). Returns one row dict per structure (complex.id + features [+ scores]); a structure that fails yields {"complex.id": id, "error": ...} so the batch stays resilient.

The two stages run in sequence and never compete for the machine.

Search is one arda call per organism plus one mmseqs MHC search, each given every core, over the whole set. Featurisation is where the time actually goes — a 100-pose probe spends 96 s there against 2.4 s of arda and 0.9 s of MHC search — and it is pure Python/numpy, so threads > 1 runs it in that many worker processes. The flag keeps its name for compatibility; it has always meant “how much of this machine may I use”.

It used to mean concurrent threads over chunk-sized batches, which was the wrong shape twice over: the GIL serialised the 94 % of the work that dominates, and each batch spawned its own mmseqs, so N batches asked for N x cores. Sharding the same work across independent subprocesses was measured 8x faster, which is what this now does directly.

chunk is retained for signature compatibility and is no longer used.

autodetect_species searches organism and mouse so a mis-declared cohort is still typed correctly. That doubles the annotation cost, so pass False when the organism is known — it halves the mmseqs work and changes nothing else.

mechanics appends the tcren.mechanics koff proxies (n_spring, S_tot, K_tens, K_shear, aniso, rupture_force, rupture_work, couple_*) to the same rows. They need the same annotated structure the descriptors do, so computing them here costs only their own arithmetic — running tcren mechanics separately repeats the whole parse and both mmseqs searches, and returns a second table keyed differently.

Parameters:

organism (str)
full (bool)
scores (bool)
with_p_real (bool)
threads (int)
chunk (int)
autodetect_species (bool)
mechanics (bool)
_mmseqs_threads (int)
_cohort_scores (bool)

Return type:

list[dict]

tcren.recognition.frozen_recognizers()[source]#

Load the shipped real-vs-shuffled recognizers (logistic, bn) from tcren.data (cached).

logistic is the headline distribution-aware BayesianLogisticRecognizer (shuffle_logistic.json.gz); bn is the GaussianBNClassifier (shuffle_bn.json.gz). Feed rows from recognition_features() to real_probability().

tcren.recognition.real_probability(rows, *, recognizers=None)[source]#

P(real) for feature rows from recognition_features().

rows is a dict or a list of dicts keyed by RECOGNITION_FEATURES. Returns {"logistic": p, "bn": p} — the headline logistic recognizer and the Gaussian BN, each an array of P(genuine TCR–pMHC interface). NaN features are imputed to the training mean by each model.

Both models were fitted before the duplicate columns were dropped, so they still ask for names like e_cdr12; _FROZEN_ALIASES points those at the column that now carries the value.

Return type:: dict[str, ndarray]

tcren.recognition.forced_pose_score(feats)[source]#

P(forced) — probability a pose is an AF-forced interface, from FORCED_POSE_MODEL.

feats is a row from recognition_features() with full=True (it needs the CDR3-frame cdr3b_reach/cdr3b_topep/cdr3a_ext plus core dock_d/extent/n_contacts_tp/ chain_balance). extent_per_ct is derived as extent / n_contacts_tp. Returns NaN if any required feature is missing/undefined. High = “too good to be true” (see FORCED_POSE_MODEL).

Parameters:: feats (dict[str, float])
Return type:: float

tcren.recognition.kit_score(p_bind, iptm)[source]#

Synergistic AF × tcren binder score: z(p_bind) + z(iptm) over the scored cohort.

Combines the intrinsic tcren binder score (tcren.binder.binder_score(), from recognize --scores) with the AlphaFold/TCRmodel2 ipTM that ships free with every model. On the TCRvdb raw-label benchmark this fixed no-fit combination beats either alone at precision (macro-PR 0.847 vs ipTM 0.782 / p_bind 0.804; precision 0.969 at 10% recall vs ipTM 0.861; Δ macro-PR vs ipTM +0.065, 95% CI [+0.022, +0.100], P(Δ>0)=1.00). A leave-epitope-out logistic on the same two inputs gives the more conservative +0.041 [+0.005, +0.076] — a different estimator, not this score. Higher = more binder-like.

Cohort-relative: z standardizes over the input arrays, so pass the whole set of AF models you are ranking (not one structure). NaNs are ignored by the mean/sd and propagate to their own entries.

Parameters:

p_bind – tcren binder scores for the cohort (recognize --scores p_bind column).
iptm – the matching AlphaFold ipTM values.

Returns:

The combined ranking score, one per structure.

Return type:

ndarray

tcren.mechanics module#

Interface mechanics: the TCR↔pMHC contact map as a network of breakable springs.

Each inter-body residue contact is a Hookean spring (stiffness from atomic-contact multiplicity), anchored at the two Cα atoms. Two quantities are exposed:

stiffness_tensor() — the linear-response stiffness tensor K = Σ kᵢ ûᵢ⊗ûᵢ of the interface, split into a tensile component along the docking axis and an in-plane shear component.
rupture() — a static steered-unbinding cartoon: rigidly displace the pMHC body along a pull direction, letting springs break past a strain threshold, and record the peak resisting force and the cumulative work.

Rationale (validated on ATLAS, 2026): these mechanical measures track the kinetic stability of the complex — the dissociation off-rate koff — far better than the equilibrium ΔG/Kd (Bell–Evans: rupture resistance reflects the dissociation barrier, not the well depth). This is the physically apt axis for the TCR, a mechanosensor whose pMHC bonds are catch bonds. Pure-numpy, single structure, no MD.

tcren.mechanics.WEIGHTS = ('unit', 'count', 'invdist2')#: Spring-stiffness models for an interface contact. "unit" = 1 per contact (pure topology); "count" = heavy-atom-pair multiplicity; "invdist2" = multiplicity / dist² (Hookean-ish, the validated default).

class tcren.mechanics.InterfaceSprings(a, b, k, rest, axis)[source]#

Bases: object

The TCR↔pMHC spring network of one structure.

a/b are the (n, 3) Cα anchor coordinates on the TCR and pMHC sides; k the (n,) spring stiffnesses; rest the (n,) rest lengths |b − a|; axis the unit docking axis (stiffness- weighted separation of the two interface centroids, pointing TCR→pMHC).

Parameters:

a (ndarray)
b (ndarray)
k (ndarray)
rest (ndarray)
axis (ndarray)

a: ndarray#

b: ndarray#

k: ndarray#

rest: ndarray#

axis: ndarray#

tcren.mechanics.interface_springs(structure, *, cutoff=8.0, weight='invdist2')[source]#

Build the TCR↔pMHC interface spring network from residue contacts.

Parameters:

structure (Structure) – An annotated structure (chains typed; peptide present).
cutoff (float) – Heavy-atom contact cutoff (Å) defining a spring.
weight (str) – Spring-stiffness model, one of WEIGHTS.

Returns:

The InterfaceSprings (empty arrays if no inter-body contact is found).

Return type:

InterfaceSprings

tcren.mechanics.stiffness_tensor(structure, *, cutoff=8.0, weight='invdist2')[source]#

Linear-response stiffness descriptors of the TCR↔pMHC interface.

Forms K = Σ kᵢ ûᵢ⊗ûᵢ over interface springs and resolves it along the docking axis.

Returns a dict with:: S_tot (trace K, total stiffness), K_tens (tensile, along the docking axis), K_shear (in-plane, S_tot − K_tens), aniso (K_shear / K_tens), lam_max/lam_min (extreme eigenvalues), n_spring. All nan if < 3 springs.

Parameters:

structure (Structure)
cutoff (float)
weight (str)

Return type:

dict[str, float]

tcren.mechanics.rupture(structure, *, direction='tensile', cutoff=8.0, weight='invdist2', break_strain=0.5, steps=80)[source]#

Steered-unbinding cartoon: pull the pMHC body off the TCR and measure rupture resistance.

The pMHC anchors are rigidly translated along the pull direction; each spring resists in tension (Hookean) and is removed once its strain exceeds break_strain. Integrates until all springs break.

Parameters:

direction (str) – "tensile" (docking axis; the default and the right choice for affinity/koff ranking), "shear" (stiffest in-plane), or "auto" (the minimum-force of tensile and shear — the easiest rupture path). For affinity prediction prefer "tensile" or stiffness_tensor() K_tens: on ATLAS the "auto"/"shear" force is size-confounded (its koff correlation collapses under a weighted-size control) while tensile survives it.
break_strain (float) – fractional extension at which a spring breaks (the tuning knob; 0.5 = 50 %).
steps (int) – displacement increments.
structure (Structure)
cutoff (float)
weight (str)

Return type:

dict[str, float]

Returns a dict: rupture_force (peak resisting force along the pull), rupture_work: (∫ force · displacement), n_spring, break_strain. nan if < 3 springs.

tcren.mechanics.coupling_residues(structure, *, cutoff=5.0)[source]#

Residues that couple the pre-formed internal scaffold to the TCR↔pMHC interface.

Counts residues sitting in both an intra-body contact and the binding interface (computed from the bound complex; the intra-body scaffold — Vα↔Vβ pairing, peptide-in-groove — is pre-formed, so this approximates the internal-contact residues the interface recruits):

couple_pep — peptide residues contacting both MHC (groove-anchored) and TCR (dual-role); couple_mhc — MHC residues contacting both peptide and TCR (groove-rim presenting); couple_tcr — TCR residues in the Vα↔Vβ interface that also contact pMHC (combining-site apex).

Also returns couple_total and n_interface (interface residue count, a size denominator).

Note: as a binding estimator these are weak/underpowered on current data — couple_pep tracks the dissociation off-rate (koff) at r≈−0.34 (class I, borderline; partially survives interface-size control), but the TCR/MHC sets are near-null. Useful primarily as interpretable structural descriptors.

Parameters:

structure (Structure)
cutoff (float)

Return type:

dict[str, int]

tcren.mechanics.interface_mechanics(structure, *, cutoff=8.0, weight='invdist2', direction='tensile', break_strain=0.5)[source]#

Every mechanics column for one annotated structure, in one dict.

The union of stiffness_tensor(), rupture() and coupling_residues() under their shipped defaults. This is the one definition of “the mechanics row”: both tcren mechanics and tcren recognize --mechanics call it, so the two agree by construction rather than by two call sites being kept in step.

The structure must already be annotated — chain typing and the MHC call are batch operations that belong to the caller’s single mmseqs search, not to a per-structure helper. An unannotated MHC chain silently empties the TCR:MHC half of the spring network.

Parameters:

structure (Structure) – An annotated structure (chains typed, MHC called, peptide present).
cutoff (float) – Heavy-atom contact cutoff (Å) defining a spring.
weight (str) – Spring-stiffness model, one of WEIGHTS.
direction (str) – Rupture pull direction — tensile, shear or auto.
break_strain (float) – Fractional extension at which a spring breaks.

Returns:

n_spring, S_tot, K_tens, K_shear, aniso, rupture_force, rupture_work, couple_pep/couple_mhc/couple_tcr/couple_total and n_interface.

Return type:

dict[str, float]

Example

>>> interface_mechanics(s)["K_tens"]
41.7

tcren.clashes module#

Steric-clash detection at the peptide interface.

interface_clashes() counts heavy-atom van der Waals overlaps between the peptide chain and its TCR/MHC partners — the signature of a forced or wrong-register pose, e.g. an AlphaFold / TCRmodel peptide swap that seats the peptide non-physically. A clash is a non-bonded heavy-atom pair whose separation is shorter than the sum of their Bondi vdW radii by more than a tolerance (Molprobity uses 0.4 Å for a “bad” clash; 0.6 Å here marks a severe one).

This is a structure-quality check: a generated complex with a heavy clash burden is geometrically non-physical, so its contact energy is read off a distorted interface (see tcren.refine.register for the register-specific diagnostic and correction). The pairwise overlap scan is a native _geom kernel; the numpy implementation behind it (_clash_pairs_numpy()) is kept as the reference and as a fallback where the extension is unavailable.

tcren.clashes.BONDI_RADII: dict[str, float] = {'BR': 1.85, 'C': 1.7, 'CL': 1.75, 'F': 1.47, 'H': 1.2, 'I': 1.98, 'N': 1.55, 'O': 1.52, 'P': 1.8, 'S': 1.8, 'SE': 1.9}#: Bondi van der Waals radii (Å) by element symbol; _DEFAULT_RADIUS covers anything unlisted.

class tcren.clashes.ClashPair(peptide_residue, partner_residue, partner_chain_type, overlap)[source]#

Bases: object

A single peptide↔partner clashing residue pair.

Parameters:

peptide_residue (str)
partner_residue (str)
partner_chain_type (str)
overlap (float)

peptide_residue: str#

partner_residue: str#

partner_chain_type: str#

overlap: float#

class tcren.clashes.ClashReport(n_clashes, n_severe, max_overlap, clash_score, by_partner, worst, n_peptide_atoms)[source]#

Bases: object

Interface steric-clash summary for a chain-typed complex.

Variables:

n_clashes (int) – Heavy-atom pairs overlapping by more than tolerance.
n_severe (int) – Subset overlapping by more than severe.
max_overlap (float) – Largest single overlap (Å); 0.0 if clash-free.
clash_score (float) – Sum of all overlaps > tolerance (Å) — a total-burden measure.
by_partner (dict[str, int]) – Clash count per partner chain_type (TRA/TRB/MHCa/…).
worst (tuple[tcren.clashes.ClashPair, ...]) – Up to top worst clashing residue pairs, largest overlap first.
n_peptide_atoms (int) – Heavy atoms on the peptide chain (denominator context).

Parameters:

n_clashes (int)
n_severe (int)
max_overlap (float)
clash_score (float)
by_partner (dict[str, int])
worst (tuple[ClashPair, ...])
n_peptide_atoms (int)

n_clashes: int#

n_severe: int#

max_overlap: float#

clash_score: float#

by_partner: dict[str, int]#

worst: tuple[ClashPair, ...]#

n_peptide_atoms: int#

property clashing: bool#: True if any pair overlaps by more than the tolerance.

tcren.clashes.interface_clashes(structure, *, tolerance=0.4, severe=0.6, top=8)[source]#

Detect steric clashes between the peptide chain and its partners.

Two heavy atoms on different chains clash when their separation is shorter than the sum of their Bondi vdW radii by more than tolerance Å. Only peptide↔partner pairs are examined (partners = every non-peptide chain); intra-chain and peptide-internal pairs are ignored, so no bonded pair is ever counted.

Parameters:

structure (Structure) – A chain-typed complex with a peptide chain (chain_type == 'PEPTIDE').
tolerance (float) – vdW overlap (Å) above which a pair counts as a clash (Molprobity 0.4).
severe (float) – overlap (Å) above which a clash is also counted as severe.
top (int) – How many worst residue pairs to return.

Returns:

A ClashReport.

Raises:

ValueError – If the structure has no peptide chain.

Return type:

ClashReport

tcren.clashes.has_clashes(structure, *, tolerance=0.4)[source]#

Convenience predicate: does the peptide interface have any clash beyond tolerance?

Parameters:

structure (Structure)
tolerance (float)

Return type:

bool

tcren.stability module#

Contact stability / fragility at the TCR:peptide interface.

The 5 Å contact cutoff is a hard threshold on a continuous distance: a contact at 4.9 Å has almost no positional slack (a ~1 Å shift, or tightening the cutoff, kills it), while one at 3.5 Å is robust. contact_stability() reads that slack directly off the contact map. For each TCR:peptide contact — a receptor-residue / peptide-residue pair whose closest heavy-atom pair is within cutoff (the same contact definition as tcren.contacts.geometry.all_atom_contacts()) — the margin m = cutoff - dmin is the contact’s positional tolerance, and a rigid isotropic shift of size delta loses it with probability clip((delta - m) / (2*delta), 0, 1). Aggregated over the interface these give mean_margin (how deep the contacts sit), frac_robust (fraction with m >= delta) and exp_lost (expected contacts lost under a delta shift) — a coordinate-only read of interface positional confidence, the physical analogue of an interface PAE.

The per-residue-pair minimum-distance scan is a native _geom kernel; the numpy implementation behind it (_contact_stability_numpy()) is the reference and a fallback where the extension is unavailable.

class tcren.stability.StabilityReport(n_contacts, mean_margin, frac_marg_lt1, frac_robust, exp_lost)[source]#

Bases: object

TCR:peptide contact-stability summary.

Variables:

n_contacts (int) – Receptor-residue / peptide-residue contacts within cutoff.
mean_margin (float) – Mean cutoff - dmin over contacts (Å); larger = contacts sit deeper.
frac_marg_lt1 (float) – Fraction of contacts with margin below delta (fragile).
frac_robust (float) – Fraction of contacts with margin at least delta (robust to a delta shift).
exp_lost (float) – Expected number of contacts lost under a delta Å isotropic shift.

Parameters:

n_contacts (int)
mean_margin (float)
frac_marg_lt1 (float)
frac_robust (float)
exp_lost (float)

n_contacts: int#

mean_margin: float#

frac_marg_lt1: float#

frac_robust: float#

exp_lost: float#

tcren.stability.contact_stability(structure, *, cutoff=5.0, delta=1.0)[source]#

Contact stability / fragility of the TCR:peptide interface.

Parameters:

structure (Structure) – A chain-typed complex with a peptide chain and at least one receptor chain.
cutoff (float) – Contact distance cutoff (Å); the closest heavy-atom pair of a residue pair must be within this to count as a contact (matches all_atom_contacts).
delta (float) – Positional shift (Å) used for the fragility metrics (frac_marg_lt1, exp_lost).

Returns:

A StabilityReport.

Raises:

ValueError – If the structure has no peptide chain or no receptor chain.

Return type:

StabilityReport

tcren.pipeline module#

End-to-end TCRen pipeline: structure → annotation → orientation → contacts → score.

One call takes a TCR-pMHC structure all the way through the tcren workflow:

import the structure (C-gene trimmed);
annotate chains — TCR loci/CDRs via arda, MHC allele/class/role + groove regions;
superimpose onto the canonical database (canonical Cα frame; optional);
markup + contacts — the per-residue region table and the 5 Å contact map;
score each interface with its potential: TCRen for TCR↔peptide, MJ for TCR↔MHC and peptide↔MHC, plus the total.

The interface energy is the sum of the residue-pair potential over the observed contacts of that interface (the closest-atom contact per residue pair, as everywhere in tcren).

class tcren.pipeline.PipelineResult(pdb_id, mhc_calls, markup, contacts, scores, oriented=None, rmsd=None, extra=<factory>)[source]#

Bases: object

Everything the pipeline produces for one structure.

extra carries the interface-sanity flag when the complex was oriented (superimpose=True): real_interface (bool — False marks assay noise / a failed dock; see tcren.binder.is_real_interface()), and the raw descriptors it was computed from (n_contacts, scanning_angle, pitch_angle). With superimpose=False the docking angles are unavailable, so real_interface is None (scanning_angle/pitch_angle None) while n_contacts is still set.

Parameters:

pdb_id (str)
mhc_calls (list[MhcCall])
markup (DataFrame)
contacts (DataFrame)
scores (dict[str, float])
oriented (Structure | None)
rmsd (float | None)
extra (dict)

pdb_id: str#

mhc_calls: list[MhcCall]#

markup: DataFrame#

contacts: DataFrame#

scores: dict[str, float]#

oriented: Structure | None#

rmsd: float | None#

extra: dict#

tcren.pipeline.run(structure, organism='human', superimpose=True, db_dir=None, cutoff=5.0, potentials=None, tcr_regions='all', contact_weight='residue', reference_aa=None)[source]#

Run the full pipeline on one structure (path or parsed Structure).

Parameters:

structure (str | Path | Structure) – a structure file (any tcren-readable format) or an already-parsed structure.
organism (str) – organism for TCR annotation.
superimpose (bool) – also orient onto the canonical database (sets oriented + rmsd).
db_dir (str | Path | None) – canonical database for superimpose (default data/Canonical2026).
cutoff (float) – contact distance threshold (Å).
potentials (dict[str, str | Potential | None] | None) – optional per-interface potential override mapping an interface name ("tcr_peptide", "tcr_mhc", "peptide_mhc") to a Potential, a bundled name ("tcren"/"mj"/"keskin"), a CSV path, or None. None (or a missing entry) keeps the default family for that interface, so the default output is unchanged.
tcr_regions (str) – which TCR regions to keep on the TCR side of the TCR-containing interfaces ("all" default = no filter = legacy behaviour; "cdr" or "cdr+fr" to restrict).
contact_weight (str) – "residue" (default, legacy) weights each contacting residue pair by 1 on all three interfaces; "atomic" weights each pair by its n_atom_contacts heavy-atom-pair count (the contact map is then built with count_atoms=True). Applies to tcr_peptide, tcr_mhc and peptide_mhc alike.
reference_aa (str | None) – if set (typically "A"), also report the reference-normalised energies delta_<interface> and delta_total — each interface’s tcren.ddg.reference_delta(), i.e. its energy minus the energy of a poly-reference_aa peptide threaded onto the same contact map. Off by default, so the default scores dict is unchanged.

Returns:

A PipelineResult with the markup, contacts, per-interface scores and (if requested) the canonical-frame oriented structure.

Return type:

PipelineResult

tcren.pipeline.score_row(result)[source]#

Flatten a PipelineResult to a one-row scores dict (for a CSV table).

The d_* reference-normalised columns are present only when the pipeline was run with reference_aa set.

Parameters:: result (PipelineResult)
Return type:: dict

tcren.oracle module#

One-call facade over the tcren pipeline for the paper notebooks.

summarize_structure() turns a single TCR-pMHC structure into a bundle of ready-to-tabulate polars frames by composing the existing milestones:

S1+S2 — tcren.pipeline.run() (annotate → orient → contacts → per-interface scores), giving the scores, markup and contacts frames;
S3 — tcren.scoring_rank.percentile_rank() of the structure’s native peptide against a random pMHC background, giving the rank frame;
S4 — tcren.ddg.alanine_scan() of the native peptide, giving the ddg frame.

Nothing is re-derived here: the facade only orchestrates the milestone functions and collects their outputs. The scores frame is therefore byte-identical to what run produces for the same structure and arguments.

tcren.oracle.summarize_structure(structure, *, organism='human', superimpose=True, potentials=None, tcr_regions='all', background=1000, seed=0, alanine=False, contact_weight='residue')[source]#

Summarise one TCR-pMHC structure into a bundle of tables.

Composes the tcren milestones on a single structure: the full pipeline (S1+S2) for per-interface energies, region markup and contacts; the percentile rank of the native peptide against a random background (S3); and, optionally, the per-position alanine scan (S4).

Parameters:

structure (str | Path | Structure) – A structure file (any tcren-readable format) or a parsed Structure.
organism (str) – Organism for TCR annotation.
superimpose (bool) – Also orient onto the canonical database (adds rmsd to scores).
potentials (dict[str, str | Potential | None] | None) – Optional per-interface potential override, forwarded to tcren.pipeline.run(). The TCR↔peptide potential is also used for the rank and ddg frames.
tcr_regions (str) – Which TCR regions to keep on the TCR side ("all" default, "cdr" or "cdr+fr"); forwarded to every milestone.
background (int) – Number of random background peptides for the percentile rank (S3).
seed (int) – Random seed for the background generation.
alanine (bool) – When True, compute the per-position alanine scan (S4) for the ddg frame. When False (default) the scan is skipped and ddg is an empty frame with the same schema.
contact_weight (str) – "residue" (default, legacy: each contacting residue pair contributes potential[a, b] x 1) or "atomic" (each pair contributes potential[a, b] x n_atom_contacts). Applies to all three interfaces of the scores frame and to the rank/ddg frames. "residue" keeps every output byte-identical to the legacy facade.

Returns:

scores — one row of per-interface energies (and rmsd if superimposed), identical to tcren.pipeline.run()’s scores.
rank — one row: the native peptide’s energy and its rank_pct against the background.
ddg — the alanine scan (columns pos/wt_aa/ddG); empty unless alanine=True.
markup — the per-residue region-markup table.
contacts — the annotated residue-contact table.

Return type:

Mapping with five polars frames

tcren.refine package#

Peptide substitution + potential-guided refinement.

substitute_peptide() threads a new sequence onto the peptide backbone; refine_peptide() runs a knowledge-based rigid-body Monte-Carlo refinement of the peptide pose via the compiled tcren._refine kernel. The refinement energy is the DOPE atom-level distance-dependent statistical potential (Shen & Sali, Protein Science 2006), used here independently of the TCRen/MJ potentials tcren scores epitopes with — so the pose is not optimised against the same quantity it is later scored with. This is a lightweight, knowledge-based refine, NOT physics relaxation (that is Rosetta FlexPepDock, as a subprocess).

tcren.refine.substitute_peptide(structure, new_peptide, chain_type='PEPTIDE')[source]#

Return a copy of structure with the peptide chain threaded to new_peptide.

The peptide backbone (and Cβ) is preserved; side-chain atoms beyond Cβ are dropped (and Cβ too for any position mutated to glycine). new_peptide must equal the peptide length and use the 20 standard one-letter amino acids.

Raises:

ValueError – if there is no peptide chain, the length differs, or a code is non-standard.

Parameters:

structure (Structure)
new_peptide (str)
chain_type (str)

Return type:

Structure

tcren.refine.refine_peptide(structure, *, shell=12.0, restraint_w=0.5, n_steps=2000, trans_sigma=0.2, rot_sigma=0.05, temp0=1.0, temp1=0.05, seed=0)[source]#

Rigid-body refine the peptide pose against its TCR+MHC partners; (structure, energy).

The energy is the DOPE atom-level distance-dependent statistical potential summed over all peptide$leftrightarrow$partner heavy-atom pairs within DOPE’s range (its short-range bins are repulsive, so it provides its own clash term), plus a harmonic restraint to the input pose (restraint_w) that keeps the search local. Only partner atoms within shell Å of the peptide are considered. The structure must be chain-typed (peptide = chain of chain_type == 'PEPTIDE'). Requires the compiled _refine ext + the bundled DOPE table.

Parameters:

structure (Structure)
shell (float)
restraint_w (float)
n_steps (int)
trans_sigma (float)
rot_sigma (float)
temp0 (float)
temp1 (float)
seed (int)

tcren.refine.interface_energy(structure, *, shell=12.0)[source]#

DOPE interaction energy across the peptide↔partner interface (lower = more favourable).

The structure must be chain-typed (peptide = chain_type == 'PEPTIDE'). Only partner heavy atoms within shell Å of the peptide are considered (beyond DOPE’s range they contribute nothing). Requires the compiled tcren._relax ext + the bundled DOPE table.

Parameters:

structure (Structure)
shell (float)

Return type:

float

tcren.refine.predict_anchors(peptide, structure=None)[source]#

Predict peptide anchors, preferring the structure’s MHC-class call over the length heuristic.

Pass a chain-typed, MHC-annotated structure (see tcren.mhc.annotate_mhc()) to use its real class assignment; otherwise the class is inferred from peptide length.

Parameters:

peptide (str)
structure (Structure | None)

Return type:

Decomposition

tcren.refine.native_peptide(structure)[source]#

One-letter sequence of the structure’s PEPTIDE chain (raises if absent).

Parameters:: structure (Structure)
Return type:: str

class tcren.refine.Decomposition(peptide, mhc_class, anchors, tcr_facing, presentation)[source]#

Bases: object

Anchor vs TCR-facing split of a peptide.

Parameters:

peptide (str)
mhc_class (str)
anchors (tuple[int, ...])
tcr_facing (str)
presentation (str)

peptide: str#

mhc_class: str#

anchors: tuple[int, ...]#

tcr_facing: str#

presentation: str#

tcren.refine.peptide_rmsd(model, reference, anchors=())[source]#

Peptide backbone / Cα / anchor-Cα RMSD of model vs reference in the MHC-groove frame.

Both structures must be chain-typed and MHC-annotated. model is superposed onto reference by their shared groove Cα; the peptide RMSDs are then computed over atoms (and residues) present in both. anchors are 0-based peptide residue indices for the anchor-Cα RMSD.

Parameters:

model (Structure)
reference (Structure)
anchors (tuple[int, ...])

Return type:

PeptideRMSD

class tcren.refine.PeptideRMSD(backbone_rmsd, ca_rmsd, anchor_ca_rmsd, n_backbone, n_ca, groove_rmsd)[source]#

Bases: object

Peptide RMSD of a model against a reference, after MHC-groove superposition.

Parameters:

backbone_rmsd (float)
ca_rmsd (float)
anchor_ca_rmsd (float)
n_backbone (int)
n_ca (int)
groove_rmsd (float)

backbone_rmsd: float#

ca_rmsd: float#

anchor_ca_rmsd: float#

n_backbone: int#

n_ca: int#

groove_rmsd: float#

tcren.refine.check_register(model, reference=None, *, anchor_rmsd_cut=2.0, tolerance=0.4, top=8)[source]#

Assess whether model’s peptide is clashing and/or in the wrong register.

Parameters:

model (Structure) – The generated complex to check (chain-typed; peptide chain_type == 'PEPTIDE').
reference (Structure | None) – A correctly-registered pose of the same complex (crystal / trusted model). When given, the anchor-Cα RMSD in the MHC-groove frame decides wrong_register.
anchor_rmsd_cut (float) – Anchor-Cα RMSD (Å) above which the register is called wrong.
tolerance (float) – vdW-overlap tolerance passed to tcren.clashes.interface_clashes().
top (int) – Worst clashing pairs to keep in the report.

Returns:

A RegisterReport.

Return type:

RegisterReport

tcren.refine.fix_register(model, template, *, engine='ccd', seed=0, **engine_kwargs)[source]#

Correct a wrong-register model by re-threading its peptide onto template’s register.

The model’s peptide sequence is threaded onto the template’s correctly-registered backbone (backbone + Cβ preserved) and re-refined by the chosen open-source engine, so the returned complex carries the model’s peptide in the template’s groove pose. This is the FlexPepDock-functional correction path (tcren.refine.model_peptide()); no PyRosetta.

Parameters:

model (Structure) – The wrong-register complex (source of the peptide sequence).
template (Structure) – A correctly-registered complex of equal peptide length (source of the pose), e.g. the crystal structure of the same clone.
engine (str) – Modelling engine (ccd/dope run out of the box; openmm/promod3 optional).
seed (int) – RNG seed for engines that sample.
**engine_kwargs – Forwarded to the engine (e.g. n_steps for dope).

Returns:

A tcren.refine.engines.base.ModelResult with the corrected structure, engine score, and anchors used.

Raises:

ValueError – If the model and template peptides differ in length (equal length is required for backbone-preserving substitution).

class tcren.refine.RegisterReport(clashes, anchors, backbone_rmsd, anchor_rmsd, groove_rmsd, wrong_register, reason)[source]#

Bases: object

Variables:

clashes (tcren.clashes.ClashReport) – Interface steric-clash burden of the model (always computed).
anchors (tuple[int, ...]) – 0-based peptide anchor indices used for the anchor RMSD.
backbone_rmsd (float) – Peptide backbone RMSD vs the reference (Å); nan without a reference.
anchor_rmsd (float) – Anchor-Cα RMSD vs the reference (Å) — the register signal; nan without one.
groove_rmsd (float) – MHC-groove superposition residual (Å); nan without a reference.
wrong_register (bool | None) – True/False when a reference is given (anchor RMSD over the cut), None when it cannot be determined (no reference).
reason (str) – Human-readable explanation of the verdict.

Parameters:

clashes (ClashReport)
anchors (tuple[int, ...])
backbone_rmsd (float)
anchor_rmsd (float)
groove_rmsd (float)
wrong_register (bool | None)
reason (str)

clashes: ClashReport#

anchors: tuple[int, ...]#

backbone_rmsd: float#

anchor_rmsd: float#

groove_rmsd: float#

wrong_register: bool | None#

reason: str#

property suspect: bool#: True if the pose is a wrong register, or (absent a reference) carries a severe clash.

Backbone-preserving peptide substitution.

score_peptides scores a candidate peptide virtually (it re-indexes the potential matrix over the native contact map — no atoms move). When you want to actually re-dock / refine a candidate you first need its coordinates: substitute_peptide() threads an equal-length sequence onto the existing peptide backbone, keeping N/Cα/C/O(+Cβ) and dropping the old side-chain atoms beyond Cβ (a refiner / rotamer repack rebuilds them). Pure data-model manipulation; returns a new structure.

tcren.refine.substitute.substitute_peptide(structure, new_peptide, chain_type='PEPTIDE')[source]#

Return a copy of structure with the peptide chain threaded to new_peptide.

Raises:

ValueError – if there is no peptide chain, the length differs, or a code is non-standard.

Parameters:

structure (Structure)
new_peptide (str)
chain_type (str)

Return type:

Structure

Peptide register / forced-pose detection and correction.

A generated (AlphaFold / TCRmodel) TCR:pMHC complex can seat the peptide in the wrong register — anchors not in their MHC pockets, the peptide shifted along the groove — which corrupts the TCR-facing contacts the score reads. On the CPL benchmark this is the ila1 anomaly: the forced AF pose drives the raw TCRen ranking below chance (ROC 0.35), and re-seating the peptide on the correctly registered crystal recovers it (≈0.77).

check_register() flags such poses. It always reports the interface clash burden (via tcren.clashes.interface_clashes()); when a correctly-registered reference (a crystal, or any trusted pose of the same complex) is supplied it also measures the anchor-Cα RMSD in the MHC-groove frame (tcren.refine.peptide_rmsd()) — the reliable register signal, since a mis-registered peptide’s anchors land far from where they should sit. Note a heavy clash burden alone does not prove a register error: AlphaFold peptide-swap models are routinely clashy, so register needs the reference.

fix_register() corrects a wrong-register model by re-threading its peptide sequence onto the correctly-registered template backbone and re-refining through the open-source modelling path (tcren.refine.model_peptide()) — the FlexPepDock-functional replacement, no PyRosetta.

class tcren.refine.register.RegisterReport(clashes, anchors, backbone_rmsd, anchor_rmsd, groove_rmsd, wrong_register, reason)[source]#

Bases: object

Variables:

clashes (tcren.clashes.ClashReport) – Interface steric-clash burden of the model (always computed).
anchors (tuple[int, ...]) – 0-based peptide anchor indices used for the anchor RMSD.
backbone_rmsd (float) – Peptide backbone RMSD vs the reference (Å); nan without a reference.
anchor_rmsd (float) – Anchor-Cα RMSD vs the reference (Å) — the register signal; nan without one.
groove_rmsd (float) – MHC-groove superposition residual (Å); nan without a reference.
wrong_register (bool | None) – True/False when a reference is given (anchor RMSD over the cut), None when it cannot be determined (no reference).
reason (str) – Human-readable explanation of the verdict.

Parameters:

clashes (ClashReport)
anchors (tuple[int, ...])
backbone_rmsd (float)
anchor_rmsd (float)
groove_rmsd (float)
wrong_register (bool | None)
reason (str)

clashes: ClashReport#

anchors: tuple[int, ...]#

backbone_rmsd: float#

anchor_rmsd: float#

groove_rmsd: float#

wrong_register: bool | None#

reason: str#

property suspect: bool#: True if the pose is a wrong register, or (absent a reference) carries a severe clash.

tcren.refine.register.check_register(model, reference=None, *, anchor_rmsd_cut=2.0, tolerance=0.4, top=8)[source]#

Assess whether model’s peptide is clashing and/or in the wrong register.

Parameters:

model (Structure) – The generated complex to check (chain-typed; peptide chain_type == 'PEPTIDE').
reference (Structure | None) – A correctly-registered pose of the same complex (crystal / trusted model). When given, the anchor-Cα RMSD in the MHC-groove frame decides wrong_register.
anchor_rmsd_cut (float) – Anchor-Cα RMSD (Å) above which the register is called wrong.
tolerance (float) – vdW-overlap tolerance passed to tcren.clashes.interface_clashes().
top (int) – Worst clashing pairs to keep in the report.

Returns:

A RegisterReport.

Return type:

RegisterReport

tcren.refine.register.fix_register(model, template, *, engine='ccd', seed=0, **engine_kwargs)[source]#

Correct a wrong-register model by re-threading its peptide onto template’s register.

Parameters:

model (Structure) – The wrong-register complex (source of the peptide sequence).
template (Structure) – A correctly-registered complex of equal peptide length (source of the pose), e.g. the crystal structure of the same clone.
engine (str) – Modelling engine (ccd/dope run out of the box; openmm/promod3 optional).
seed (int) – RNG seed for engines that sample.
**engine_kwargs – Forwarded to the engine (e.g. n_steps for dope).

Returns:

A tcren.refine.engines.base.ModelResult with the corrected structure, engine score, and anchors used.

Raises:

ValueError – If the model and template peptides differ in length (equal length is required for backbone-preserving substitution).

Peptide anchor prediction (which residues bury into the MHC groove).

Class I uses fixed positions counted from both termini (P2 + the C-terminal PΩ); class II has no fixed register, so a one-pass heuristic slides a 9-mer window and scores each register by a P1-hydrophobic / P4,P6,P9-not-Pro/Gly rule (a cheap proxy for GibbsCluster/NNAlign register inference), then reports P1/P4/P6/P9 of the best core. The anchor residues are the groove-facing (“presentation”) side; the rest of the peptide is TCR-facing.

This logic is portable and depends only on the standard library. It is lifted from the antigenomics mhcmatch/seqtree decomposition primitives (which themselves carry no tcren dependency at runtime); predict_anchors() adds a thin tcren wrapper that prefers the structure’s own MHC-class call (from tcren.mhc) over the length heuristic when a typed structure is available.

tcren.refine.anchors.infer_class(peptide)[source]#

Heuristic MHC class from peptide length (≤ 11 → class I, else class II).

Parameters:: peptide (str)
Return type:: str

tcren.refine.anchors.anchor_indices(peptide, mhc_class)[source]#

0-based anchor positions: class-I P2/PΩ, class-II core P1/P4/P6/P9. mhc_class in {MHCI, MHCII}.

Parameters:

peptide (str)
mhc_class (str)

Return type:

tuple[int, …]

class tcren.refine.anchors.Decomposition(peptide, mhc_class, anchors, tcr_facing, presentation)[source]#

Bases: object

Anchor vs TCR-facing split of a peptide.

Parameters:

peptide (str)
mhc_class (str)
anchors (tuple[int, ...])
tcr_facing (str)
presentation (str)

peptide: str#

mhc_class: str#

anchors: tuple[int, ...]#

tcr_facing: str#

presentation: str#

tcren.refine.anchors.decompose(peptide, mhc_class=None)[source]#

Split peptide into anchor (groove-facing) and TCR-facing parts.

Parameters:

peptide (str)
mhc_class (str | None)

Return type:

Decomposition

tcren.refine.anchors.predict_anchors(peptide, structure=None)[source]#

Predict peptide anchors, preferring the structure’s MHC-class call over the length heuristic.

Pass a chain-typed, MHC-annotated structure (see tcren.mhc.annotate_mhc()) to use its real class assignment; otherwise the class is inferred from peptide length.

Parameters:

peptide (str)
structure (Structure | None)

Return type:

Decomposition

tcren.refine.anchors.native_peptide(structure)[source]#

One-letter sequence of the structure’s PEPTIDE chain (raises if absent).

Parameters:: structure (Structure)
Return type:: str

Peptide RMSD between two poses of the same complex.

The benchmark question is “how close is a re-modelled peptide to its native crystal pose?” We answer it the way the rest of the codebase measures geometry: superpose the two complexes on the conserved MHC groove Cα (tcren.orient.align._matched_anchors() + Biopython’s SVDSuperimposer), then compute the peptide backbone RMSD in that common MHC frame. Superposing on the MHC — not on the peptide — means the peptide RMSD reflects how well the pose sits in the groove, which is the quantity that matters for downstream TCR-contact scoring.

class tcren.refine.rmsd.PeptideRMSD(backbone_rmsd, ca_rmsd, anchor_ca_rmsd, n_backbone, n_ca, groove_rmsd)[source]#

Bases: object

Peptide RMSD of a model against a reference, after MHC-groove superposition.

Parameters:

backbone_rmsd (float)
ca_rmsd (float)
anchor_ca_rmsd (float)
n_backbone (int)
n_ca (int)
groove_rmsd (float)

backbone_rmsd: float#

ca_rmsd: float#

anchor_ca_rmsd: float#

n_backbone: int#

n_ca: int#

groove_rmsd: float#

tcren.refine.rmsd.peptide_rmsd(model, reference, anchors=())[source]#

Peptide backbone / Cα / anchor-Cα RMSD of model vs reference in the MHC-groove frame.

Parameters:

model (Structure)
reference (Structure)
anchors (tuple[int, ...])

Return type:

PeptideRMSD

Interface interaction energy — the ΔΔG core (native, DOPE, no PyRosetta).

interface_energy(structure) sums the DOPE atom-level statistical potential over peptide↔partner heavy-atom pairs — the interaction energy across the TCR/MHC:peptide interface. Because only the cross (peptide↔partner) terms are summed, it is already E_bound − E_separated for the interaction, i.e. the interface ΔG contribution of the peptide. This is the quantity a sampling ΔΔG differences between a mutant and the native (each first repacked/relaxed); see the _relax kernel roadmap in CPP_REWRITE.md.

Reuses the bundled DOPE table + atom-class map from tcren.refine (_dope), and the compiled tcren._relax kernel.

tcren.refine.interface.interface_energy(structure, *, shell=12.0)[source]#

DOPE interaction energy across the peptide↔partner interface (lower = more favourable).

Parameters:

structure (Structure)
shell (float)

Return type:

float

Open-source peptide (re)modelling: substitute → predict anchors → model in the groove.

model_peptide() is the single entry point that replaces the licensed FlexPepDock/MODELLER path. It threads a candidate peptide onto the groove backbone (tcren.refine.substitute_peptide()), predicts which residues anchor into the MHC (tcren.refine.anchors.predict_anchors()), and hands off to a modelling engine (dope and ccd run out of the box; openmm and promod3 are optional). The returned ModelResult carries the refined structure, an engine score, and the anchor set used.

tcren.refine.model.model_peptide(structure, new_peptide=None, *, engine='ccd', seed=0, **engine_kwargs)[source]#

Model new_peptide into structure’s groove with the chosen engine.

Parameters:

structure (Structure) – A chain-typed, MHC-annotated complex (peptide chain chain_type == 'PEPTIDE').
new_peptide (str | None) – Candidate one-letter sequence (equal length to the native peptide). None re-models the native sequence in place — the self-reconstruction case used by the benchmark.
engine (str) – One of ENGINES (dope, ccd, openmm, promod3).
seed (int) – RNG seed (engines that sample).
**engine_kwargs – Forwarded to the engine’s run (e.g. perturb/anchor_targets for ccd, n_steps for dope).

Returns:

A ModelResult with the refined structure, engine score, and anchors used.

Return type:

ModelResult

class tcren.refine.model.ModelResult(structure, energy, engine, anchors, iterations=0, info=<factory>)[source]#

Bases: object

Outcome of (re)modelling a peptide into the groove.

Parameters:

structure (Structure)
energy (float)
engine (str)
anchors (tuple[int, ...])
iterations (int)
info (dict)

structure: Structure#

energy: float#

engine: str#

anchors: tuple[int, ...]#

iterations: int#

info: dict#

FlexPepDock oracle — the accuracy ceiling the open engines are measured against.

Rosetta FlexPepDock is the reference-standard peptide refiner. It is a protocol inside Rosetta, so there are two ways to reach it:

PyRosetta API (preferred here) — pyrosetta.rosetta.protocols.flexpep_docking. Installed via pyrosetta-installer (academic license). No external binary needed.
External binary — a licensed FlexPepDocking executable via $ROSETTA_BIN / rosetta_bin=.

Either way this is an oracle, never a shipped tcren dependency: the RMSD FlexPepDock achieves from a displaced start is the lower bound on error the license-free engines (ccd/dope/openmm) should approach (see CPP_REWRITE.md). Refinement is slow (tens of seconds to minutes per structure), so the oracle column in fold_benchmark.py is meant for a subset, not the full sweep.

KNOWN ISSUE (2026-07): via the PyRosetta path below, FlexPepDockingProtocol().apply(pose) is a no-op on our 5-chain TCR-pMHC complexes — the dumped pose equals the input (peptide RMSD Δ=0.000 vs the un-refined baseline across every test pair). FlexPepDock refine needs the receptor↔peptide jump / FoldTree set up (the command-line app infers it from the last-chain-peptide convention; bare protocol construction here does not). The oracle is therefore currently NON-FUNCTIONAL and its numbers must NOT be reported as FlexPepDock accuracy. Fix by configuring the FlexPepDock FoldTree + -flexPepDocking flags (or use the external binary path) before trusting the column.

tcren.refine.oracle_flexpep.pyrosetta_available()[source]#

Return type:: bool

tcren.refine.oracle_flexpep.flexpep_available(rosetta_bin=None)[source]#

True if FlexPepDock is reachable (PyRosetta installed, or a binary resolvable).

Parameters:: rosetta_bin (str | None)
Return type:: bool

tcren.refine.oracle_flexpep.flexpep_refine(structure, *, rosetta_bin=None, seed=0, nstruct=1, extra_flags=())[source]#

Refine the peptide pose with FlexPepDock; return the refined Structure.

Uses the PyRosetta API if installed, else a resolvable FlexPepDocking binary. The structure should carry full chains (imported with keep_c_gene=True). Raises RuntimeError if neither route is available.

Parameters:

structure (Structure)
rosetta_bin (str | None)
seed (int)
nstruct (int)
extra_flags (tuple[str, ...])

Return type:

Structure

Common types for peptide-modelling engines.

Every engine takes a chain-typed Structure whose peptide has already been threaded to the candidate sequence (via tcren.refine.substitute_peptide()) plus the predicted anchors, and returns a ModelResult: the refined structure, an engine-specific energy/score, and bookkeeping. Engines that need an unavailable dependency raise EngineUnavailable at call time (never at import), so import tcren always succeeds.

exception tcren.refine.engines.base.EngineUnavailable[source]#

Bases: RuntimeError

Raised when an engine’s backend (OpenMM, ProMod3, …) is not installed.

class tcren.refine.engines.base.ModelResult(structure, energy, engine, anchors, iterations=0, info=<factory>)[source]#

Bases: object

Outcome of (re)modelling a peptide into the groove.

Parameters:

structure (Structure)
energy (float)
engine (str)
anchors (tuple[int, ...])
iterations (int)
info (dict)

structure: Structure#

energy: float#

engine: str#

anchors: tuple[int, ...]#

iterations: int#

info: dict#

class tcren.refine.engines.base.Engine(*args, **kwargs)[source]#

Bases: Protocol

A peptide-modelling backend.

name: str#

available()[source]#

True if the backend can run in this environment (cheap import probe, no side effects).

Return type:: bool

run(structure, decomp, *, seed=0)[source]#

Model the (already substituted) peptide; return the refined structure + score.

Parameters:

structure (Structure)
decomp (Decomposition)
seed (int)

Return type:

ModelResult

DOPE rigid-body Monte-Carlo engine (wraps the existing tcren._refine kernel).

This is the engine that already ships: a knowledge-based rigid-body refine of the peptide pose against its TCR+MHC partners under the DOPE statistical potential. It does not use the anchors (the whole peptide moves as a rigid body), but it is the always-available baseline every other engine is compared against. See tcren.refine.refine_peptide().

class tcren.refine.engines.dope.DopeEngine[source]#

Bases: object

name = 'dope'#

available()[source]#

Return type:: bool

run(structure, decomp, *, seed=0, n_steps=2000)[source]#

Parameters:

structure (Structure)
decomp (Decomposition)
seed (int)
n_steps (int)

Return type:

ModelResult

CCD anchor-restrained closure engine (wraps the tcren._fold C++ kernel).

Cyclic Coordinate Descent drives the peptide’s anchor Cα onto target positions (predicted MHC-groove pocket centroids) while the rest of the backbone follows as a kinematic linkage. This is the license-free geometric path: no Rosetta, no MODELLER, only the bundled stdlib-only _fold kernel.

Draft scope (ponytail: marked so the simplification reads as intent, not ignorance):

Operates on the Cα trace with consecutive-Cα rotatable bonds, and writes the closed pose back by rigid per-residue translation (each residue’s atoms shift by its Cα displacement). The kernel preserves Cα–Cα distances exactly, but because adjacent residues receive different translations the inter-residue peptide-bond geometry (C(i)–N(i+1) ≈ 1.33 Å, φ/ψ) is only approximate — intra- residue geometry is intact, the chain as a whole is a Cα-trace model, NOT a physically valid all-atom backbone. The output therefore MUST be followed by an energy refine (DOPE / OpenMM) to regularise the peptide bonds; the upgrade path is a full N–Cα–C kinematic chain + rotamer repack (OpenMM/ProMod3 do this). Do not treat the ccd output, on its own, as a finished structure.
Anchor targets must be supplied. In the self-reconstruction benchmark they are the native anchor Cα; predicting pocket centroids de novo (groove-pocket geometry) is the open piece flagged in STATUS.md and left to the scoring/orient layer.

class tcren.refine.engines.ccd.CcdEngine[source]#

Bases: object

name = 'ccd'#

available()[source]#

Return type:: bool

run(structure, decomp, *, seed=0, anchor_targets=None, perturb=0.0, max_iter=1000, tol=0.08)[source]#

Parameters:

structure (Structure)
decomp (Decomposition)
seed (int)
anchor_targets (ndarray | None)
perturb (float)
max_iter (int)
tol (float)

Return type:

ModelResult

OpenMM anchor-restrained relaxation engine (physics, MIT-licensed, optional).

The open replacement for FlexPepDock’s physics relaxation: build an OpenMM system for the complex, freeze the receptor (TCR + MHC) by zeroing masses, pin the anchor Cα to their target pockets with harmonic restraints (the physics analog of MODELLER’s forms.gaussian distance restraints), and run local energy minimisation. PDBFixer rebuilds the side chains that substitute_peptide() stripped, so the output is a full-atom relaxed peptide (an improvement over the backbone-only threading).

Optional dependency: conda install -c conda-forge openmm pdbfixer. Raises EngineUnavailable (never an ImportError at module import) when absent.

This engine is a reference/accuracy oracle for the native C++ rewrite (see CPP_REWRITE.md): its relaxed pose and relative energies validate the compact _relax minimiser that will replace it in the dependency-free path — OpenMM’s force field itself is not reimplemented.

class tcren.refine.engines.openmm_engine.OpenMMEngine[source]#

Bases: object

name = 'openmm'#

available()[source]#

Return type:: bool

run(structure, decomp, *, seed=0, anchor_targets=None, restraint_k=500.0, forcefield='amber14-all.xml', max_iter=0, tolerance=5.0)[source]#

Parameters:

structure (Structure)
decomp (Decomposition)
seed (int)
anchor_targets (ndarray | None)
restraint_k (float)
forcefield (str)
max_iter (int)
tolerance (float)

Return type:

ModelResult

ProMod3 side-chain reconstruction engine (Apache-2.0, optional).

ProMod3 is the open (Apache-2.0, license-key-free) replacement for MODELLER. Its cleanest, most directly useful capability for a groove peptide whose backbone is already placed (threaded or CCD-closed) is rotamer side-chain reconstruction (modelling.ReconstructSidechains): it rebuilds the side chains that substitute_peptide() strips, packed against the real receptor context with ProMod3’s rotamer library — the SCWRL/packer role. Full backbone loop building (promod3.loop CCD/KIC + fragments, the MODELLER-loopmodel analog) is the documented upgrade; repack is the piece that works via a stable API and complements the geometric (ccd) and physics (openmm) engines.

Optional dependency: conda install -c bioconda openstructure promod3. Raises EngineUnavailable (never an ImportError at import) when absent. Reference oracle for the native C++ rewrite (see CPP_REWRITE.md): validates a future native rotamer packer.

class tcren.refine.engines.promod3_engine.ProMod3Engine[source]#

Bases: object

name = 'promod3'#

available()[source]#

Return type:: bool

run(structure, decomp, *, seed=0, anchor_targets=None)[source]#

Parameters:

structure (Structure)
decomp (Decomposition)
seed (int)
anchor_targets (ndarray | None)

Return type:

ModelResult

Data paths#

tcren.paths module#

Filesystem locations for tcren’s reference data.

The library’s runtime dataset lives in the repo data/ directory (or $TCREN_DATA_DIR): the canonical Native2026 structure set (HF isalgo/tcren_structures, gitignored), PDB_date.tsv and orient_metadata.json. Structures are fetched lazily; nothing here is bundled into the installed package.

tcren.paths.data_dir()[source]#

Root of the runtime dataset: $TCREN_DATA_DIR or the repo data/ directory.

Return type:: Path

tcren.paths.native_dir()[source]#

Directory holding the canonical Native2026 structures (data/Native2026).

Return type:: Path

tcren.paths.reference_structure_path(pdb_id)[source]#

Resolve a canonical reference structure by id (plain/gzipped PDB/mmCIF).

Looks under data/Native2026 first; if absent (e.g. a pip-installed library with no repo data/), lazily downloads it from the HF dataset into the HF cache. This makes orienting a new, non-canonical structure work out of the box for both the library and CLI.

Raises FileNotFoundError if it is neither local nor fetchable.

Parameters:: pdb_id (str)
Return type:: Path

Analysis#

tcren.analysis module#

Dataset-level analyses for the TCRen contact statistics and potentials.

Helpers for the analysis notebook / benchmarks, following the TCRen manuscript logic: potential heatmaps and comparisons, the distribution of TCR↔peptide contacts per structure and per region, and how contacts distribute over peptide / CDR3 positions as a function of peptide / CDR3 length. They take the manuscript contact + summary tables as explicit paths (the committed oracle lives under tests/assets/oracle/).

tcren.analysis.load_interface_contacts(contact_maps, summary)[source]#

Load and enrich the manuscript TCR↔peptide contact table.

Adds, per contact: peptide_pos (0-based peptide position = residue.index.to), peptide_len, cdr3_len (CDR3α/β length for the contacting TCR chain), cdr3_rel_pos (residue index relative to the chain’s first contacting CDR3 residue — a relative position, since the committed table carries no region start), and the nonred flag.

Parameters:

contact_maps (str | Path)
summary (str | Path)

Return type:

DataFrame

tcren.analysis.contacts_per_structure(df, nonred_only=True)[source]#

Number of TCR↔peptide contacts per structure (with TRA/TRB split).

Parameters:

df (DataFrame)
nonred_only (bool)

Return type:

DataFrame

tcren.analysis.region_contact_counts(df, nonred_only=True)[source]#

Total contacts by TCR region (CDR1/2/3, FR) and chain (TRA/TRB).

Parameters:

df (DataFrame)
nonred_only (bool)

Return type:

DataFrame

tcren.analysis.position_distribution(df, side='peptide', nonred_only=True)[source]#

Contact counts by position, stratified by chain/molecule length.

Parameters:

df (DataFrame) – enriched contacts (see load_interface_contacts()).
side (str) – "peptide" (peptide position vs peptide length) or "cdr3a" / "cdr3b" (relative CDR3 position vs CDR3 length, for that TCR chain).
nonred_only (bool) – restrict to non-redundant structures.

Returns:

length, position, n_contacts.

Return type:

Long counts

tcren.analysis.potential_long(potential)[source]#

Heatmap-ready long form of a potential: residue.aa.from, residue.aa.to, value.

Parameters:: potential (Potential)
Return type:: DataFrame

tcren.analysis.compare_potentials(a, b)[source]#

Join two potentials on their amino-acid pairs and add their difference.

Returns residue.aa.from, residue.aa.to, value_a, value_b, diff over shared pairs.

Parameters:

a (Potential)
b (Potential)

Return type:

DataFrame

tcren.analysis.potential_matrix(potential)[source]#

Dense matrix + (from_labels, to_labels) for plotting a potential heatmap.

Parameters:: potential (Potential)
Return type:: tuple[ndarray, list[str], list[str]]

Orientation#

tcren.orient.align module#

Bring a structure into a canonical reference frame by MHC superposition.

A query complex is oriented onto a native reference by superposing the conserved MHC groove Cα atoms (the helix/floor residues from tcren.mhc.regions). Because every structure is aligned to the same reference, all oriented complexes share one frame — the basis for overlaying structures and for 2D interface projection. Correspondence between query and reference groove residues is established by sequence alignment, so different alleles/numbering are handled.

class tcren.orient.align.OrientationResult(rotation, translation, rmsd, n_anchor_atoms, reference_id)[source]#

Bases: object

Rigid transform that maps a structure onto the canonical reference frame.

Parameters:

rotation (ndarray)
translation (ndarray)
rmsd (float)
n_anchor_atoms (int)
reference_id (str)

rotation: ndarray#

translation: ndarray#

rmsd: float#

n_anchor_atoms: int#

reference_id: str#

tcren.orient.align.align_to_native(structure, reference_id=None)[source]#

Compute the transform orienting structure onto a native reference by MHC.

structure must already be chain-typed and MHC-annotated (see tcren.mhc.annotate_mhc()). The reference (default a canonical complex for the structure’s MHC class) is loaded from the Native2026 dataset (tcren.paths).

Parameters:

structure (Structure)
reference_id (str | None)

Return type:

OrientationResult

tcren.orient.align.apply_transform(structure, result)[source]#

Return a copy of structure with the orientation transform applied to all atoms.

Parameters:

structure (Structure)
result (OrientationResult)

Return type:

Structure

tcren.orient.superimpose module#

Superimpose query structures onto a canonical database by MHC.

Unlike tcren.orient.run_folder() (which builds a canonical set from native complexes using the per-class derived frame), superimpose() brings a new structure into the canonical frame defined by an existing database (data/Canonical2026 by default).

How it works, per query structure:

Chain-type + MHC-annotate the query; read its MHC class (MHCI/MHCII) and species.
Select every database structure with the same class and species (from the database’s orient_metadata.json).
Superpose the query’s conserved groove Cα onto each selected database structure (sequence alignment establishes the residue correspondence, so alleles/numbering differ freely).
Average the resulting rigid transforms — translations by mean, rotations by the chordal (SVD-orthonormalised) mean — into one consensus placement, and apply it.

Because every database member already sits in the same canonical frame, each superposition independently yields a canonical placement; averaging over the whole matching ensemble cancels the per-structure groove variation.

tcren.orient.superimpose.superimpose(structure, db_dir=None, organism='human', annotate=True)[source]#

Superimpose structure onto a canonical database by MHC (see module docstring).

structure is chain-typed + MHC-annotated here unless annotate=False (used by the threaded batch driver, which annotates the whole input set in one mmseqs pass first). The ensemble alignment itself is mmseqs-free, so it is the part safe to run on a thread pool. db_dir defaults to data/Canonical2026. Returns the oriented, A–E renamed structure and the consensus OrientationResult (averaged over the ensemble).

Parameters:

structure (Structure)
db_dir (str | Path | None)
organism (str)
annotate (bool)

Return type:

tuple[Structure, OrientationResult]

tcren.orient.frame module#

Canonical TCR-pMHC frame by PCA: z ≈ PC1 (MHC→TCR), y ≈ PC2 (peptide), x ≈ PC3.

Every structure is first superposed onto a per-class native reference by its MHC groove Cα (tcren.orient.align.align_to_native()); a fixed per-class rotation R_canon then maps that reference frame into the canonical axes. R_canon is obtained by centring the reference complex’s Cα cloud at its centre of mass and taking its principal axes (PCA):

z = PC1 (largest variance, the MHC→TCR long axis), signed +z toward the TCR so the MHC sits at −z;
y = PC2 (the groove/peptide axis), signed +y toward the peptide C-terminus;
x = PC3 (the thin axis), signed for a right-handed frame.

R_canon + the variance fractions are cached in the bundled tcren/data/canonical_frame.json so orientation is reproducible and inspectable. When no native database is available the same PCA axes are fit directly from the query (the PCA fallback).

class tcren.orient.frame.CanonResult(rotation, translation, rmsd, n_anchor_atoms, reference_id, frame, reversed_dock=None, chain_map=<factory>)[source]#

Bases: object

Composed rigid transform that maps a structure into the canonical frame.

Parameters:

rotation (ndarray)
translation (ndarray)
rmsd (float)
n_anchor_atoms (int)
reference_id (str | None)
frame (Literal['native', 'pca'])
reversed_dock (bool | None)
chain_map (dict[str, str])

rotation: ndarray#

translation: ndarray#

rmsd: float#

n_anchor_atoms: int#

reference_id: str | None#

frame: Literal['native', 'pca']#

reversed_dock: bool | None#

chain_map: dict[str, str]#

tcren.orient.frame.canonical_frame(structure, reference_id=None, force_pca=False)[source]#

Compose the MHC superposition with the per-class R_canon (native), or fit the canonical axes directly from the query (PCA fallback when no DB / too few anchors).

Parameters:

structure (Structure)
reference_id (str | None)
force_pca (bool)

Return type:

CanonResult

tcren.orient.frame.build_canonical_frame()[source]#

(Re)compute R_canon for each class reference and return the artifact dict.

Writes nothing; the caller serialises to tcren/data/canonical_frame.json.

Return type:: dict

tcren.orient.exceptions module#

Detect reverse-docked TCR-pMHC complexes (a biological exception, flagged not flipped).

The canonical frame is fixed by peptide polarity (+y = peptide C-terminus). The conserved diagonal docking then places the VDJ chain (TRB/TRD, Vβ) on the peptide-C side (+y) and the VJ chain (TRA/TRG, Vα) on the peptide-N side (−y) — consistent with the CDR footprint CDR1α·CDR2α·CDR3α·CDR3β·CDR2β·CDR1β laid out N→C. A genuinely reverse-docked TCR lands with the α/β sides mirrored. We report it; we never force-flip, because the orientation is meaningful.

tcren.orient.exceptions.detect_reverse_dock(structure, rotation, translation, margin=2.0)[source]#

Apply the canonical transform and check the TCR α/β handedness.

Canonical: VDJ (TRB/TRD) at +y (peptide-C side) and VJ (TRA/TRG) at −y. Returns True when the VJ chain is on the +y side of the VDJ chain by more than margin Å (reverse dock), False for a canonical dock, and None when a TCR side is missing.

Parameters:

structure (Structure)
rotation (ndarray)
translation (ndarray)
margin (float)

Return type:

bool | None

tcren.orient.chains module#

Select a single TCR-pMHC complex and rename its chains to the canonical A–E scheme.

tcren.orient.chains.select_primary_complex(structure)[source]#

Keep one mutually-contacting chain per canonical role (one TCR-pMHC complex).

Chosen as a connected unit so chains that do not touch the peptide (notably β2m) are not grabbed from another copy: the primary peptide is the one most embedded in an MHC-α groove (then most TCR contacts, then shortest); the TCR α/β and MHC-α are taken by contacts to that peptide; β2m / MHC-β by contacts to the chosen MHC-α. No-op for single complexes.

Parameters:: structure (Structure)
Return type:: Structure

tcren.orient.chains.rename_chains(structure)[source]#

Return a copy with only the canonical complex, chain ids remapped per CHAIN_RENAME, plus the old→new map.

Chains with no canonical role (tags, additives, unrelated proteins) are dropped so the output is exactly the A–E TCR-pMHC complex. Raises ValueError if two source chains map to the same canonical id (unresolved multi-copy — run select_primary_complex() first).

Parameters:: structure (Structure)
Return type:: tuple[Structure, dict[str, str]]

tcren.orient.pipeline module#

Orchestrate canonicalization of TCR-pMHC structures into the common MHC frame.

tcren.orient.pipeline.canonicalize_structure(structure, reference_id=None, force_pca=False, select_primary=True)[source]#

Orient an (already chain-typed + MHC-annotated) structure into the canonical frame.

Returns the oriented, A–E renamed structure and the populated CanonResult (transform, frame, rmsd, reverse-dock flag, chain map). Coordinates are transformed; the chain roles drive the rename, so order matters (frame + reverse-dock are read before the transform clears region markup).

Parameters:

structure (Structure)
reference_id (str | None)
force_pca (bool)
select_primary (bool)

Return type:

tuple[Structure, CanonResult]

tcren.orient.pipeline.align_to_canonical(structure, reference_id=None, organism='human', force_pca=False)[source]#

Align a NEW (parsed) structure onto the Native2026 canonical frame.

Runs chain typing + MHC annotation, then canonicalize_structure(). The stored per-class R_canon is reused, so the result is in the same frame as the dataset and the composed transform in the returned CanonResult replays the placement exactly.

Parameters:

structure (Structure)
reference_id (str | None)
organism (str)
force_pca (bool)

Return type:

tuple[Structure, CanonResult]

tcren.orient.pipeline.check_oriented_complex(structure, max_peptide_len=25, max_offset=25.0, max_tcr_gap=15.0, max_orphan=70.0)[source]#

Geometric sanity check on an oriented A–E complex; (ok, reason).

Rejects structures whose canonical placement is inconsistent: missing / overlong peptide, peptide not at the groove centre (≈ origin), the TCR not engaging the peptide, or any chain stranded far from the complex (an orphan copy that survived primary-complex selection).

Parameters:

max_peptide_len (int)
max_offset (float)
max_tcr_gap (float)
max_orphan (float)

tcren.orient.pipeline.run_folder(structures, out, metadata=None, organism='human', reference_id=None, force_pca=False, threads=None, mmcif=False, compress=False)[source]#

Canonicalize a file or folder of structures; write oriented structures + a metadata table.

Output format follows mmcif (.cif vs .pdb) and compress (trailing .gz); plain PDB by default (pass compress=True to rebuild the gzipped Canonical2026 set).

Annotation is BATCHED — one mmseqs search for all TCR chains (per organism) and one for all MHC chains across the whole set (mmseqs parallelises internally; never per-structure, never Python-threaded). Only the embarrassingly-parallel, mmseqs-free stages — parsing and the structural alignment + write — use a thread pool (threads worker threads, default os.cpu_count()).

Parameters:

structures (str | Path)
out (str | Path)
metadata (str | Path | None)
organism (str)
reference_id (str | None)
force_pca (bool)
threads (int | None)
mmcif (bool)
compress (bool)

Return type:

DataFrame

tcren.orient.pipeline.run_superimpose(structures, out, db_dir=None, organism='human', mmcif=False, compress=False, threads=None)[source]#

Superimpose input structure(s) onto a canonical database; write oriented structures.

structures is a file, directory, .tar.gz, or a shell glob. out is an output directory, or — for a single input — a structure file whose extension must match mmcif/compress. Annotation is BATCHED (one mmseqs pass over all inputs); only the mmseqs-free ensemble alignment + write runs on the thread pool (threads workers, default all cores). See tcren.orient.superimpose() for the MHC-ensemble method.

Parameters:

structures (str | Path)
out (str | Path)
db_dir (str | Path | None)
organism (str)
mmcif (bool)
compress (bool)
threads (int | None)

Return type:

DataFrame

tcren.orient.docking module#

TCR docking geometry: crossing angle and incident (tilt) angle.

The TCR “docking angle” (Rudolph, Stanfield & Wilson 2006; Garcia et al.) describes how the αβ (or γδ) TCR sits on top of the peptide-MHC groove. It is computed here directly from the canonical frame (tcren.orient.frame.canonical_frame()), so no external package (TCRdock / STCRpy) is required:

the crossing angle is the angle between the Vα→Vβ pseudo-axis projected into the MHC groove plane and the groove long axis (peptide N→C, canonical +y — collinear with the MHC α1 helix to within a few degrees). Reported on [0, 180); canonical αβ TCRs cluster around ~20–70°. A signed variant ([-180, 180)) carries the handedness of the docking.
the incident (tilt) angle is the elevation of the same Vα→Vβ vector out of the groove plane (canonical z is the MHC→TCR normal): positive when Vβ rides higher above the groove than Vα.

The Vα/Vβ landmarks are the centroids of the variable-domain Cα atoms of the two receptor chains (TRA/TRB for αβ, TRG/TRD for γδ). The frame is fit from the query itself (PCA), so the calculation needs neither the native database nor mmseqs once the structure is chain-typed.

class tcren.orient.docking.DockingAngles(crossing_angle, crossing_angle_signed, incident_angle, cell_type, n_va, n_vb)[source]#

Bases: object

TCR docking geometry relative to the MHC groove (all angles in degrees).

Parameters:

crossing_angle (float)
crossing_angle_signed (float)
incident_angle (float)
cell_type (str)
n_va (int)
n_vb (int)

crossing_angle: float#

crossing_angle_signed: float#

incident_angle: float#

cell_type: str#

n_va: int#

n_vb: int#

tcren.orient.docking.crossing_incident_from_vector(v_canon)[source]#

(crossing, crossing_signed, incident) degrees from a Vα→Vβ vector in canonical axes.

v_canon is [vx, vy, vz] along canonical x (groove width), y (groove long axis, peptide N→C) and z (MHC→TCR normal). The crossing angle is measured in the groove plane (xy) from the long axis; the incident angle is the elevation out of that plane.

Parameters:: v_canon (ndarray)
Return type:: tuple[float, float, float]

tcren.orient.docking.docking_angles(structure)[source]#

Crossing + incident angle of a chain-typed TCR-pMHC complex.

The structure must already be chain-typed (classify_chains) and MHC-annotated (annotate_mhc) so the canonical frame can be fit. The frame is taken from the query’s own Cα cloud (PCA), so the result needs no native database.

Parameters:: structure (Structure) – a chain-typed, MHC-annotated TCR-pMHC structure.
Returns:: A DockingAngles with the crossing and incident angles.
Raises:: ValueError – if a receptor chain pair (TRA/TRB or TRG/TRD) is missing, or the canonical frame is degenerate.
Return type:: DockingAngles

tcren.orient.tcrdock_geometry module#

TCR:pMHC docking geometry — native reimplementation of TCRdock’s rigid-body parameterisation.

Reimplemented from the Bradley lab’s TCRdock (phbradley/TCRdock, MIT license), commit c5a7af42eeb0c2a4492a4d4fe803f1f9aafb6193 (2024-03-04), specifically tcrdock/docking_geometry.py, tcrdock/mhc_util.py, tcrdock/tcr_util.py and tcrdock/superimpose.py. No TCRdock code is imported; the geometry is ported to tcren’s own Structure and annotation.

The docking geometry describes how the TCR sits on the peptide-MHC groove as a rigid-body transform between two coordinate frames (“stubs”):

the MHC stub — from the ~180° pseudo-symmetry of the class-I α1α2 (or class-II α1β1) β-sheet floor: x points toward the peptide, z from one half of the sheet to the other, origin at the two halves’ midpoint;
the TCR stub — from the ~180° pseudo-symmetry relating Vα and Vβ: x toward the CDR loops, z from Vα to Vβ, origin at the two domains’ midpoint.

Six numbers (DockingGeometry) fix the relative pose: d (frame separation), torsion (dihedral about the MHC–TCR line), and tcr_unit_y/z + mhc_unit_y/z (each frame’s direction to the other, in the other’s local axes). For interpretable in-plane / tilt angles use tcren.orient.docking.docking_angles() (crossing_angle = the groove-plane “scanning” angle; incident_angle = the tilt); this module adds the full rigid-body pose that those two scalars do not capture.

MHC-I core positions are mapped by BLOSUM-aligning the α chain to TCRdock’s class-I template; TCR core positions are the conserved IMGT framework positions, located from tcren’s arda region markup.

Provenance note (validated on 618 TCRvdb TCRmodel2 models, 2026-07-05): tcren’s crossing_angle reproduces the “scanning_angle” reported by upstream AF/TCRmodel2 annotation tables (r≈0.88), so that quantity is a genuine, reproducible interface geometry. The upstream “pitch_angle”, however, matches no clean geometric angle (best correlate is d, r≈0.42) and discriminates TCRvdb binders better (macro-PR≈0.72) than any clean docking feature computed here (d≈0.64, torsion≈0.62, tilt≈0.58) — i.e. its extra signal is not reproducible from coordinates and is likely AlphaFold-confidence contamination. Prefer these documented, crystal-computable descriptors over the opaque upstream pitch.

class tcren.orient.tcrdock_geometry.DockingGeometry(d, torsion, tcr_unit_y, tcr_unit_z, mhc_unit_y, mhc_unit_z)[source]#

Bases: object

TCR:pMHC rigid-body docking geometry — TCRdock’s 6-parameter form.

Parameters:

d (float)
torsion (float)
tcr_unit_y (float)
tcr_unit_z (float)
mhc_unit_y (float)
mhc_unit_z (float)

d: float#

torsion: float#

tcr_unit_y: float#

tcr_unit_z: float#

mhc_unit_y: float#

mhc_unit_z: float#

to_dict()[source]#

Return type:: dict

tcren.orient.tcrdock_geometry.docking_geometry(structure)[source]#

Compute the TCRdock docking geometry of a chain-typed, MHC-annotated TCR-pMHC structure.

The structure must already be chain-typed (tcren.annotation.classify_chains()) and MHC-annotated (tcren.mhc.annotate_mhc()) with arda region markup on the TCR chains.

Parameters:: structure – a chain-typed, annotated TCR-pMHC Structure.
Returns:: The DockingGeometry.
Raises:: ValueError – if the MHC β-sheet core or a complete TCR Vα/Vβ core cannot be located.
Return type:: DockingGeometry

tcren.orient.graft module#

Build chimeric TCR:pMHC complexes by grafting one complex’s TCR onto another’s pMHC.

Given a host and a donor TCR:pMHC complex, substitute_tcr() produces a new complex that keeps the host peptide + MHC and the donor TCR, rigidly positioned by one of two anchors:

by="mhc" — superpose the donor MHC-groove Cα onto the host MHC-groove Cα, then carry the donor TCR into the host frame. The donor TCR keeps its native docking geometry relative to MHC (the pose it adopts on its own pMHC, transferred onto the host groove).
by="tcr" — superpose the donor TCR Cα onto the host TCR Cα, then drop the host TCR. The donor TCR inherits the host’s docking pose (it lands where the host TCR sat).

Both yield host pMHC + donor TCR; only the superposition anchor differs. Residue correspondence between the two MHCs (or the two TCRs) is by sequence alignment, so different alleles / V-genes are handled. The MHC path needs both inputs MHC-annotated (tcren.mhc.annotate_mhc()); both paths need chain typing (tcren.annotation.classify_chains()).

tcren.orient.graft.substitute_tcr(host, donor, *, by='mhc')[source]#

Graft the donor TCR onto the host pMHC → a chimeric TCR:pMHC Structure.

The result keeps the host peptide + MHC chains and the donor TCR chains (relabelled to avoid id collisions), with the donor TCR rigidly placed by the by anchor:

"mhc" — superpose the donor MHC groove onto the host MHC groove (both inputs must be tcren.mhc.annotate_mhc()-annotated); the donor TCR keeps its native docking geometry.
"tcr" — superpose the donor TCR onto the host TCR; the donor TCR inherits the host pose.

Both inputs must be chain-typed (tcren.annotation.classify_chains()). Raises ValueError if by is invalid, the host lacks a peptide or MHC chain, the donor lacks a TCR chain, or too few matched Cα anchors are found to superpose.

Parameters:

host (Structure)
donor (Structure)
by (str)

Return type:

Structure

2D projection & visualization#

tcren.project2d.frame module#

Project an oriented TCR-peptide-MHC interface onto a single 2D plane.

The “optimal plane” is the MHC groove plane. Two routes produce it:

native (default) — orient the structure onto a canonical reference via tcren.orient.align.align_to_native(); the canonical groove plane is xy and its normal is z, so projecting is just dropping z. The transform is applied only to the extracted Cα coordinates (apply_transform would clear region annotations).
pca (fallback) — fit a plane to the groove-floor Cα by SVD when no native database is available. The normal is the lowest-variance axis; z is sign-oriented toward the peptide.

The residues projected are the CDR1-3 loops (TRA/TRB), the peptide, and the MHC groove helices/floor — the residues that line the interface.

class tcren.project2d.frame.ProjectionResult(keys, coords3d, frame, reference_id=None, rmsd=None)[source]#

Bases: object

2D projection of selected interface residues onto the groove plane.

Parameters:

keys (list[tuple[str, int]])
coords3d (ndarray)
frame (Literal['native', 'pca'])
reference_id (str | None)
rmsd (float | None)

keys: list[tuple[str, int]]#

coords3d: ndarray#

frame: Literal['native', 'pca']#

reference_id: str | None#

rmsd: float | None#

property uv: ndarray#: In-plane (u, v) coordinates.

property height: ndarray#: Signed distance above the groove plane (z).

tcren.project2d.frame.project_structure(structure, reference_id=None, force_pca=False)[source]#

Project the interface residues of an annotated structure onto the groove plane.

Parameters:

structure (Structure) – A chain-typed, MHC-annotated structure.
reference_id (str | None) – Canonical reference complex id (defaults per MHC class).
force_pca (bool) – Skip the native frame and fit the plane by PCA.

Returns:

A ProjectionResult.

Return type:

ProjectionResult

tcren.project2d.tables module#

Canonical polars tables for complementarity maps: residue markup + contacts.

Both tables are tidy and keyed for joins on (structure_id, structure_chain, aa_index) where aa_index is the 0-based position in the chain’s protein sequence (the legacy mir residue.index); residue_index is the author/PDB numbering (a display label). The contacts table wraps the parity-preserving all_atom_contacts(), so it is compatible with the original TCRen contact calculation.

tcren.project2d.tables.residue_markup_table(structure, projection=None)[source]#

Per-residue markup with chain/region annotation and groove-frame coordinates.

Parameters:

structure (Structure) – A chain-typed, annotated structure.
projection (ProjectionResult | None) – Optional ProjectionResult; its groove-frame coordinates are joined on (chain_id, seq_index) (x/y/z and in-plane u/v/height).

Returns:

structure_id, structure_chain, complex_chain, complex_region, residue_index, aa_index, aa_len, aa, x, y, z, u, v, height.

Return type:

Columns

tcren.project2d.tables.ca_contacts_table(structure, threshold=8.0)[source]#

Inter-chain Cα–Cα “chain contacts” within threshold Å (bold edges in the map).

Columns: structure_id, structure_chain_1, structure_chain_2, aa_index_1, aa_index_2, ca_dist. Complements the all-atom contacts_table() (which gives the dashed inter-residue edges).

Parameters:

structure (Structure)
threshold (float)

Return type:

DataFrame

tcren.project2d.tables.classify_contact(aa1, aa2, atom1, atom2, dist)[source]#

Heuristically classify a residue–residue contact from its closest atom pair.

Returns one of salt_bridge, hydrogen_bond, aromatic, hydrophobic, polar, other. Cutoffs are pragmatic, not a force field — documented and kept in this pure function for easy tuning.

Parameters:

aa1 (str)
aa2 (str)
atom1 (str)
atom2 (str)
dist (float)

Return type:

str

tcren.project2d.tables.contacts_table(structure, threshold=5.0)[source]#

Inter-chain residue contacts within threshold Å, classified by bond type.

Parameters:

structure (Structure) – The structure.
threshold (float) – Distance cutoff in Å (3 ≤ threshold ≤ 6, default 5 — the original TCRen contact distance).

Returns:

structure_id, structure_chain_1, structure_chain_2, residue_index_1, residue_index_2, aa_index_1, aa_index_2, min_dist, contact_type, backbone_1, backbone_2. One row per unordered inter-chain residue pair (TCRen-compatible).

Return type:

Columns

tcren.project2d.tables.region_pair_contacts(structure, kind='closest', cutoff=None)[source]#

Inter-chain residue contacts annotated with the region pair they bridge.

kind selects the contact definition: "closest" (closest heavy-atom pair, all_atom_contacts(), default 5 Å — the original TCRen definition, and the only kind that carries a contact_type bond classification), "cb" (Cβ representative atom, default 8 Å), or "ca" (Cα representative atom, default 12 Å). Region labels are ordered canonically per row so a pair is direction-independent (cdr3``↔``peptide == peptide``↔``cdr3).

Returns columns: structure_id, complex_chain_1, region_1, complex_chain_2, region_2, aa_index_1, aa_index_2, min_dist plus, for kind="closest", contact_type.

Parameters:

structure (Structure)
kind (str)
cutoff (float | None)

Return type:

DataFrame

tcren.project2d.tables.region_pair_summary(structure, kind='closest', cutoff=None)[source]#

Per region-pair contact counts (+ bond-type breakdown for kind="closest").

Aggregates region_pair_contacts() to one row per (region_1, region_2) pair with n_contacts. For kind="closest" it adds a column per bond type (n_hydrogen_bond, n_salt_bridge, n_aromatic, n_hydrophobic, n_polar, n_other), so hydrogen bonds (and the rest) are reported for every region pair in the complex, not just the MHC interface.

Parameters:

structure (Structure)
kind (str)
cutoff (float | None)

Return type:

DataFrame

tcren.project2d.pockets module#

Approximate MHC peptide-binding pocket (A–F) markers along the peptide track.

The six class-I pockets (A–F) line the groove from the peptide N-terminus (A ≈ P1) to its C-terminus (F ≈ P-Ω). Without explicit pocket geometry we place A–F evenly along the projected peptide Cα track — a labelling aid for the 2D map and 3D view, not a structural pocket definition.

tcren.project2d.pockets.pocket_markers(markup)[source]#

Place A–F markers along the projected peptide (u, v) track.

Returns columns pocket, peptide_pos, u, v (empty if the peptide is not projected).

Parameters:: markup (DataFrame)
Return type:: DataFrame

tcren.viz.svg2d module#

Direct SVG builder for 2D complementarity maps.

Renders projected interface residues as squares (AA + number), Cα–Cα “chain contacts” as bold lines, and closest-atom inter-residue contacts as dashed lines. Every element carries its data as data-* attributes plus a <title> tooltip, so the SVG is both a figure and a queryable, metadata-bearing artifact. Pure string building — no dependencies.

tcren.viz.svg2d.render_complementarity_map(markup, contacts=None, ca_contacts=None, pockets=None, show_chains=None, draw_backbone=True, width=900, height=700, margin=60.0)[source]#

Render a complementarity map to an SVG string.

Parameters:

markup (DataFrame) – residue markup table (needs u, v; non-null rows are drawn).
contacts (DataFrame | None) – closest-atom inter-residue contacts (dashed) with structure_chain_1/2, aa_index_1/2, min_dist, contact_type.
ca_contacts (DataFrame | None) – Cα–Cα chain contacts (bold) with structure_chain_1/2, aa_index_1/2, ca_dist.
pockets (DataFrame | None) – optional A–F pocket markers with pocket, u, v.
show_chains (list[str] | None) – if given, only draw residues whose complex_chain is in this set (e.g. ["tra", "trb", "peptide"] to hide MHC). Contacts/edges are kept only between shown residues.
draw_backbone (bool) – connect consecutive residues within each chain (sequence-adjacent aa_index) with a thin backbone trace.
width (int) – canvas geometry.
height (int) – canvas geometry.
margin (float) – canvas geometry.

Returns:

SVG markup (string).

Return type:

str

tcren.viz.pocket3d module#

Interactive 3D peptide-pocket view with CDR1-3 overlay (py3Dmol / 3Dmol.js).

Renders the MHC groove (cartoon + optional translucent surface, histo.fyi style), the peptide as sticks, and the TCR CDR1-3 loops as Cα traces coloured by the shared palette, all in the canonical (oriented) frame. py3Dmol is imported lazily so the rest of the package has no hard 3D dependency.

tcren.viz.pocket3d.view_pocket_cdr(structure, reference_id=None, surface=True, width=700, height=500)[source]#

Build a py3Dmol view of the groove + peptide + CDR loops (oriented frame).

Parameters:

structure (Structure) – a chain-typed, MHC-annotated structure.
reference_id (str | None) – native reference for orientation (identity frame if unavailable).
surface (bool) – draw a translucent groove surface.
width (int) – viewer size.
height (int) – viewer size.

Returns:

A py3Dmol.view (call .show() in a notebook).

tcren.viz.pymol module#

Publication-ready PyMOL renders of canonically-oriented TCR–pMHC complexes.

A render of an oriented structure is only interpretable if the reader can tell which way the canonical frame points, so the centrepiece here is the axis gizmo — a thin, arrow-headed triad in a corner of the image, turning with the camera, naming what each direction means.

The frame it draws. tcren.orient.frame puts every structure into one frame by PCA of the reference complex, and CANONICAL_AXES is that frame written out in words:

axis	definition in code	figure label	equivalent in the literature
`x`	PC3, the thin axis	`groove width`	`PC2_MHC`, groove width
`y`	PC2, +y to peptide C-term	`peptide N→C`	`PC1_MHC`, groove long axis
`z`	PC1, +z toward the TCR	`pMHC→TCR`	`PC3_MHC`, groove normal

The principal-component numbers differ from the docking-geometry literature (SwiftTCR, TCR3d) because those fit the MHC groove alone while tcren.orient.frame fits the whole complex, in which the MHC→TCR direction carries the most variance. The three directions are the same three; only their ranking differs. Naming them for what they are is the point — pMHC→TCR is readable in a figure, z is not.

How the gizmo is placed. In its own render pass, not by projecting a world position into the corner. PyMOL’s orthoscopic viewport does not span the world height that field_of_view and the camera distance imply — measured on a real scene it is out by about a quarter — so a gizmo placed by that arithmetic lands off-frame. Rendering the triad alone on a transparent background and compositing it at pixel coordinates makes its size and position exact by construction, and leaves nothing to reverse-engineer. It also guarantees the molecule can never occlude it.

Why the geometry lives here. PyMOL runs under its own interpreter and cannot import tcren, so the scripts this module emits are lists of literals and every decision behind them is ordinary testable Python.

Example

>>> from tcren.viz.pymol import render, overlay_scene
>>> render(overlay_scene(["1ao7"], "data/Canonical2026"), "overlay.png")

class tcren.viz.pymol.Axis(letter, label, short, definition)[source]#

Bases: object

One canonical axis: its letter, what it means, and how to label it in a figure.

Parameters:

letter (str)
label (str)
short (str)
definition (str)

letter: str#

label: str#

short: str#

definition: str#

tcren.viz.pymol.CANONICAL_AXES: tuple[Axis, Axis, Axis] = (Axis(letter='x', label='groove width', short='width', definition='PC3, the thin axis; across the cleft, α1↔α2 helix separation'), Axis(letter='y', label='peptide N→C', short='N→C', definition='PC2, the groove/peptide axis, signed toward the peptide C-terminus'), Axis(letter='z', label='pMHC→TCR', short='TCR', definition='PC1, the MHC→TCR long axis, signed toward the TCR; the MHC sits at −z'))#: The canonical frame of tcren.orient.frame, in figure-ready words. short is for the corner gizmo, where anything longer than a few characters collides with the arrows.

tcren.viz.pymol.CORNERS: dict[str, tuple[int, int]] = {'bottom-left': (-1, -1), 'bottom-right': (1, -1), 'top-left': (-1, 1), 'top-right': (1, 1)}#: Corner anchors, as the sign of (x, y) in image space with y up.

tcren.viz.pymol.PALETTES: dict[str, tuple[tuple[float, float, float], ...]] = {'mono': ((0.2, 0.2, 0.2), (0.2, 0.2, 0.2), (0.2, 0.2, 0.2)), 'okabe-ito': ((0.835, 0.369, 0.0), (0.0, 0.62, 0.451), (0.0, 0.447, 0.698))}#: mono — one restrained grey for all three arrows, told apart by their labels. The default, because the structure already spends colour on chains and an orientation gizmo that competes with the molecule makes a worse figure. okabe-ito is the colourblind-safe triad for when the axes themselves are the subject.

tcren.viz.pymol.CHAIN_COLOURS: dict[str, tuple[str, str]] = {'A': ('Vα', 'marine'), 'B': ('Vβ', 'orange'), 'C': ('peptide', 'yellow'), 'D': ('MHC α', 'grey70'), 'E': ('MHC β / β2m', 'grey60')}#: Chain roles after tcren orient, and the colours used throughout these scenes.

tcren.viz.pymol.gizmo_cgo(*, arm=10.0, radius=0.3, head_length=0.3, head_radius=2.8, palette='mono')[source]#

The axis triad as CGO, at the world origin, in world units.

Absolute size does not matter: the triad is rendered on its own and scaled to the gizmo tile, so only the proportions here are visible. What they control is how the figure reads — a thin shaft with a distinct arrowhead, rather than the fat default axes that dominate a panel.

Parameters:

arm (float) – Arrow length.
radius (float) – Shaft radius. Thin relative to arm is the intent.
head_length (float) – Arrowhead length as a fraction of arm.
head_radius (float) – Arrowhead base radius as a multiple of radius.
palette (str) – A key of PALETTES.

Returns:

(cgo, tips) — the flat CGO float list, and the three arrow tips, where labels go.

Raises:

ValueError – If palette is unknown, or a size is non-positive.

Return type:

tuple[list[float], list[tuple[float, float, float]]]

Example

>>> cgo, tips = gizmo_cgo()
>>> tips[1]
(0.0, 10.0, 0.0)

tcren.viz.pymol.label_points(rotation, tips, *, offset=1.42, head_on=0.25)[source]#

Where each axis label goes, pushed clear of its arrow in screen space.

Anchoring a label to its tip in world coordinates fails exactly where it matters: an axis pointing at the viewer foreshortens to a dot, and its label lands on top of the origin and the other two labels. Pushing outward along the axis’s projected direction keeps every label clear whatever the camera does, and an axis too head-on to have a projected direction is pushed down-left instead, where it reads as belonging to the dot at the origin.

Parameters:

rotation – The nine floats of the camera rotation, column-major as PyMOL stores them.
tips – The three arrow tips from gizmo_cgo().
offset (float) – How far past the tip to sit, as a multiple of the arm length.
head_on (float) – Projected length below which an axis counts as pointing at the viewer.

Returns:

Three world-space points, rounded for embedding in a script.

Example

>>> pts = label_points([1,0,0, 0,1,0, 0,0,1], [(10,0,0), (0,10,0), (0,0,10)])
>>> pts[0][0] > 10          # pushed out past the tip
True

tcren.viz.pymol.gizmo_scene(rotation, *, axes=(Axis(letter='x', label='groove width', short='width', definition='PC3, the thin axis; across the cleft, α1↔α2 helix separation'), Axis(letter='y', label='peptide N→C', short='N→C', definition='PC2, the groove/peptide axis, signed toward the peptide C-terminus'), Axis(letter='z', label='pMHC→TCR', short='TCR', definition='PC1, the MHC→TCR long axis, signed toward the TCR; the MHC sits at −z')), short_labels=True, label_size=26.0, label_colour='gray20', label_offset=1.42, **cgo_kwargs)[source]#

A PyMOL scene holding only the triad, seen under rotation.

Parameters:

rotation (list[float]) – The first nine floats of the main scene’s cmd.get_view() — the same camera rotation, so the triad reports the orientation the molecule is actually drawn in.
axes (tuple[Axis, Axis, Axis]) – The axes to draw, in x, y, z order.
short_labels (bool) – Use Axis.short (fits a small tile) rather than Axis.label.
label_size (float) – PyMOL label_size. Large, because the tile is later scaled down.
label_colour (str) – Any PyMOL colour name.
label_offset (float) – Where the label sits along the arm, as a multiple of its length.
**cgo_kwargs – Forwarded to gizmo_cgo().

Returns:

A PyMOL script body. Labels are pseudoatoms because CGO has no text primitive; they ray-trace like any other label.

Raises:

ValueError – If rotation is not nine floats.

Return type:

str

tcren.viz.pymol.probe_rotation(scene, *, pymol_bin=None)[source]#

Run scene and report the camera rotation it ends up with.

The gizmo must be drawn under the same rotation as the molecule, and that is not known until the scene has loaded, turned and zoomed — cmd.zoom leaves the rotation alone but cmd.turn and cmd.orient do not, so reading it back beats assuming it.

Parameters:

scene (str) – A PyMOL script body that loads structures and sets the view.
pymol_bin (str | None) – Override the pymol executable.

Returns:

The first nine floats of cmd.get_view().

Return type:

list[float]

tcren.viz.pymol.render(scene, png, *, size=(1200, 1200), dpi=300, gizmo=True, corner='bottom-left', gizmo_scale=0.19, gizmo_margin=0.025, pymol_bin=None, **gizmo_kwargs)[source]#

Ray-trace scene to png, with the canonical-frame gizmo in a corner.

Parameters:

scene (str) – A PyMOL script body: load structures, style them, set and zoom the view.
png (str | Path) – Output path.
size (tuple[int, int]) – (width, height) in pixels. 1200 px at 300 dpi is a 4-inch figure panel.
dpi (int) – Written into the PNG so a document places it at a known physical size.
gizmo (bool) – Draw the axis triad. Turn it off where the frame is not the point.
corner (str) – Which corner the gizmo goes in; a key of CORNERS.
gizmo_scale (float) – Gizmo tile width as a fraction of the image width.
gizmo_margin (float) – Gizmo inset from the edges, as a fraction of the image width.
pymol_bin (str | None) – Override the pymol executable.
**gizmo_kwargs – Forwarded to gizmo_scene() / gizmo_cgo().

Returns:

The path written.

Return type:

Path

Example

>>> render('cmd.load("1ao7.pdb")\ncmd.show("cartoon")', "fig.png")

tcren.viz.pymol.composite(base_png, tile_png, out_png, *, corner='bottom-left', scale=0.19, margin=0.025)[source]#

Paste the gizmo tile into a corner of the render, preserving alpha.

Parameters:

base_png – The molecule render.
tile_png – The gizmo render, transparent outside the arrows.
out_png – Where to write. May be base_png.
corner (str) – A key of CORNERS.
scale (float) – Tile width as a fraction of the base image width.
margin (float) – Inset from the edges, as a fraction of the base image width.

Returns:

The path written.

Raises:

ValueError – If corner is unknown or scale/margin leave no room.

tcren.viz.pymol.overlay_scene(ids, canon_dir, *, limit=8, transparency=0.55)[source]#

Superpose a set of oriented structures, seen side-on.

Parameters:

ids – PDB ids present in canon_dir as <id>.pdb.gz.
canon_dir – The canonical (oriented) structure directory.
limit (int) – Draw at most this many; past roughly eight the overlay stops being readable.
transparency (float) – Cartoon transparency, so the spread of the ensemble shows through.

Returns:

A PyMOL scene body for render().

Return type:

str

tcren.viz.pymol.groove_scene(pid, canon_dir, *, surface=False)[source]#

One complex from above the groove: peptide as sticks in the MHC cleft.

The layout histo.fyi uses for its structure pages — a pale MHC with the peptide threaded along the cleft, which is the most legible way to show what is presented.

Parameters:

pid – PDB id.
canon_dir – The canonical (oriented) structure directory.
surface (bool) – Add a translucent molecular surface over the MHC ribbon.

Returns:

A PyMOL scene body for render().

Return type:

str

tcren.viz.pymol.interface_scene(pid, canon_dir, cdr_resi)[source]#

The recognition interface: peptide plus the CDR loops that touch it.

Parameters:

pid – PDB id.
canon_dir – The canonical (oriented) structure directory.
cdr_resi – {"TRA": [...], "TRB": [...]} PDB residue numbers, as from annotating the pre-orientation structure — orientation preserves numbering.

Returns:

A PyMOL scene body for render().

Return type:

str

tcren.viz.pymol.residue_importance(structure, *, interface='tcr_peptide', cutoff=5.0, potential=None, tcr_regions='all')[source]#

Per-residue share of the interface, on both the physics and the geometry axis.

The interface energy Φ is a sum over residue–residue contacts, so it decomposes exactly: each residue’s share is the sum of φ(aa_i, aa_j) over the contacts it makes across the interface. That is the quantity to colour a figure by — it says which residues the score is actually made of, where the total only says how large it is.

Two columns come back because they answer different questions. phi is the energy share, negative being favourable; n_contacts is the geometric share, how much of the interface the residue physically occupies. A residue can be large on one and small on the other, and that difference is usually the interesting part.

Parameters:

structure – An annotated structure (chains typed, MHC called).
interface (str) – Which interface to decompose — tcr_peptide, tcr_mhc or peptide_mhc.
cutoff (float) – Heavy-atom contact distance (Å).
potential – A tcren.potential.Potential; defaults to the bundled TCRen for tcr_peptide and MJ for the MHC interfaces, matching tcren.pipeline.
tcr_regions (str) – Which TCR regions count — all, cdr or cdr+fr.

Returns:

chain.id, residue.index, residue.aa, region.type, n_contacts, phi. One row per residue that touches the interface, either side of it, sorted most-favourable first.

Return type:

A polars frame

Example

>>> residue_importance(s).head(3)

tcren.viz.pymol.importance_scene(pid, canon_dir, importance, *, by='phi', regions=('CDR3', 'PEPTIDE'), spectrum='blue_white_red')[source]#

Colour the recognition interface by each residue’s share of it.

The residues that carry the score are shown as sticks on a colour ramp; everything else stays pale so the eye goes to the ramp. Values ride in on the B-factor column, which is what PyMOL’s spectrum reads.

The ramp is centred on zero when colouring by phi, so blue and red mean favourable and unfavourable rather than merely “less” and “more” — a ramp fitted to the observed range would paint the least-favourable residue red even in an interface where every contact is stabilising.

Parameters:

pid – PDB id.
canon_dir – The canonical (oriented) structure directory.
importance – The frame from residue_importance().
by (str) – Which column to colour by — phi (energy share) or n_contacts (geometric share).
regions – Which region types to draw as coloured sticks. PEPTIDE matches the peptide chain, whose residues carry no CDR region label.
spectrum (str) – Any PyMOL spectrum name.

Returns:

A PyMOL scene body for render().

Raises:

ValueError – If by is not a column of importance.

Return type:

str

tcren.viz.palette module#

Colorblind-safe palette for complementarity maps (Okabe-Ito), shared by 2D and 3D.

Hue encodes the complex chain (tra/trb/peptide/mhca/mhcb); CDR loops get distinct shades.

tcren.viz.palette.color_for(complex_chain, complex_region)[source]#

Resolve a fill color: region shade if known, else chain hue, else grey.

Parameters:

complex_chain (str | None)
complex_region (str | None)

Return type:

str

Reference data & reproduction#

tcren.recent module#

Fetch recent TCR-pMHC structures from the RCSB PDB into data/pdb_recent (gitignored).

Two entry points, both gzipping validated mmCIF into the destination:

fetch_ids() — download specific PDB ids (e.g. the Native2026 set) straight from RCSB.
discover_similar() — RCSB full-text search for TCR:peptide:MHC entries (optionally released after a date / excluding ids we already have), to surface new structures.

Robustness notes baked in per the PDB’s current state:

IDs may be longer than 4 characters (extended pdb_0000XXXX accessions) — handled.
The PDB is deprecating split ``.pdb`` files for large structures, so we always pull mmCIF (.cif.gz), which tcren reads natively.

Every kept structure is annotated (batched, one mmseqs pass) and must have all 5 required chains — MHCα, β2m or MHCβ, peptide, and the two TCR chains (TRA/TRB or TRG/TRD) — else it is dropped. huggingface_hub/network are not involved; this uses requests + the RCSB APIs.

tcren.recent.recent_dir()[source]#

The gitignored destination for fetched structures (data/pdb_recent).

Return type:: Path

tcren.recent.discover_similar(after_date=None, limit=200, timeout=30.0)[source]#

RCSB full-text search for TCR:peptide:MHC structures; return candidate PDB ids.

after_date (YYYY-MM-DD) restricts to entries released on/after it (find new structures). This is the keyword-driven discovery step; the returned ids are downloaded and then strictly validated by _has_required_chains() (an agent can further curate the keyword set or the returned list before a fetch).

Parameters:

after_date (str | None)
limit (int)
timeout (float)

Return type:

list[str]

tcren.recent.fetch_ids(ids, dest=None, organism='human')[source]#

Download ids from RCSB into dest as .cif.gz, keep only complete complexes.

Incomplete (missing one of the 5 required chains) or unparseable downloads are removed. Returns a summary dict {requested, downloaded, complete, kept}.

Parameters:

ids (list[str])
dest (Path | None)
organism (str)

Return type:

dict

tcren.recent.native2026_ids()[source]#

PDB ids of the local Native2026 set (the seed for a pdb_recent refresh).

Return type:: list[str]

tcren.paper.bootstrap module#

Bootstrap data for the Nat Comput Sci reproduction notebooks.

Fetches the structure sets from the Hugging Face dataset isalgo/tcren_structures and the pinned vdjdb release into the shared notebooks/data/ dir, and stages the paper’s non-structure data there plus the legacy (mir/R) text outputs used as the regression oracle under notebooks/natcompsci2022/results_legacy/. Network access uses requests only (no huggingface_hub dependency). All copied csv/tsv/txt are gzipped; the downloaded HF structures are gitignored.

tcren.paper.bootstrap.fetch_hf_structures(data_dir, folders=('Native2022', 'Native2026', 'PolyV2022', 'Bobisse', 'Bigot'), force=False, timeout=120.0)[source]#

Download the HF structure folders into <data_dir>/structures/ (gitignored).

Uses huggingface_hub.snapshot_download when available (robust LFS download with resume + retries); otherwise falls back to a retrying requests loop. Neither is a hard package dependency — this is an optional reproduction tool.

Parameters:

data_dir (Path)
folders (tuple[str, ...])
force (bool)
timeout (float)

Return type:

dict[str, int]

tcren.paper.bootstrap.fetch_vdjdb(data_dir, date='2022-03-30', timeout=300.0)[source]#

Download the pinned vdjdb release and gzip its slim table into <data_dir>/vdjdb/.

Parameters:

data_dir (Path)
date (str)
timeout (float)

Return type:

Path

tcren.paper.bootstrap.fetch_pdb_dates(data_dir, pdb_ids, timeout=30.0)[source]#

Fetch PDB initial-release dates from RCSB into <data_dir>/PDB_date.csv.gz.

RCSB release dates are external published metadata (needed for the holdout date-split). Queried per entry via the RCSB Data API; cached, so re-runs only fetch missing ids.

Parameters:

data_dir (Path)
pdb_ids (list[str])
timeout (float)

Return type:

Path

tcren.paper.bootstrap.copy_external_inputs(data_dir, repo_data=None)[source]#

Stage the allowed externally published inputs (gzipped) into <data_dir>/.

These are the only non-structure inputs the reproduced pipeline may consume: the MJ/Keskin potentials, the Birnbaum set, IEDB, and the Bobisse/Bigot candidate lists.

Parameters:

data_dir (Path)
repo_data (Path | None)

Return type:

int

tcren.paper.bootstrap.copy_legacy_results(paper_dir=PosixPath('/home/runner/work/tcren/tcren/notebooks/natcompsci2022'), repo_data=None)[source]#

Stage the legacy mir/R outputs into data_legacy/ (gzipped) — comparison only.

Holds the 2022 baselines the new tcren results are measured against: the published TCRen matrix, mir contacts, the old non-redundancy summary, the legacy MHC annotation, the mir annotation oracle, paper source_data/ and the other-tools (TITAN/ERGO-II) outputs. Never consumed as a pipeline input.

Parameters:

paper_dir (Path)
repo_data (Path | None)

Return type:

int

tcren.paper.bootstrap.bootstrap(paper_dir=PosixPath('/home/runner/work/tcren/tcren/notebooks/natcompsci2022'), data_dir=PosixPath('/home/runner/work/tcren/tcren/notebooks/data'), structures=True, canonical=False, **_legacy_flags)[source]#

Fetch the HF structure sets into data_dir (notebooks/data by default).

Each set lands in its own folder directly under data_dir (e.g. notebooks/data/ Native2026/) — no structures/ wrapper — and is gitignored. The non-structure inputs (vdjdb, Birnbaum, MJ/Keskin, IEDB, epitope lists, PDB dates) and the legacy comparison baselines are already committed under <paper_dir>/data_legacy/, so they are not re-fetched here. Pass canonical=True to also fetch the re-oriented Canonical2026 set.

Parameters:

paper_dir (Path)
data_dir (Path)
structures (bool)
canonical (bool)

Return type:

dict

tcren.paper.helpers module#

Helpers for the Nat Comput Sci 2022 reproduction notebooks.

contact_table replaces the legacy mir extract_contact_map (it returns the same TCR↔peptide contact columns the R analyses consume, computed through the tcren pipeline). compare is the small regression utility behind 07_compare_legacy.ipynb.

tcren.paper.helpers.contact_table(structure, cutoff=5.0, count_atoms=False)[source]#

TCR↔peptide contact table for an annotated structure (the mir-replacement).

The structure must already be chain-typed (classify_chains) and MHC-annotated (annotate_mhc). Returns the columns the R benchmarks use: pdb.id, chain.type.from, region.type.from, residue.index.from, residue.index.to, pos.from, pos.to, residue.aa.from, residue.aa.to.

When count_atoms is set, an extra n_atom_contacts column (the heavy-atom-pair count per residue pair) is carried through for atomic-weighted scoring. Default False keeps the schema byte-identical to the legacy output.

Parameters:

structure (Structure)
cutoff (float)
count_atoms (bool)

Return type:

DataFrame

tcren.paper.helpers.annotate_structure_set(struct_dir, on_error='skip', count_atoms=False)[source]#

Run the tcren pipeline over a folder of PDBs → (contacts, markup) tables.

Replaces the legacy mir batch annotation. contacts is the stacked TCR↔peptide contact_table(); markup is one row per structure with the CDR3α/CDR3β/peptide sequences + species (the inputs to non-redundancy clustering and the benchmarks). Species is auto-detected per structure by alignment score (human vs mouse). All chains across the whole folder are annotated in a single mmseqs call per organism (the per-call process overhead dominates, so dataset-level batching is far faster than per-structure annotation).

When count_atoms is set, each contact row carries an n_atom_contacts heavy-atom-pair count (needed for atomic-weighted scoring).

Parameters:

struct_dir (str | Path)
on_error (str)
count_atoms (bool)

Return type:

tuple[DataFrame, DataFrame]

tcren.paper.helpers.mhc_annotation(struct_dir, ids=None, organism='human', on_error='skip')[source]#

Per-structure MHC allele + class for a folder (tcren mapper) — fully batched.

Replaces the legacy PDB_MHC_annotation table. ids restricts to those PDB ids. Every chain is TCR-typed in one batched arda call (so MHC candidates can be found), then every candidate MHC chain across the whole folder is searched against the MHC reference in a single mmseqs easy_search (mmseqs parallelises internally — no Python process/thread pool, which would either deadlock on fork or re-pay the fixed mmseqs startup cost per structure). Returns pdb.id, mhc.class, mhc.allele, status.

Parameters:

organism (str)
on_error (str)

Return type:

DataFrame

tcren.paper.helpers.annotate_batch(structures, arda, organisms=('human', 'mouse'), threads=0)[source]#

Annotate every chain of every structure with one mmseqs call per organism.

Public since 2.3.0. Batching matters: arda/mmseqs costs seconds per call, so annotating a 1,000-structure cohort one at a time is minutes of pure index rebuild. Every benchmark that scores a cohort needs this, which is why four downstream scripts were reaching into the private name.

threads caps the mmseqs thread count for THIS call. It matters whenever batches are annotated concurrently: arda.annotate_sequences does not forward a thread count, so annotate_records falls back to its threads=0 default, which mmseqs reads as all cores. Twelve concurrent batches on a 16-core machine then ask for 192 threads and spend their time context-switching instead of searching. Pass max(1, cpu_count // concurrency). 0 keeps the all-cores default, which is right for a single non-concurrent call.

Returns records[struct_idx][organism][chain_id] — the per-structure slices fed to classify_chains() as precomputed_records.

Parameters:: threads (int)
Return type:: list[dict[str, dict[str, dict]]]

tcren.paper.helpers.compare(old_path, new_path, keys, value_cols=None, tol=1e-06)[source]#

Compare two tables on keys and report row-set + max numeric differences.

Returns {rows_old, rows_new, matched, only_old, only_new, max_abs_diff, status} where status is "pass" when the key sets agree and every shared numeric column differs by ≤ tol.

Parameters:

old_path (str | Path)
new_path (str | Path)
keys (list[str])
value_cols (list[str] | None)
tol (float)

Return type:

dict

Command line#

tcren.cli module#

Command-line interface for tcren.

Commands are grouped in tcren --help:

Scoring & prediction

tcren score — end-to-end candidate-epitope scoring (drop-in for run_TCRen.R).
tcren rank — percentile-rank a peptide’s energy against a random pMHC background.
tcren ddg — ΔΔG of peptide mutations (fast virtual-matrix path; alanine scan / neoantigens).
tcren binder — TCR binder vs non-binder from AF-orthogonal interface geometry.
tcren energy — DOPE atom-level interface interaction energy (the ΔΔG e_native scorer).
tcren mechanics — interface mechanics (stiffness / rupture / coupling) — the koff proxies.
tcren scoring — per-interface contact energies Φ (--delta for ΔΦ, --geometry for Q).

Annotation & contacts

tcren annotate — chain typing + region markup (TCR CDR/FR, MHC groove, peptide; --pseudo).
tcren contacts — annotated residue-pair contact table for an interface.

Orientation & refinement

tcren superimpose — orient structure(s) onto the canonical database by MHC.
tcren refine — potential-guided peptide-pose refinement (DOPE MC; optional --substitute).

Reference data & potentials

tcren orient — build a canonical database from native complexes.
tcren derive-potential — derive a TCRen potential from a contact-map table.
tcren fetch-data / fetch-recent — fetch reference sets / recent RCSB TCR-pMHC entries.
tcren build-mhc-ref — build the IMGT/HLA + mouse MHC allele reference.

Info

tcren info — version + dependency check.
tcren paper … — Nat Comput Sci 2022 reproduction helpers.

tcren.cli.paper_bootstrap(structures=<typer.models.OptionInfo object>, canonical=<typer.models.OptionInfo object>)[source]#

Fetch HF structure sets into notebooks/data/<Set>/ (gitignored; non-structure inputs are already committed under natcompsci2022/data_legacy/).

Parameters:

structures (bool)
canonical (bool)

Return type:

None

tcren.cli.info()[source]#

Show version and dependency availability.

Return type:: None

tcren.cli.annotate(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, pseudo=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>)[source]#

Annotate chains and emit a per-residue region-markup table.

Covers TCR (CDR/FR), MHC groove (helices/floor) and peptide in one pass — --regions restricts the output to one chain class. --pseudo additionally marks the NetMHCpan MHC pseudosequence residues (region MPS). MHC groove + MPS require MHC annotation, which runs automatically when needed.

Parameters:

structures (Path)
out (Path)
regions (str)
pseudo (bool)
organism (str)

Return type:

None

tcren.cli.contacts(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>, interface=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>)[source]#

Compute and emit an annotated contact table.

Parameters:

structures (Path)
out (Path)
cutoff (float)
interface (str)
regions (str)
organism (str)

Return type:

None

tcren.cli.orient(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, metadata=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, reference_id=<typer.models.OptionInfo object>, force_pca=<typer.models.OptionInfo object>, threads=<typer.models.OptionInfo object>, mmcif=<typer.models.OptionInfo object>, compress=<typer.models.OptionInfo object>)[source]#

Build a canonical database: orient native TCR-pMHC complexes into the common MHC frame.

Derives the per-class canonical frame and writes every complex into it (A–E chains). This is how the bundled Canonical2026 set is produced; use superimpose to bring a new structure into an existing canonical database.

Parameters:

structures (Path)
out (Path)
metadata (Path)
organism (str)
reference_id (str)
force_pca (bool)
threads (int)
mmcif (bool)
compress (bool)

Return type:

None

tcren.cli.shuffle(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, n=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, within_class=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, compress=<typer.models.OptionInfo object>)[source]#

Generate wrong-TCR-on-real-pMHC decoys (a Shuffled set) for recognition models.

Keeps each oriented complex’s pMHC intact and grafts on n different complexes’ TCRs (a within-MHC-class derangement, so no decoy reproduces a real pairing). Real (label 1) vs these decoys (label 0) trains a label-free TCR-recognition classifier. Inputs must be co-framed — run tcren orient first.

Parameters:

structures (Path)
out (Path)
n (int)
seed (int)
within_class (bool)
organism (str)
compress (bool)

Return type:

None

tcren.cli.superimpose(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, db=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, threads=<typer.models.OptionInfo object>, mmcif=<typer.models.OptionInfo object>, compress=<typer.models.OptionInfo object>)[source]#

Superimpose structure(s) onto a canonical database by MHC.

Detects each input’s MHC chains, class, and species, then superposes its conserved groove Cα onto every database structure of the same class and species and averages the transforms into one consensus placement. The database defaults to data/Canonical2026 (populated at install).

-s accepts a file, directory, .tar.gz, or a shell glob. -o is an output directory, or — for a single input — a structure file whose extension must match --mmCIF/--compress. Annotation is one batched mmseqs call; -t threads the alignment + write.

Parameters:

structures (str)
out (Path)
db (Path)
organism (str)
threads (int)
mmcif (bool)
compress (bool)

Return type:

None

tcren.cli.derive_potential(contact_maps=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, summary=<typer.models.OptionInfo object>, nonred=<typer.models.OptionInfo object>, structure_dir=<typer.models.OptionInfo object>, redundancy_t=<typer.models.OptionInfo object>, variant=<typer.models.OptionInfo object>, pseudocount=<typer.models.OptionInfo object>, loo=<typer.models.OptionInfo object>)[source]#

Derive a TCRen potential from observed contacts.

Provide contacts either as a precomputed -i CSV or as a --structure-dir of PDBs (assembled via annotate_structure_set); pass exactly one. With a structure directory, --redundancy-t additionally restricts derivation to one representative per non-redundant cluster of αβ complexes (PDBs→contacts→cluster→derive in one call).

Parameters:

contact_maps (Path | None)
out (Path)
summary (Path | None)
nonred (bool)
structure_dir (Path | None)
redundancy_t (float | None)
variant (str)
pseudocount (int)
loo (bool)

Return type:

None

tcren.cli.fetch_data(canonical=<typer.models.OptionInfo object>)[source]#

Populate data/ with the reference structure sets from the HF dataset.

Run once at install (setup.sh does this). Fetches Native2026 (orientation references) and, by default, Canonical2026 (the default superimpose database) into $TCREN_DATA_DIR / repo data/. Skips folders already present.

Parameters:: canonical (bool)
Return type:: None

tcren.cli.build_mhc_ref(species=<typer.models.OptionInfo object>, force_download=<typer.models.OptionInfo object>)[source]#

Download and curate the MHC allele reference (IMGT/HLA + UniProt mouse).

Parameters:

species (str)
force_download (bool)

Return type:

None

tcren.cli.fetch_recent(dest=<typer.models.OptionInfo object>, discover=<typer.models.OptionInfo object>, after=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>)[source]#

Download recent TCR-pMHC structures from RCSB into data/pdb_recent.

Seeds with the Native2026 ids; with –discover also full-text-searches RCSB for new entries. Each is pulled as mmCIF (.cif.gz; handles extended PDB ids), annotated, and kept only if it has all 5 required chains (MHCa + b2m/MHCb + peptide + TCR pair).

Parameters:

dest (Path)
discover (bool)
after (str)
organism (str)

Return type:

None

tcren.cli.score(structures=<typer.models.OptionInfo object>, candidates=<typer.models.OptionInfo object>, potential=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, interface=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>)[source]#

Score candidate epitopes against input structures (end-to-end pipeline).

Parameters:

structures (Path)
candidates (Path)
potential (str | None)
out (Path)
interface (str)
regions (str)
organism (str)
cutoff (float)

Return type:

None

tcren.cli.ddg_cmd(structures=<typer.models.OptionInfo object>, native=<typer.models.OptionInfo object>, alanine_scan=<typer.models.OptionInfo object>, mutant=<typer.models.OptionInfo object>, potential=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, interface=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>)[source]#

ΔΔG of peptide mutations (fast virtual-matrix path; no atoms move).

Re-scores the mutant sequence on the native contact map; ddG = E(native) - E(mutant) (positive => STABILISING: the mutant scores lower, which is the better binder). Use --alanine-scan for a per-position scan, or one or more --mutant for specific neoantigen substitutions.

Parameters:

structures (Path)
native (str)
alanine_scan (bool)
mutant (list[str])
potential (str | None)
out (Path)
interface (str)
regions (str)
organism (str)
cutoff (float)

Return type:

None

tcren.cli.cpl_cmd(structures=<typer.models.OptionInfo object>, peptide=<typer.models.OptionInfo object>, position=<typer.models.OptionInfo object>, mutation=<typer.models.OptionInfo object>, to_mixture=<typer.models.OptionInfo object>, reference=<typer.models.OptionInfo object>, potential=<typer.models.OptionInfo object>, mhc_potential=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>)[source]#

Predict a combinatorial-peptide-library response matrix from a template TCR:pMHC structure.

One row per (peptide position, amino acid) cell. Every cell carries BOTH peptide-bearing interfaces summed – TCRen over TCR:peptide plus Miyazawa-Jernigan over peptide:MHC – because the assay reads activation, which needs the peptide presented as well as the receptor engaged.

Two reference states are reported, and a cell is only meaningful against one of them: effect_equimolar scores a residue against the 1/20 mixture, which is the CPL background and the right axis to compare with a measured matrix; effect_wild_type scores it against the residue the template carries, which is the mutation-scan / neoantigen question. Positive is favourable on both.

Three narrower questions, all from the same matrix:

: –position 5 every substitution at position 5, best first –position 5 –mutation W just that one cell –position 5 –to-mixture the cost of giving position 5 up to the 1/20 mixture

Parameters:

structures (Path)
peptide (str)
position (int)
mutation (str)
to_mixture (bool)
reference (str)
potential (str | None)
mhc_potential (str | None)
out (Path)
regions (str)
organism (str)
cutoff (float)

Return type:

None

tcren.cli.rank(structures=<typer.models.OptionInfo object>, candidates=<typer.models.OptionInfo object>, potential=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, interface=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, background=<typer.models.OptionInfo object>, background_source=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>)[source]#

Percentile-rank peptides’ TCRen energy against a random pMHC background.

For each structure, scores the supplied candidate peptides (or the structure’s own peptide when -c is omitted) together with --background random peptides of the same length and reports rank_pct — the fraction of background scoring at least as well (lower energy = better binder, so a small rank_pct means a strong binder).

Parameters:

structures (Path)
candidates (Path)
potential (str | None)
out (Path)
interface (str)
regions (str)
background (int)
background_source (Path)
seed (int)
organism (str)
cutoff (float)

Return type:

None

tcren.cli.binder(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>, features_only=<typer.models.OptionInfo object>)[source]#

Predict TCR binder/non-binder from AF-orthogonal interface geometry (native _geom + frozen model).

Scores each complex from interface size, dual-chain balance, H-bonds, buried ΔSASA and the CDR1/2-vs-CDR3α TCRen potential — signal that beats AlphaFold/TCRmodel2 confidence for ranking candidate TCRs against a fixed pMHC. All descriptors are computed natively (no PyRosetta/Biopython SASA/sklearn). Low p_bind = unlikely binder.

Parameters:

structures (Path)
out (Path)
organism (str)
cutoff (float)
features_only (bool)

Return type:

None

tcren.cli.recognize(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, features_only=<typer.models.OptionInfo object>, full=<typer.models.OptionInfo object>, scores=<typer.models.OptionInfo object>, mechanics=<typer.models.OptionInfo object>, cohort=<typer.models.OptionInfo object>, iptm=<typer.models.OptionInfo object>, invert_f_thresh=<typer.models.OptionInfo object>, threads=<typer.models.OptionInfo object>, autodetect_species=<typer.models.OptionInfo object>)[source]#

Full interface descriptor table + joint P(real) for each TCR-pMHC complex (one TSV row per PDB).

One row per structure with the complete recognition feature set (tcren.recognition.RECOGNITION_FEATURES): docking geometry (pitch, crossing, the 6 TCRdock rigid-body params), per-interface energies F_{tcr_pep,tcr_mhc,pep_mhc} and poly-alanine dF, CDR-loop energies F_{cdr12,cdr3a,cdr3b}, contact-type tallies, ΔSASA burial and the MHC-class indicator — plus p_real (the distribution-aware Bayesian logistic) and p_real_bn (the Gaussian BN): the joint probability the complex is a genuine recognition interface rather than a wrong-TCR shuffle. --full also emits the 18 CDR3-local frame descriptors (the FramePose strain layer). --scores adds the recommended fit-free q_bind (binder-ID; the directional-decorrelated interface-quality score, calibrated on the native crystal reference so it is defined per structure and transfers) and s_strain (forced-pose), alongside the fitted p_bind / p_forced. --features-only skips the models. Output is TSV.

--mechanics appends the koff proxies tcren mechanics reports — stiffness tensor, steered rupture, coupling residues — to these same rows. Prefer it to running the two commands: they need the identical annotated structure, so the second command repeats the parse and both mmseqs searches to produce a second table in a different format (CSV) under a different key (pdb.id) that then has to be joined. One flag costs about a sixth of the descriptor pass and returns one table.

Complementary scorer on the same inputs: tcren ddg (per-mutation alanine/neoantigen ΔΔF).

Examples:

tcren recognize -s models/ -o out.tsv                        # descriptors + P(real)
tcren recognize -s models/ --scores -o out.tsv               # + fit-free q_bind/s_strain + fitted p_bind
tcren recognize -s models/ --cohort -o out.tsv               # Q_geom + F_score + z(Q)±z(F) (no ipTM needed)
tcren recognize -s models.tar.gz --iptm meta.tsv -o out.tsv  # + AF synergy: z(ipTM)+z(Q_geom), z(ipTM)+z(Q)+z(F)
tcren recognize -s models/ --scores --mechanics -t 0 -o out.tsv   # every descriptor this tool reports, one table

Synergy with AlphaFold, made automatic. The contact energy F reads real binding chemistry but inverts on forced poses (benchmark ledger C27/C42), and ipTM is AlphaFold’s own pose-confidence signal — so with --iptm the command flags each low-ipTM (forced) pose in F_invert and emits z(Q)+z(F|iptm), which applies +z(F) to confident poses and -z(F) to forced ones (threshold --invert-f-thresh, default 0.5). It also prints how many poses are forced. Without --iptm it tells you F is being trusted unconditionally and how to gate it (--iptm or s_strain). z(ipTM)+z(Q_geom) is the geometry-only channel that is robust to the inversion without needing F.

Parameters:

structures (str)
out (Path)
organism (str)
features_only (bool)
full (bool)
scores (bool)
mechanics (bool)
cohort (bool)
iptm (Path | None)
invert_f_thresh (float)
threads (int)
autodetect_species (bool)

Return type:

None

tcren.cli.scoring(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, no_superimpose=<typer.models.OptionInfo object>, db=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>, tcr_peptide_potential=<typer.models.OptionInfo object>, tcr_mhc_potential=<typer.models.OptionInfo object>, peptide_mhc_potential=<typer.models.OptionInfo object>, regions=<typer.models.OptionInfo object>, contact_weight=<typer.models.OptionInfo object>, delta=<typer.models.OptionInfo object>, reference_aa=<typer.models.OptionInfo object>, geometry=<typer.models.OptionInfo object>, skip_errors=<typer.models.OptionInfo object>, threads=<typer.models.OptionInfo object>)[source]#

Score structures: per-interface contact energies Φ (and ΔΦ, and interface geometry).

This is scoring only — it reads structures and writes numbers. The preparation steps (canonicalisation, region mapping, Cα / contact / atom-distance matrices) are separate commands: tcren annotate, tcren superimpose, tcren contacts.

Columns F_tcr_pep, F_tcr_mhc, F_pep_mhc are the three interface terms Φ_TP, Φ_TM, Φ_PM; F_total is their sum Φ. With --delta each also gets its poly-alanine-referenced counterpart dF_* (ΔΦ_TP, ΔΦ_TM≡0, ΔΦ_PM) and dF_total = ΔΦ. The names match tcren recognize, so the two tables join on pdb.id. ΔΦ is the score to use across candidates that each carry their own generated pose, where raw Φ partly reads the pose geometry rather than the peptide sequence.

--geometry appends the interface descriptors (buried surface burial, peptide coverage n_pep_contacted, chain_balance, n_hbond, docking pitch/crossing) and Q — the directional, decorrelated interface-quality score, standardised against the native-crystal reference so it is defined for a single structure (tcren.q_score()). For the complete descriptor catalogue plus P(real), use tcren recognize.

Each interface’s potential can be overridden with a bundled name (tcren/mj/ keskin) or a CSV path; an unset option keeps the default family for that interface.

Examples:

tcren scoring -s complex.pdb.gz -o scores.csv
tcren scoring -s a.pdb.gz -s b.pdb.gz --delta          # repeat -s, or comma-separate
tcren scoring -s 'models/*.pdb.gz' --delta --geometry  # quote the glob
tcren scoring -s models/ --delta -t 8                  # a directory, 8 workers
tcren scoring -s models.txt --delta                    # one path per line
tcren scoring -s models.tar.gz --regions cdr           # CDR contacts only

Scoring a cohort is embarrassingly parallel and dominated by the per-structure mmseqs annotation, so -t is worth setting for anything above a handful of structures (-t 0 uses every core).

Parameters:

structures (list[str])
out (Path)
no_superimpose (bool)
db (Path)
organism (str)
cutoff (float)
tcr_peptide_potential (str)
tcr_mhc_potential (str)
peptide_mhc_potential (str)
regions (str)
contact_weight (str)
delta (bool)
reference_aa (str)
geometry (bool)
skip_errors (bool)
threads (int)

Return type:

None

tcren.cli.pipeline(ctx)[source]#

Deprecated alias for tcren scoring (this command never ran the full pipeline).

Parameters:: ctx (Context)
Return type:: None

tcren.cli.energy(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, relax=<typer.models.OptionInfo object>, shell=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>)[source]#

DOPE atom-level interaction energy across the peptide↔partner interface (the _relax kernel).

Sums the DOPE potential over peptide↔partner heavy-atom pairs — the interface ΔG contribution of the peptide (lower = more favourable). With --relax it also reports the energy after a rigid-body DOPE refinement (tcren.refine_peptide()) and the relaxation gap = e_native − e_relax. This is the single-structure scorer behind the ΔΔG benchmark (e_native/e_relax).

Parameters:

structures (str)
out (Path)
relax (bool)
shell (float)
organism (str)
seed (int)

Return type:

None

tcren.cli.mechanics(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, cutoff=<typer.models.OptionInfo object>, weight=<typer.models.OptionInfo object>, direction=<typer.models.OptionInfo object>, break_strain=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, threads=<typer.models.OptionInfo object>, autodetect_species=<typer.models.OptionInfo object>)[source]#

Interface mechanics — the koff proxies: stiffness tensor + steered rupture + coupling residues.

Treats the TCR↔pMHC contact map as a network of breakable springs and reports, per structure: n_spring, S_tot/K_tens/K_shear/aniso (stiffness tensor), rupture_force/ rupture_work (steered unbinding), and couple_pep/couple_total (coupling residues). Validated on ATLAS: the tensile stiffness / rupture resistance track the dissociation off-rate (koff) far better than the equilibrium ΔG/Kd (Bell–Evans; the TCR is a mechanosensor).

Parameters:

structures (str)
out (Path)
cutoff (float)
weight (str)
direction (str)
break_strain (float)
organism (str)
threads (int)
autodetect_species (bool)

Return type:

None

tcren.cli.refine(structures=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, substitute=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>, n_steps=<typer.models.OptionInfo object>, restraint_w=<typer.models.OptionInfo object>, seed=<typer.models.OptionInfo object>, mmcif=<typer.models.OptionInfo object>, compress=<typer.models.OptionInfo object>)[source]#

Potential-guided rigid-body refinement of the peptide pose (knowledge-based, not physics).

Optionally --substitute a new equal-length peptide first, then run a Monte-Carlo refinement scored by the DOPE atom-level statistical potential (restrained to the input pose; independent of the TCRen/MJ scoring potentials). Writes one structure per input and prints the final DOPE energy. (For physics-grade relaxation use Rosetta FlexPepDock externally.)

Parameters:

structures (str)
out (Path)
substitute (str)
organism (str)
n_steps (int)
restraint_w (float)
seed (int)
mmcif (bool)
compress (bool)

Return type:

None

tcren.cli.substitute_tcr_cmd(host=<typer.models.OptionInfo object>, donor=<typer.models.OptionInfo object>, out=<typer.models.OptionInfo object>, by=<typer.models.OptionInfo object>, organism=<typer.models.OptionInfo object>)[source]#

Graft the donor TCR onto the host pMHC → a chimeric TCR:pMHC complex.

Keeps the host peptide + MHC and the donor TCR. --by mhc superposes the donor MHC groove onto the host groove (the donor TCR keeps its native docking geometry); --by tcr superposes the donor TCR onto the host TCR (the donor TCR inherits the host’s docking pose). Both inputs are chain-typed automatically (and, for --by mhc, MHC-annotated).

Parameters:

host (Path)
donor (Path)
out (Path)
by (str)
organism (str)

Return type:

None