Skip to content

fasta_lake.tune#

fasta_lake.tune #

FastaLake tune — parameter tuning from existing SAGE results.

Takes a directory of SAGE entrapment searches, applies a grid of filters post-hoc, and recommends parameters that achieve a user-specified target FDP.

No re-search required — operates on existing results.sage.tsv files.

CLI: fasta-lake tune --results-dir PATH --target-fdp 0.02

Output: tuning_report.csv — per-filter FDP, target count, pareto verdict recommended.json — the winning filter as a SAGE config diff

FilterSpec dataclass #

One parameter combination to test.

Source code in fasta_lake/tune.py 
36
37
38
39
40
41
42
43
@dataclass
class FilterSpec:
    """One parameter combination to test."""

    name: str
    q_threshold: float
    min_length: int
    min_score: float

default_filter_grid #

default_filter_grid() -> list[FilterSpec]

Standard grid. len × q × score cartesian product.

Source code in fasta_lake/tune.py 
46
47
48
49
50
51
52
53
54
def default_filter_grid() -> list[FilterSpec]:
    """Standard grid. len × q × score cartesian product."""
    grid = []
    for q in [0.01, 0.005, 0.001]:
        for L in [7, 8, 9, 10]:
            for s in [0.0, 0.5, 0.7]:
                name = f"q{q:.3f}_len{L}_score{s:.1f}"
                grid.append(FilterSpec(name, q, L, s))
    return grid

apply_filter #

apply_filter(rows: list[dict], spec: FilterSpec, entrap_prefix: str) -> dict

Count target / decoy / entrapment PSMs under a filter.

Source code in fasta_lake/tune.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def apply_filter(rows: list[dict], spec: FilterSpec, entrap_prefix: str) -> dict:
    """Count target / decoy / entrapment PSMs under a filter."""
    target = decoy = entrap = 0
    for r in rows:
        if r["q"] > spec.q_threshold:
            continue
        if r["len"] < spec.min_length:
            continue
        if r["score"] < spec.min_score:
            continue
        if r["is_decoy"]:
            decoy += 1
        elif entrap_prefix in r["proteins"]:
            entrap += 1
        else:
            target += 1
    return {"target": target, "decoy": decoy, "entrap": entrap}

parse_tsv #

parse_tsv(path: Path) -> list[dict]

Parse required Sage fields for the historical spectrum-q tuning diagnostic.

Return q, length, discriminant score, protein text and the legacy decoy flag for each row. Missing columns and invalid numeric fields raise errors.

Source code in fasta_lake/tune.py 
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
def parse_tsv(path: Path) -> list[dict]:
    """Parse required Sage fields for the historical spectrum-q tuning diagnostic.

    Return q, length, discriminant score, protein text and the legacy decoy flag
    for each row. Missing columns and invalid numeric fields raise errors.
    """
    rows = []
    with open(path) as f:
        r = csv.DictReader(f, delimiter="\t")
        require_columns(
            r.fieldnames, ("spectrum_q", "peptide_len", "sage_discriminant_score", "proteins"), path
        )
        for row in r:
            p = row["proteins"]
            rows.append(
                {
                    "q": float(row["spectrum_q"]),
                    "len": int(row["peptide_len"]),
                    "score": float(row["sage_discriminant_score"]),
                    "proteins": p,
                    "is_decoy": "rev_" in p,
                }
            )
    return rows

compute_fdp_combined #

compute_fdp_combined(entrap: int, target_plus_decoy: int, r: float) -> float | None

Noble 2025 combined-method FDP as percentage; None when undefined.

A sample with no accepted PSMs has no FDP, not an FDP of zero; returning 0.0 pulled the median that picks the winning filter towards zero.

