Skip to content

Re-referencing

High-level entry point: ChemicalShifts.reref. The routines below operate on raw DataFrames.

Dispatch

compute_offsets

compute_offsets(df, method, n_std=_N_STD_OUTLIER)

Fit per-atom re-referencing offsets for a long-format shift table.

Parameters:

Name Type Description Default
df DataFrame

Long-format shifts with at least Comp_ID, Seq_ID, Atom_ID, Val.

required
method (lacs, panav)
'lacs'
n_std int

LACS only — BMRB-statistics outlier threshold (default 4).

_N_STD_OUTLIER

Returns:

Name Type Description
offsets dict {atom: float | None}

Correction per atom such that corrected = Val - offset; None where fitting failed. LACS atoms: CA, CB, C, N, H. PANAV atoms: N, CA, CB, C.

check dict {atom: bool}

Whether re-referencing succeeded for each atom.

Returns ``(None, None)`` if no backbone shifts are present.
Source code in makeshift/reref/__init__.py
def compute_offsets(df, method, n_std=_N_STD_OUTLIER):
    """
    Fit per-atom re-referencing offsets for a long-format shift table.

    Parameters
    ----------
    df : DataFrame
        Long-format shifts with at least Comp_ID, Seq_ID, Atom_ID, Val.
    method : {'lacs', 'panav'}
    n_std : int
        LACS only — BMRB-statistics outlier threshold (default 4).

    Returns
    -------
    offsets : dict {atom: float | None}
        Correction per atom such that ``corrected = Val - offset``; None where
        fitting failed. LACS atoms: CA, CB, C, N, H. PANAV atoms: N, CA, CB, C.
    check : dict {atom: bool}
        Whether re-referencing succeeded for each atom.

    Returns ``(None, None)`` if no backbone shifts are present.
    """
    if method not in ("lacs", "panav"):
        raise ValueError(f"method must be 'lacs' or 'panav', got {method!r}")

    work = df.copy()
    work = work.loc[
        work.Atom_ID.isin(_BACKBONE) | work.Atom_ID.str.contains("^HA", na=False)
    ]

    # average GLY HA/HA2 into one HA per residue
    gly = (work["Comp_ID"] == "GLY") & work["Atom_ID"].str.contains("HA", na=False)
    for seq_id, group in work[gly].groupby("Seq_ID"):
        mask = (work["Seq_ID"] == seq_id) & gly
        work.loc[mask, "Val"] = group["Val"].mean()
        idx = work.loc[mask].index
        work.loc[idx[0], "Atom_ID"] = "HA"
        if len(idx) > 1:
            work.loc[idx[1:], "Atom_ID"] = "HA2"

    work = work.loc[work.Atom_ID.isin(_BACKBONE)].copy()
    if work.empty:
        return None, None

    if method == "lacs":
        return reref_lacs(work, n_std=n_std)

    return reref_panav(work)

apply_offsets

apply_offsets(df, offsets)

Return a copy of df with Val corrected by per-atom offsets (Val - offset).

Atoms absent from offsets (or whose offset is None) are left unchanged, so sidechain shifts pass through untouched.

Source code in makeshift/reref/base.py
def apply_offsets(df, offsets):
    """
    Return a copy of ``df`` with ``Val`` corrected by per-atom ``offsets``
    (``Val - offset``).

    Atoms absent from ``offsets`` (or whose offset is None) are left
    unchanged, so sidechain shifts pass through untouched.
    """
    applied = {a: v for a, v in (offsets or {}).items() if v is not None}
    out = df.copy()
    if applied:
        out["Val"] = out["Val"] - out["Atom_ID"].map(applied).fillna(0.0)
    return out

LACS

reref_lacs

reref_lacs(df, n_std=_N_STD_OUTLIER)

Compute per-atom LACS offsets for a backbone shift table.

