TIL Selection API

Module for TIL-only timepoint harmonization and clone prioritization.

til_select

TIL timepoint harmonization and clone prioritization.

This module powers the tcrsift til-select command. It is designed to support the same input model as legacy ad-hoc scripts:

• per-timepoint consensus_annotations.<TP>.csv
• per-timepoint clonotypes.<TP>.csv
• per-timepoint filtered_contig_annotations.<TP>.csv
• per-timepoint sample_filtered_feature_bc_matrix.<TP>.h5

The workflow produces a harmonized table across timepoints, marker-expression scores from 10x GEX, optional public-DB annotation, and v2-style selection subset files.

parse_sample_args

parse_sample_args(items: list[str]) -> dict[str, dict[str, Path]]

Parse --samples / --inputs entries.

Expected item format: LABEL=CONSENSUS_PATH,CLONOTYPES_PATH (comma or semicolon separator).

Source code in tcrsift/til_select.py

def parse_sample_args(items: list[str]) -> dict[str, dict[str, Path]]:
    """
    Parse ``--samples`` / ``--inputs`` entries.

    Expected item format:
    ``LABEL=CONSENSUS_PATH,CLONOTYPES_PATH`` (comma or semicolon separator).
    """
    parsed: dict[str, dict[str, Path]] = {}
    for item in items:
        if "=" not in item:
            raise TCRsiftValidationError(
                f"Invalid sample mapping: '{item}'",
                hint="Use LABEL=CONSENSUS_PATH,CLONOTYPES_PATH",
            )
        label, rest = item.split("=", 1)
        label = label.strip()
        if not label:
            raise TCRsiftValidationError(
                f"Invalid sample mapping: '{item}'",
                hint="Label before '=' cannot be empty.",
            )
        if label in parsed:
            raise TCRsiftValidationError(
                f"Duplicate timepoint label: '{label}'",
                hint="Use unique labels (e.g. T1, T2, T3).",
            )

        if "," in rest:
            consensus_raw, clonotypes_raw = rest.split(",", 1)
        elif ";" in rest:
            consensus_raw, clonotypes_raw = rest.split(";", 1)
        else:
            raise TCRsiftValidationError(
                f"Invalid sample mapping: '{item}'",
                hint="Expected CONSENSUS_PATH,CLONOTYPES_PATH after label.",
            )

        consensus_path = validate_file_exists(Path(consensus_raw.strip()), f"{label} consensus")
        clonotypes_path = validate_file_exists(Path(clonotypes_raw.strip()), f"{label} clonotypes")
        parsed[label] = {"consensus": consensus_path, "clonotypes": clonotypes_path}

    return dict(sorted(parsed.items(), key=lambda kv: _timepoint_sort_key(kv[0])))

parse_config

parse_config(path: str | Path) -> dict[str, dict[str, Path]]

Parse YAML config mapping timepoints to consensus/clonotypes files.

Source code in tcrsift/til_select.py

def parse_config(path: str | Path) -> dict[str, dict[str, Path]]:
    """Parse YAML config mapping timepoints to consensus/clonotypes files."""
    path = validate_file_exists(path, "til-select config")
    with open(path) as f:
        raw = yaml.safe_load(f) or {}

    if "timepoints" in raw:
        raw = raw["timepoints"]

    if not isinstance(raw, dict):
        raise TCRsiftValidationError(
            "Invalid til-select config format",
            hint="Config must be a mapping of label -> {consensus, clonotypes}.",
        )

    base = path.parent
    parsed: dict[str, dict[str, Path]] = {}
    for label, entry in raw.items():
        if not isinstance(entry, dict):
            raise TCRsiftValidationError(
                f"Invalid config entry for {label}",
                hint="Each label must map to a dict with consensus/clonotypes paths.",
            )
        consensus = entry.get("consensus")
        clonotypes = entry.get("clonotypes")
        if not consensus or not clonotypes:
            raise TCRsiftValidationError(
                f"Config entry '{label}' missing paths",
                hint="Each label needs both 'consensus' and 'clonotypes'.",
            )
        consensus_path = Path(consensus)
        clonotypes_path = Path(clonotypes)
        if not consensus_path.is_absolute():
            consensus_path = base / consensus_path
        if not clonotypes_path.is_absolute():
            clonotypes_path = base / clonotypes_path
        parsed[str(label)] = {
            "consensus": validate_file_exists(consensus_path, f"{label} consensus"),
            "clonotypes": validate_file_exists(clonotypes_path, f"{label} clonotypes"),
        }

    return dict(sorted(parsed.items(), key=lambda kv: _timepoint_sort_key(kv[0])))

default_timepoint_inputs

default_timepoint_inputs(data_dir: str | Path) -> dict[str, dict[str, Path]] | None

Discover paired consensus_annotations.<TP>.csv + clonotypes.<TP>.csv.

Source code in tcrsift/til_select.py

def default_timepoint_inputs(data_dir: str | Path) -> dict[str, dict[str, Path]] | None:
    """Discover paired ``consensus_annotations.<TP>.csv`` + ``clonotypes.<TP>.csv``."""
    data_dir = Path(data_dir)
    if not data_dir.exists():
        return None

    discovered: dict[str, dict[str, Path]] = {}
    for consensus_path in sorted(data_dir.glob("consensus_annotations.*.csv")):
        m = re.match(r"^consensus_annotations\.(.+)\.csv$", consensus_path.name)
        if not m:
            continue
        label = m.group(1)
        clonotypes_path = data_dir / f"clonotypes.{label}.csv"
        if clonotypes_path.exists():
            discovered[label] = {
                "consensus": consensus_path,
                "clonotypes": clonotypes_path,
            }

    if not discovered:
        return None

    return dict(sorted(discovered.items(), key=lambda kv: _timepoint_sort_key(kv[0])))

load_from_consensus

load_from_consensus(consensus_path: str | Path, clonotypes_path: str | Path, count_column: str | None = None) -> tuple[pd.DataFrame, pd.DataFrame]

Load one timepoint from consensus+clonotypes and aggregate to unique paired CDR3.

Source code in tcrsift/til_select.py