Source code in fasta_lake/tune.py 
102
103
104
105
106
107
108
109
110
111
def compute_fdp_combined(entrap: int, target_plus_decoy: int, r: float) -> float | None:
    """Noble 2025 combined-method FDP as percentage; ``None`` when undefined.

    A sample with no accepted PSMs has no FDP, not an FDP of zero; returning
    0.0 pulled the median that picks the winning filter towards zero.
    """
    denom = target_plus_decoy + entrap
    if denom == 0:
        return None
    return entrap * (1 + 1 / r) / denom * 100

find_sample_results #

find_sample_results(results_dir: Path) -> Iterator[tuple[str, str, Path]]

Yield (sample_id, entrap_type, results.sage.tsv path) for each sample.

Source code in fasta_lake/tune.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
def find_sample_results(results_dir: Path) -> Iterator[tuple[str, str, Path]]:
    """Yield (sample_id, entrap_type, results.sage.tsv path) for each sample."""
    # Expected layout: {results_dir}/{sample_id}/{entrap_type}/results.sage.tsv
    # Where entrap_type in {entrapment, entrapment_plant, entrapment_archaea}
    subdir_map = {
        "entrapment": "shuffled",
        "entrapment_plant": "plant",
        "entrapment_archaea": "archaea",
    }
    for sample_dir in sorted(results_dir.iterdir()):
        if not sample_dir.is_dir():
            continue
        for subdir_name, entrap_type in subdir_map.items():
            tsv = sample_dir / subdir_name / "results.sage.tsv"
            if tsv.exists():
                yield sample_dir.name, entrap_type, tsv

tune #

tune(results_dir: Path, target_fdp_pct: float = 2.0, entrap_db_sizes: dict[str, int] | None = None, target_db_size_fn=None, grid: list[FilterSpec] | None = None) -> dict

Tune parameters by post-hoc filter sweep.

Parameters:

Name Type Description Default
results_dir Path

Directory with per-sample subdirs (each containing entrapment/, entrapment_plant/, entrapment_archaea/ with results.sage.tsv).

 required
target_fdp_pct float

Maximum acceptable median shuffled FDP (percent). Winner is the filter with highest target PSM count that satisfies this threshold.

 2.0
entrap_db_sizes dict[str, int] | None

Map entrap_type -> protein count. Defaults to FastaLake standard sizes.

 None
target_db_size_fn

Callable (sample_id) -> int returning the target FASTA protein count. If None, assumes 75000 (typical MP metaproteomics).

 None
grid list[FilterSpec] | None

List of FilterSpec to evaluate. Defaults to default_filter_grid().

 None

Returns:

