fasta_lake.study

fasta_lake.study

Group accepted study peptides and assign LFQ features directly by peptide.

The greedy cover follows the reviewed publication procedure. A representative labels a group and explains its assigned peptides; it is not proof of their unique biological origin. Pre-search selection and search confidence are unchanged.

canonical_peptide

canonical_peptide(peptide: str) -> str

Remove Sage modifications and use the declared I/L equivalence.

Sage writes residue mods in brackets (PEPT[+15.9949]LDEK) and N-terminal mods as a bracket followed by - ([+42.0106]-PEPTLDEK); both go.

Source code in fasta_lake/study.py

def canonical_peptide(peptide: str) -> str:
    """Remove Sage modifications and use the declared I/L equivalence.

    Sage writes residue mods in brackets (``PEPT[+15.9949]LDEK``) and N-terminal
    mods as a bracket followed by ``-`` (``[+42.0106]-PEPTLDEK``); both go.
    """
    return _MOD.sub("", peptide).replace("-", "").replace("I", "L")

member_set_key

member_set_key(proteins: str) -> str

Preserve the full distinct accession membership, without source-tag stripping.

Source code in fasta_lake/study.py

def member_set_key(proteins: str) -> str:
    """Preserve the full distinct accession membership, without source-tag stripping."""
    return ";".join(sorted({p.strip() for p in proteins.split(";") if p.strip()}))

greedy_set_cover

greedy_set_cover(protein_peptides: dict[str, set[str]]) -> tuple[list, set[str]]

Choose study representatives by maximum uncovered peptide gain.

Parameters:

Name	Type	Description	Default
protein_peptides	dict[str, set[str]]	Observed canonical peptides for each compatible protein accession.	required

Returns:

Type	Description
tuple	Selection-ordered (accession, newly covered peptide count) pairs and the set of covered peptides. Equal gains resolve by lexical accession.

Notes

Gains are recalculated as evidence is covered. This is post-search reporting, distinct from pre-search razor's fixed evidence counts. Greedy cover does not guarantee the smallest possible protein set; inputs are not modified.

Source code in fasta_lake/study.py

def greedy_set_cover(protein_peptides: dict[str, set[str]]) -> tuple[list, set[str]]:
    """Choose study representatives by maximum uncovered peptide gain.

    Parameters
    ----------
    protein_peptides : dict[str, set[str]]
        Observed canonical peptides for each compatible protein accession.

    Returns
    -------
    tuple
        Selection-ordered (accession, newly covered peptide count) pairs and the
        set of covered peptides. Equal gains resolve by lexical accession.

    Notes
    -----
    Gains are recalculated as evidence is covered. This is post-search reporting,
    distinct from pre-search razor's fixed evidence counts. Greedy cover does
    not guarantee the smallest possible protein set; inputs are not modified.
    """
    heap = [(-len(peptides), protein) for protein, peptides in protein_peptides.items()]
    heapq.heapify(heap)
    covered, selected = set(), []
    while heap:
        negative_gain, protein = heapq.heappop(heap)
        gain = len(protein_peptides[protein] - covered)
        if not gain:
            continue
        if gain != -negative_gain:
            heapq.heappush(heap, (-gain, protein))
            continue
        selected.append((protein, gain))
        covered.update(protein_peptides[protein])
    return selected, covered

study_dictionary

study_dictionary(protein_peptides: dict[str, set[str]], searched: set[str]) -> dict

Assign accepted peptides to a fixed set of study representatives.

Parameters:

Name	Type	Description	Default
protein_peptides	dict[str, set[str]]	Accepted canonical peptides associated with observed target accessions.	required
searched	set[str]	Accessions in the union of FASTAs actually searched for the study.	required

Returns:

Type	Description
dict	Selected representatives, covered peptides, selection ranks, peptide_to_group assignments and unobserved searched accessions.

Raises:

Type	Description
ValueError	If an observed accession is outside the searched universe.

Notes

Each peptide belongs to the first selected representative explaining it. A protein may be compatible with peptides assigned to several groups. This function checks graph membership; group_searches additionally verifies peptide containment in the actual representative sequences.

Source code in fasta_lake/study.py

def study_dictionary(protein_peptides: dict[str, set[str]], searched: set[str]) -> dict:
    """Assign accepted peptides to a fixed set of study representatives.

    Parameters
    ----------
    protein_peptides : dict[str, set[str]]
        Accepted canonical peptides associated with observed target accessions.
    searched : set[str]
        Accessions in the union of FASTAs actually searched for the study.

    Returns
    -------
    dict
        Selected representatives, covered peptides, selection ranks,
        peptide_to_group assignments and unobserved searched accessions.

    Raises
    ------
    ValueError
        If an observed accession is outside the searched universe.

    Notes
    -----
    Each peptide belongs to the first selected representative explaining it.
    A protein may be compatible with peptides assigned to several groups.
    This function checks graph membership; group_searches additionally verifies
    peptide containment in the actual representative sequences.
    """
    outside = set(protein_peptides) - searched
    if outside:
        raise ValueError(f"{len(outside)} observed proteins are outside the searched universe")
    selected, covered = greedy_set_cover(protein_peptides)
    peptide_to_group = {}
    for protein, _ in selected:
        for peptide in sorted(protein_peptides[protein]):
            peptide_to_group.setdefault(peptide, protein)
    return dict(
        selected=selected,
        covered=covered,
        rank={protein: i for i, (protein, _) in enumerate(selected)},
        peptide_to_group=peptide_to_group,
        unobserved=searched - set(protein_peptides),
    )

