adesomics.mapping.connection

 1import ast
 2import heapq
 3import numpy as np
 4from cobra.util.solver import linear_reaction_coefficients
 5
 6from adesomics.io.gsmm import Model
 7from adesomics.io.count_matrix import CountMatrix
 8
 9
10class Connector:
11    """
12    Binds an expression matrix to a GSMM via GPR rules. -> 
13    integration with E-Flux
14    """
15
16    def __init__(self, model: Model, counts: CountMatrix, max_bound: float = 1000.0):
17        self.model = model
18        self.counts = counts
19        self.tpm = counts.tpm()
20        self.max_bound = max_bound
21        self.gene_mapping = self._map_genes()
22
23    def _map_genes(self) -> dict:
24        model_ids = {g.id for g in self.model.model.genes}
25        count_ids = set(self.counts.counts.index)
26        mapping = {gid: gid for gid in count_ids if gid in model_ids}
27        print(f"{len(mapping)} / {len(model_ids)} model genes mapped to count matrix.")
28        return mapping
29
30    def _gpr_score(self, node, expr):
31        """Missing genes are unknown (NaN), not absent (0)."""
32        NA = float("nan")
33        if isinstance(node, ast.Expression): #usually not accessed
34            return self._gpr_score(node.body, expr)
35        if isinstance(node, ast.Name):
36            return expr.get(node.id, NA)
37        if isinstance(node, ast.BoolOp):
38            vals = [self._gpr_score(v, expr) for v in node.values]
39            known = [v for v in vals if not np.isnan(v)]
40            if not known:
41                return NA
42            return min(known) if isinstance(node.op, ast.And) else sum(known)
43        return NA
44
45    def _reaction_expression(self, sample):
46        expr = self.tpm[sample].to_dict()
47        scores = {}
48        for rxn in self.model.model.reactions:
49            rule = rxn.gene_reaction_rule.strip()   
50            if not rule:                            
51                continue
52            tree = ast.parse(rule, mode="eval").body
53            scores[rxn.id] = self._gpr_score(tree, expr)
54        return scores
55
56    def _compute_scale(self, samples, percentile: float = 95) -> float:
57        pooled = []
58        for s in samples:
59            pooled.extend(v for v in self._reaction_expression(s).values()
60                        if not np.isnan(v))
61        if not pooled:
62            raise ValueError("No reaction could be scored — check gene ID mapping")
63        return float(np.percentile(pooled, percentile))
64
65    def eflux(self, sample, scale=None):
66        if not np.isfinite(scale) or scale <= 0:
67            raise ValueError(f"Invalid scale: {scale}")
68        if scale is None:
69            scale = self._compute_scale(self.tpm.columns)
70        if scale == 0:
71            raise ValueError("Expression scale is zero -> cannot apply E-Flux because all reactions would be unconstrained. Check gene ID mapping/overlap.")
72        scores = {k: v for k, v in self._reaction_expression(sample).items()
73                if not np.isnan(v)}                 # NaN -> unconstrained
74
75        m = self.model.model.copy()
76        for rxn in m.reactions:
77            if rxn.id not in scores:
78                continue                              
79            ratio = min(scores[rxn.id] / scale, 1.0)  
80            lb, ub = rxn.lower_bound, rxn.upper_bound
81            rxn.bounds = (min(lb, 0) * ratio if lb < 0 else lb,
82                        max(ub, 0) * ratio if ub > 0 else ub)
83        return m.optimize()
84
85    def compare(self, samples, percentile=95):
86        scale = self._compute_scale(samples, percentile)
87        eflux_results = {}
88        for sample in samples:
89            if sample not in self.tpm.columns:
90                raise ValueError(f"Sample '{sample}' not found in count matrix.")
91            eflux_results[sample] = self.eflux(sample, scale)
92        return eflux_results