Skip to content

fasta_lake.multi_omics.exact#

fasta_lake.multi_omics.exact #

Strict sequence and TPM inputs for optional post-selection molecular unions.

Protein identity follows FastaLake hash version 5.3; peptide I/L equivalence is deliberately absent from protein hashing. No cluster representative substitution.

sha256 #

sha256(path)

Hash file bytes in 8 MiB blocks and return the SHA-256 hex digest.

Source code in fasta_lake/multi_omics/exact.py 
20
21
22
23
24
25
26
def sha256(path):
    """Hash file bytes in 8 MiB blocks and return the SHA-256 hex digest."""
    digest = hashlib.sha256()
    with Path(path).open("rb") as stream:
        for block in iter(lambda: stream.read(8 << 20), b""):
            digest.update(block)
    return digest.hexdigest()

read_fasta #

read_fasta(path, *, max_records=None)

Return accession -> (full header, normalized sequence, SHA-256).

Source code in fasta_lake/multi_omics/exact.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def read_fasta(path, *, max_records=None):
    """Return accession -> (full header, normalized sequence, SHA-256)."""
    if max_records is not None and (type(max_records) is not int or max_records < 1):
        raise ValueError("max_records must be a positive integer")
    path = Path(path)
    opener = gzip.open if path.suffix == ".gz" else open
    result, header, parts = {}, None, []

    def save():
        """Store a normalized record with its exact-sequence hash, enforcing ID limits."""
        if header is None:
            return
        accession = header.split()[0] if header.split() else ""
        sequence = normalize_sequence("".join(parts)).decode("ascii")
        if not accession or not sequence or accession in result:
            raise ValueError(f"Empty/duplicate FASTA accession or empty sequence: {accession!r}")
        result[accession] = (header, sequence, hashlib.sha256(sequence.encode("ascii")).hexdigest())
        if max_records is not None and len(result) > max_records:
            raise ValueError(
                f"FASTA exceeds the {max_records} record limit; select a smaller subset"
            )

    with opener(path, "rt", encoding="utf-8") as stream:
        for line in stream:
            if line.startswith(">"):
                save()
                header, parts = line[1:].rstrip("\r\n"), []
            elif line.strip():
                if header is None:
                    raise ValueError("Sequence before first FASTA header")
                parts.append(line.strip())
        save()
    if not result:
        raise ValueError("Empty FASTA")
    return result

numeric #

numeric(value, *, nan_is_missing=False)

Parse a nonnegative finite TPM while preserving explicit missing tokens.

Blank, NA, N/A, null and none return None. NaN is missing only when requested. A missing TSV field, negative value or nonfinite number raises ValueError.

Source code in fasta_lake/multi_omics/exact.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def numeric(value, *, nan_is_missing=False):
    """Parse a nonnegative finite TPM while preserving explicit missing tokens.

    Blank, NA, N/A, null and none return None. NaN is missing only when requested.
    A missing TSV field, negative value or nonfinite number raises ValueError.
    """
    if value is None:
        raise ValueError("Malformed TSV row: missing field")
    if value.strip().lower() in {"", "na", "n/a", "null", "none"} or (
        nan_is_missing and value.strip().lower() == "nan"
    ):
        return None
    number = float(value)
    if not math.isfinite(number) or number < 0:
        raise ValueError(f"Invalid nonnegative finite TPM: {value!r}")
    return number

read_tpm #

read_tpm(path, *, nan_is_missing=False)

Read plain or gzipped gene TPM rows without converting missing values to zero.

Require unique prodigal_protein IDs and metaG_tpm/metaT_tpm columns. Return a gene-keyed dictionary with parsed dna/rna and original raw strings; malformed rows, duplicate genes and invalid numeric values raise ValueError.

Source code in fasta_lake/multi_omics/exact.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def read_tpm(path, *, nan_is_missing=False):
    """Read plain or gzipped gene TPM rows without converting missing values to zero.

    Require unique prodigal_protein IDs and metaG_tpm/metaT_tpm columns.
    Return a gene-keyed dictionary with parsed dna/rna and original raw strings;
    malformed rows, duplicate genes and invalid numeric values raise ValueError.
    """
    required = {"prodigal_protein", "metaG_tpm", "metaT_tpm"}
    rows = {}
    opener = gzip.open if Path(path).suffix == ".gz" else open
    with opener(path, "rt", newline="") as stream:
        reader = csv.DictReader(stream, delimiter="\t")
        fields = reader.fieldnames or []
        if len(fields) != len(set(fields)) or not required <= set(fields):
            raise ValueError("TPM requires unique prodigal_protein, metaG_tpm, metaT_tpm columns")
        for row in reader:
            gene = row["prodigal_protein"]
            if None in row or not gene or gene != gene.strip() or gene in rows:
                raise ValueError(f"Malformed or duplicate TPM gene: {gene!r}")
            rows[gene] = dict(
                dna=numeric(row["metaG_tpm"], nan_is_missing=nan_is_missing),
                rna=numeric(row["metaT_tpm"], nan_is_missing=nan_is_missing),
                raw_dna=row["metaG_tpm"],
                raw_rna=row["metaT_tpm"],
            )
    if not rows:
        raise ValueError("TPM table has no genes")
    return rows