Type Description
dict with keys: all_results (per-filter per-entrap summary), winner (best filter).
Source code in fasta_lake/tune.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
def tune(
    results_dir: Path,
    target_fdp_pct: float = 2.0,
    entrap_db_sizes: dict[str, int] | None = None,
    target_db_size_fn=None,
    grid: list[FilterSpec] | None = None,
) -> dict:
    """
    Tune parameters by post-hoc filter sweep.

    Parameters
    ----------
    results_dir
        Directory with per-sample subdirs (each containing entrapment/, entrapment_plant/,
        entrapment_archaea/ with results.sage.tsv).
    target_fdp_pct
        Maximum acceptable median shuffled FDP (percent). Winner is the filter
        with highest target PSM count that satisfies this threshold.
    entrap_db_sizes
        Map entrap_type -> protein count. Defaults to FastaLake standard sizes.
    target_db_size_fn
        Callable (sample_id) -> int returning the target FASTA protein count.
        If None, assumes 75000 (typical MP metaproteomics).
    grid
        List of FilterSpec to evaluate. Defaults to `default_filter_grid()`.

    Returns
    -------
    dict with keys: all_results (per-filter per-entrap summary), winner (best filter).
    """
    if grid is None:
        grid = default_filter_grid()
    if entrap_db_sizes is None:
        entrap_db_sizes = {"shuffled": 6500, "plant": 16418, "archaea": 505}
    if target_db_size_fn is None:

        def target_db_size_fn(_sid):
            """Supply the historical 75,000-protein fallback when no size callback exists."""
            return 75000

    # Collect per-sample rows
    per_sample: dict[tuple[str, str], list[dict]] = {}
    for sample_id, entrap_type, tsv in find_sample_results(results_dir):
        per_sample[(sample_id, entrap_type)] = parse_tsv(tsv)

    if not per_sample:
        raise RuntimeError(f"No entrapment result TSVs found under {results_dir}")

    # Evaluate each filter spec × each entrapment type across all samples
    summary = []
    for spec in grid:
        for entrap_type, entrap_db in entrap_db_sizes.items():
            prefix = ENTRAP_PREFIXES[entrap_type]
            target_counts = []
            fdps = []
            for (sid, etype), rows in per_sample.items():
                if etype != entrap_type:
                    continue
                target_db = target_db_size_fn(sid)
                r = entrap_db / target_db if target_db else 0.0
                if r <= 0:
                    continue
                c = apply_filter(rows, spec, prefix)
                fdp = compute_fdp_combined(c["entrap"], c["target"] + c["decoy"], r)
                if fdp is None:
                    continue  # undefined for this sample: not a zero
                target_counts.append(c["target"])
                fdps.append(fdp)

            if not fdps:
                continue
            summary.append(
                {
                    "filter": spec.name,
                    "q_threshold": spec.q_threshold,
                    "min_length": spec.min_length,
                    "min_score": spec.min_score,
                    "entrap_type": entrap_type,
                    "n_samples": len(fdps),
                    "median_target": statistics.median(target_counts),
                    "median_fdp_pct": statistics.median(fdps),
                    "mean_fdp_pct": statistics.mean(fdps),
                    "p90_fdp_pct": sorted(fdps)[int(0.9 * len(fdps))],
                }
            )

    # Pick winner: highest target count where shuffled median FDP ≤ target_fdp_pct
    shuffled_rows = [r for r in summary if r["entrap_type"] == "shuffled"]
    passing = [r for r in shuffled_rows if r["median_fdp_pct"] <= target_fdp_pct]

    if not passing:
        # No filter passes → recommend the tightest
        winner = min(shuffled_rows, key=lambda r: r["median_fdp_pct"]) if shuffled_rows else None
        winner_verdict = "NO_FILTER_PASSES — showing tightest available"
    else:
        winner = max(passing, key=lambda r: r["median_target"])
        winner_verdict = f"PASSES @ target_fdp_pct={target_fdp_pct}"

    return {
        "all_results": summary,
        "winner": winner,
        "winner_verdict": winner_verdict,
        "target_fdp_pct": target_fdp_pct,
        "n_samples": len(set(sid for sid, _ in per_sample.keys())),
        "n_entrap_types": len(set(et for _, et in per_sample.keys())),
    }

write_reports #

write_reports(report: dict, out_dir: Path) -> None

Write tuning CSV, recommended settings and supporting reports.

Create the directory as needed and overwrite report files. An empty result list produces no new tuning_report.csv.

Source code in fasta_lake/tune.py 
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
def write_reports(report: dict, out_dir: Path) -> None:
    """Write tuning CSV, recommended settings and supporting reports.

    Create the directory as needed and overwrite report files. An empty result
    list produces no new tuning_report.csv.
    """
    out_dir.mkdir(parents=True, exist_ok=True)

    # CSV report
    csv_path = out_dir / "tuning_report.csv"
    if report["all_results"]:
        cols = list(report["all_results"][0].keys())
        with open(csv_path, "w") as f:
            w = csv.DictWriter(f, fieldnames=cols)
            w.writeheader()
            w.writerows(report["all_results"])

    # Recommended preset / filter
    rec_path = out_dir / "recommended.json"
    recommended = {
        "winner_verdict": report["winner_verdict"],
        "target_fdp_pct": report["target_fdp_pct"],
        "winner": report["winner"],
    }
    rec_path.write_text(json.dumps(recommended, indent=2, default=float))