determine linear pathways with net stoichiometry

modelpruner/graph/linear_pathway.py

@@ -18,38 +18,135 @@ class LinearPathway:
         init_sid: initial pathway node
         ng: Class Network Graph to work on
         """
-
         self.init_sid = init_sid
         self.balance_c2_sids = set(balance_c2_sids)
         self.ng = network_graph
         self.net_reactants = {}
         self.net_products = {}
+        self.segments = []
+        self.c2_sids_pwy = []
+        self.inconsistant = False
 
-        # construct the linear pathway
-        self.segments = [self.get_initial_segment()]
-        self.add_downstream_segments()
-        self.add_upstream_segments()
+        # construct the main linear pathway
+        main_segments = self.get_lin_pwy_containing(init_sid)
+        self.segments.extend(main_segments)
 
-        # collect pathway data
-        self.sids = []
-        self.rids = [self.segments[0].producing.rid]
-        self.reversible = self.segments[0].reversible
-        scale = 1.0
-        direction = 1.0 if self.segments[0].producing.forward is True else -1.0
-        self.rids_scale = {self.segments[0].producing.rid: scale * direction}
-        for pwy_seg in self.segments:
-            self.sids.append(pwy_seg.sid)
-            self.rids.append(pwy_seg.consuming.rid)
-            scale *= pwy_seg.producing.stoic / pwy_seg.consuming.stoic
-            direction = 1.0 if pwy_seg.consuming.forward is True else -1.0
-            self.rids_scale[pwy_seg.consuming.rid] = scale * direction
+        main_data = self.collect_pwy_data(main_segments)
+        self.c2_sids_pwy.extend(main_data['sid_seq'])
+        self.rids_scale = main_data['rids_scale']
 
+        # get_net_stoichiometry for current sub_network
         net_stoic = self.get_net_stoichiometry()
         self.net_reactants = {sid: -stoic for sid, stoic in net_stoic.items() if stoic < 0.0}
         self.net_products = {sid: stoic for sid, stoic in net_stoic.items() if stoic > 0.0}
+        c2_reactants = set(self.net_reactants.keys()).intersection(self.balance_c2_sids)
+        c2_products = set(self.net_products.keys()).intersection(self.balance_c2_sids)
+
+        # remove from balance sids the consumed sids
+        self.balance_c2_sids -= set(main_data['sid_seq'])
+
+        # adding branches to incorporate connected c2 nodes
+        connected_c2_sids = c2_reactants | c2_products
+        while len(connected_c2_sids) > 0:
+            branch_sid = connected_c2_sids.pop()
+
+            branch_segments = self.get_lin_pwy_containing(branch_sid)
+            self.segments.extend(branch_segments)
+
+            branch_data = self.collect_pwy_data(branch_segments)
+            self.c2_sids_pwy.extend(branch_data['sid_seq'])
+            self.rids_scale |= self.rescale_branch(branch_sid, branch_segments, branch_data)
+
+            net_stoic = self.get_net_stoichiometry()
+            self.net_reactants = {sid: -stoic for sid, stoic in net_stoic.items() if stoic < 0.0}
+            self.net_products = {sid: stoic for sid, stoic in net_stoic.items() if stoic > 0.0}
+            c2_reactants = set(self.net_reactants.keys()).intersection(self.balance_c2_sids)
+            c2_products = set(self.net_products.keys()).intersection(self.balance_c2_sids)
+            self.balance_c2_sids -= set(branch_data['sid_seq'])
+
+            connected_c2_sids = c2_reactants | c2_products
+
+    def rescale_branch(self, branch_sid, branch_segments, branch_data):
+        """Rescale branch reactions, considering joining reaction.
+
+        :param branch_sid: c2 metabolite used as reactant or product in main pathway
+        :type branch_sid: str
+        :param branch_segments: pathway segments of the branch
+        :type branch_segments: list of PathwaySegments
+        :param branch_data: data collected from the branch
+        :type branch_data: dict
+        :return: rescaled reactions of the branch, excluding the joint reaction
+        :rtype: dict (keys: reaction ids, values: reaction scale)
+        """
+        pwy_idx = branch_data['sid_seq'].index(branch_sid)
+        joining_seg = branch_segments[pwy_idx]
+
+        # determine joining reaction with main pathway
+        if joining_seg.producing.rid in self.rids_scale:
+            joint_rid = joining_seg.producing.rid
+        else:
+            joint_rid = joining_seg.consuming.rid
+
+        # rescale reactions in the branch wrt to the joint reaction with the main pathway
+        branch_rids_scale = branch_data['rids_scale']
+        main_scale = self.rids_scale[joint_rid] / branch_rids_scale[joint_rid]
+        rids_scale = {}
+        for branch_rid, branch_scale in branch_rids_scale.items():
+            if branch_rid not in self.rids_scale:
+                rids_scale[branch_rid] = branch_scale * main_scale
+        return rids_scale
+
+    def collect_pwy_data(self, pwy_segs):
+        """Collect date from linear pathway
+
+        Determine sequence of metabolites and sequence of reactions
+        Determine reversivility of the pathway
+        Determine relative scale of the reaction in the pathway
+          Initial reaction set at a scale of 1.0
+          Subsequent reactions scaling depend on input vs.
+            output stoichiometry of a pathway segment as well
+            as of reaction directionality used in the pathway
+
+        :param pwy_segs: pathway segments of the pathway
+        :type pwy_segs: list of PathwaySegments
+        :return: information from linear pathway
+        :rtype: dict
+        """
+        sid_seq = []
+        rid_seq = [pwy_segs[0].producing.rid]
+        reversible = pwy_segs[0].reversible
+        scale = 1.0
+        direction = 1.0 if pwy_segs[0].producing.forward is True else -1.0
+        rids_scale = {pwy_segs[0].producing.rid: scale * direction}
+        for pwy_seg in pwy_segs:
+            sid_seq.append(pwy_seg.sid)
+            rid_seq.append(pwy_seg.consuming.rid)
+            scale *= pwy_seg.producing.stoic / pwy_seg.consuming.stoic
+            direction = 1.0 if pwy_seg.consuming.forward is True else -1.0
+            rids_scale[pwy_seg.consuming.rid] = scale * direction
+        return {'init_sid': self.init_sid, 'sid_seq': sid_seq, 'rid_seq': rid_seq,
+                'rids_scale': rids_scale, 'reversible': reversible}
+
+    def get_lin_pwy_containing(self, c2_sid):
+        """Determine the linear pathway containing a specific c2 node.
+
+        :param c2_sid: metabolite with connectivity 2
+        :type c2_sid: str
+        :return: pathway segments building the linear pathway
+        :rtype: list of PathwaySegments
+        """
+        segments = [self.get_initial_segment(c2_sid)]
+        segments = self.add_downstream_segments(segments)
+        segments = self.add_upstream_segments(segments)
+        return segments
 
     def get_net_stoichiometry(self):
+        """determine net stoichiometry of the pathway
 
