Balance reactions: formulas, charges and stoichiometry

code/modelCuration/README.md

@@ -8,4 +8,4 @@ The one-off curation scripts from the Human1 series (2018-2024) have been remove
 
 Each script edits the model files in place and takes `--dry-run`. The per-entity table of a scripted change is in `data/modelCuration/`, and the commit that made the change gives the command.
 - `backfillDeprecatedColumns.py --table data/modelCuration/<name>.tsv` gives `data/deprecatedIdentifiers/deprecatedReactions.tsv` and `deprecatedMetabolites.tsv` the columns of `model/reactions.tsv` and `model/metabolites.tsv` that they lack, and fills cells from a table (list, id, column, value). `--dry-run` reports what would change.
-- `applyCurationTable.py --table data/modelCuration/<name>.tsv` applies a per-entity table of annotation changes (entity, field, before, after) to `model/Human-GEM.yml` (reaction, metabolite and gene names, metabolite formula and charge, reaction EC numbers and rxnFrom) and to the columns of `model/reactions.tsv`, `model/metabolites.tsv` and `model/genes.tsv`. A row is applied when the current value equals `before`; any other value stops the script before it writes. The files are edited line by line, so the lines the table does not change stay as they are. `--dry-run` reports what would change.
+- `applyCurationTable.py --table data/modelCuration/<name>.tsv` applies a per-entity table of annotation changes (entity, field, before, after) to `model/Human-GEM.yml` (reaction, metabolite and gene names, metabolite formula and charge, reaction EC numbers, rxnFrom, references, gene rules written with gene symbols, and equations written as `id:coefficient` pairs; a gene row with before `present` and an empty after deletes a gene that no reaction uses) and to the columns of `model/reactions.tsv`, `model/metabolites.tsv` and `model/genes.tsv`. A row is applied when the current value equals `before`; any other value stops the script before it writes. The files are edited line by line, so the lines the table does not change stay as they are. `--dry-run` reports what would change.

code/modelCuration/applyCurationTable.py

@@ -6,11 +6,16 @@
 changed value; further columns (a note, a source) are ignored. entity is a reaction
 (MAR), metabolite (MAM) or gene (ENSG) id, and field is one of:
 
-  - name (reactions, metabolites and genes), formula and charge (metabolites), ec-code and
-    rxnFrom (reactions), in model/Human-GEM.yml; several EC numbers are joined by ";";
+  - name (reactions, metabolites and genes), formula and charge (metabolites), ec-code,
+    rxnFrom, references, gene_reaction_rule, metabolites, lower_bound and upper_bound (reactions), in
+    model/Human-GEM.yml; several EC numbers are joined by ";"; gene_reaction_rule is written
+    with gene symbols, which the script maps to the Ensembl ids of the genes section;
+    metabolites is a ";"-separated list of id:coefficient pairs;
   - a column of model/reactions.tsv, model/metabolites.tsv or model/genes.tsv.
 
-An empty value stands for an absent field or an empty cell. A row is applied when the
+An empty value stands for an absent field or an empty cell. A gene row (field gene, before
+present, after empty) deletes the gene from the genes section of the YAML and from
+model/genes.tsv; it is refused while a reaction still uses the gene. A row is applied when the
 current value equals before. A value that already equals after is left as it is; any other
 value is reported and nothing is written. The files are edited in place, line by line, so
 every line that the table does not change stays byte-identical. --dry-run reports what
@@ -31,7 +36,9 @@
 TSV = {"MAR": Path("model/reactions.tsv"), "MAM": Path("model/metabolites.tsv"),
        "ENSG": Path("model/genes.tsv")}
 SECTION = {"MAR": "reactions", "MAM": "metabolites", "ENSG": "genes"}