group_searches

group_searches(search_root: str | Path, output: str | Path, q_threshold: float = 0.01, *, config_name: str = 'config.json', acquisition_names: list[str] | None = None) -> dict

Build study groups and sum quantities from completed acquisition searches.

Parameters:

Name	Type	Description	Default
search_root	str or Path	Directory containing one completed Sage search subdirectory per acquisition. Each contains a configuration, results.sage.tsv and single-acquisition lfq.tsv.	required
output	str or Path	New directory outside search_root. Existing paths are refused.	required
q_threshold	float	Inclusive target peptide-q threshold in (0, 1], applied within each search.	0.01
config_name	('config.json', 'results.json')	Filename recording the exact searched FASTA and one mzML path.	"config.json"
acquisition_names	list[str] or None	Explicit acquisition directory names. None uses all direct subdirectories; the final roster is sorted and must contain unique, nonempty names.	None

Returns:

Type	Description
dict	The exported summary: roster, grouping rule, counts, intensity accounting and checks. Tables, representatives, peptide evidence and input hashes are written beside summary.json.

Raises:

Type	Description
ValueError	For inconsistent acquisitions, malformed evidence, accession/sequence conflicts, invalid confidence settings or incompatible representatives.
FileExistsError	If the requested output already exists.

Notes

Greedy cover uses only accepted evidence and actually searched proteins. Positive finite LFQ features require an accepted peptide in that acquisition. Each feature is assigned once; missing cells stay blank. Local members and representative presence are retained separately. This adds no study-wide or protein-group FDR estimate. Invalid inputs leave no completed bundle.

Source code in fasta_lake/study.py