+        considering direction scale, including directionality
+        :return: overall metabolite stoichiometry
+        :rtype: dict (keys: metabolites ids; values: stoichiometry)
+        """
         substrates = set()
         for rid in self.rids_scale:
             substrates |= set(self.ng.edges[rid].reactants)
@@ -62,37 +159,39 @@ class LinearPathway:
             for sid, stoic in self.ng.edges[rid].products.items():
                 net_stoic[sid] += stoic * scale
 
-        return net_stoic
+        return {sid: stoic for sid, stoic in net_stoic.items() if stoic != 0.0}
 
-    def get_initial_segment(self):
+    def get_initial_segment(self, c2_sid):
         """get initial pathway segment for inital node.
 
         Node with upstream and downstream connectivity,
         from where connect subsequently other downstream
         and upstream pathway segments
 
+        :param c2_sid: metabolite with connectivity 2
+        :type c2_sid: str
         :return: inital pathway segment
         :rtype: PathwaySegment
         """
-        self.balance_c2_sids.remove(self.init_sid)
-        connectivity = self.get_us_ds_neighbors(self.init_sid)
+        self.balance_c2_sids.remove(c2_sid)
+        connectivity = self.get_us_ds_neighbors(c2_sid)
         return PathwaySegment(connectivity)
 
-    def add_downstream_segments(self):
+    def add_downstream_segments(self, segments):
         """add downstream segments """
 
         lcount = 1
-        upstream_seg = self.segments[-1]
+        upstream_seg = segments[-1]
         while ((upstream_seg.consuming is not None) and
                (upstream_seg.consuming.nb is not None) and
-               (lcount < 50)):
+               (lcount < 100)):
             lcount += 1
             pwy_rev = upstream_seg.reversible
             self.balance_c2_sids.remove(upstream_seg.consuming.nb)
             connectivity = self.get_ds_neighbor(upstream_seg)
             pwy_seg = PathwaySegment(connectivity)
             pwy_seg.set_reversible(pwy_rev and connectivity['reversible'])
-            self.segments.append(pwy_seg)
+            segments.append(pwy_seg)
 
             if connectivity['consuming'] is None:
                 # termination of the pathway at a c1 metabolite (blocked node)
@@ -103,45 +202,47 @@ class LinearPathway:
                 if pwy_rev is True:
                     if connectivity['reversible'] is False:
                         # irreversible termination in a reversible pathway
-                        for pwy_seg in self.segments:
+                        for pwy_seg in segments:
                             pwy_seg.set_reversible(False)
                     else:
                         # reversible termination in a reversible pathway
-                        for pwy_seg in self.segments:
+                        for pwy_seg in segments:
                             pwy_seg.reverse_direction()
-                        self.segments.reverse()
+                        segments.reverse()
 
             if connectivity['consuming'].nb is not upstream_seg.sid:
                 # downstream connectivity to a c2 node
                 if (pwy_rev is True) and (connectivity['reversible'] is False):
-                    for pwy_seg in self.segments:
+                    for pwy_seg in segments:
                         pwy_seg.set_reversible(False)
             elif pwy_rev is True:
                 # pathway reversal situation
                 assert connectivity['reversible'] is False
-                for pwy_seg in self.segments[:-1]:
+                for pwy_seg in segments[:-1]:
                     pwy_seg.set_reversible(False)
                     pwy_seg.reverse_direction()
-                self.segments.reverse()
+                segments.reverse()
             else:
+                self.inconsistant = True
                 print(f'inconsistent configuration at {pwy_seg.sid}, starting from {self.init_sid}')
                 break
-            upstream_seg = self.segments[-1]
+            upstream_seg = segments[-1]
+        return segments
 
-    def add_upstream_segments(self):
+    def add_upstream_segments(self, segments):
 
         lcount = 1
-        downstream_seg = self.segments[0]
+        downstream_seg = segments[0]
         while ((downstream_seg.producing is not None) and
                (downstream_seg.producing.nb in self.balance_c2_sids) and
-               (lcount < 50)):
+               (lcount < 100)):
             lcount += 1
             pwy_rev = downstream_seg.reversible
             self.balance_c2_sids.remove(downstream_seg.producing.nb)
             connectivity = self.get_us_neighbor(downstream_seg)
             pwy_seg = PathwaySegment(connectivity)
             pwy_seg.set_reversible(pwy_rev and connectivity['reversible'])
-            self.segments.insert(0, pwy_seg)
+            segments.insert(0, pwy_seg)
 
             if connectivity['producing'] is None:
                 # termination of the pathway at a c1 metabolite (blocked node)
@@ -150,26 +251,28 @@ class LinearPathway:
             if connectivity['producing'].nb is None:
                 # termination of the pathway with connectivity towards the network
                 if (pwy_rev is True) and (connectivity['reversible'] is False):
-                    for pwy_seg in self.segments:
+                    for pwy_seg in segments:
                         pwy_seg.set_reversible(False)
                 break
 
             if connectivity['producing'].nb is not downstream_seg.sid:
                 # upstream connectivity to a c2 node
                 if (pwy_rev is True) and (connectivity['reversible'] is False):
-                    for pwy_seg in self.segments:
+                    for pwy_seg in segments:
                         pwy_seg.set_reversible(False)
             elif pwy_rev is True:
                 # pathway reversal situation
                 assert connectivity['reversible'] is False
-                for pwy_seg in self.segments[:-1]:
+                for pwy_seg in segments[:-1]:
                     pwy_seg.set_reversible(False)
                     pwy_seg.reverse_direction()
-                self.segments.reverse()
+                segments.reverse()
             else:
+                self.inconsistant = True
                 print(f'inconsistent configuration at {pwy_seg.sid}, starting from {self.init_sid}')
                 break
-            downstream_seg = self.segments[0]
+            downstream_seg = segments[0]
+        return segments
 
     def get_neighbor_node(self, substrates):
         """Get a c2 node from reactants/product sids