Source code in makeshift/reref/lacs.py
def reref_lacs(df, n_std=_N_STD_OUTLIER):
    """Compute per-atom LACS offsets for a backbone shift table."""

    from ..chemshift import ChemicalShifts  # reuse secondary-shift / CSI logic

    df = df.loc[:, ~df.columns.duplicated()].copy()
    df["reref_mask"] = True
    outliers = df.apply(
        lambda r: _is_outlier(r["Comp_ID"], r["Atom_ID"], r["Val"], n_std), axis=1
    )
    df.loc[outliers, "reref_mask"] = False

    # secondary shift: CA/CB/N/H from the random-coil table; C' from its own
    df["secondary_shift"] = df.apply(ChemicalShifts._secondary_shift, axis=1)
    c_mask = df["Atom_ID"] == "C"
    df.loc[c_mask, "secondary_shift"] = df.loc[c_mask].apply(
        lambda r: _c_prime_secondary_shift(r["Comp_ID"], r["Val"]), axis=1
    )

    # CSI (CA - CB secondary shift) per residue, strict so a CA-only value
    # never enters the fit; csi_prev is the previous residue's CSI
    work = ChemicalShifts(df)
    df["csi"] = df.apply(lambda r: work._csi_raw(r, strict=True), axis=1)
    csi_by_seq = df.drop_duplicates("Seq_ID").set_index("Seq_ID")["csi"]
    df["csi_prev"] = df["Seq_ID"].map(lambda s: csi_by_seq.get(s - 1, np.nan))

    raw = {}
    for atom in ("CA", "CB"):
        x, y = _prepare_fit_data(df[df["Atom_ID"] == atom], "csi")
        raw[atom] = _piecewise_offset(x, y)

    x, y = _prepare_fit_data(df[df["Atom_ID"] == "C"], "csi")
    raw["C"] = _single_line_offset(x, y)

    for atom in ("N", "H"):
        expected_slope, slope_tight, threshold = _ATOM_PARAMS[atom]
        x, y = _prepare_fit_data(df[df["Atom_ID"] == atom], "csi_prev")
        raw[atom] = _lacs_offset_linear(
            x, y, expected_slope, slope_tight, threshold, _SLOPE_TOL[atom]
        )

    check = {atom: _fitted(raw.get(atom)) for atom in _LACS_ATOMS}
    offsets = {atom: (raw[atom] if check[atom] else None) for atom in _LACS_ATOMS}
    return offsets, check

PANAV

reref_panav

reref_panav(df)

Compute cumulative per-atom PANAV offsets over two rounds. Returns (offsets, check), offsets such that corrected = Val - offset.

Source code in makeshift/reref/panav.py
def reref_panav(df):
    """Compute cumulative per-atom PANAV offsets over two rounds.
    Returns (offsets, check), offsets such that corrected = Val - offset."""
    df = df.copy()
    df["Atom_ID"] = df["Atom_ID"].replace("HA2", "HA")
    df["outlier_1"] = False
    df["outlier_2"] = False

    check = {atom: True for atom in _PANAV_ATOMS}
    cumulative = {atom: 0.0 for atom in _PANAV_ATOMS}

    for round_i in range(2):
        df[["ss_probs", "ss_max"]] = df.apply(
            _panav_ss_probs, axis=1, result_type="expand"
        )
        ha = df.loc[df.Atom_ID == "HA"]
        df[["offset", "ss_max"]] = df.apply(
            lambda r: _panav_get_offset(r, ha.loc[ha.Seq_ID == r["Seq_ID"]]),
            axis=1, result_type="expand",
        )

        outlier_col = f"outlier_{round_i + 1}"
        df[outlier_col] = df.apply(_panav_judge_outlier, axis=1)

        round_offsets = {}
        for atom in _PANAV_ATOMS:
            mask = (df.Atom_ID == atom) & (~df["outlier_1"])
            if round_i == 1:
                mask &= ~df["outlier_2"]
            vals = np.array(df.loc[mask, "offset"], dtype=float)

            if np.isnan(vals).all():
                round_offsets[atom] = None
                check[atom] = False
                continue

            m, s = np.nanmean(vals), np.nanstd(vals)
            vals[(vals > m + 3 * s) | (vals < m - 3 * s)] = np.nan
            if (~np.isnan(vals)).sum() < 25:
                round_offsets[atom] = None
                check[atom] = False
                continue

            round_offsets[atom] = float(np.nanmean(vals))

        apply_now = {a: (0.0 if v is None else v) for a, v in round_offsets.items()}
        df["Val"] = df.apply(
            lambda r: r["Val"] - apply_now.get(r["Atom_ID"], 0.0), axis=1
        )
        for atom in _PANAV_ATOMS:
            if round_offsets.get(atom) is not None:
                cumulative[atom] += round_offsets[atom]

        failed = [a for a, ok in check.items() if not ok]
        df.loc[df["Atom_ID"].isin(failed), outlier_col] = True

    offsets = {atom: (cumulative[atom] if check[atom] else None)
               for atom in _PANAV_ATOMS}
    return offsets, check