def group_searches(
    search_root: str | Path,
    output: str | Path,
    q_threshold: float = 0.01,
    *,
    config_name: str = "config.json",
    acquisition_names: list[str] | None = None,
) -> dict:
    """Build study groups and sum quantities from completed acquisition searches.

    Parameters
    ----------
    search_root : str or Path
        Directory containing one completed Sage search subdirectory per acquisition.
        Each contains a configuration, results.sage.tsv and single-acquisition lfq.tsv.
    output : str or Path
        New directory outside search_root. Existing paths are refused.
    q_threshold : float
        Inclusive target peptide-q threshold in (0, 1], applied within each search.
    config_name : {"config.json", "results.json"}
        Filename recording the exact searched FASTA and one mzML path.
    acquisition_names : list[str] or None
        Explicit acquisition directory names. None uses all direct subdirectories;
        the final roster is sorted and must contain unique, nonempty names.

    Returns
    -------
    dict
        The exported summary: roster, grouping rule, counts, intensity accounting
        and checks. Tables, representatives, peptide evidence and input hashes are
        written beside summary.json.

    Raises
    ------
    ValueError
        For inconsistent acquisitions, malformed evidence, accession/sequence
        conflicts, invalid confidence settings or incompatible representatives.
    FileExistsError
        If the requested output already exists.

    Notes
    -----
    Greedy cover uses only accepted evidence and actually searched proteins.
    Positive finite LFQ features require an accepted peptide in that acquisition.
    Each feature is assigned once; missing cells stay blank. Local members and
    representative presence are retained separately. This adds no study-wide
    or protein-group FDR estimate. Invalid inputs leave no completed bundle.
    """
    if not math.isfinite(q_threshold) or not 0 < q_threshold <= 1:
        raise ValueError("q_threshold must be in (0, 1]")
    if config_name not in {"config.json", "results.json"}:
        raise ValueError("config_name must be config.json or results.json")
    requested = Path(output).absolute()
    if requested.exists() or requested.is_symlink():
        raise FileExistsError(f"Output exists: {requested}; choose a new directory")
    root, out = Path(search_root).resolve(), requested.resolve()
    if root == out or root in out.parents:
        raise ValueError("Output must be outside the search input directory")
    if acquisition_names is None:
        samples = sorted(p for p in root.iterdir() if p.is_dir())
    else:
        if not acquisition_names or len(acquisition_names) != len(set(acquisition_names)):
            raise ValueError("Acquisition names must be nonempty and unique")
        if any(Path(name).name != name or name in {".", ".."} for name in acquisition_names):
            raise ValueError("Acquisition names must be directory names, not paths")
        samples = [root / name for name in sorted(acquisition_names)]
        if any(not sample.is_dir() for sample in samples):
            raise ValueError("A declared acquisition directory is missing")
    if not samples:
        raise ValueError("No acquisition directories found")
    graph, sequences, gates, searched_by_sample = defaultdict(set), {}, {}, {}
    registry = []
    acquisitions = {}
    for sample in samples:
        config = sample / config_name
        raw = config.read_bytes()
        registry.append({"path": str(config), "sha256": hashlib.sha256(raw).hexdigest()})
        cfg = json.loads(raw)
        paths = cfg.get("mzml_paths")
        if not isinstance(paths, list) or len(paths) != 1 or not isinstance(paths[0], str):
            raise ValueError(f"{sample}: config must record exactly one mzml_paths entry")
        acquisition = _sage_filename(paths[0])
        if not acquisition:
            raise ValueError(f"{sample}: empty acquisition identity")
        acquisitions[sample.name] = acquisition
        fasta = Path(cfg["database"]["fasta"])
        if not fasta.is_absolute():
            raise ValueError(f"{sample}: config must record an absolute searched FASTA path")
        searched = _read_fasta(fasta, sequences, registry)
        searched_by_sample[sample.name] = searched
        gate = set()
        for line, row, _ in _rows(
            sample / "results.sage.tsv",
            {"peptide", "proteins", "label", "peptide_q", "filename", "rank"},
            registry=registry,
            acquisition=acquisition,
        ):
            if not _rank_one(row, sample / "results.sage.tsv", line):
                continue
            if row["label"] != "1" or not _passing(row["peptide_q"], q_threshold):
                continue
            peptide = canonical_peptide(row["peptide"])
            if not peptide:
                continue
            members = _members(row["proteins"])
            if members - searched:
                raise ValueError(f"{sample.name}: accepted proteins are outside its searched FASTA")
            gate.add(peptide)
            for protein in members:
                graph[protein].add(peptide)
        gates[sample.name] = gate
    dictionary = study_dictionary(graph, set(sequences))
    protein_counts = Counter(p for peptides in graph.values() for p in peptides)
    if not dictionary["covered"]:
        raise ValueError("No target peptide evidence passed the declared threshold")
    for peptide, group in dictionary["peptide_to_group"].items():
        if peptide not in sequences[group].replace("I", "L"):
            raise ValueError(f"Representative {group} does not contain assigned peptide {peptide}")

    # Build in a private temporary directory. Invalid input leaves no finished output.
    out.parent.mkdir(parents=True, exist_ok=True)
    with tempfile.TemporaryDirectory(prefix=f".{out.name}-", dir=out.parent) as temporary:
        temp = Path(temporary)
        summary = _quantify(
            samples, dictionary, gates, searched_by_sample, temp, registry, acquisitions
        )
        _peptide_evidence(samples, dictionary, gates, protein_counts, temp, registry, q_threshold)
        summary.update(
            schema_version=1,
            acquisition_names=[sample.name for sample in samples],
            assignment_rule="direct canonical peptide to first selected explaining representative",
            grouping_rule="greedy uncovered peptide cover; lexical accession gain ties",
            inferred_study_groups=len(dictionary["selected"]),
            searched_universe=len(sequences),
            proteins_with_observed_evidence=len(graph),
            unobserved_searched_proteins=len(dictionary["unobserved"]),
            pooled_canonical_peptides=len(dictionary["covered"]),
            representative_sequence_compatibility_checked=True,
            acquisition_identity_checked=True,
            identification_rank_policy="require explicit positive integer rank; accept rank 1 only",
            lfq_feature_policy="unique modified peptide and charge within each acquisition",
            q_threshold=q_threshold,
            config_name=config_name,
            confidence="post-search reporting; no new protein-group or study-wide FDR estimate",
        )
        assigned = Counter(dictionary["peptide_to_group"].values())
        _write_table(
            temp / "study_groups.tsv",
            [
                "study_group",
                "selection_rank",
                "observed_support_peptides",
                "assigned_peptides",
            ],
            ((p, i, len(graph[p]), assigned[p]) for i, (p, _) in enumerate(dictionary["selected"])),
        )
        _write_table(
            temp / "peptide_to_study_group.tsv",
            ["canonical_peptide", "study_group"],
            sorted(dictionary["peptide_to_group"].items()),
        )
        # Export exactly the reporting representatives for optional annotation.
        # This does not collapse the retained peptide/member ambiguity ledgers.
        with (temp / "representatives.fasta").open("x") as stream:
            for protein, _ in dictionary["selected"]:
                stream.write(f">{protein}\n{sequences[protein]}\n")
        summary["annotation_targets"] = {
            "fasta": "representatives.fasta",
            "sha256": hashlib.sha256((temp / "representatives.fasta").read_bytes()).hexdigest(),
            "scope": "Selected reporting representatives; annotation does not prove unique origin",
        }
        module = Path(__file__)
        registry.append(
            {"path": str(module), "sha256": hashlib.sha256(module.read_bytes()).hexdigest()}
        )
        (temp / "input_manifest.json").write_text(json.dumps(registry, indent=2) + "\n")
        (temp / "summary.json").write_text(json.dumps(summary, indent=2) + "\n")
        # Exclusive creation also protects an output created concurrently after preflight.
        out.mkdir(exist_ok=False)
        for path in sorted(temp.iterdir(), key=lambda p: (p.name == "summary.json", p.name)):
            path.rename(out / path.name)
    return summary