Generation Probability (Pgen) API

Lightweight generation-probability estimator for TCR β/α CDR3 sequences. Used to flag likely-public CDR3s so DB matches against them can be appropriately discounted (#58).

Why a proxy, not OLGA

The gold-standard tool is OLGA (Sethna 2019), which fits explicit V/D/J/insertion models and computes Pgen by dynamic programming. OLGA gives the right answer but ships ~30 MB of pre-trained model files plus a runtime dependency chain.

Most tcrsift use cases (publicness discounting, ranking convergent vs. private sequences) just need the rank order, not numerically calibrated Pgens. This module provides that proxy in ~250 lines with no runtime dependency beyond numpy.

Decomposition

log Pgen(CDR3, V, J, chain)
  ≈ log P(V)             # V-gene usage marginal
  + log P(J)             # J-gene usage marginal
  + log P(length)        # CDR3 AA length, per chain
  + log P(n_inserted)    # total non-templated nt
  + log P(CDR3 | length) # uniform-AA composition baseline

The N-insertion term is the largest contributor to publicness: sequences with 0 non-templated nucleotides are convergently generatable across donors and dominate public-match noise. We model N-insertion length as geometric per junction (β has two junctions, α has one); the total N-insertion count is estimated from observed CDR3 nt length minus a typical templated V+J+D contribution.

Calibration

The estimator runs ~10 log units lower than real OLGA Pgens because of the uniform-over-20 AA composition baseline. Default publicness_score cutoffs are calibrated to this estimator's scale (-30 / -18); pass OLGA cutoffs (-20 / -8) if feeding in real OLGA values, or use auto_quantile=True for distribution-based auto-calibration.

Usage

from tcrsift.pgen import annotate_publicness

clones = annotate_publicness(clones)
public = clones[clones["publicness"] >= 0.5]
print(f"{len(public)} likely-public clones; DB matches should be discounted")

Wired into annotate_clonotypes:

from tcrsift import annotate_clonotypes

annotated = annotate_clonotypes(
    clones, vdjdb_path="vdjdb.tsv", add_publicness=True,
)

pgen

Lightweight Pgen estimator for TCR β / α CDR3 sequences (#58).

The gold-standard generation-probability tool is OLGA (Sethna et al. 2019), which fits V/D/J/insertion models and computes Pgen by dynamic programming. OLGA gives the right answer but ships ~30 MB of pre-trained model files plus a runtime dependency chain.

Most tcrsift use cases (publicness discounting in DB annotation, ranking convergent vs. private sequences) just need a coarse Pgen proxy — values that rank sequences from "highly generatable" to "clearly private", not necessarily numerically calibrated. This module provides that proxy in <200 lines, with no runtime dependency beyond numpy.

Decomposition:

log Pgen(CDR3, V, J, chain)
  ≈ log P(V)             # gene-usage marginal
  + log P(J)             # gene-usage marginal
  + log P(length)        # CDR3 AA length, per chain
  + log P(n_inserted)    # estimated total non-templated nt
  + log P(CDR3 | length) # length-normalized AA-composition term

For an observed sequence, n_inserted is estimated from the CDR3 nucleotide length minus a typical templated boundary contribution (see :data:_TEMPLATED_NT_TYPICAL). When the caller doesn't have nucleotide-level data, we approximate from AA length.

Marginal parameters bundled inline (not as external data files) — they're small (<3 KB total), and inlining keeps the module self-contained and easy to audit. Sources are documented per-table.

For real Pgens use OLGA; for ranking/discounting :func:compute_pgen is sufficient.

pgen_single

pgen_single(cdr3_aa: str, v_gene: str | None, j_gene: str | None, *, chain: str = 'beta', n_inserted: int | None = None) -> float

Return log10 Pgen estimate for one sequence.

Source code in tcrsift/pgen.py

def pgen_single(
    cdr3_aa: str,
    v_gene: str | None,
    j_gene: str | None,
    *,
    chain: str = "beta",
    n_inserted: int | None = None,
) -> float:
    """Return log10 Pgen estimate for one sequence."""
    parts = pgen_components(
        cdr3_aa, v_gene, j_gene, chain=chain, n_inserted=n_inserted
    )
    return float(sum(parts.values()))

pgen_components

pgen_components(cdr3_aa: str, v_gene: str | None, j_gene: str | None, *, chain: str = 'beta', n_inserted: int | None = None) -> dict[str, float]

Decompose log10 Pgen into its five terms.

Useful for debugging: shows which term is driving a sequence toward "public" or "private".

Source code in tcrsift/pgen.py

def pgen_components(
    cdr3_aa: str,
    v_gene: str | None,
    j_gene: str | None,
    *,
    chain: str = "beta",
    n_inserted: int | None = None,
) -> dict[str, float]:
    """Decompose log10 Pgen into its five terms.

    Useful for debugging: shows which term is driving a sequence
    toward "public" or "private".
    """
    chain = chain.lower()
    if chain not in ("alpha", "beta"):
        raise ValueError(f"chain must be 'alpha' or 'beta', got {chain!r}")
    if n_inserted is None:
        n_inserted = _estimate_n_inserted(cdr3_aa, chain)
    L = len(cdr3_aa) if isinstance(cdr3_aa, str) else 0
    LN10 = math.log(10.0)
    return {
        "log10_p_v":          _v_usage_log(v_gene or "", chain) / LN10,
        "log10_p_j":          _j_usage_log(j_gene or "", chain) / LN10,
        "log10_p_length":     _length_log(L, chain) / LN10,
        "log10_p_n_inserted": _n_inserted_log(n_inserted, chain) / LN10,
        "log10_p_composition": _composition_log(cdr3_aa) / LN10,
    }

compute_pgen

compute_pgen(df: DataFrame, *, cdr3_col: str = 'CDR3_beta', v_gene_col: str = 'beta_v_gene', j_gene_col: str = 'beta_j_gene', chain: str = 'beta', n_inserted_col: str | None = None, output_col: str = 'log10_pgen') -> pd.Series

Compute log10 Pgen estimate for each row in df.

Parameters:

Name	Type	Description	Default
df	DataFrame	Per-clone frame with at least the cdr3_col.	required
cdr3_col	str	Column with CDR3 amino-acid sequence.	'CDR3_beta'
v_gene_col	str	Columns with V / J gene calls (allele suffix tolerated).	'beta_v_gene'
j_gene_col	str	Columns with V / J gene calls (allele suffix tolerated).	'beta_v_gene'
chain	str	"alpha" or "beta".	'beta'
n_inserted_col	str | None	Optional column with pre-computed total N-insertion length (nt). When None, estimated from CDR3 AA length and a typical templated contribution from :data:_TEMPLATED_NT_TYPICAL.	None
output_col	str	Name applied to the returned Series.	'log10_pgen'

Returns:

Type	Description
Series	log10 Pgen estimate per row, indexed like df and named output_col. Rows whose cdr3_col is empty or non-string get NaN.

Source code in tcrsift/pgen.py

def compute_pgen(
    df: pd.DataFrame,
    *,
    cdr3_col: str = "CDR3_beta",
    v_gene_col: str = "beta_v_gene",
    j_gene_col: str = "beta_j_gene",
    chain: str = "beta",
    n_inserted_col: str | None = None,
    output_col: str = "log10_pgen",
) -> pd.Series:
    """Compute log10 Pgen estimate for each row in ``df``.

    Parameters
    ----------
    df : pd.DataFrame
        Per-clone frame with at least the ``cdr3_col``.
    cdr3_col : str
        Column with CDR3 amino-acid sequence.
    v_gene_col, j_gene_col : str
        Columns with V / J gene calls (allele suffix tolerated).
    chain : str
        ``"alpha"`` or ``"beta"``.
    n_inserted_col : str | None
        Optional column with pre-computed total N-insertion length
        (nt). When None, estimated from CDR3 AA length and a typical
        templated contribution from :data:`_TEMPLATED_NT_TYPICAL`.
    output_col : str
        Name applied to the returned Series.

    Returns
    -------
    pd.Series
        log10 Pgen estimate per row, indexed like ``df`` and named
        ``output_col``. Rows whose ``cdr3_col`` is empty or non-string
        get NaN.
    """
    chain = chain.lower()
    if chain not in ("alpha", "beta"):
        raise ValueError(f"chain must be 'alpha' or 'beta', got {chain!r}")
    if cdr3_col not in df.columns:
        raise ValueError(f"compute_pgen: missing {cdr3_col!r} column")

    cdr3 = df[cdr3_col]
    valid_mask = cdr3.apply(lambda s: isinstance(s, str) and bool(s))

    # Vectorize the four CDR3-length-derived terms (which don't depend
    # on V/J at all) and the V/J-only terms (which are dict lookups,
    # naturally per-row).
    LN10 = math.log(10.0)
    aa_lens = cdr3.where(valid_mask, "").apply(
        lambda s: len(s) if isinstance(s, str) else 0
    )

    # Length component — gaussian density at the AA length.
    mu, sigma = _CDR3_LENGTH_PARAMS[chain]
    z = (aa_lens.astype(float) - mu) / sigma
    pdf = np.exp(-0.5 * z * z) / (sigma * math.sqrt(2 * math.pi))
    log10_p_length = np.log(np.maximum(pdf, 1e-30)) / LN10

    # N-insertion component.
    if n_inserted_col and n_inserted_col in df.columns:
        n_ins_arr = pd.to_numeric(df[n_inserted_col], errors="coerce")
        n_ins_arr = n_ins_arr.where(
            n_ins_arr.notna(),
            other=(3 * aa_lens - _TEMPLATED_NT_TYPICAL[chain]),
        )
    else:
        n_ins_arr = 3 * aa_lens - _TEMPLATED_NT_TYPICAL[chain]
    n_ins_arr = n_ins_arr.clip(lower=0).astype(int)
    p_ins = _N_INSERT_GEOM_P[chain]
    if chain == "beta":
        # Negative binomial(r=2, p): P(N=n) = (n+1) (1-p)^n p^2
        log10_p_n = (
            np.log(n_ins_arr.astype(float) + 1.0)
            + n_ins_arr.astype(float) * math.log(1 - p_ins)
            + 2 * math.log(p_ins)
        ) / LN10
    else:
        log10_p_n = (
            n_ins_arr.astype(float) * math.log(1 - p_ins) + math.log(p_ins)
        ) / LN10

    # AA composition baseline — length × log(1/20).
    log10_p_composition = aa_lens.astype(float) * (_LOG_AA_BASELINE / LN10)

    # V / J gene-usage logs — dict-lookup per row.
    v_series = df[v_gene_col] if v_gene_col in df.columns else pd.Series([None] * len(df), index=df.index)
    j_series = df[j_gene_col] if j_gene_col in df.columns else pd.Series([None] * len(df), index=df.index)
    log10_p_v = v_series.apply(lambda g: _v_usage_log(g or "", chain) / LN10)
    log10_p_j = j_series.apply(lambda g: _j_usage_log(g or "", chain) / LN10)

    total = (
        log10_p_v + log10_p_j + log10_p_length
        + log10_p_n + log10_p_composition
    )
    # Empty/missing CDR3 → NaN, regardless of V/J availability.
    total = total.where(valid_mask, other=np.nan)
    total.name = output_col
    return total

publicness_score

publicness_score(log10_pgen: Series | Iterable[float], *, low_pgen_cutoff: float | None = -30.0, high_pgen_cutoff: float | None = -18.0, auto_quantile: bool = False, quantile_low: float = 0.1, quantile_high: float = 0.9) -> pd.Series

Map log10 Pgen → [0, 1] publicness in a monotone-decreasing way.

• log10 Pgen ≥ high_pgen_cutoff (very generatable) → 1.0
• log10 Pgen ≤ low_pgen_cutoff (very private) → 0.0
• linear interpolation in between.

A high publicness score means the sequence is likely to be public and DB matches against it should be discounted.

Defaults are calibrated to this module's coarse Pgen estimator, not OLGA. Typical β CDR3s land around log10 ≈ −25; very public short sequences with canonical V/J around log10 ≈ −18. The estimator is uniformly more negative than OLGA because its AA composition term is uniform-over-20 rather than position-specific. Override the cutoffs if you're feeding in real OLGA values (use ~−20 / −8 for OLGA).

Parameters:

Name	Type	Description	Default
log10_pgen	Series | Iterable[float]	Per-row log10 Pgen estimates. When a Series, its index is preserved in the returned Series.	required
low_pgen_cutoff	float | None	Fixed cutoffs. Ignored when auto_quantile=True.	-30.0
high_pgen_cutoff	float | None	Fixed cutoffs. Ignored when auto_quantile=True.	-30.0
auto_quantile	bool	When True, derive cutoffs from the input distribution rather than a fixed pair. Useful when the input's Pgen scale is unknown or when the fixed cutoffs would saturate the score for most rows (e.g. a narrow distribution all within ~5 log units).	False
quantile_low	float	Quantile cutoffs used when auto_quantile=True. Defaults (0.10, 0.90) clip to the central 80% of the input range.	0.1
quantile_high	float	Quantile cutoffs used when auto_quantile=True. Defaults (0.10, 0.90) clip to the central 80% of the input range.	0.1

Source code in tcrsift/pgen.py

def publicness_score(
    log10_pgen: pd.Series | Iterable[float],
    *,
    low_pgen_cutoff: float | None = -30.0,
    high_pgen_cutoff: float | None = -18.0,
    auto_quantile: bool = False,
    quantile_low: float = 0.10,
    quantile_high: float = 0.90,
) -> pd.Series:
    """Map log10 Pgen → [0, 1] publicness in a monotone-decreasing way.

    - ``log10 Pgen ≥ high_pgen_cutoff`` (very generatable) → 1.0
    - ``log10 Pgen ≤ low_pgen_cutoff`` (very private) → 0.0
    - linear interpolation in between.

    A high publicness score means the sequence is likely to be public
    and DB matches against it should be discounted.

    Defaults are calibrated to *this module's* coarse Pgen estimator,
    not OLGA. Typical β CDR3s land around log10 ≈ −25; very public
    short sequences with canonical V/J around log10 ≈ −18. The
    estimator is uniformly more negative than OLGA because its AA
    composition term is uniform-over-20 rather than position-specific.
    Override the cutoffs if you're feeding in real OLGA values
    (use ~−20 / −8 for OLGA).

    Parameters
    ----------
    log10_pgen : pd.Series | Iterable[float]
        Per-row log10 Pgen estimates. When a Series, its index is
        preserved in the returned Series.
    low_pgen_cutoff, high_pgen_cutoff : float | None
        Fixed cutoffs. Ignored when ``auto_quantile=True``.
    auto_quantile : bool
        When True, derive cutoffs from the input distribution rather
        than a fixed pair. Useful when the input's Pgen scale is
        unknown or when the fixed cutoffs would saturate the score for
        most rows (e.g. a narrow distribution all within ~5 log units).
    quantile_low, quantile_high : float
        Quantile cutoffs used when ``auto_quantile=True``. Defaults
        (0.10, 0.90) clip to the central 80% of the input range.
    """
    if isinstance(log10_pgen, pd.Series):
        arr = log10_pgen.astype(float)
    else:
        arr = pd.Series(np.asarray(list(log10_pgen), dtype=float))

    if auto_quantile:
        non_nan = arr.dropna()
        if non_nan.empty:
            # Nothing to calibrate against — return all-NaN preserving index.
            return pd.Series(np.nan, index=arr.index)
        low_pgen_cutoff = float(non_nan.quantile(quantile_low))
        high_pgen_cutoff = float(non_nan.quantile(quantile_high))
        if high_pgen_cutoff <= low_pgen_cutoff:
            # Degenerate distribution (constant or near-constant) —
            # publicness is undefined; return 0.5 for everything.
            return pd.Series(
                np.where(arr.notna(), 0.5, np.nan), index=arr.index
            )

    if low_pgen_cutoff is None or high_pgen_cutoff is None:
        raise ValueError(
            "low_pgen_cutoff/high_pgen_cutoff required when auto_quantile=False"
        )
    width = high_pgen_cutoff - low_pgen_cutoff
    if width <= 0:
        raise ValueError("high_pgen_cutoff must be > low_pgen_cutoff")
    clipped = arr.clip(lower=low_pgen_cutoff, upper=high_pgen_cutoff)
    score = (clipped - low_pgen_cutoff) / width
    return score.where(arr.notna(), other=np.nan)

annotate_publicness

annotate_publicness(df: DataFrame, *, cdr3_col: str = 'CDR3_beta', v_gene_col: str = 'beta_v_gene', j_gene_col: str = 'beta_j_gene', chain: str = 'beta', pgen_col: str = 'log10_pgen', publicness_col: str = 'publicness') -> pd.DataFrame

Augment df with log10_pgen + publicness columns.

Convenience wrapper combining :func:compute_pgen and :func:publicness_score. Does not modify the input frame.

Source code in tcrsift/pgen.py

def annotate_publicness(
    df: pd.DataFrame,
    *,
    cdr3_col: str = "CDR3_beta",
    v_gene_col: str = "beta_v_gene",
    j_gene_col: str = "beta_j_gene",
    chain: str = "beta",
    pgen_col: str = "log10_pgen",
    publicness_col: str = "publicness",
) -> pd.DataFrame:
    """Augment ``df`` with ``log10_pgen`` + ``publicness`` columns.

    Convenience wrapper combining :func:`compute_pgen` and
    :func:`publicness_score`. Does not modify the input frame.
    """
    out = df.copy()
    out[pgen_col] = compute_pgen(
        out,
        cdr3_col=cdr3_col,
        v_gene_col=v_gene_col,
        j_gene_col=j_gene_col,
        chain=chain,
    )
    out[publicness_col] = publicness_score(out[pgen_col])
    return out