top_cutoff #

top_cutoff(values, percent, rounding='floor')

Top percentage among positive GENES; zero disables, ties all included.

Source code in fasta_lake/multi_omics/exact.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def top_cutoff(values, percent, rounding="floor"):
    """Top percentage among positive GENES; zero disables, ties all included."""
    if not math.isfinite(percent) or not 0 <= percent <= 100:
        raise ValueError("DNA top percent must be finite and in [0, 100]")
    if rounding not in {"floor", "ceil"}:
        raise ValueError("Rank rounding must be floor or ceil")
    values = list(values)
    if any(v is not None and (not math.isfinite(v) or v < 0) for v in values):
        raise ValueError("Invalid TPM values")
    positive = sorted(v for v in values if v is not None and v > 0)
    if not positive or percent == 0:
        return None
    rank = max(1, (math.floor if rounding == "floor" else math.ceil)(len(positive) * percent / 100))
    return positive[-rank]

write_table #

write_table(path, fields, rows)

Write dictionary rows as a TSV, using deterministic gzip for a .gz path.

The supplied field order defines the columns. Existing destination files are overwritten; callers enforce new output directories where required.

Source code in fasta_lake/multi_omics/exact.py 
130
131
132
133
134
135
136
137
138
139
140
141
def write_table(path, fields, rows):
    """Write dictionary rows as a TSV, using deterministic gzip for a .gz path.

    The supplied field order defines the columns. Existing destination files
    are overwritten; callers enforce new output directories where required.
    """
    path = Path(path)
    opener = open_gzip_text if path.suffix == ".gz" else (lambda p: open(p, "w", newline=""))
    with opener(path) as stream:
        writer = csv.DictWriter(stream, fieldnames=fields, delimiter="\t", lineterminator="\n")
        writer.writeheader()
        writer.writerows(rows)

load_source #

load_source(manifest_path)

Load one explicitly specimen-qualified quantification reference.

Checksums and all-gene agreement detect changed/mismatched inputs. The provenance record must establish how the original quantification reference was linked to these proteins; the caller's specimen assertion is not inferred from filenames. A hash alone cannot establish biological specimen identity.

Source code in fasta_lake/multi_omics/exact.py 
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
def load_source(manifest_path):
    """Load one explicitly specimen-qualified quantification reference.

    Checksums and all-gene agreement detect changed/mismatched inputs. The
    provenance record must establish how the original quantification reference
    was linked to these proteins; the caller's specimen assertion is not inferred
    from filenames. A hash alone cannot establish biological specimen identity.
    """
    manifest_path = Path(manifest_path).resolve()
    obj = json.loads(manifest_path.read_text())
    if obj.get("schema") != "fastalake.molecular-source.v1":
        raise ValueError("Unsupported molecular source schema")
    for field in ("specimen", "reference_id"):
        if not isinstance(obj.get(field), str) or not obj[field].strip():
            raise ValueError(f"Missing {field}")
    inputs = {}
    for field in ("tpm", "proteins", "provenance"):
        record = obj[field]
        path = Path(record["path"])
        if not path.is_absolute():
            path = manifest_path.parent / path
        path = path.resolve(strict=True)
        if sha256(path) != record["sha256"]:
            raise ValueError(f"{field} checksum mismatch")
        inputs[field] = path
    if obj.get("reference_status") != "verified_original_quantification_reference":
        raise ValueError("Original quantification reference is not verified")
    genes = read_tpm(inputs["tpm"])
    proteins = read_fasta(inputs["proteins"])
    if set(genes) != set(proteins):
        missing = sorted(set(genes) - set(proteins))
        extra = sorted(set(proteins) - set(genes))
        raise ValueError(
            f"All-gene reference mismatch: missing={len(missing)}, extra={len(extra)}; "
            f"examples={missing[:5]}"
        )
    return obj, inputs, genes, proteins

build_union #

build_union(baseline, source_manifest, specimen, output, *, dna_percent=1.0, dna_absolute=None, rna_absolute=None, rounding='floor', rna_percent=0, background=None, metabolites=False)

Keep every baseline sequence and append independently selected molecular proteins.

Source code in fasta_lake/multi_omics/exact.py 
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
def build_union(
    baseline,
    source_manifest,
    specimen,
    output,
    *,
    dna_percent=1.0,
    dna_absolute=None,
    rna_absolute=None,
    rounding="floor",
    rna_percent=0,
    background=None,
    metabolites=False,
):
    """Keep every baseline sequence and append independently selected molecular proteins."""
    from .selection import build_selection

    return build_selection(
        source_manifest,
        specimen,
        output,
        baseline=baseline,
        background=background,
        dna_percent=dna_percent,
        rna_percent=rna_percent,
        dna_absolute=dna_absolute,
        rna_absolute=rna_absolute,
        rounding=rounding,
        metabolites=metabolites,
    )