@@ -201,31 +304,12 @@ class LinearPathway:
         :rtype: dict
         """
         cur_node = self.ng.nodes[init_sid]
+        assert cur_node.connectivity == 2
 
         rids = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids)]
         revs = [val['rev'] for val in (cur_node.consuming_rids | cur_node.producing_rids).values()]
         rev = bool(np.all(np.array(revs)))
-        producing_leg = None
-        consuming_leg = None
 
-        if cur_node.connectivity == 1:
-            # c1 nodes are actually blocked nodes where a pathway ends
-            rid = rids[0]
-
-            if rid in cur_node.consuming_rids:
-                stoic = cur_node.consuming_rids[rid]['stoic']
-                nb = self.get_neighbor_node(self.ng.edges[rid].products)
-                consuming_leg = ReactionLeg(rid, stoic, nb, True)
-            else:
-                stoic = cur_node.producing_rids[rid]['stoic']
-                nb = self.get_neighbor_node(self.ng.edges[rid].reactants)
-                if rev is False:
-                    producing_leg = ReactionLeg(rid, stoic, nb, True)
-                else:
-                    # in case of a reversible producing reaction, configure consuming path
-                    consuming_leg = ReactionLeg(rid, stoic, nb, False)
-
-        if cur_node.connectivity == 2:
         if rev is True:
             # process a node with reversible fluxes on both sides
             #  one being considered as producing, one as consuming reaction
@@ -310,12 +394,7 @@ class LinearPathway:
     def get_ds_neighbor(self, us_pwy_seg):
         """get downstream neighbor at us_sid coming from p_rid.
 