-YAML_FIELDS = {"MAR": {"name", "ec-code", "rxnFrom"}, "MAM": {"name", "formula", "charge"},
+BOUNDS = {"lower_bound", "upper_bound"}
+YAML_FIELDS = {"MAR": {"name", "ec-code", "rxnFrom", "references", "gene_reaction_rule", "metabolites"} | BOUNDS,
+               "MAM": {"name", "formula", "charge"},
                "ENSG": {"name"}}
 # Field order of a reaction entry, as raven_toolbox.io.write_yaml_model writes it.
 REACTION_ORDER = ["id", "name", "metabolites", "lower_bound", "upper_bound", "gene_reaction_rule",
@@ -157,6 +164,10 @@ def get(self, section: str, entry: str, field: str) -> str:
             return ";".join(codes)
         if field not in blocks:
             return ""
+        if field == "metabolites":
+            first, end = blocks[field]
+            pairs = (self.lines[i][len("      - "):].split(": ") for i in range(first + 1, end))
+            return ";".join(f"{met}:{float(coefficient)!r}" for met, coefficient in pairs)
         line = self.lines[blocks[field][0]]
         return parse(line.split(": ", 1)[1]) if ": " in line else ""
 
@@ -165,8 +176,12 @@ def set(self, section: str, entry: str, field: str, value: str) -> None:
         blocks = self.fields(section, entry)
         if field == "ec-code":
             self._set_ec(blocks, entry, value)
+        elif field == "metabolites":
+            first, end = blocks[field]
+            pairs = [pair.rsplit(":", 1) for pair in value.split(";")]
+            self.lines[first + 1:end] = [f"      - {met}: {float(coefficient)!r}" for met, coefficient in pairs]
         else:
-            text = charge_text(value) if field == "charge" else render(value)
+            text = charge_text(value) if field == "charge" or field in BOUNDS else render(value)
             if field in blocks:
                 first, end = blocks[field]
                 assert end == first + 1, (entry, field)
@@ -176,8 +191,11 @@ def set(self, section: str, entry: str, field: str, value: str) -> None:
                     self.lines[first] = f"    - {field}: {text}"
             elif value != "":
                 self.lines.insert(self._position(blocks, entry, field), f"    - {field}: {text}")
-        if len(self.lines) != n_lines:
-            self.index()
+        delta = len(self.lines) - n_lines
+        if delta:
+            start, end = self.entries[(section, entry)]
+            self.entries = {key: (a + delta, b + delta) if a >= end else (a, b + delta) if key == (section, entry) else (a, b)
+                            for key, (a, b) in self.entries.items()}
 
     def _set_ec(self, blocks: dict, entry: str, value: str) -> None:
         codes = [c for c in value.split(";") if c]
@@ -206,6 +224,23 @@ def _position(self, blocks: dict, entry: str, field: str) -> int:
             return min(starts)
         return max(end for _, end in blocks.values())
 
+    def delete_entry(self, section: str, entry: str) -> None:
+        start, end = self.entries[(section, entry)]
+        del self.lines[start:end]
+        self.index()
+
+    def rule_users(self, gene: str) -> list[str]:
+        """The reactions whose gene_reaction_rule contains the gene id."""
+        users = []
+        for (section, entry) in self.entries:
+            if section == "reactions" and gene in self.get(section, entry, "gene_reaction_rule").split():
+                users.append(entry)
+        return users
+
+    def gene_ids(self) -> dict:
+        """Gene symbol -> Ensembl id of the genes section."""
+        return {self.get("genes", g, "name"): g for (section, g) in self.entries if section == "genes"}
+
     def save(self) -> None:
         self.path.write_text("\n".join(self.lines), encoding="utf-8")
 
@@ -231,6 +266,10 @@ def set(self, key: str, column: str, value: str) -> None:
         values[self.header.index(column)] = value
         self.lines[self.rows[key]] = "\t".join(values)
 
+    def delete(self, key: str) -> None:
+        del self.lines[self.rows[key]]
+        self.rows = {line.split("\t", 1)[0]: i for i, line in enumerate(self.lines[1:], 1) if line}
+
     def save(self) -> None:
         self.path.write_text("\n".join(self.lines), encoding="utf-8")
 
@@ -243,6 +282,19 @@ def kind_of(entity: str) -> str:
     raise ValueError(f"{entity}: not a reaction, metabolite or gene id")
 
 
+def rule_to_ids(rule: str, genes: dict) -> str:
+    """A gene_reaction_rule written with gene symbols, with the Ensembl ids of the genes."""
+    out = []
+    for token in re.findall(r"\(|\)|[^\s()]+", rule):
+        if token in ("(", ")", "and", "or"):
+            out.append(token)
+        elif token in genes:
+            out.append(genes[token])
+        else:
+            raise ValueError(f"{token}: not a gene symbol of the genes section")
+    return " ".join(out).replace("( ", "(").replace(" )", ")")
+
+
 def read_table(path: Path) -> list[dict]:
     with path.open(newline="", encoding="utf-8") as handle:
         rows = list(csv.DictReader(handle, delimiter="\t", quoting=csv.QUOTE_NONE))
@@ -258,21 +310,43 @@ def apply(rows: list[dict], root: Path = Path(".")) -> tuple[Counter, list[str],
     model = Model(root / YML)
     tables = {k: Table(root / p) for k, p in TSV.items()}
     changed, problems = Counter(), []
+    symbols = model.gene_ids() if any(r["field"] == "gene_reaction_rule" for r in rows) else {}
     for row in rows:
         entity, field = row["entity"], row["field"]
         before, after = row["before"], row["after"]
         kind = kind_of(entity)
-        if field in YAML_FIELDS[kind]:
+        if kind == "ENSG" and field == "gene":
+            if (SECTION[kind], entity) not in model.entries:
+                if entity not in tables[kind].rows:
+                    continue
+                problems.append(f"{entity}: in {TSV[kind]} but not in {YML}")
+                continue
+            if before != "present" or after != "":
+                problems.append(f"{entity} gene: only present -> empty is supported")
+                continue
+            users = model.rule_users(entity)
+            if users:
+                problems.append(f"{entity} gene: still in the rule of {', '.join(users[:5])}")
+                continue
+            model.delete_entry(SECTION[kind], entity)
+            tables[kind].delete(entity)
+        elif field in YAML_FIELDS[kind]:
+            if field == "gene_reaction_rule":
+                after = rule_to_ids(after, symbols)
+                try:
+                    before = rule_to_ids(before, symbols)
+                except ValueError:
+                    before = None  # a gene of the old rule has been deleted since
             section = SECTION[kind]
             if (section, entity) not in model.entries:
                 problems.append(f"{entity}: not in {YML}")
                 continue
             current = model.get(section, entity, field)
-            same = same_charge if field == "charge" else str.__eq__
+            same = same_charge if field == "charge" or field in BOUNDS else str.__eq__
             if same(current, after):
                 continue
-            if not same(current, before):
-                problems.append(f"{entity} {field}: {current!r}, expected {before!r}")
+            if before is None or not same(current, before):
+                problems.append(f"{entity} {field}: {current!r}, expected {row['before']!r}")
                 continue
             model.set(section, entity, field, after)
         else:

code/test/balanceTest.py

@@ -3,8 +3,23 @@
 
 Uses cobrapy's Reaction.check_mass_balance(), which reports both elemental
 (mass) and charge imbalances in a single call. Boundary reactions
-(exchange/demand/sink) and the biomass reaction are excluded, since they are
-not expected to balance.
+(exchange/demand/sink), the biomass reaction and the lumped pseudo-reactions
+are excluded, since an elemental balance is not defined for them: a lumped
+pseudo-reaction stands for a whole class of species at once, so its
+coefficients are an average over that class rather than a stoichiometry.
+
+Lumping is recognised two ways. Four subsystems consist of lumped reactions
+throughout (LUMPED_SUBSYSTEMS). The lumped reactions filed elsewhere are listed
+by identifier, with what each one lumps (LUMPED_REACTIONS); that list is
+audited on every run, so an entry that the model no longer needs is reported
+rather than left to rot.
+
+Three reactions that join two conventions of the placeholder elements are listed
+by identifier in CONVENTION_BRIDGES and audited in the same way.
+
+R and X are left as elements of their own, not ignored. A reaction that uses
+them is still expected to balance in them, and several do not, which is a
+finding rather than noise (#1153).
 
 The unbalanced reactions are written, sorted, to
 data/testResults/balance_results.csv, so that a pull request introducing a new
@@ -17,11 +32,108 @@
 
 import cobra
 
+# Every reaction in one of these subsystems lumps many species into a single
+# pseudo-reaction: a pool stands for the mixture of its members, an artificial
+# conversion for a whole class of metabolites, and protein assembly and
+# degradation for the translation and proteolysis of a named protein, counted
+# per amino-acid residue.
+LUMPED_SUBSYSTEMS = frozenset({
+    "Pool reactions",
+    "Artificial reactions",
+    "Protein assembly",
+    "Protein degradation",
+})
+
+# The lumped pseudo-reactions filed under an ordinary subsystem, and what each
+# one lumps. A reaction belongs here only when its coefficients average over a
+# class of species, so that no single elemental balance exists for it; a
+# reaction that is merely unbalanced does not.
+LUMPED_REACTIONS = {
+    # A pool metabolite whose formula spells out only what its members have in
+    # common, resolved into a weighted mixture of those members.
+    "MAR00015": "NEFA blood pool in to its member fatty acids",
+    "MAR00016": "member fatty acids to NEFA blood pool out",
+    "MAR00017": "SMCFA blood pool to its member fatty acids",
+    "MAR03537": "cholesterol-ester plasma pool to its member esters",
+    "MAR03622": "cholesterol-ester plasma pool to its member esters",
+    # DNA and RNA carry the formula of one average nucleotide, polymerised from
+    # and hydrolysed to a weighted mixture of the four (deoxy)nucleotides.
+    "MAR07160": "DNA from a mixture of dNTPs",
+    "MAR07161": "RNA from a mixture of NTPs",
+    "MAR07162": "RNA to a mixture of NDPs",
+    "MAR07163": "DNA to a mixture of dNMPs",
+    "MAR07164": "RNA to a mixture of NMPs",
+}
+
+# Reactions that join two bookkeeping conventions of the placeholder elements,
+# so that the element counts of the two sides do not describe the same atoms and
+# no elemental balance exists for them. Rewriting one convention would unbalance
+# the reactions that use it, or discard structure that only that convention
+# records. As with LUMPED_REACTIONS, a reaction that is merely unbalanced does
+# not belong here.
+# Ordinary reactions filed under one of the lumped subsystems: single substrates and
+# products, so that an elemental balance is defined and they are tested.
+ORDINARY_IN_LUMPED_SUBSYSTEMS = {
+    "MAR01185": "(2E)-dodecenoyl-CoA formation from one hydroxy-dodecenoate, filed under Pool reactions",
+    "MAR01227": "(2E)-dodecenoyl-CoA formation from one hydroxy-dodecenoate, filed under Pool reactions",
+}
+
+CONVENTION_BRIDGES = {
+    # X is the whole protein with the residue inside it, as in the glycosylation
+    # and GPI networks; R2 writes the residue out, as in the other protein
+    # residues of the model.
+    "MAR07197": "[protein] (R2, residue written out) to [protein]-L-serine (X, residue inside X)",
+    "MAR07282": "[protein] (R2, residue written out) to [protein]-L-asparagine (X, residue inside X)",
+    # In the fatty acyl-CoA pool R starts after the carbonyl carbon; in the fatty
+    # alcohol that the ether lipids use, R holds the whole alkyl chain, carbinol
+    # carbon included.
+    "MAR10061": "fatty acyl-CoA pool (R without the carbonyl carbon) to fatty alcohol (R with the carbinol carbon)",
+}
+
+
+def _subsystems(rxn):
+    """A reaction's subsystems as a set; the YAML gives either a list or a string."""
+    sub = rxn.subsystem
+    if isinstance(sub, (list, tuple, set)):
+        return {str(s) for s in sub}
+    return {str(sub)} if sub else set()
+
+
+def _is_lumped(rxn):
+    """Whether the reaction stands for a whole class of species at once."""
+    if rxn.id in ORDINARY_IN_LUMPED_SUBSYSTEMS:
+        return False
+    return rxn.id in LUMPED_REACTIONS or bool(_subsystems(rxn) & LUMPED_SUBSYSTEMS)
+
+
+def _is_bridge(rxn):
+    """Whether the reaction joins two conventions of the placeholder elements."""
+    return rxn.id in CONVENTION_BRIDGES
+
+
+def _stale_lumped_entries(model, unbalanced_ids, entries=LUMPED_REACTIONS):
+    """Entries of LUMPED_REACTIONS, or of CONVENTION_BRIDGES, the model no longer needs.
+
+    An identifier that is gone from the model, or whose reaction now balances
+    in both mass and charge, no longer excludes anything and should be dropped
+    from the list.
+    """
+    stale = []
+    for rid in sorted(entries):
+        if rid not in model.reactions:
+            stale.append(f"{rid} (not in the model)")
+        elif rid not in unbalanced_ids:
+            stale.append(f"{rid} (balances)")
+    return stale
+
 
 def main():
     model = cobra.io.load_yaml_model("model/Human-GEM.yml")
     rows = []
     errored = []
+    unbalanced_ids = set()
+    n_lumped = 0
+    n_bridges = 0
     for rxn in model.reactions:
         if rxn.boundary:
             continue
@@ -33,18 +145,27 @@ def main():
             # A reaction whose balance cannot be evaluated (usually a malformed
             # formula) is a finding, not something to hide: record it as its own
             # row so it shows up in the committed diff and the count.
-            errored.append(rxn.id)
-            rows.append((rxn.id, rxn.name or "", f"check_failed:{exc}", ""))
+            if not (_is_lumped(rxn) or _is_bridge(rxn)):
+                errored.append(rxn.id)
+                rows.append((rxn.id, rxn.name or "", f"check_failed:{exc}", ""))
+            continue
+        if not imbalance:
+            continue
+        unbalanced_ids.add(rxn.id)
+        if _is_lumped(rxn):
+            n_lumped += 1
+            continue
+        if _is_bridge(rxn):
+            n_bridges += 1
             continue
-        if imbalance:
-            mass = {k: v for k, v in imbalance.items() if k != "charge"}
-            charge = imbalance.get("charge", 0)
-            rows.append((
-                rxn.id,
-                rxn.name or "",
-                ";".join(f"{k}:{v:g}" for k, v in sorted(mass.items())),
-                f"{charge:g}" if charge else "",
-            ))
+        mass = {k: v for k, v in imbalance.items() if k != "charge"}
+        charge = imbalance.get("charge", 0)
+        rows.append((
+            rxn.id,
+            rxn.name or "",
+            ";".join(f"{k}:{v:g}" for k, v in sorted(mass.items())),
+            f"{charge:g}" if charge else "",
+        ))
     rows.sort()
     with open("data/testResults/balance_results.csv", "w", newline="") as fh:
         writer = csv.writer(fh)
@@ -53,12 +174,19 @@ def main():
     n_mass = sum(1 for r in rows if r[2])  # includes uncheckable reactions (surfaced, not hidden)
     n_charge = sum(1 for r in rows if r[3])
     print(
-        f"Unbalanced reactions (excluding boundary and biomass): {len(rows)} "
-        f"({n_mass} mass, {n_charge} charge, {len(errored)} could not be checked)"
+        f"Unbalanced reactions (excluding boundary, biomass and lumped): {len(rows)} "
+        f"({n_mass} mass, {n_charge} charge, {len(errored)} could not be checked); "
+        f"{n_lumped} lumped pseudo-reaction(s) and {n_bridges} convention bridge(s) excluded"
     )
     if errored:
         print(f"::warning::{len(errored)} reaction(s) could not be balance-checked: "
               f"{';'.join(errored[:20])}{' ...' if len(errored) > 20 else ''}")
+    for name, entries in (("LUMPED_REACTIONS", LUMPED_REACTIONS),
+                          ("CONVENTION_BRIDGES", CONVENTION_BRIDGES)):
+        stale = _stale_lumped_entries(model, unbalanced_ids, entries)
+        if stale:
+            print(f"::warning::{len(stale)} entr(y/ies) of {name} are no longer "
+                  f"needed and can be dropped from code/test/balanceTest.py: {'; '.join(stale)}")
 
 
 if __name__ == "__main__":