(* -------------------------------------------------------------------------- *)
(*           Analysis of Protein Signalling in MAPK Cascades                  *) 
(*        based on computing change-of-mind equilibria (Eq^{CoM})             *)
(*                  in Coversion/Preference games                             *)
(* -------------------------------------------------------------------------- *)
(*     Feature      : Provide 3 main functions for                            *)
(*                    1. Construct C/P game (as a graph of situations showing *) 
(*                       a group of proteins are ready to change)             *)
(*                    2. Compute CoM Equilibria in the graph by finding       *) 
(*                       terminal Strongly Connected Components               *)
(*                    3. Analyze paths in CoM graph as pathways               *)
(*     Input        : List of reactions consisting of 3 sets of proteins      *)
(*                      <substrates, products, singleton set of an enzyme>    *)
(*                    and List of protein names to focus on (others are left  *)
(*                    as auxiliaries, e.g., ATP, ADP, Pi, H20)                *)
(*     Output       : List of CoMEq graphs for those focused proteins         *)
(*     Post-process : By given a CoMEq graph with given source & taget names  *)
(*                    Pathway Analyzer compute pathways for those points and  *)
(*                    show the pathways in the graph with green lines         *)
(* ---------------------------------------------------------------------------*)
(*     Developed by : Jittisak Senachak, JAIST                                *)
(*                    Cooperating with Rene Vestergaard, JAIST                *)
(*     Developed in : Wolfram Mathematica(tm) 5.2                             *)
(*     Date         : May 22, 2006                                            *)
(*     All rights reserved by the developers.                                 *)
(* -------------------------------------------------------------------------- *)

<< DiscreteMath`Combinatorica`