-        Only for c1 and c2 nodes
-
-        c1 node: this is a terminal node (blocked node),
-          as only exising reaction is consumed by p_rid
-
-        c2 node: downstream is either a consuming reaction or
+        downstream is either a consuming reaction or
         a reversible producing reaction
         downstream node can be another c1/c2 node, or the network connection
         reversible flag is set according to flow through this node
@@ -327,10 +406,10 @@ class LinearPathway:
         """
         sid = us_pwy_seg.consuming.nb
         cur_node = self.ng.nodes[sid]
+        assert cur_node.connectivity == 2
+
         revs = [val['rev'] for val in (cur_node.consuming_rids | cur_node.producing_rids).values()]
         rev = bool(np.all(np.array(revs)))
-        producing_leg = None
-        consuming_leg = None
 
         # collect information for upstream leg
         p_rid = us_pwy_seg.consuming.rid
@@ -343,10 +422,6 @@ class LinearPathway:
             p_stoic = cur_node.consuming_rids[p_rid]['stoic']
             p_fwd = False
 
-        if cur_node.connectivity == 1:
-            producing_leg = ReactionLeg(p_rid, p_stoic, us_nb, p_fwd)
-
-        if cur_node.connectivity == 2:
         c_rid = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids) if rid != p_rid][0]
         if c_rid in cur_node.consuming_rids:
             c_stoic = cur_node.consuming_rids[c_rid]['stoic']
@@ -372,10 +447,10 @@ class LinearPathway:
 
         sid = ds_pwy_seg.producing.nb
         cur_node = self.ng.nodes[sid]
+        assert cur_node.connectivity == 2
+
         revs = [val['rev'] for val in (cur_node.consuming_rids | cur_node.producing_rids).values()]
         rev = bool(np.all(np.array(revs)))
-        producing_leg = None
-        consuming_leg = None
 
         # collect information for downstream leg
         c_rid = ds_pwy_seg.producing.rid
@@ -388,10 +463,6 @@ class LinearPathway:
             c_stoic = cur_node.producing_rids[c_rid]['stoic']
             c_fwd = False
 
-        if cur_node.connectivity == 1:
-            consuming_leg = ReactionLeg(c_rid, c_stoic, ds_nb, c_fwd)
-
-        if cur_node.connectivity == 2:
         p_rid = [rid for rid in (cur_node.consuming_rids | cur_node.producing_rids) if rid != c_rid][0]
         if p_rid in cur_node.producing_rids:
             p_stoic = cur_node.producing_rids[p_rid]['stoic']
@@ -411,10 +482,3 @@ class LinearPathway:
             consuming_leg = ReactionLeg(p_rid, p_stoic, us_nb, not p_fwd)
 
         return {'sid': sid, 'producing': producing_leg, 'consuming': consuming_leg, 'reversible': rev}
-
-    def get_reactions(self):
-
-        rids = set()
-        for pwy_seg in self.segments:
-            rids |= pwy_seg.get_reactions()
-        return rids

modelpruner/graph/linear_pathways.py

@@ -8,57 +8,64 @@ from modelpruner.graph.linear_pathway import LinearPathway
 
 class LinearPathways:
 
-    # TODO: rework again to first get lin_pathways from each of the nodes (excluding biomass precursors)
-    #  idea to get as long pathways as possible
-    #  sort for longest pathways
-    #  remove duplicated pathways
-    #  accept longest pathways first and drop pathways that contain already consumed reactions
-    #  check if c2_nodes remain, and for these request again pathways
+    def __init__(self, network_graph):
+        """Initialize LinearPathways.
 