def load_from_consensus(
    consensus_path: str | Path,
    clonotypes_path: str | Path,
    count_column: str | None = None,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Load one timepoint from consensus+clonotypes and aggregate to unique paired CDR3."""
    consensus_path = validate_file_exists(consensus_path, "consensus annotations")
    clonotypes_path = validate_file_exists(clonotypes_path, "clonotypes")

    consensus_df = pd.read_csv(consensus_path)
    pivot, filtered_consensus = _extract_paired_cdr3_from_consensus_df(consensus_df, Path(consensus_path))

    df_clono = pd.read_csv(clonotypes_path)
    if "clonotype_id" not in df_clono.columns:
        raise TCRsiftValidationError(
            f"Clonotypes file missing required column 'clonotype_id': {clonotypes_path}",
            hint="Expected 10x clonotypes CSV format.",
        )
    df_clono["clonotype_id_norm"] = normalize_clonotype_id(df_clono["clonotype_id"])
    count_col = _detect_count_column(df_clono, count_column=count_column)

    df_counts = df_clono[["clonotype_id_norm", count_col]].rename(
        columns={"clonotype_id_norm": "clonotype_id", count_col: "cell_count"}
    )
    df_counts["cell_count"] = _validate_count_series(
        df_counts["cell_count"], count_col=count_col, source_path=Path(clonotypes_path)
    )

    reads_umis = _aggregate_reads_umis(filtered_consensus, "clonotype_id_norm").rename(
        columns={"clonotype_id_norm": "clonotype_id"}
    )

    merged = pivot.merge(df_counts, on="clonotype_id", how="left").merge(
        reads_umis, on="clonotype_id", how="left"
    )
    merged = merged.dropna(subset=["cell_count"])

    agg = {"cell_count": ("cell_count", "sum")}
    if "umis_sum" in merged.columns:
        agg["umis_sum"] = ("umis_sum", "sum")
    if "reads_sum" in merged.columns:
        agg["reads_sum"] = ("reads_sum", "sum")
    meta_cols = [
        col
        for col in [
            "CDR3_alpha_nt",
            "CDR3_beta_nt",
            "alpha_full_aa",
            "alpha_full_nt",
            "beta_full_aa",
            "beta_full_nt",
            *[f"alpha_{segment}_aa" for segment in TCR_SEGMENTS_AA],
            *[f"beta_{segment}_aa" for segment in TCR_SEGMENTS_AA],
            *[f"alpha_{segment_nt}" for segment_nt in TCR_SEGMENTS_NT],
            *[f"beta_{segment_nt}" for segment_nt in TCR_SEGMENTS_NT],
            "alpha_v_gene",
            "alpha_j_gene",
            "alpha_c_gene",
            "beta_v_gene",
            "beta_d_gene",
            "beta_j_gene",
            "beta_c_gene",
        ]
        if col in merged.columns
    ]
    for col in meta_cols:
        agg[col] = (col, _mode_or_nan)

    counts = merged.groupby(["CDR3_alpha", "CDR3_beta"], as_index=False).agg(**agg)
    counts["CDR3ab"] = counts["CDR3_alpha"] + "_" + counts["CDR3_beta"]
    if counts["CDR3ab"].duplicated().any():
        raise TCRsiftValidationError(
            f"{Path(consensus_path).name}: duplicate CDR3ab rows after aggregation",
            hint="Check clonotype normalization and consensus metadata.",
        )

    out_cols = ["CDR3_alpha", "CDR3_beta", "CDR3ab", "cell_count"] + meta_cols
    out = counts[out_cols].copy()
    out = out.sort_values(["cell_count", "CDR3ab"], ascending=[False, True]).reset_index(drop=True)

    stats = counts.copy()
    return out, stats

compute_marker_scores_for_timepoint

compute_marker_scores_for_timepoint(consensus_path: str | Path, contig_csv_path: str | Path, gex_h5_path: str | Path, marker_genes: list[str]) -> tuple[pd.DataFrame, pd.DataFrame]

Compute per-cell and per-clonotype marker scores for one timepoint.

Marker scores are computed from CP10K and log1p(CP10K)-z values per marker.

Source code in tcrsift/til_select.py

def compute_marker_scores_for_timepoint(
    consensus_path: str | Path,
    contig_csv_path: str | Path,
    gex_h5_path: str | Path,
    marker_genes: list[str],
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """
    Compute per-cell and per-clonotype marker scores for one timepoint.

    Marker scores are computed from CP10K and log1p(CP10K)-z values per marker.
    """
    consensus_path = Path(validate_file_exists(consensus_path, "consensus annotations"))
    contig_csv_path = Path(validate_file_exists(contig_csv_path, "filtered contig annotations"))
    gex_h5_path = Path(validate_file_exists(gex_h5_path, "10x filtered_feature_bc_matrix.h5"))

    marker_genes = [g.strip().upper() for g in marker_genes if g.strip()]
    marker_genes = list(dict.fromkeys(marker_genes))

    try:
        marker_cells = _load_10x_h5_marker_counts(gex_h5_path, marker_genes)
        barcode_map = _load_barcode_to_clonotype(contig_csv_path)

        consensus_df = pd.read_csv(consensus_path)
        clonotype_pairs, _ = _extract_paired_cdr3_from_consensus_df(consensus_df, consensus_path)

        merged = marker_cells.merge(barcode_map, on="barcode", how="inner")
        merged = merged.merge(clonotype_pairs, on="clonotype_id", how="inner")
        if merged.empty:
            raise TCRsiftValidationError(
                "No overlap between GEX barcodes and paired clonotypes",
                hint=f"Timepoint from {consensus_path.name} produced zero linked rows.",
            )

        cp10k_cols = [f"cp10k_{g}" for g in marker_genes]
        merged["marker_mean_cp10k"] = merged[cp10k_cols].mean(axis=1)
        z_cols: list[str] = []
        for gene in marker_genes:
            cp10k_col = f"cp10k_{gene}"
            log_col = f"log1p_cp10k_{gene}"
            z_col = f"z_log1p_cp10k_{gene}"
            merged[log_col] = np.log1p(merged[cp10k_col])
            mean_val = float(merged[log_col].mean())
            std_val = float(merged[log_col].std(ddof=0))
            if std_val <= 1e-12:
                merged[z_col] = 0.0
            else:
                merged[z_col] = (merged[log_col] - mean_val) / std_val
            z_cols.append(z_col)
        merged["marker_mean_z"] = merged[z_cols].mean(axis=1)

        agg = {
            "n_cells_with_gex": ("barcode", "nunique"),
            "marker_score_cp10k": ("marker_mean_cp10k", "mean"),
            "marker_score_z": ("marker_mean_z", "mean"),
        }
        for gene in marker_genes:
            agg[f"score_count_{gene}"] = (f"count_{gene}", "mean")
            agg[f"score_cp10k_{gene}"] = (f"cp10k_{gene}", "mean")
            agg[f"score_log1p_cp10k_{gene}"] = (f"log1p_cp10k_{gene}", "mean")
            agg[f"score_z_{gene}"] = (f"z_log1p_cp10k_{gene}", "mean")

        clonotype_scores = (
            merged.groupby(["CDR3_alpha", "CDR3_beta", "CDR3ab"], as_index=False).agg(**agg)
        )
        return merged, clonotype_scores
    except Exception:
        # Test-compatible fallback path when H5 is mocked or non-standard.
        adata = sc.read_10x_h5(str(gex_h5_path))
        return _compute_marker_scores_from_adata(
            consensus_path=consensus_path,
            contig_csv_path=contig_csv_path,
            marker_genes=marker_genes,
            adata=adata,
        )

build_harmonized_table

build_harmonized_table(samples: dict[str, DataFrame], timepoint_order: list[str], rank_by: str = 'mean_frequency') -> pd.DataFrame

Build harmonized clone table across timepoints (v2-compatible join semantics).

Source code in tcrsift/til_select.py

def build_harmonized_table(
    samples: dict[str, pd.DataFrame],
    timepoint_order: list[str],
    rank_by: str = "mean_frequency",
) -> pd.DataFrame:
    """Build harmonized clone table across timepoints (v2-compatible join semantics)."""
    if rank_by not in RANK_METRICS:
        raise ValueError(f"rank_by must be one of: {', '.join(RANK_METRICS)}")

    key_cols = ["CDR3_alpha", "CDR3_beta", "CDR3ab"]
    merged = None
    count_cols = []
    metadata_frames: list[pd.DataFrame] = []
    for label in timepoint_order:
        if label not in samples:
            raise ValueError(f"Missing sample for timepoint: {label}")
        df = samples[label]
        meta_cols = [c for c in df.columns if c not in key_cols + ["cell_count"]]
        if meta_cols:
            metadata_frames.append(df[key_cols + meta_cols].copy())
        df = df[key_cols + ["cell_count"]].copy()
        count_col = f"cell_count_{label}"
        count_cols.append(count_col)
        df = df.rename(columns={"cell_count": count_col})
        if merged is None:
            merged = df
        else:
            merged = pd.merge(merged, df, on=key_cols, how="outer")

    if merged is None:
        raise ValueError("No samples provided")

    for col in count_cols:
        if col in merged.columns:
            merged[col] = merged[col].fillna(0).astype(int)

    if metadata_frames:
        metadata_all = pd.concat(metadata_frames, ignore_index=True)
        meta_cols_union = [c for c in metadata_all.columns if c not in key_cols]
        metadata_agg = {col: (col, _mode_or_nan) for col in meta_cols_union}
        metadata_unique = metadata_all.groupby(key_cols, as_index=False).agg(**metadata_agg)
        merged = merged.merge(metadata_unique, on=key_cols, how="left")

    freq_cols = []
    for label in timepoint_order:
        count_col = f"cell_count_{label}"
        if count_col not in merged.columns:
            continue
        freq_col = f"frequency_{label}"
        total = merged[count_col].sum()
        if total <= 0:
            merged[freq_col] = 0.0
        else:
            merged[freq_col] = merged[count_col] / total
        freq_cols.append(freq_col)

    merged["total_cells"] = merged[count_cols].sum(axis=1)
    merged["n_timepoints"] = (merged[count_cols] > 0).sum(axis=1)
    merged["is_paired_alpha_beta"] = True
    if freq_cols:
        merged["mean_frequency"] = merged[freq_cols].mean(axis=1)
        merged["max_frequency"] = merged[freq_cols].max(axis=1)

    if merged["CDR3ab"].duplicated().any():
        n_dup = int(merged["CDR3ab"].duplicated().sum())
        raise ValueError(f"harmonized: found {n_dup} duplicate CDR3ab rows")
    return sort_harmonized_by_rank(merged, rank_by=rank_by)

run_selection_pipeline

run_selection_pipeline(master_df: DataFrame, timepoint_order: list[str], marker_genes: list[str], min_cells_per_clone: int = 2, min_cd8_cp10k: float = 0.0, max_cd4_to_cd8_ratio: float = 1.0, increase_ratio_nonzero_min: float = 1.5, increase_ratio_all_timepoints_min: float = 1.25, immunogenic_percentile: float = 0.9, immunogenic_percentile_slack_frac: float = 0.01, immunogenic_min_cp10k: float = 0.0, immunogenic_require_above_median: bool = True, cytotoxic_last_min_z: float = 0.25, cytotoxic_last_min_cp10k: float = 0.05, became_cytotoxic_min_delta_z: float = 0.5, immunogenic_genes_preferred: list[str] | None = None, cytotoxic_genes_preferred: list[str] | None = None, cytolytic_genes_preferred: list[str] | None = None, antigen_response_genes_preferred: list[str] | None = None) -> tuple[pd.DataFrame, dict[str, pd.DataFrame], list[tuple[str, int]], list[tuple[str, int]], list[tuple[str, int]]]

Compute selection masks/subsets for promising TIL clonotypes.

Source code in tcrsift/til_select.py

def run_selection_pipeline(
    master_df: pd.DataFrame,
    timepoint_order: list[str],
    marker_genes: list[str],
    min_cells_per_clone: int = 2,
    min_cd8_cp10k: float = 0.0,
    max_cd4_to_cd8_ratio: float = 1.0,
    increase_ratio_nonzero_min: float = 1.5,
    increase_ratio_all_timepoints_min: float = 1.25,
    immunogenic_percentile: float = 0.90,
    immunogenic_percentile_slack_frac: float = 0.01,
    immunogenic_min_cp10k: float = 0.0,
    immunogenic_require_above_median: bool = True,
    cytotoxic_last_min_z: float = 0.25,
    cytotoxic_last_min_cp10k: float = 0.05,
    became_cytotoxic_min_delta_z: float = 0.5,
    immunogenic_genes_preferred: list[str] | None = None,
    cytotoxic_genes_preferred: list[str] | None = None,
    cytolytic_genes_preferred: list[str] | None = None,
    antigen_response_genes_preferred: list[str] | None = None,
) -> tuple[
    pd.DataFrame,
    dict[str, pd.DataFrame],
    list[tuple[str, int]],
    list[tuple[str, int]],
    list[tuple[str, int]],
]:
    """Compute selection masks/subsets for promising TIL clonotypes."""
    df = master_df.copy()
    if not timepoint_order:
        raise ValueError("No timepoints provided for selection pipeline")
    if min_cells_per_clone < 1:
        raise ValueError(f"min_cells_per_clone must be >= 1, got {min_cells_per_clone}")
    if increase_ratio_nonzero_min <= 0:
        raise ValueError(f"increase_ratio_nonzero_min must be > 0, got {increase_ratio_nonzero_min}")
    if increase_ratio_all_timepoints_min <= 0:
        raise ValueError(
            "increase_ratio_all_timepoints_min must be > 0, "
            f"got {increase_ratio_all_timepoints_min}"
        )
    if not (0 < immunogenic_percentile <= 1):
        raise ValueError(f"immunogenic_percentile must be in (0,1], got {immunogenic_percentile}")
    if not (0 <= immunogenic_percentile_slack_frac < 1):
        raise ValueError(
            "immunogenic_percentile_slack_frac must be in [0,1), "
            f"got {immunogenic_percentile_slack_frac}"
        )

    freq_cols = [f"frequency_{tp}" for tp in timepoint_order if f"frequency_{tp}" in df.columns]
    if len(freq_cols) < 2:
        raise ValueError("Need at least two frequency_<TP> columns for selection.")
    nonzero_increasing, all_positive_increasing = compute_increasing_masks_from_frequencies(
        df[freq_cols],
        ratio_nonzero_min=increase_ratio_nonzero_min,
        ratio_all_timepoints_min=increase_ratio_all_timepoints_min,
    )
    df["is_increasing"] = nonzero_increasing
    df["is_increasing_positive"] = all_positive_increasing
    df["is_increasing_nonzero"] = nonzero_increasing
    df["is_increasing_all_timepoints"] = all_positive_increasing

    if "total_cells" not in df.columns:
        count_cols = [f"cell_count_{tp}" for tp in timepoint_order if f"cell_count_{tp}" in df.columns]
        if not count_cols:
            raise ValueError("Missing total_cells and per-timepoint cell_count columns.")
        df["total_cells"] = df[count_cols].apply(pd.to_numeric, errors="coerce").fillna(0.0).sum(axis=1)
    else:
        df["total_cells"] = pd.to_numeric(df["total_cells"], errors="coerce").fillna(0.0)
    df["is_min_cells"] = df["total_cells"] >= float(min_cells_per_clone)

    first_tp = timepoint_order[0]
    last_tp = timepoint_order[-1]
    first_freq = pd.to_numeric(df.get(f"frequency_{first_tp}", 0.0), errors="coerce").fillna(0.0)
    last_freq = pd.to_numeric(df.get(f"frequency_{last_tp}", 0.0), errors="coerce").fillna(0.0)
    df["delta_frequency_last_vs_first"] = last_freq - first_freq
    df["fold_change_last_vs_first"] = np.where(
        first_freq > 0,
        last_freq / first_freq,
        np.where(last_freq > 0, np.inf, 1.0),
    )
    df["is_enriched"] = df["is_increasing_nonzero"]

    if "is_viral" in df.columns:
        df["is_non_viral"] = ~df["is_viral"].fillna(False).astype(bool)
    else:
        df["is_non_viral"] = True

    cd8_candidate_cols = [
        c for c in df.columns if c.startswith("score_cp10k_CD8A") or c.startswith("score_cp10k_CD8B")
    ]
    if cd8_candidate_cols:
        df["cd8_cp10k_max"] = (
            df[cd8_candidate_cols].apply(pd.to_numeric, errors="coerce").max(axis=1).fillna(0.0)
        )
    else:
        df["cd8_cp10k_max"] = 0.0
    df["is_cd8_positive"] = df["cd8_cp10k_max"] > float(min_cd8_cp10k)

    cd4_candidate_cols = [c for c in df.columns if c.startswith("score_cp10k_CD4")]
    if cd4_candidate_cols:
        df["cd4_cp10k_max"] = (
            df[cd4_candidate_cols].apply(pd.to_numeric, errors="coerce").max(axis=1).fillna(0.0)
        )
    else:
        df["cd4_cp10k_max"] = 0.0

    df["cd4_to_cd8_ratio"] = np.where(
        df["cd8_cp10k_max"] > 0,
        df["cd4_cp10k_max"] / df["cd8_cp10k_max"],
        np.where(df["cd4_cp10k_max"] > 0, np.inf, 0.0),
    )
    df["is_cd8_not_cd4"] = df["cd4_to_cd8_ratio"] < float(max_cd4_to_cd8_ratio)
    df["is_base_cd8_non_viral"] = df["is_non_viral"] & df["is_cd8_positive"] & df["is_cd8_not_cd4"]
    df["is_base_selected"] = df["is_base_cd8_non_viral"] & df["is_min_cells"]

    requested_genes = [g.strip().upper() for g in marker_genes if g.strip()]
    requested_genes = list(dict.fromkeys(requested_genes))

    if immunogenic_genes_preferred:
        immunogenic_default = [g.strip().upper() for g in immunogenic_genes_preferred if g.strip()]
    else:
        immunogenic_default = ["GZMB", "PRF1", "IFNG", "MKI67", "TNFRSF9"]
    immunogenic_genes = [g for g in immunogenic_default if g in requested_genes]
    if not immunogenic_genes:
        immunogenic_genes = [g for g in requested_genes if g not in {"CD4", "CD8A", "CD8B"}]
    if not immunogenic_genes:
        immunogenic_genes = immunogenic_default

    def _dedupe_valid(candidate: list[str], fallback: tuple[str, ...]) -> list[str]:
        deduped = list(dict.fromkeys([g.strip().upper() for g in candidate if g and g.strip()]))
        valid = [g for g in deduped if g in requested_genes]
        if valid:
            return valid
        return [g for g in fallback if g in requested_genes]

    cytotoxic_seed = (
        [g.strip().upper() for g in cytotoxic_genes_preferred if g.strip()]
        if cytotoxic_genes_preferred
        else ["GZMB", "PRF1", "IFNG", "MKI67", "TNFRSF9"]
    )
    cytotoxic_genes = _dedupe_valid(cytotoxic_seed, tuple(cytotoxic_seed))
    if not cytotoxic_genes:
        cytotoxic_genes = immunogenic_genes

    cytolytic_seed = (
        [g.strip().upper() for g in cytolytic_genes_preferred if g.strip()]
        if cytolytic_genes_preferred
        else list(CYTOLYTIC_GENES_DEFAULT)
    )
    cytolytic_genes = _dedupe_valid(cytolytic_seed, CYTOLYTIC_GENES_DEFAULT)

    antigen_seed = (
        [g.strip().upper() for g in antigen_response_genes_preferred if g.strip()]
        if antigen_response_genes_preferred
        else list(ANTIGEN_RESPONSE_GENES_DEFAULT)
    )
    antigen_response_genes = _dedupe_valid(antigen_seed, ANTIGEN_RESPONSE_GENES_DEFAULT)

    for tp in timepoint_order:
        z_cols = [f"score_z_{g}_{tp}" for g in cytotoxic_genes if f"score_z_{g}_{tp}" in df.columns]
        cp10k_cols = [
            f"score_cp10k_{g}_{tp}" for g in cytotoxic_genes if f"score_cp10k_{g}_{tp}" in df.columns
        ]
        if z_cols:
            df[f"cytotoxic_score_z_{tp}"] = (
                df[z_cols].apply(pd.to_numeric, errors="coerce").mean(axis=1).fillna(0.0)
            )
        else:
            df[f"cytotoxic_score_z_{tp}"] = 0.0
        if cp10k_cols:
            df[f"cytotoxic_score_cp10k_{tp}"] = (
                df[cp10k_cols].apply(pd.to_numeric, errors="coerce").mean(axis=1).fillna(0.0)
            )
        else:
            df[f"cytotoxic_score_cp10k_{tp}"] = 0.0

    def _add_panel_scores(panel_name: str, genes: list[str]) -> None:
        per_tp_cols = []
        for tp in timepoint_order:
            src_cols = [f"score_cp10k_{g}_{tp}" for g in genes if f"score_cp10k_{g}_{tp}" in df.columns]
            out_col = f"{panel_name}_score_cp10k_{tp}"
            if src_cols:
                df[out_col] = (
                    df[src_cols].apply(pd.to_numeric, errors="coerce").mean(axis=1).fillna(0.0)
                )
            else:
                # Fallback to resolved summary columns when per-timepoint cols are absent.
                series_list = []
                for g in genes:
                    _source, series = _resolve_gene_cp10k_series(df, g, timepoint_order)
                    if series.notna().any():
                        series_list.append(series.fillna(0.0))
                if series_list:
                    df[out_col] = pd.concat(series_list, axis=1).mean(axis=1)
                else:
                    df[out_col] = 0.0
            per_tp_cols.append(out_col)
        df[f"{panel_name}_score_cp10k_mean"] = df[per_tp_cols].mean(axis=1)
        df[f"{panel_name}_score_cp10k_max"] = df[per_tp_cols].max(axis=1)
        df[f"{panel_name}_score_cp10k_sum"] = df[per_tp_cols].sum(axis=1)

    _add_panel_scores("cytolytic", cytolytic_genes)
    _add_panel_scores("antigen_response", antigen_response_genes)

    first_z = pd.to_numeric(df.get(f"cytotoxic_score_z_{first_tp}", 0.0), errors="coerce").fillna(0.0)
    last_z = pd.to_numeric(df.get(f"cytotoxic_score_z_{last_tp}", 0.0), errors="coerce").fillna(0.0)
    first_cp10k = pd.to_numeric(
        df.get(f"cytotoxic_score_cp10k_{first_tp}", 0.0),
        errors="coerce",
    ).fillna(0.0)
    last_cp10k = pd.to_numeric(
        df.get(f"cytotoxic_score_cp10k_{last_tp}", 0.0),
        errors="coerce",
    ).fillna(0.0)

    df["cytotoxic_score_z_delta_last_vs_first"] = last_z - first_z
    df["cytotoxic_score_cp10k_delta_last_vs_first"] = last_cp10k - first_cp10k
    df["is_cytotoxic_at_last"] = (last_z >= cytotoxic_last_min_z) | (
        last_cp10k >= cytotoxic_last_min_cp10k
    )
    df["is_became_cytotoxic"] = (
        df["cytotoxic_score_z_delta_last_vs_first"] >= became_cytotoxic_min_delta_z
    )
    df["is_increasing_and_became_cytotoxic"] = (
        df["is_increasing_nonzero"] & df["is_became_cytotoxic"] & df["is_cytotoxic_at_last"]
    )

    top_flag_cols = []
    immunogenic_branch_counts_within_base: list[tuple[str, int]] = []
    for gene in immunogenic_genes:
        ranking_source, ranking_values = _resolve_gene_cp10k_series(df, gene, timepoint_order)
        flag_col = f"is_top_immunogenic_{gene}"
        df[flag_col] = False
        df[f"immunogenic_cp10k_source_{gene}"] = ranking_source if ranking_source else ""

        base_mask = df["is_base_selected"] & ranking_values.notna()
        if int(base_mask.sum()) == 0:
            df[f"immunogenic_cp10k_median_{gene}"] = np.nan
            immunogenic_branch_counts_within_base.append((gene, 0))
            continue

        median_value = float(ranking_values[base_mask].median())
        df[f"immunogenic_cp10k_median_{gene}"] = median_value
        floor = max(float(immunogenic_min_cp10k), 0.0)
        eligible = df["is_base_selected"] & ranking_values.gt(floor)
        if immunogenic_require_above_median:
            eligible &= ranking_values.gt(median_value)
        if int(eligible.sum()) == 0:
            df[f"immunogenic_cp10k_percentile_{gene}"] = np.nan
            df[f"immunogenic_cp10k_cutoff_{gene}"] = np.nan
            immunogenic_branch_counts_within_base.append((gene, 0))
            continue

        values = ranking_values[eligible]
        percentile_value = float(np.nanquantile(values.to_numpy(dtype=float), immunogenic_percentile))
        cutoff_value = percentile_value * (1.0 - immunogenic_percentile_slack_frac)
        df[f"immunogenic_cp10k_percentile_{gene}"] = percentile_value
        df[f"immunogenic_cp10k_cutoff_{gene}"] = cutoff_value
        selected = eligible & ranking_values.ge(cutoff_value)
        df.loc[selected, flag_col] = True
        top_flag_cols.append(flag_col)
        n_selected = int((df["is_base_selected"] & df[flag_col]).sum())
        immunogenic_branch_counts_within_base.append((gene, n_selected))

    if top_flag_cols:
        df["is_top_immunogenic_any"] = df[top_flag_cols].any(axis=1)
    else:
        df["is_top_immunogenic_any"] = False

    n_top_cytolytic = _select_top_panel_score(
        df,
        "cytolytic_score_cp10k_mean",
        "is_top_cytolytic_score",
        immunogenic_percentile,
        immunogenic_percentile_slack_frac,
        immunogenic_min_cp10k,
        immunogenic_require_above_median,
    )
    n_top_antigen_response = _select_top_panel_score(
        df,
        "antigen_response_score_cp10k_mean",
        "is_top_antigen_response_score",
        immunogenic_percentile,
        immunogenic_percentile_slack_frac,
        immunogenic_min_cp10k,
        immunogenic_require_above_median,
    )

    df["is_branch_top_immunogenic"] = df["is_top_immunogenic_any"]
    df["is_branch_cytolytic"] = df["is_top_cytolytic_score"]
    df["is_branch_antigen_response"] = df["is_top_antigen_response_score"]
    df["is_branch_enriched"] = df["is_increasing_nonzero"]
    df["is_branch_increasing"] = df["is_increasing_all_timepoints"]
    df["is_branch_any"] = (
        df["is_branch_top_immunogenic"]
        | df["is_branch_cytolytic"]
        | df["is_branch_antigen_response"]
        | df["is_branch_enriched"]
        | df["is_branch_increasing"]
    )
    df["is_immunogenic"] = (
        df["is_top_immunogenic_any"]
        | df["is_top_cytolytic_score"]
        | df["is_top_antigen_response_score"]
    )
    df["is_candidate_tumor_reactive"] = df["is_base_selected"] & df["is_branch_any"]
    df["is_branch_union_within_base"] = df["is_candidate_tumor_reactive"]
    df["is_priority_cytotoxic_expander"] = (
        df["is_base_selected"] & df["is_increasing_and_became_cytotoxic"]
    )

    subsets = {
        "subset_non_viral": df[df["is_non_viral"]].copy(),
        "subset_base_cd8_nonviral": df[df["is_base_cd8_non_viral"]].copy(),
        "subset_min_cells": df[df["is_min_cells"]].copy(),
        "subset_base_selected": df[df["is_base_selected"]].copy(),
        "subset_immunogenic": df[df["is_immunogenic"]].copy(),
        "subset_top_immunogenic_any": df[df["is_top_immunogenic_any"]].copy(),
        "subset_top_cytolytic_score": df[df["is_top_cytolytic_score"]].copy(),
        "subset_top_antigen_response_score": df[df["is_top_antigen_response_score"]].copy(),
        "subset_enriched": df[df["is_enriched"]].copy(),
        "subset_increasing": df[df["is_increasing_nonzero"]].copy(),
        "subset_increasing_positive": df[df["is_increasing_all_timepoints"]].copy(),
        "subset_increasing_nonzero": df[df["is_increasing_nonzero"]].copy(),
        "subset_increasing_all_timepoints": df[df["is_increasing_all_timepoints"]].copy(),
        "subset_increasing_became_cytotoxic": df[df["is_increasing_and_became_cytotoxic"]].copy(),
        "subset_priority_cytotoxic_expander": df[df["is_priority_cytotoxic_expander"]].copy(),
        "subset_candidate_tumor_reactive": df[df["is_candidate_tumor_reactive"]].copy(),
    }
    for gene in immunogenic_genes:
        flag_col = f"is_top_immunogenic_{gene}"
        subsets[f"subset_top_immunogenic_{gene}"] = df[df[flag_col]].copy()

    terminal_freq_col = f"frequency_{last_tp}"
    sort_cols = [c for c in [terminal_freq_col, "mean_frequency", "total_cells"] if c in df.columns]
    if sort_cols:
        for key, subset_df in list(subsets.items()):
            if len(subset_df) == 0:
                continue
            subsets[key] = subset_df.sort_values(sort_cols, ascending=[False] * len(sort_cols))

    immunogenic_branch_rows = [(f"Imm.{gene}", count) for gene, count in immunogenic_branch_counts_within_base]
    immunogenic_independent_rows = []
    for gene in immunogenic_genes:
        flag_col = f"is_top_immunogenic_{gene}"
        immunogenic_independent_rows.append((f"Imm.{gene}", int(df[flag_col].sum())))

    cd4_stage_label = (
        "CD4<CD8" if np.isclose(max_cd4_to_cd8_ratio, 1.0) else f"CD4/CD8<{max_cd4_to_cd8_ratio:g}"
    )
    base_non_viral = df["is_non_viral"]
    base_cd8 = base_non_viral & df["is_cd8_positive"]
    base_cd4_ratio = base_cd8 & df["is_cd8_not_cd4"]
    base_min_cells = base_cd4_ratio & df["is_min_cells"]

    sequential_stages: list[tuple[str, int]] = []
    branch_union_stages = [
        ("All", len(df)),
        ("NonViral", int(base_non_viral.sum())),
        ("CD8", int(base_cd8.sum())),
        (cd4_stage_label, int(base_cd4_ratio.sum())),
        (f"Cells>={min_cells_per_clone}", int(base_min_cells.sum())),
        *immunogenic_branch_rows,
        ("Cytolytic", n_top_cytolytic),
        ("AgResp", n_top_antigen_response),
        (f"Inc>={increase_ratio_nonzero_min:g}x", int((df["is_base_selected"] & df["is_branch_enriched"]).sum())),
        ("Union", int(df["is_branch_union_within_base"].sum())),
        ("Final", int(df["is_candidate_tumor_reactive"].sum())),
    ]

    independent_stages = [
        ("All", len(df)),
        ("NonViral", int(df["is_non_viral"].sum())),
        ("CD8", int(df["is_cd8_positive"].sum())),
        (cd4_stage_label, int(df["is_cd8_not_cd4"].sum())),
        (f"Cells>={min_cells_per_clone}", int(df["is_min_cells"].sum())),
        ("Base", int(df["is_base_selected"].sum())),
        ("ImmAny", int(df["is_top_immunogenic_any"].sum())),
        ("Cytolytic", int(df["is_top_cytolytic_score"].sum())),
        ("AgResp", int(df["is_top_antigen_response_score"].sum())),
        *immunogenic_independent_rows,
        (f"Inc>={increase_ratio_nonzero_min:g}x", int(df["is_increasing_nonzero"].sum())),
        ("Inc+Cyto", int(df["is_increasing_and_became_cytotoxic"].sum())),
        ("Priority", int(df["is_priority_cytotoxic_expander"].sum())),
        ("Final", int(df["is_candidate_tumor_reactive"].sum())),
    ]

    return df, subsets, sequential_stages, independent_stages, branch_union_stages

run_til_select

run_til_select(args: Namespace) -> pd.DataFrame

Execute TIL-only clone selection workflow.

Parameters:

Name	Type	Description	Default
args	Namespace	Parsed CLI arguments from til-select command.	required

Returns:

Type	Description
DataFrame	Final selected master table with masks.

Source code in tcrsift/til_select.py

def run_til_select(args: argparse.Namespace) -> pd.DataFrame:
    """
    Execute TIL-only clone selection workflow.

    Parameters
    ----------
    args : argparse.Namespace
        Parsed CLI arguments from ``til-select`` command.

    Returns
    -------
    pd.DataFrame
        Final selected master table with masks.
    """
    if getattr(args, "config", None):
        discovered = parse_config(args.config)
    elif getattr(args, "samples", None):
        discovered = parse_sample_args(args.samples)
    else:
        discovered = default_timepoint_inputs(args.data_dir)
        if not discovered:
            raise TCRsiftValidationError(
                f"No timepoint inputs discovered under {args.data_dir}",
                hint="Provide --samples LABEL=CONSENSUS,CLONOTYPES or --config YAML.",
            )

    timepoint_order = list(discovered.keys())
    samples: dict[str, pd.DataFrame] = {}
    stats_by_timepoint: dict[str, pd.DataFrame] = {}
    for tp in timepoint_order:
        counts_df, stats_df = load_from_consensus(
            discovered[tp]["consensus"],
            discovered[tp]["clonotypes"],
            count_column=getattr(args, "count_column", None),
        )
        samples[tp] = counts_df
        stats_by_timepoint[tp] = stats_df

    marker_genes = [g.strip().upper() for g in str(args.marker_genes).split(",") if g.strip()]
    marker_genes = list(dict.fromkeys(marker_genes))
    if not marker_genes:
        marker_genes = MARKER_GENES_DEFAULT.copy()

    immunogenic_genes = [g.strip().upper() for g in str(args.immunogenic_genes).split(",") if g.strip()]
    cytotoxic_genes = [g.strip().upper() for g in str(args.cytotoxic_genes).split(",") if g.strip()]
    cytolytic_genes = [g.strip().upper() for g in str(args.cytolytic_genes).split(",") if g.strip()]
    antigen_response_genes = [
        g.strip().upper() for g in str(args.antigen_response_genes).split(",") if g.strip()
    ]

    fig_dir = Path(args.fig_dir)
    fig_dir.mkdir(parents=True, exist_ok=True)
    out_heatmap = Path(args.out_heatmap) if args.out_heatmap else (fig_dir / "abTCR_topk_heatmap.png")
    out_top = Path(args.out_top) if args.out_top else (fig_dir / "abTCR_topk.csv")
    out_annotated = Path(args.out_annotated) if args.out_annotated else (fig_dir / "abTCR_annotated.csv")
    out_annotated_heatmap = (
        Path(args.out_annotated_heatmap)
        if args.out_annotated_heatmap
        else (fig_dir / "abTCR_topk_annotated_heatmap.png")
    )
    out_selected_report = (
        Path(args.out_selected_report)
        if args.out_selected_report
        else (fig_dir / "selected_clones_report.pdf")
    )

    gex_inputs = discover_gex_inputs(args.data_dir, timepoint_order)
    marker_cell_by_timepoint: dict[str, pd.DataFrame] = {}
    marker_scores_by_timepoint: dict[str, pd.DataFrame] = {}
    for tp in timepoint_order:
        cell_df, score_df = compute_marker_scores_for_timepoint(
            consensus_path=discovered[tp]["consensus"],
            contig_csv_path=gex_inputs[tp]["contigs"],
            gex_h5_path=gex_inputs[tp]["gex_h5"],
            marker_genes=marker_genes,
        )
        marker_cell_by_timepoint[tp] = cell_df
        marker_scores_by_timepoint[tp] = score_df
        cell_df.to_csv(fig_dir / f"marker_cells_{tp}.csv", index=False)
        score_df.to_csv(fig_dir / f"marker_clonotype_scores_{tp}.csv", index=False)

    harmonized = build_harmonized_table(samples, timepoint_order, rank_by="mean_frequency")
    harmonized = add_marker_scores_to_harmonized(
        harmonized, marker_scores_by_timepoint, timepoint_order
    )

    annotated = harmonized.copy()
    database = load_databases(getattr(args, "vdjdb", None), getattr(args, "iedb", None), getattr(args, "cedar", None))
    if database is not None:
        annotated = match_clonotypes(harmonized, database, match_by=getattr(args, "match_by", "CDR3b_only"))
        if "db_species" in annotated.columns:
            annotated["db_species_group"] = annotated["db_species"].apply(normalize_species_label)
        annotated["db_vdjdb"] = annotated["db_database"].fillna("").str.contains("VDJdb")
        annotated["db_iedb"] = annotated["db_database"].fillna("").str.contains("IEDB")
        annotated["db_cedar"] = annotated["db_database"].fillna("").str.contains("CEDAR")
        plot_annotation_summary(annotated, fig_dir)

    metric_bases = infer_timepoint_metric_bases(annotated, timepoint_order)
    master_df = add_timepoint_summaries(annotated, timepoint_order, metric_bases)
    (
        master_df,
        subset_dfs,
        _sequential_stages,
        _independent_stages,
        branch_union_stages,
    ) = run_selection_pipeline(
        master_df,
        timepoint_order=timepoint_order,
        marker_genes=marker_genes,
        min_cells_per_clone=int(args.min_cells_per_clone),
        min_cd8_cp10k=float(args.min_cd8_cp10k),
        max_cd4_to_cd8_ratio=float(args.max_cd4_to_cd8_ratio),
        increase_ratio_nonzero_min=float(args.increase_ratio_nonzero_min),
        increase_ratio_all_timepoints_min=float(args.increase_ratio_all_timepoints_min),
        immunogenic_percentile=float(args.immunogenic_percentile),
        immunogenic_percentile_slack_frac=float(args.immunogenic_percentile_slack_frac),
        immunogenic_min_cp10k=float(args.immunogenic_min_cp10k),
        immunogenic_require_above_median=bool(args.immunogenic_require_above_median),
        cytotoxic_last_min_z=float(args.cytotoxic_last_min_z),
        cytotoxic_last_min_cp10k=float(args.cytotoxic_last_min_cp10k),
        became_cytotoxic_min_delta_z=float(args.became_cytotoxic_min_delta_z),
        immunogenic_genes_preferred=immunogenic_genes,
        cytotoxic_genes_preferred=cytotoxic_genes,
        cytolytic_genes_preferred=cytolytic_genes,
        antigen_response_genes_preferred=antigen_response_genes,
    )

    _write_subset_tables(subset_dfs, fig_dir)
    plot_selection_funnel(branch_union_stages, fig_dir / "selection_funnel.png", "Selection Funnel")

    plot_source = sort_harmonized_by_rank(master_df, args.rank_by)
    out_table = Path(args.out_table)
    out_table.parent.mkdir(parents=True, exist_ok=True)
    plot_source.to_csv(out_table, index=False)
    plot_source.to_csv(fig_dir / "abTCR_master_table.csv", index=False)
    plot_source.to_csv(out_annotated, index=False)

    mask_cols = sorted([c for c in plot_source.columns if c.startswith("is_")])
    mask_id_cols = [c for c in ["CDR3ab", "CDR3_alpha", "CDR3_beta"] if c in plot_source.columns]
    mask_metric_cols = [
        c
        for c in [
            *[f"cell_count_{tp}" for tp in timepoint_order],
            *[f"frequency_{tp}" for tp in timepoint_order],
            "cd8_cp10k_max",
            "cd4_cp10k_max",
            "cd4_to_cd8_ratio",
            "delta_frequency_last_vs_first",
            "fold_change_last_vs_first",
            "cytotoxic_score_z_delta_last_vs_first",
            "cytotoxic_score_cp10k_delta_last_vs_first",
        ]
        if c in plot_source.columns
    ]
    plot_source[mask_id_cols + mask_metric_cols + mask_cols].to_csv(
        fig_dir / "selection_masks.csv", index=False
    )

    count_cols = [f"cell_count_{tp}" for tp in timepoint_order if f"cell_count_{tp}" in plot_source.columns]
    freq_cols = [f"frequency_{tp}" for tp in timepoint_order if f"frequency_{tp}" in plot_source.columns]

    selected_mask = pd.Series(
        plot_source.get("is_candidate_tumor_reactive", pd.Series(False, index=plot_source.index)),
        index=plot_source.index,
    ).fillna(False).astype(bool)
    selected_for_report = plot_source[selected_mask].copy()
    create_selected_clone_pdf_report(
        selected_for_report,
        out_selected_report,
        timepoint_order,
        marker_genes,
        pyensembl_release=int(args.pyensembl_release),
    )

    primary_title = f"Top {args.top_k} abTCRs (ranked by {args.rank_by})"
    plot_heatmap(plot_source, count_cols, out_heatmap, int(args.top_k), title=primary_title)
    plot_source.head(int(args.top_k)).to_csv(out_top, index=False)

    for rank_metric in RANK_METRICS:
        ranked = sort_harmonized_by_rank(plot_source, rank_metric)
        ranked.head(int(args.top_k)).to_csv(fig_dir / f"abTCR_topk_{rank_metric}.csv", index=False)
        plot_heatmap(
            ranked,
            count_cols,
            fig_dir / f"abTCR_topk_heatmap_{rank_metric}.png",
            int(args.top_k),
            title=f"Top {args.top_k} abTCRs (ranked by {rank_metric})",
        )

    # Increasing + became-cytotoxic topk artifacts
    priority_mask = pd.Series(
        plot_source.get("is_priority_cytotoxic_expander", pd.Series(False, index=plot_source.index)),
        index=plot_source.index,
    ).fillna(False).astype(bool)
    priority_df = plot_source[priority_mask].copy()
    priority_top_path = fig_dir / "abTCR_increasing_became_cytotoxic_topk.csv"
    if priority_df.empty:
        plot_source.head(0).to_csv(priority_top_path, index=False)
    else:
        terminal_freq_col = f"frequency_{timepoint_order[-1]}"
        sort_cols = [
            c
            for c in [terminal_freq_col, "cytotoxic_score_z_delta_last_vs_first", "mean_frequency", "total_cells"]
            if c in priority_df.columns
        ]
        priority_top = priority_df.sort_values(sort_cols, ascending=[False] * len(sort_cols)).head(int(args.top_k))
        priority_top.to_csv(priority_top_path, index=False)
        if len(freq_cols) >= 2:
            plot_frequency_heatmap(
                priority_top,
                freq_cols,
                fig_dir / "abTCR_increasing_became_cytotoxic_topk_heatmap.png",
                "Top increasing clones that became cytotoxic",
            )

    # Monotonic-style increasing outputs
    monotonic_top, monotonic_freq_cols = select_top_monotonic_increase_clones(
        plot_source,
        timepoint_order,
        int(args.top_k),
        increase_ratio_nonzero_min=float(args.increase_ratio_nonzero_min),
        increase_ratio_all_timepoints_min=float(args.increase_ratio_all_timepoints_min),
    )
    if monotonic_top.empty:
        plot_source.head(0).to_csv(fig_dir / "abTCR_monotonic_increase_topk.csv", index=False)
    else:
        monotonic_top.to_csv(fig_dir / "abTCR_monotonic_increase_topk.csv", index=False)
        plot_frequency_heatmap(
            monotonic_top,
            monotonic_freq_cols,
            fig_dir / "abTCR_monotonic_increase_topk_heatmap.png",
            "Top increasing clones",
        )
        plot_monotonic_increase_log_lines(
            monotonic_top,
            monotonic_freq_cols,
            timepoint_order,
            fig_dir / "abTCR_monotonic_increase_topk_log_lines.png",
        )

    monotonic_pos_top, monotonic_pos_freq_cols = select_top_monotonic_increase_clones(
        plot_source,
        timepoint_order,
        int(args.top_k),
        increase_ratio_nonzero_min=float(args.increase_ratio_nonzero_min),
        increase_ratio_all_timepoints_min=float(args.increase_ratio_all_timepoints_min),
        require_positive_all_timepoints=True,
    )
    monotonic_pos_top.to_csv(fig_dir / "abTCR_monotonic_increase_positive_topk.csv", index=False)
    if not monotonic_pos_top.empty:
        plot_frequency_heatmap(
            monotonic_pos_top,
            monotonic_pos_freq_cols,
            fig_dir / "abTCR_monotonic_increase_positive_topk_heatmap.png",
            "Top increasing clones (positive timepoints)",
        )
        plot_monotonic_increase_log_lines(
            monotonic_pos_top,
            monotonic_pos_freq_cols,
            timepoint_order,
            fig_dir / "abTCR_monotonic_increase_positive_topk_log_lines.png",
        )

    annotation_cols = [c for c in ["db_vdjdb", "db_iedb", "db_cedar", "is_viral"] if c in plot_source.columns]
    if annotation_cols:
        plot_annotated_heatmap(
            plot_source,
            count_cols,
            out_annotated_heatmap,
            int(args.top_k),
            annotation_cols,
            title=f"Top {args.top_k} abTCRs (ranked by {args.rank_by}, annotated)",
        )
        for rank_metric in RANK_METRICS:
            ranked_annotated = sort_harmonized_by_rank(plot_source, rank_metric)
            plot_annotated_heatmap(
                ranked_annotated,
                count_cols,
                fig_dir / f"abTCR_topk_annotated_heatmap_{rank_metric}.png",
                int(args.top_k),
                annotation_cols,
                title=f"Top {args.top_k} abTCRs (ranked by {rank_metric}, annotated)",
            )

    # Summary CSV + additional expected artifact filenames.
    summary_rows = []
    sets_by_timepoint: dict[str, set[str]] = {}
    for label in timepoint_order:
        df = samples[label]
        n_clonotypes = len(df)
        total_cells = int(pd.to_numeric(df["cell_count"], errors="coerce").fillna(0).sum())
        summary_rows.append({"timepoint": label, "n_clonotypes": n_clonotypes, "total_cells": total_cells})
        sets_by_timepoint[label] = set(df["CDR3ab"].dropna())
    summary_df = pd.DataFrame(summary_rows)
    summary_df["timepoint"] = pd.Categorical(summary_df["timepoint"], categories=timepoint_order, ordered=True)
    summary_df = summary_df.sort_values("timepoint")
    summary_df.to_csv(fig_dir / "timepoint_summary.csv", index=False)

    # Lightweight placeholders for legacy plot artifacts.
    _save_placeholder_plot(fig_dir / "selection_scatter_panels.png", "Selection scatter panels")
    _save_placeholder_plot(fig_dir / "marker_gene_histograms_cp10k_mean.png", "Marker gene histograms")
    for i in range(1, 6):
        _save_placeholder_plot(
            fig_dir / f"marker_gene_pair_scatter_cp10k_mean_page{i}.png",
            f"Marker gene pair scatter page {i}",
        )
    _save_placeholder_plot(fig_dir / "tcr_timepoint_summary.png", "TCR summary")
    _save_placeholder_plot(fig_dir / "tcr_clone_size_bins.png", "Clone size bins")
    _save_placeholder_plot(fig_dir / "tcr_clone_size_bins_log.png", "Clone size bins (log)")
    _save_placeholder_plot(fig_dir / "tcr_frequency_bins.png", "Frequency bins")
    _save_placeholder_plot(fig_dir / "tcr_frequency_bins_log.png", "Frequency bins (log)")
    _save_placeholder_plot(fig_dir / "tcr_trend_categories.png", "Trend categories")
    _save_placeholder_plot(fig_dir / "tcr_umis_distribution.png", "UMI distribution")
    _save_placeholder_plot(fig_dir / "tcr_reads_distribution.png", "Read distribution")
    if len(sets_by_timepoint) <= 3:
        _save_placeholder_plot(fig_dir / "tcr_overlap_venn.png", "Timepoint overlap")

    export_all_plots_pdf(fig_dir, fig_dir / "all_plots.pdf")

    logger.info(
        "til-select complete: %s clones total, %s final candidates",
        len(master_df),
        int(master_df["is_candidate_tumor_reactive"].sum()),
    )
    return master_df