(* - Auxillary Functions - *)
SelectNode[e_, vertices_List] := 
   Map[ #[[1]]&, Position[vertices, e]]
                
GetIndexedCoM[{x_, y_, z_}, vertices_List] := 
   Block[ {selectX = selectY = Range[Length[vertices]]},
          Do[ selectX = Intersection[selectX, SelectNode[ x[[i]], vertices ]],
              {i, Length[x]}
          ];
          Do[ selectY = Intersection[selectY, SelectNode[ y[[i]], vertices ]],
              {i, Length[y]}
          ];
          Map[List[z, #]&, Tuples[ List[ selectX, selectY ]]]
   ]

(* Remove reactions that is no important proteins as a substrate or a product *)
FilterReaction[spList_List, important_List, auxiliary_List] :=
    Block[ {filter = {}, sp1 = {}, sp2 = {}},
           Do[ sp1  = Complement[spList[[i,1]], auxiliary];
               sp2  = Complement[spList[[i,2]], auxiliary];

               If[ (Length[sp1] > 0) && (Length[sp2] > 0) &&
                   (Apply[And, Map[MemberQ[important,#]&, Union[sp1,sp2]]]),
                   AppendTo[filter, i]
               ],
               {i, Length[spList]}
           ];
           Return[filter]
   ]

(* Remove all reactions having enzymes in the enzyme_list *)
Inhibit[reactions_List, enzymes_List] :=
    Block[ {filter = {}},
           Do[
               If[ !MemberQ[enzymes, First[Flatten[reactions[[i,3]]]]],
(*         AppendTo[filter, reactions[[i]]]; *)
           AppendTo[filter, reactions[[i]]];
               ],
               {i, Length[reactions]}
           ];
           Return[filter]
   ]

MakeEnzymeEdgeLists[enzPro_List, enzymeCoM_List] :=
   Block[ {j = 1, temp = {}, enzEdgeList = {}},
          Do
          [ temp = {};
              While[ (j <= Length[enzymeCoM]) && 
                     (List[enzPro[[i]]] == enzymeCoM[[j,1]]),
                     AppendTo[temp, enzymeCoM[[j,2]]];
                     j++
               ];
               AppendTo[enzEdgeList, temp],
              {i, Length[enzPro]}
          ];
          Return[enzEdgeList]
   ]
   

MakeTopologicalRootList[enzPro_List, enzEdgeList_List] :=
   Block[ {tVertices = {}, tVertexId = {}, tArrows = {}, g1 = {},
           tSortedList = {}, tRootList = {}, n},
          Do[
              tVertices = Union[Flatten[enzEdgeList[[i]]]];
              tVertexId = Table[ tVertices[[j]] -> j, {j, Length[tVertices]} ];
              tArrows   = Replace[enzEdgeList[[i]], tVertexId, 2];

              g1 = SetGraphOptions[ EmptyGraph[Length[tVertices]],
                                    EdgeDirection -> On 
                   ];
              g1 = AddEdges[g1, tArrows ];
              tSortedList = Map[ tVertices[[#]]&, TopologicalSort[g1]];
  
             (* Take roots of n tological sorted lists *)
              n = Length[ConnectedComponents[g1]];
              AppendTo[tRootList, List[List[enzPro[[i]]], Take[tSortedList, n]]]
              ,{i, Length[enzEdgeList]}
          ];
          Return[tRootList]
   ]

(* Ideal case (no duplication in data): get three tiers of activation *)
GenCollapseInfoTTL[topo_List] :=
   Block[ {cInfo = {}, pos = {}, p},
          Do[
              If[ (cInfo == {}), AppendTo[cInfo, Flatten[topo[[i]]]],
                  pos = Position[cInfo, topo[[i,2,1]]];
                  If[ (pos != {}),
                      p = First[First[pos]];
                      cInfo[[p]] = Union[cInfo[[p]],Flatten[topo[[i]]]],
                                   AppendTo[cInfo, Flatten[topo[[i]]]]
                  ]
              ],
              {i, Length[topo]}
          ];
          Return[Map[Union[Flatten[#]]&,cInfo]]
   ]

(* Show FACT : all enzymes in the same node can activate any proteins inside *)
GenCollapseInfoTTS[topo_List] :=
   Block[ {cInfo = {}, pos = {}, p},
          Do[
              If[ (cInfo == {}), AppendTo[cInfo, topo[[i]]],
                  pos = Position[cInfo, topo[[i,2]]];
                  If[ (pos != {}),
                      p = First[First[pos]];
                      cInfo[[p]] = Union[cInfo[[p]], topo[[i]]],
                                   AppendTo[cInfo,topo[[i]]]
                  ]
              ],
              {i, Length[topo]}
          ];
          Return[Map[Union[Flatten[#]]&,cInfo]]
    ]

(* Show what enzyme a protein can be activated *) 
(*      (a protein surrounding with enzymes) *)
GenCollapseInfoTFS[topo_List] :=
   Block[ {cInfo = {}, node = {}},
          node  = Apply[Union, Map[#[[2]]&, topo]];
          cInfo = Apply[Union, Table[Map[List[List[#],topo[[i,1]]]&,topo[[i,2]]]
                                     ,{i, Length[topo]}]];
          Return[MapThread[Prepend[Flatten[#2],#1]&, 
                                   List[node,MakeEnzymeEdgeLists[node, cInfo]]]]
   ]

(* Show what protein an enzyme can activate *)
(*      (an enzyme surrounding with proteins) *)
GenCollapseInfoFTS[topo_List] :=
   Block[ {cInfo = {}, node = {}},
          node  = Apply[Union, Map[#[[1]]&, topo]];
          cInfo = Apply[Union, Table[Map[List[topo[[i,1]],List[#]]&,topo[[i,2]]]
                                     ,{i, Length[topo]}]];
          Return[MapThread[Prepend[Flatten[#2],#1]&,
                                   List[node,MakeEnzymeEdgeLists[node, cInfo]]]]
   ]

(* ?? *)
GenCollapseInfoFFS[topo_List] :=
   Block[ {cInfo = {}, node = {}},
          cInfo = Apply[Union, Table[Map[Prepend[topo[[i,1]],#]&, topo[[i,2]]]
                                     ,{i, Length[topo]}]];
          Return[cInfo]
   ]

GetCollapse[enzPro_List, enzymeCoM_List, verticesNoAux_List] :=
   Block[ {enzEdgeLists = {}, topoRootList = {}, 
           nodesToRemove = {}, collapseInfo = {}, vertexId = {},
           inCol = False, outCol = False, strict = True},

          Print["Find a new graph with Collapse/Duplicate", Date[]];

      (* - 1. Group arrows to each enzyme -- *)
          enzEdgeList = MakeEnzymeEdgeLists[enzPro, enzymeCoM];

      (* - 2. Topological Sort for each enzyme in enzEdgeList - *)
          topoRootList = MakeTopologicalRootList[enzPro, enzEdgeList];

      (* - 3. Generate Collapsed Info - *)
          vertexId = Table[verticesNoAux[[i]] -> i, {i, Length[verticesNoAux]}];
          nodesToRemove = Union[Flatten[Replace[Union[topoRootList],
                                                vertexId,2]]];

          vertexId = Table[i -> verticesNoAux[[i]], {i, Length[verticesNoAux]}];
          If[(inCol && outCol && !strict),
             collapseInfo = Replace[GenCollapseInfoTTL[topoRootList],vertexId,2]
          ];
          If[(inCol && outCol && strict),
             collapseInfo = Replace[GenCollapseInfoTTS[topoRootList],vertexId,2]
          ];
          If[(inCol && !outCol && strict),
             collapseInfo = Replace[GenCollapseInfoTFS[topoRootList],vertexId,2]
          ];
          If[(!inCol && outCol && strict),
             collapseInfo = Replace[GenCollapseInfoFTS[topoRootList],vertexId,2]
          ];
          If[(!inCol && !outCol && strict),
             collapseInfo = Replace[GenCollapseInfoFFS[topoRootList],vertexId,2]
          ];

      (* --- *)
          Return[Union[Map[Flatten, collapseInfo],
                      verticesNoAux[[Complement[Range[Length[verticesNoAux]], 
                                           nodesToRemove]]]]];
   ]

(* --- AUX: Compute Change-of-Mind Equilibria --- *)
(* ---      by SCCs and theirs terminals      --- *)

ExistNode[g_Graph, n_Integer] :=
   MemberQ[Range[V[g]], n];
IsNoEdge[g_Graph, n_Integer] :=
   (ExistNode[g, n] && (OutDegree[g,n] == 0) && (InDegree[g,n] == 0))
IsOneSelf[g_Graph, n_Integer] :=
   Block[ {adjEdge = ToAdjacencyLists[g]},
          Return[ ExistNode[g, n] && 
                  (OutDegree[g,n] == 1) && (InDegree[g,n] == 1) && 
                  adjEdge[[n]] == List[n]
          ];
   ]
IsTerminal[g_Graph, n_Integer] := 
   Block[ {adjEdge = ToAdjacencyLists[g]},
      Return[ MemberQ[Range[V[g]], n] &&
                  ((adjEdge[[n]] == List[]) || (adjEdge[[n]] == List[n]))
          ];
   ]
IsSingleton[terminal_List, n_Integer] :=
   (MemberQ[Range[Length[terminal]], n]) && (Length[terminal[[n]]] == 1);

Shrunken[g_Graph, sccC_List] :=
   Block[ {sg, edge ={}},

      (* - 1. Assign sccIds for all nodes in sccC list - *)
          Do[  
              Do[ sccIds[[ sccC[[i,j]] ]] = i, 
                  {j, Length[sccC[[i]]]}
              ],{i, Length[sccC]}
          ];

     (* - 2. Create Shrunken Graph with following conditions - *)
     (* -    - Nodes : all distinct SCCIds                   - *)
     (* -    - Edges : all edges that connects among Nodes   - *)
     (* -            (no edges connecting inside same SCCId) - *)
          edgeList = ToOrderedPairs[g];
          sg       = SetGraphOptions[ EmptyGraph[Length[sccC] ], 
                                     EdgeDirection -> On ];
          sg       = SetGraphOptions[ sg, VertexColor -> Red ];
          sg       = SetGraphOptions[ sg, VertexStyle -> Disk[Large] ];
          sg       = SetVertexLabels[ sg, sccC];

          Do[ edge = List[ sccIds[[ edgeList[[i,1]] ]], 
                           sccIds[[ edgeList[[i,2]] ]]]; 
              If[ (edge[[1]] != edge[[2]]),
                  sg = AddEdge[ sg, edge]
              ],
              {i, Length[edgeList]}
          ];
     (*   sg      = RemoveMultipleEdges[sg]; *)
          sg
   ]

CreateTGraph[tNodeList_List, nodeLabels_List] :=
   Block[ { tNamesMap = {}, tIdsMap = {}, tNodesMap = {}, tNodeNames = {}, 
            tg, edge = {}},

          tNamesMap  = Table[i -> nodeLabels[[i]], {i,Length[nodeLabels]}];
          tNodeNames = Map[Replace[#, tNamesMap]&, tNodeList];

      (* - 1. Get a mapping of tg\'s edges to g\'s edges - *) 
          tIdsMap = Table[ tNodeList[[j]] -> j, {j, Length[tNodeList]} ];

      (* - 2. Create Nodes: all nodeIds in terminal list- *)
          tg  = SetGraphOptions[ EmptyGraph[Length[tNodeList]], 
                                 EdgeDirection -> On
                ];
          tg  = SetGraphOptions[tg, VertexColor -> Green];
          tg  = SetGraphOptions[tg, VertexStyle -> Disk[Large]];
          tg  = SetVertexLabels[tg, tNodeNames];

     (* - 3. Add edges: For any 2 g\'edges which are in terminal list and - *)
     (* -                                        share SCCId              - *)
          Do[ edge = edgeList[[j]];
              If[ (MemberQ[tNodeList, edge[[1]]]) && 
                  (MemberQ[tNodeList, edge[[2]]]) &&
                  (sccIds[[ edge[[1]] ]] == sccIds[[ edge[[2]] ]]),
                  tg = AddEdge[tg, List[ Replace[edge[[1]], tIdsMap], 
                                         Replace[edge[[2]], tIdsMap] ] 
                       ]
              ],
              {j, Length[edgeList]}
          ];
          tg   = RemoveMultipleEdges[tg];
          Return[tg];
   ]

(* --- AUX: Direct Path Analyzer (OwnBT) --- *)     
ExistsTagged[v_Integer] :=
   Do[ If[ (table[[v, i+1]] == "TAGGED"), Return[True]],
       {i, Length[table[[v]]]-1}
   ]

OwnBTRec[cn_Integer] :=
   Block[ {},
          If[ (MemberQ[drain, cn]), Return[True]];
          If[ (table[[cn,1]] == "VISITED"), Return[False]];

          table[[cn, 1]] = "VISITED";
          Do[ If[ (OwnBTRec[adjEdge[[cn,j]]]), 
                  table[[cn,j+1]] = "TAGGED"
              ],
              {j, Length[adjEdge[[cn]]]}
          ];
          table[[cn, 1]] = "NOTVISIT";

          Return[ExistsTagged[cn]]
  ]

OwnBT[g_Graph, source_List, drain_List] :=
   Block[ {adjEdge = ToAdjacencyLists[g], table = {}, row = path = {} },

     (* - 1. Prepare table for keeping VISTED/TAGGED info - *)
          Do[ AppendTo[table, Join[ List["NOTVISIT"],
                                    Table["NOTTAG", {j, Length[adjEdge[[i]]]}]
                              ]
              ],
              {i, V[g]}
          ];

     (* - 2. Do backtrack for search all direct paths - *)
          Do[ OwnBTRec[source[[i]]],
              {i, Length[source]}
          ];

     (* - 3. Interpret computed information (TAGGED info) - *)
          Do[
              Do[ If[ (table[[i,j+1]] == "TAGGED"),
                      AppendTo[path, {i, adjEdge[[i,j]]}]
                  ],
                  {j, Length[table[[i]]]-1}
              ],
              {i, Length[table]}
          ];

          Return[path]
   ]

(* ====================== *)
(* === Main Functions === *)

constructCPG[allReactions_List, important_List] :=
   Block[ {spList = {}, sp = {}, auxiliary = {}, reactions = {}, enzymes = {}, 
           proteins = {}, vertices = {}, enzymeCoM = {}, CoM = {}, enzPro = {},
           g, gList ={}, nodesToRemove = {}},

     (* - 1. Extract relevant information from reaction - *)
          spList   = Map[List[#[[1]], #[[2]]]&, allReactions];
          sp       = Union[Flatten[spList]];

         (* cut off information in reactions, but not in focus *)
          auxiliary= Union[Complement[sp, important]];
          proteins = Union[Intersection[sp, important]];
         (* reactions : allReaction filtered by important *)
          reactions= allReactions[[FilterReaction[spList,important,auxiliary]]];
          enzymes  = Apply[Union, Map[#[[3]]&, reactions]];

     (* - 2. Set a vertex as a protein having aux. around - *)
          Print["Creating vertices : ", Date[]]; 
          vertices = Map[ Union[auxiliary,List[#]]&, proteins ];
          Print[Length[vertices], " vertices ... ", Date[], " Done."];

     (* - 3. Get edges(CoM) from reaction information (reaction as CoM) -- *)
          Print["Creating Change-of-Mind relations : ", Date[]]; 
          enzymeCoM = Apply[Union, Map[ GetIndexedCoM[#,vertices]&, reactions]];
          CoM       = Union[Map[#[[2]]&, enzymeCoM]];
          Print[ Length[CoM], " (", Length[enzymeCoM], " enzyme-specified) ",
                 "Change-of-Mind ... ", Date[], " Done."];

     (* - 4. Regenerate a collapsed graph g -- *)
          Print["Creating collapsed vertices : ", Date[]];
          enzPro = Intersection[enzymes,proteins];
          vertices = GetCollapse[ enzPro, 
                                  Select[enzymeCoM,MemberQ[enzPro,#[[1,1]]]&], 
                                  Map[Complement[#,auxiliary]&, vertices]
                     ];
          vertices = Map[ Union[auxiliary,#]&, vertices ];
          Print[Length[vertices], " collapsed vertices... ", Date[], " Done."];

          Print["Creating collapsed Change-of-Mind : ", Date[]];
          enzymeCoM = Apply[Union, Map[ GetIndexedCoM[#,vertices]&, reactions]];
          CoM       = Union[Map[#[[2]]&, enzymeCoM]];
          Print[Length[CoM], " (", Length[enzymeCoM], " enzyme-specified) ",
                "Change-of-Mind ... ", Date[], " Done."];

          g = SetGraphOptions[ EmptyGraph[Length[vertices]],
                               EdgeDirection -> On
              ];
          g = AddEdges[g, CoM];
     (* - Cut auxiliary info in g - *)
          g = SetVertexLabels[g,Map[Complement[#,auxiliary]&, vertices]] ;
          Print[Date[], " Done."];

          gList = List[ Union[ Flatten[CoM] , {VertexColor -> Blue } ]];
          g = SetGraphOptions[ g, gList];
          Return[g];
   ]

findCoMEq[g_Graph] :=
   Block[ {sccC = {}, edgeList = {}, sccIds = {}, nodeLabels = {},
           sg, terminal = {}, singTerm = {}, tgList = {}},

          sccIds = Table[-1, {V[g]}];
          edgeList = ToOrderedPairs[g];
          nodeLabels = GetVertexLabels[g];

     (* - 1. Compute SCC (list of list(group) of strongly connected nodes)- *)
          sccC = StronglyConnectedComponents[g];

     (* - 2. Compute Shrunken SCC Graph - *)
          sg = Shrunken[ g, sccC ];

     (* - 3. Identify Terminal SCC (from sccList) - *)
          terminal = sccC[[ Select[Range[Length[sccC]], (IsTerminal[sg,#])& ]]];
      singTerms = terminal[[ Select[Range[Length[terminal]], (IsSingleton[terminal,#])&]]];
          Print["There are ", Length[terminal], " CoM equilibria. (of which ", 
                              Length[singTerms], " are pure Nash equilibria)"];

     (* - 4. Make a list of CoM graph - *)
          Do[
             AppendTo[tgList, CreateTGraph[ terminal[[i]], nodeLabels ]],
             {i, Length[terminal]}
          ];
          Return[tgList]
  ]

doDPA[tg_Graph, sourceNames_List, drainNames_List] :=
   Block[ {gList = {}, path = {}, source = {}, drain = {}, 
           nodeMap = {}, vertexLabel = {}},

     (* - 1. Get nodeIds for source & drain - *)
          nodeLabels = GetVertexLabels[tg];
          Do[ 
              Do[ AppendTo[nodeMap, nodeLabels[[i,j]] -> i],
                  {j, Length[nodeLabels[[i]]]}
              ],
              {i, Length[nodeLabels]}
          ];
          source = Map[Replace[#, nodeMap]&, sourceNames];
          drain  = Map[Replace[#, nodeMap]&, drainNames];

     (* - 2. Find Direct Path (by Backtracking technique) - *)
          path = OwnBT[ tg, source, drain ];
          Map[ Replace[#, Table[i->nodeLabels[[i]], {i,Length[nodeLabels]}],2]&, 
               path]
   ]

comEq[reactions_List, focus_List] :=
   findCoMEq[constructCPG[reactions, focus]]
   
PathAnalyse[label_String, g_Graph, source_List, sink_List] :=
  Block[ {dpa = doDPA[g, source, sink], nodeLabels = {}, nodeMap = {}},
    
         Print[ "The direct paths from ", source, " to ", sink, " in ", label, 
                " involve ", Length[Union[Flatten[dpa]]], 
                " unique representative compounds: "];
         Print[Union[Flatten[dpa]]];
         
         Print[ "The direct paths from ", source, " to ", sink, " in ", label, 
                " involve ", Length[Union[dpa]], " unique reactions:"];
         Print[ColumnForm[Union[dpa]]];
         
         nodeLabels = GetVertexLabels[g];
         nodeMap = Table[nodeLabels[[i]] -> i, {i, Length[nodeLabels]}];
         ShowGraph[g, Union[ Map[Replace[#, nodeMap]&, dpa, 2], 
                             {EdgeColor -> Green, EdgeStyle -> Thick}]]; 
                             
        Return[List[Union[Flatten[dpa]],Union[dpa]]];
  ]