-    def __init__(self, network_graph, protected_rids=None):
-        """
-
-        Construct linear pathways from metabolites with a connectivity of 2
+        LinearPathways contain c2 metaboltes (metabolites with connectivity of 2).
+        c2 metabolites can be protected using protect_reactions() method,
+         e.g. protect c2 metabolites connected to biomass function
+              protect c2 metabolites connected ot exchange functions
 
-        C1 metabolites (with connectivity of 1) are also included
-        Each such C1/C2 metabolite will only be consumed in once pathway
+        :param network_graph: network graph for fba model
+        :type network_graph: NetworkGraph
         """
-
         self.ng = network_graph
+        self.c2_sids = {node.sid: None for node in self.ng.nodes.values() if node.is_c2_node()}
+        self.pwy_rids = {}
+        self.pathways = []
+        self.inconsistant_pathways = []
+
+    def protect_reactions(self, protected_rids):
+        """Protect metabolites connected to selected reactions.
 
-        # TODO: biomass protected, external metabolites protected
-        # metabolites connected to protected reactions get protected
+        :param protected_rids: list of reaction ids to be protected
+        :type protected_rids: list of str
+        """
         protected_sids = set()
-        if protected_rids is not None:
         for rid in protected_rids:
             if rid in self.ng.edges:
                 protected_sids |= set(self.ng.edges[rid].reactants.keys())
                 protected_sids |= set(self.ng.edges[rid].products.keys())
+        for sid in protected_sids:
+            if sid in self.c2_sids:
+                del self.c2_sids[sid]
 
-        self.c2_sids = {node.sid: None for node in self.ng.nodes.values()
-                        if node.sid not in protected_sids and node.is_c2_node()}
-        self.pwy_rids = {}
-        self.pathways = []
+    def construct(self):
+        """Construct linear pathways for the network.
+
+        We iterate through all c2 metabolites. At each step we select one c2 metabolite
+        that has not already been consumed on one of the linear pathways. From this
+        c2 metabolite we construct a complete linear pathways, including branches that
+        consist of c2 metabolites.
+
+        We update references in c2_sids and pwy_rids to indicate in which pathway
+        respective metabolites and reactions are contained.
 
+        :return: number of linear pathways
+        :rtype: integer
+        """
         pwy_idx = 0
         for sid in sorted(self.c2_sids):
             if self.c2_sids[sid] is None:
                 balance_c2_sids = {sid for sid, idx in self.c2_sids.items() if idx is None}
-                lin_pwy = LinearPathway(sid, balance_c2_sids, network_graph)
+                lin_pwy = LinearPathway(sid, balance_c2_sids, self.ng)
+                if lin_pwy.inconsistant is False:
                     self.pathways.append(lin_pwy)
-
-                for subs_sid in self.pwy_substrates(lin_pwy):
-                    if (subs_sid in self.c2_sids) and (self.c2_sids[subs_sid] is None):
+                    for subs_sid in lin_pwy.c2_sids_pwy:
                         self.c2_sids[subs_sid] = pwy_idx
-                for rid in lin_pwy.get_reactions():
+                    for rid in lin_pwy.rids_scale:
                         self.pwy_rids[rid] = pwy_idx
+                else:
+                    self.inconsistant_pathways.append(lin_pwy)
                 pwy_idx += 1
-
-    def pwy_substrates(self, lin_pwy):
-        sids = set()
-        for pwy_seg in lin_pwy.segments:
-            sids |= set(self.ng.edges[pwy_seg.producing.rid].reactants.keys())
-            sids |= set(self.ng.edges[pwy_seg.producing.rid].products.keys())
-            sids |= set(self.ng.edges[pwy_seg.consuming.rid].reactants.keys())
-            sids |= set(self.ng.edges[pwy_seg.consuming.rid].products.keys())
-        return sids
+        return pwy_idx

modelpruner/graph/metabolite_node.py

@@ -32,7 +32,5 @@ class MetaboliteNode:
             if len(self.consuming_rids) == 2:
                 if bool(np.any(np.array([val['rev'] for val in self.consuming_rids.values()]))) is True:
                     return True
-        elif self.connectivity == 1:  # blocked metabolite
-            return True
         else:
             return False