[Bio] / FigKernelPackages / model.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/model.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.32 - (view) (download) (as text)

1 : olson 1.1
2 :     ## _*_ Perl _*_ ##
3 :     #
4 :     # model.pm
5 :     #
6 :     # Kevin Formsma
7 :     # Hope College
8 :     # Created: 6/1/2006
9 :     #
10 :     ##################
11 :    
12 :     ### BrainStorm/Planning ###
13 :     #
14 :     # Functions
15 :     # -Scenario access methods for subsystem
16 :     # -Access current genome models
17 :     # -Interface to/Integration of find_reaction_paths.cgi
18 :     # -Scenario Relationships:
19 :     # 1. Automanaged - map togather based on starting/ending compounds
20 :     # 2. User Defined - through some type of map or XML
21 :     # -Generate files for Scenarios on a organism
22 :     # 1. Path Picking - All or Select?
23 :     # -Filter written Scenarios based on Relationships
24 :     # -Report valid Scenarios for a genome
25 :     #
26 :     #
27 :     # Basic Process for a Control Script
28 :     # 1. Input Genomes for model creation and Subsystems
29 :     # 2. Define selected Scenario relationships
30 :     # 3. Create Scenario Paths for each genome and subsystem
31 :     # 4. Filter invalid Scenarios based on relationships
32 :     # 5. Combine Subsystem models for each genome into genome models
33 :     # 6. Report on the model results for each genome
34 :     #
35 :    
36 :     package model;
37 :    
38 :     use strict;
39 : parrello 1.32 no warnings 'redefine';
40 : formsma 1.14 use Scenario;
41 : olson 1.1 use FIG;
42 :     use Subsystem;
43 :     use File::Path;
44 :    
45 : parrello 1.32 my $fig = eval("FIG->new()"); ## So that this will compile properly on non-FIG systems!
46 : olson 1.1
47 :     #global variables that make process_paths work
48 :     #These need to be cleared and reloaded frequently for process path and flux writing
49 :     my (%reactions_to_substrate_arrays, %reactions_to_product_arrays, %all_compounds_to_main);
50 :     my %all_reactions;
51 :     my %scenario_cycles;
52 :     my @all_outputs_lists;
53 :     my %all_inputs;
54 :     my %all_outputs;
55 :    
56 : dejongh 1.12 #this variable is used to set the loop count for single scenario and then assembly runs, 100/25 by default.
57 : olson 1.1 my $loop_max = 100;
58 :     my $loop_max_assembly = 25;
59 :    
60 :     #These store what supersets/Subsystems we are using
61 :     my %superset_to_ss;
62 :     my %ss_to_superset;
63 :    
64 :     #Flip this bit to enable debugging
65 : olson 1.11 my $debug = int($ENV{HOPE_DEBUG});
66 : dejongh 1.21 $debug = 0;
67 : olson 1.11
68 :     sub set_fig
69 :     {
70 :     my($newfig) = @_;
71 :     $fig = $newfig;
72 :     }
73 : olson 1.1
74 :     sub new {
75 :     my $type = shift;
76 :     my $self = {};
77 :     return bless $self, $type;
78 :     }
79 :    
80 :     sub get_ss_scenarios
81 :     {
82 :     my ($ss_name) = @_;
83 :     # This is the scenario data structure storage
84 :     # %scenario_data{scenario_name}
85 :     #
86 :     my %scenario_data;
87 :    
88 :     my $subsystem = $fig->get_subsystem($ss_name);
89 : olson 1.11 if (!$subsystem)
90 :     {
91 :     warn "Cannot open subsystem $subsystem\n";
92 :     return;
93 :     }
94 :    
95 : olson 1.1 my @scenario_names = $subsystem->get_hope_scenario_names;
96 :     foreach my $name (@scenario_names)
97 :     {
98 :     $scenario_data{$name} = &get_scenario($subsystem,$name);
99 :     }
100 :    
101 :     return \%scenario_data;
102 :     }
103 :    
104 :    
105 :     sub get_scenario
106 :     {
107 :     my($subsystem,$name) = @_;
108 :    
109 :     my %data;
110 :    
111 : dejongh 1.17 my @inputs = $subsystem->get_hope_input_compounds($name);
112 :     my @outputs = $subsystem->get_hope_output_compounds($name);
113 :     my @map_ids = $subsystem->get_hope_map_ids($name);
114 :     my @additional_reactions = $subsystem->get_hope_additional_reactions($name);
115 :     my @ignore_reactions = $subsystem->get_hope_ignore_reactions($name);
116 : olson 1.1
117 :     $data{inputs} = \@inputs;
118 :     $data{outputs} = \@outputs;
119 :     $data{map_ids} = \@map_ids;
120 :     $data{additional_reactions} = \@additional_reactions;
121 :     $data{ignore_reactions} = \@ignore_reactions;
122 :    
123 :     return \%data;
124 :    
125 :     }
126 :    
127 :     sub process_init
128 :     {
129 : dejongh 1.19 my ($ss_name,$scenario_name,$genome,$assembly) = @_;
130 : olson 1.1 my (%sc_inputs, %sc_outputs);
131 :    
132 :     if ($genome eq "")
133 :     {
134 :     $genome = "All";
135 :     }
136 :     print STDERR "\nSubsystem : ".$ss_name." Scenario: $scenario_name \n" if $debug;
137 :     my $subsystem = $fig->get_subsystem($ss_name);
138 :     my $scenario_data = &get_scenario($subsystem,$scenario_name);
139 : dejongh 1.4
140 : olson 1.1 #load the other arrays
141 :     my %ss_reactions;
142 :    
143 : dejongh 1.4 if ($genome eq "All")
144 : olson 1.1 {
145 : dejongh 1.17 my %all_reactions = $subsystem->get_hope_reactions;
146 : dejongh 1.4 foreach my $role (keys %all_reactions)
147 : olson 1.1 {
148 : dejongh 1.4 map { $ss_reactions{$_} = 1 } @{$all_reactions{$role}};
149 : olson 1.1 }
150 :     }
151 : dejongh 1.4 else
152 : olson 1.1 {
153 : dejongh 1.24 my %reactions_for_genome = get_reactions_for_genome_in_subsystem($genome,$ss_name);
154 : parrello 1.32 map { $ss_reactions{$_} = 1 } keys %reactions_for_genome if keys %reactions_for_genome;
155 : olson 1.1 }
156 : dejongh 1.4
157 : olson 1.1 map { $sc_inputs{$_} = 1 } @{$scenario_data->{inputs}};
158 :    
159 :     foreach my $list (@{$scenario_data->{outputs}})
160 :     {
161 :     map { $sc_outputs{$_} = 1 } @$list;
162 :     push @all_outputs_lists, $list;
163 :     }
164 :    
165 :     map { $scenario_cycles{$_} = 1 if defined $sc_outputs{$_} } keys %sc_inputs;
166 :     map { $all_inputs{$_} = 1 } keys %sc_inputs;
167 :     map { $all_outputs{$_} = 1 } keys %sc_outputs;
168 :    
169 :     my @hope_additional_reactions = @{$scenario_data->{additional_reactions}};
170 :     my @hope_ignore_reactions = @{$scenario_data->{ignore_reactions}};
171 :     my %sc_reactions;
172 :     map { $sc_reactions{$_} = 1 } keys %ss_reactions;
173 :    
174 :     # flag additional reactions so we won't check if they are in a map
175 :     foreach my $rid (@hope_additional_reactions)
176 :     {
177 :     $sc_reactions{$rid} = 2;
178 :     }
179 :    
180 :     foreach my $rid (@hope_ignore_reactions)
181 :     {
182 :     delete $sc_reactions{$rid};
183 :     }
184 :    
185 :     #for now we do this outside of the if statment, but that might need to change
186 : dejongh 1.19 &load_substrate_and_product_arrays(\%sc_reactions,$scenario_data->{map_ids}, \%sc_inputs);
187 : olson 1.1
188 :    
189 :    
190 :    
191 :     }
192 :    
193 :     sub execute_paths
194 :     {
195 :     my ($assembly_paths,$find_first,$input_path, $output_path) = @_;
196 :    
197 :     my $num_paths = scalar @{$assembly_paths};
198 :    
199 :     my (%substrates_to_reactions, %products_to_reactions,
200 :     %reactions_to_substrates, %reactions_to_products);
201 :    
202 :     if($num_paths ==0 ) #Load the reactions normally if we are creating scenarios
203 :     {
204 :     create_reactions(\%substrates_to_reactions,\%products_to_reactions,
205 :     \%reactions_to_substrates,
206 :     \%reactions_to_products);
207 :     }
208 :     #create assembly reactions that are closed 'paths' from scenarios or assemblies
209 :     if($num_paths > 0)
210 :     {
211 :     create_assembly_reactions(\%substrates_to_reactions,\%products_to_reactions,
212 :     \%reactions_to_substrates,
213 :     \%reactions_to_products,$assembly_paths);
214 :     }
215 :    
216 :     #This deals with user specifed input/output paths, and uses these to generate
217 :     #what our input and output compounds should be for an assembly.
218 :     #This is only used for creating assemblies, has no effect on scenario creation
219 :    
220 :     if(scalar @{$input_path} && scalar @{$output_path})
221 :     {
222 :    
223 :     %all_inputs = ();
224 :     %all_outputs = ();
225 :     #we should have loaded these paths above, so lets just pull out the information we need
226 :     foreach my $path (@$input_path)
227 :     {
228 :     my $input_rxn = "$path->[-3]/$path->[-2]/$path->[-1]_R";
229 :     my @user_in = @{$reactions_to_substrates{$input_rxn}};
230 :     map{ $all_inputs{$_} = 1 } @user_in;
231 :     }
232 :     foreach my $path (@$output_path)
233 :     {
234 :     my $output_rxn = "$path->[-3]/$path->[-2]/$path->[-1]_R";
235 :     my @user_out = @{$reactions_to_products{$output_rxn}};
236 :     map{ $all_outputs{$_} = 1 } @user_out;
237 :     }
238 :     #map { $scenario_cycles{$_} = 1 if defined $all_outputs{$_} } keys %all_inputs;
239 :     }
240 :    
241 :     print STDERR "Inputs :\n" if $debug;
242 : olson 1.2 print STDERR map { $_."\n" } keys %all_inputs if $debug;
243 : olson 1.1 print STDERR "Outputs:\n" if $debug;
244 : olson 1.2 print STDERR map { $_."\n" } keys %all_outputs if $debug;
245 : olson 1.1
246 :     #filter the input/outputs lists, removing the intersection unless something
247 :     # is a known cycle
248 :     foreach my $input (keys %all_inputs)
249 :     {
250 :     if(defined $all_outputs{$input} && ! defined $scenario_cycles{$input})
251 :     {
252 :     print STDERR "Deleting $input from input and output lists\n" if $debug;
253 :     delete $all_inputs{$input};
254 :     delete $all_outputs{$input};
255 :     }
256 :     }
257 :    
258 :     my $create_assembly = 0;
259 :     $create_assembly = 1 if(scalar @{$assembly_paths} !=0);
260 :     return process_paths(\%all_inputs, \%all_outputs, \@all_outputs_lists,
261 :     \%reactions_to_substrates, \%reactions_to_products,
262 :     \%substrates_to_reactions,\%products_to_reactions,$create_assembly,$find_first);
263 :     }
264 :    
265 :    
266 :     sub create_reactions
267 :     {
268 :     my ($substrates_to_reactions, $products_to_reactions,
269 :     $reactions_to_substrates, $reactions_to_products) = @_;
270 : olson 1.2 print STDERR "building SS reactions\n" if $debug;
271 : olson 1.1 # use subsystem reactions
272 :     foreach my $drxn (map { ($_."_L", $_."_R") } keys %all_reactions)
273 :     {
274 :     foreach my $substrArr ($reactions_to_substrate_arrays{$drxn})
275 :     {
276 :     foreach my $cinfo (@$substrArr)
277 :     {
278 :     my $cpd = $cinfo->[0];
279 :     my $main = $cinfo->[2] || defined $all_inputs{$cpd}; # main in this reaction
280 :    
281 :     if ($main)
282 :     {
283 :     push(@{$reactions_to_substrates->{$drxn}}, $cpd);
284 :     push(@{$substrates_to_reactions->{$cpd}}, $drxn);
285 :     }
286 :     }
287 :     }
288 :    
289 :     foreach my $prodArr ($reactions_to_product_arrays{$drxn})
290 :     {
291 :     foreach my $cinfo (@$prodArr)
292 :     {
293 :     my $cpd = $cinfo->[0];
294 :     my $main = $cinfo->[2] || defined $all_outputs{$cpd}; # main in this reaction
295 :    
296 :     if ($main)
297 :     {
298 :     push(@{$reactions_to_products->{$drxn}}, $cpd);
299 :     push(@{$products_to_reactions->{$cpd}}, $drxn);
300 :     }
301 :     }
302 :     }
303 :     }
304 :     }
305 :    
306 :     sub create_assembly_reactions
307 :     {
308 :     my ($substrates_to_reactions, $products_to_reactions,
309 :     $reactions_to_substrates, $reactions_to_products,$assembly_paths) = @_;
310 :     my %intersection;
311 :    
312 :     foreach my $path (@$assembly_paths)
313 :     {
314 : dejongh 1.20 my $genome = shift @$path;
315 : dejongh 1.23 my $paths_dir = get_scenario_directory($genome) . "/" . join "/" , @$path;
316 : olson 1.1
317 :     my $drxn = "$path->[-3]/$path->[-2]/$path->[-1]_R";
318 :    
319 :     $all_reactions{"$path->[-3]/$path->[-2]/$path->[-1]"} = "R";
320 :    
321 :     print STDERR "Making reaction: $drxn\n" if $debug;
322 :    
323 :     open (M_INPUTS , $paths_dir."/inputs") or die ("Failed to open $paths_dir"."/inputs");
324 :     my @substrArr;
325 :    
326 :     print STDERR "Gathering Inputs:\n" if $debug;
327 :    
328 :     while (<M_INPUTS>)
329 :     {
330 :     my ($cpd, $stoich) = split "\t" , $_;
331 :    
332 :     #We are going to assume everything is a main...
333 :     if(!defined $all_compounds_to_main{$cpd})
334 :     {
335 :     $all_compounds_to_main{$cpd} = 1;
336 :     }
337 :    
338 :     if ($all_compounds_to_main{$cpd})
339 :     {
340 :     push(@{$reactions_to_substrates->{$drxn}}, $cpd);
341 :     push(@{$substrates_to_reactions->{$cpd}}, $drxn);
342 :     }
343 :    
344 :     if($all_compounds_to_main{$cpd} !=0 || !defined $all_compounds_to_main{$cpd})
345 :     {
346 :     $all_inputs{$cpd} = 1;
347 :     }
348 :     push @substrArr, [$cpd, $stoich, $all_compounds_to_main{$cpd}];
349 :    
350 :     my @names = $fig->names_of_compound($cpd);
351 :     print STDERR "\t$stoich\t$cpd\t$names[0]\t$all_compounds_to_main{$cpd}\n" if $debug;
352 :     }
353 :    
354 :     $reactions_to_substrate_arrays{$drxn} = \@substrArr;
355 :    
356 :     close M_INPUTS;
357 :    
358 :     open (M_OUTPUTS, $paths_dir."/outputs") or die("Failed to open $paths_dir"."/outputs");
359 :     my @prodArr;
360 :    
361 :     print STDERR "Gathering outputs:\n" if $debug;
362 :    
363 :     while (<M_OUTPUTS>)
364 :     {
365 :     my ($cpd, $stoich) = split "\t", $_;
366 :     print STDERR "Found $stoich $cpd\n" if $debug;
367 :    
368 :     #We are going to assume everything is a main...
369 :     if(!defined $all_compounds_to_main{$cpd})
370 :     {
371 :     $all_compounds_to_main{$cpd} = 1;
372 :     }
373 :    
374 :     if ($all_compounds_to_main{$cpd})
375 :     {
376 :     push(@{$reactions_to_products->{$drxn}}, $cpd);
377 :     push(@{$products_to_reactions->{$cpd}}, $drxn);
378 :     }
379 :    
380 :    
381 :     #This adds cycles from 'assemblys' because they weren't added earlier
382 :     foreach my $ele (@{$reactions_to_substrates->{$drxn}})
383 :     {
384 :     if($ele eq $cpd)
385 :     {
386 :     $scenario_cycles{$cpd} = 1;
387 :     }
388 :     }
389 :    
390 :     if($all_compounds_to_main{$cpd} !=0 || !defined $all_compounds_to_main{$cpd})
391 :     {
392 :     $all_outputs{$cpd} = 1;
393 :     }
394 :     push @prodArr, [$cpd,$stoich,$all_compounds_to_main{$cpd}];
395 :     my @names = $fig->names_of_compound($cpd);
396 :     print STDERR "\t$cpd\t$names[0]\t$all_compounds_to_main{$cpd}\n" if $debug;
397 :     }
398 :    
399 :     $reactions_to_product_arrays{$drxn} = \@prodArr;
400 :    
401 :     close M_OUTPUTS;
402 :    
403 :     }
404 :     }
405 :    
406 :    
407 :     sub load_substrate_and_product_arrays
408 :     {
409 : dejongh 1.19 my ($reactions, $map_ids, $sc_inputs) = @_;
410 : olson 1.1
411 :     # determine whether the reaction is in one of the maps, and get directionality accordingly
412 : dejongh 1.17 my (%reactions_in_maps, %reactions_not_in_maps, %reactions_not_in_any_map);
413 : olson 1.1
414 :     foreach my $rxn (keys %$reactions)
415 :     {
416 :     my $direction;
417 :    
418 :     if($fig->valid_reaction_id($rxn))
419 :     {
420 :     # get an array of triplets. The triplets are [reaction id][map id]
421 :     # [left to right - R, right to left - L, or both - B]
422 :     my @triplets = $fig->reaction_direction($rxn);
423 :    
424 :     foreach my $trip (@triplets)
425 :     {
426 :     foreach my $map_id (@$map_ids)
427 :     {
428 :     if (@{$trip}[1] eq $map_id)
429 :     {
430 :     my $this_direction = @{$trip}[2];
431 :    
432 :     # bidirectional in one map overrules unidirectional in another.
433 :     # opposite directions in two maps becomes bidirectional
434 :     if (! defined $direction)
435 :     {
436 :     $direction = $this_direction;
437 :     }
438 :     elsif ($direction ne "B" && ($this_direction eq "B" ||
439 :     $this_direction ne $direction))
440 :     {
441 :     $direction = "B";
442 :     }
443 :    
444 :     $reactions_in_maps{$rxn} = 1;
445 :     }
446 :     }
447 :     }
448 :    
449 :    
450 :     if(! $reactions_in_maps{$rxn})
451 :     {
452 :     my $found_in_other_map = 0;
453 :    
454 :     #reaction not in scenario map ids, try to get directionality from other maps
455 :     foreach my $trip (@triplets)
456 :     {
457 :     my $this_direction = @{$trip}[2];
458 :    
459 :     #bidreactional in one map overrules unidirectional in another
460 :     #opposite directions in two maps becomes bidirectional
461 :     if(! defined $direction)
462 :     {
463 :     $direction = $this_direction;
464 :     }
465 :     elsif($direction ne "B" && ($this_direction eq "B" || $this_direction ne $direction))
466 :     {
467 :     $direction = "B";
468 :     }
469 :    
470 :     $found_in_other_map = 1;
471 :     }
472 :     if (!$found_in_other_map)
473 :     {
474 : dejongh 1.17 # reaction not in any map, get directionality without reference to map
475 : olson 1.1 if($fig->reversible($rxn) eq "1")
476 :     {
477 :     $direction = "B";
478 :     }
479 :     else
480 :     {
481 :     $direction = "R";
482 :     }
483 : dejongh 1.17
484 :     $reactions_not_in_any_map{$rxn} = 1;
485 : olson 1.1 }
486 :     }
487 :    
488 :     if (! defined $all_reactions{$rxn} || $direction eq "B")
489 :     {
490 :     $all_reactions{$rxn} = $direction;
491 :     }
492 :     elsif ($all_reactions{$rxn} ne $direction)
493 :     {
494 :     $all_reactions{$rxn} = "B";
495 :     }
496 :    
497 :     my (@substrArr, @prodArr);
498 :    
499 :     if ($direction eq "L")
500 :     {
501 :     @substrArr = $fig->reaction2comp($rxn, 1, $map_ids);
502 :     @prodArr = $fig->reaction2comp($rxn, 0, $map_ids);
503 :     $reactions_to_substrate_arrays{$rxn . "_L"} = \@substrArr;
504 :     $reactions_to_product_arrays{$rxn . "_L"} = \@prodArr;
505 :    
506 :     }
507 :     else
508 :     {
509 :     if ($reactions_in_maps{$rxn})
510 :     {
511 :     @substrArr = $fig->reaction2comp($rxn, 0, $map_ids);
512 :     @prodArr = $fig->reaction2comp($rxn, 1, $map_ids);
513 :     }
514 :     else
515 :     {
516 :     @substrArr = $fig->reaction2comp($rxn, 0);
517 :     @prodArr = $fig->reaction2comp($rxn, 1);
518 :     }
519 :    
520 :     $reactions_to_substrate_arrays{$rxn."_R"} = \@substrArr;
521 :     $reactions_to_product_arrays{$rxn."_R"} = \@prodArr;
522 :    
523 :     if ($direction eq "B")
524 :     {
525 :     $reactions_to_substrate_arrays{$rxn."_L"} = \@prodArr;
526 :     $reactions_to_product_arrays{$rxn."_L"} = \@substrArr;
527 :     }
528 :     }
529 :    
530 :     print STDERR "\nFor $rxn, found substrates:\n" unless !$debug;
531 :     map { print STDERR "\t$_->[0]\t$_->[1]\t$_->[2]\n" unless !$debug } @substrArr;
532 :     print STDERR "For $rxn, found products:\n" unless !$debug;
533 :     map { print STDERR "\t$_->[0]\t$_->[1]\t$_->[2]\n" unless !$debug } @prodArr;
534 :    
535 :     # load "main" designation based on reactions that are in the maps
536 :     if ($reactions_in_maps{$rxn})
537 :     {
538 :     foreach my $cinfo ((@substrArr, @prodArr))
539 :     {
540 :     my $cpd = $cinfo->[0];
541 :     my $main = $cinfo->[2];
542 :    
543 :     if (defined $all_compounds_to_main{$cpd})
544 :     {
545 :     # only main if it's main in all reactions
546 :     $all_compounds_to_main{$cpd} &= $main;
547 :     }
548 :     else
549 :     {
550 :     $all_compounds_to_main{$cpd} = $main;
551 :     }
552 :     }
553 :     }
554 :     else
555 :     {
556 :     # save subs and prods for processing at end
557 :     $reactions_not_in_maps{$rxn} = [ (@substrArr, @prodArr) ];
558 :     }
559 :     }
560 :     }
561 :    
562 :     # now load "main" designation based on reactions not in the maps - but don't overrule
563 :     # what's already been loaded
564 :     my %additional_compounds_to_main;
565 :    
566 :     foreach my $rxn (keys %reactions_not_in_maps)
567 :     {
568 :     foreach my $cinfo (@{$reactions_not_in_maps{$rxn}})
569 :     {
570 :     my $cpd = $cinfo->[0];
571 :     my $main = $cinfo->[2];
572 :    
573 :     if (defined $additional_compounds_to_main{$cpd})
574 :     {
575 :     # main if it's main in any reactions not in map
576 :     $additional_compounds_to_main{$cpd} |= $main;
577 :     }
578 :     else
579 :     {
580 :     $additional_compounds_to_main{$cpd} = $main;
581 :     }
582 :     }
583 :     }
584 :    
585 :     foreach my $cpd (keys %additional_compounds_to_main)
586 :     {
587 :     if (! defined $all_compounds_to_main{$cpd})
588 :     {
589 :     $all_compounds_to_main{$cpd} = $additional_compounds_to_main{$cpd};
590 :     }
591 :     }
592 : dejongh 1.17
593 :     # the reactions that aren't in any map at all won't have any main compounds.
594 :     # mark those compounds that are in all_compounds_to_main as main.
595 :     foreach my $rxn (keys %reactions_not_in_any_map)
596 :     {
597 :     print STDERR "Checking reaction not in any map: $rxn\n" if $debug;
598 :     foreach my $cpd_array ($reactions_to_substrate_arrays{$rxn."_L"},
599 :     $reactions_to_substrate_arrays{$rxn."_R"},
600 :     $reactions_to_product_arrays{$rxn."_L"},
601 :     $reactions_to_product_arrays{$rxn."_R"})
602 :     {
603 :     if (defined $cpd_array)
604 :     {
605 : dejongh 1.19 my $at_least_one_is_main = 0;
606 :    
607 : dejongh 1.17 foreach my $cinfo (@{$cpd_array})
608 :     {
609 :     my $cpd = $cinfo->[0];
610 : dejongh 1.19 if ($all_compounds_to_main{$cpd} || exists $sc_inputs->{$cpd})
611 : dejongh 1.17 {
612 : dejongh 1.19 print STDERR "\t Setting $cpd to main\n" if $debug;
613 : dejongh 1.17 $cinfo->[2] = 1;
614 : dejongh 1.19 $at_least_one_is_main = 1;
615 :     }
616 :     }
617 :    
618 :     if ($at_least_one_is_main == 0)
619 :     {
620 :     foreach my $cinfo (@{$cpd_array})
621 :     {
622 :     my $cpd = $cinfo->[0];
623 :     print STDERR "\t Setting $cpd to provisional main\n" if $debug;
624 :     $cinfo->[2] = 2;
625 : dejongh 1.17 }
626 :     }
627 :     }
628 :     }
629 :     }
630 : olson 1.1 }
631 :    
632 :     sub process_paths
633 :     {
634 :     my ($input_cpds, $output_cpds, $outputs_lists, $reactions_to_substrates, $reactions_to_products, $substrates_to_reactions, $products_to_reactions, $create_assembly, $find_first ) = @_;
635 :    
636 :     my (%path_inputs, %path_outputs);
637 :     map { $path_inputs{$_} = 0 } keys %$input_cpds;
638 :     map { $path_outputs{$_} = 0 } keys %$output_cpds;
639 :    
640 :     my %data_results = ("infinite" => 0);
641 :    
642 :     # %compounds_to_tokens maps from compound ids to tokens placed on those compounds
643 :     # the tokens are organized in hashes mapping from token id to number of tokens with that id
644 :     my %compounds_to_tokens;
645 :     map { $compounds_to_tokens{$_} = {} } keys %all_compounds_to_main;
646 :    
647 :     # %tokens maps from token_ids to the token data structures
648 :     my %tokens;
649 :     my $token_id_counter = 1;
650 :    
651 :     # %compounds_borrowed_to_tokens maps from compounds to lists of token ids that borrowed
652 :     # compound in order to run a reaction
653 :     my %compounds_borrowed_to_tokens;
654 :    
655 :     print STDERR "\nIn process_paths, path_inputs are @{[ keys %path_inputs ]}, path_outputs are @{[ keys %path_outputs] }, scenario_cycles are @{[ keys %scenario_cycles ]} \n\n" unless !$debug;
656 :    
657 :     my $initial_pass = 1;
658 :     my $done = 0;
659 :     my $loop_counter = 1;
660 :     my $infinite_loop_check = $create_assembly ? $loop_max_assembly : $loop_max;
661 :    
662 :     # we may get to the point where we need to add some more path inputs into the mix
663 :     # to push stalled tokens
664 :     my $add_path_inputs = 0;
665 :    
666 :     while(!$done)
667 :     {
668 :     if ($initial_pass || $add_path_inputs)
669 :     {
670 :     foreach my $cpd (keys %path_inputs)
671 :     {
672 :     # place a token on each path input
673 :     my $new_token_id = $token_id_counter++;
674 :     my %new_token;
675 :     $new_token{visited_reactions} = {};
676 :     $new_token{visited_compounds} = { $cpd => 0 }; # 0 means supplied from "outside"
677 :     $new_token{token_path_inputs} = { $cpd => 1 }; # 1 means one was supplied
678 :     $new_token{initial_pass} = $initial_pass;
679 :     $compounds_to_tokens{$cpd}->{$new_token_id}++;
680 :     $tokens{$new_token_id} = \%new_token;
681 :    
682 :     print STDERR "\t\tCreated new token '$new_token_id' for path input $cpd\n" unless !$debug;
683 :     }
684 :    
685 :     $initial_pass = 0;
686 :     $add_path_inputs = 0;
687 :     }
688 :    
689 :     # Find the reactions that can run.
690 :     my %reactions_to_try;
691 :    
692 :     foreach my $cpd (keys %compounds_to_tokens)
693 :     {
694 :     foreach my $token_id (keys %{$compounds_to_tokens{$cpd}})
695 :     {
696 :     if ($compounds_to_tokens{$cpd}->{$token_id} > 0 &&
697 :     ! $tokens{$token_id}->{done})
698 :     {
699 :     # this compound has tokens that aren't done
700 :     map { $reactions_to_try{$_} = 1 } @{$substrates_to_reactions->{$cpd}};
701 :     last;
702 :     }
703 :     }
704 :     }
705 :    
706 :     print STDERR "\n\tIn loop, trying reactions @{[keys %reactions_to_try]}\n\n" unless !$debug;
707 :    
708 :     # Map the reactions that can run to the tokens they can use
709 :     my %reactions_to_tokens_available;
710 :     # Keep track of the main substrates defined by each reaction
711 :     my %reactions_to_main_substrates;
712 :    
713 :     # count up the total number of tokens needed for each compound to run
714 :     # every reaction that is ready to go
715 :     my %reactions_to_tokens_needed;
716 :    
717 :     rxn: foreach my $reaction (keys %reactions_to_try)
718 :     {
719 :     print STDERR "\tChecking reaction $reaction\n" unless !$debug;
720 :    
721 :     my @substrArr = @{$reactions_to_substrate_arrays{$reaction}};
722 :     my @prodArr = @{$reactions_to_product_arrays{$reaction}};
723 :     my %main_substrates;
724 :    
725 :     # Determine if this reaction has necessary inputs to run.
726 :     # There must be tokens available for at least one main substrate that isn't
727 :     # a path output (unless it's an initial scenario cycled compound).
728 :     # Also, any path input must have a token.
729 :     # Count number of tokens needed for each main substrate.
730 :     my %tokens_available;
731 :     my %tokens_needed;
732 :     my $reaction_can_run = 0;
733 :    
734 :     foreach my $substr (@substrArr)
735 :     {
736 :     my $cpd = @{$substr}[0];
737 :     my $stoich = @{$substr}[1];
738 : dejongh 1.19 my $is_it_main = @{$substr}[2];
739 : parrello 1.32 my $main = $is_it_main == 1 || $all_compounds_to_main{$cpd}; #main either way
740 : dejongh 1.19
741 :     if ($is_it_main == 2)
742 :     {
743 :     # provisionally main - need to see if there is a token that has visited
744 :     # this compound
745 :    
746 :     foreach my $token_id (keys %{$compounds_to_tokens{$cpd}})
747 :     {
748 :     my %visited_compounds = %{$tokens{$token_id}->{visited_compounds}};
749 :     $main = 1 if exists $visited_compounds{$cpd};
750 :     }
751 :     }
752 :    
753 : olson 1.1 $main_substrates{$cpd} = 1 if $main;
754 :     print STDERR "\t\tSubstrate: $cpd\tstoich: $stoich\tmain: $main\n" unless !$debug;
755 :    
756 :     if (! $main)
757 :     {
758 :     # on any pass we can take in non-main compounds
759 :     $tokens_needed{$cpd} = $stoich;
760 :     }
761 :     else
762 :     {
763 :     # if tokens are available for compound, check their history for the
764 :     # main compounds produced by the reaction so we don't loop back over
765 :     # previous main compounds. Also, don't use tokens on scenario inputs
766 :     # moved to in steps other than initial token creation.
767 :     my %ok_tokens;
768 :     my $num_ok_tokens = 0;
769 :    
770 :     foreach my $token_id (keys %{$compounds_to_tokens{$cpd}})
771 :     {
772 :     next if $tokens{$token_id}->{done} ||
773 :     $compounds_to_tokens{$cpd}->{$token_id} == 0;
774 :    
775 :     # now check that we aren't running an already visited reaction in reverse
776 :     my %visited_reactions = %{$tokens{$token_id}->{visited_reactions}};
777 :    
778 :     if (($reaction =~ /(.*)_R/ && defined $visited_reactions{$1."_L"}) ||
779 :     ($reaction =~ /(.*)_L/ && defined $visited_reactions{$1."_R"}))
780 :     {
781 :     print STDERR "\t\tToken '$token_id' has run the reverse reaction already\n" unless !$debug;
782 :     next;
783 :     }
784 :    
785 :     my %visited_compounds = %{$tokens{$token_id}->{visited_compounds}};
786 :    
787 :     print STDERR "\t\tToken '$token_id' has visited @{[map { ($_, $visited_compounds{$_} ) } sort { $visited_compounds{$a} <=> $visited_compounds{$b} } keys %visited_compounds]}\n" unless !$debug;
788 :    
789 :     # need to find at least one prod that hasn't been visited yet
790 :     # or was visited in a loop cycle not before the loop cycle in which
791 :     # it visited the substrate,
792 :     # or is a path output
793 :     my $prods_are_ok = 0;
794 :    
795 :     # check each main product
796 :     foreach my $prod (@{$reactions_to_products->{$reaction}})
797 :     {
798 :     if (! defined $visited_compounds{$prod} ||
799 :     $visited_compounds{$prod} >= $visited_compounds{$cpd} ||
800 :     defined $path_outputs{$prod} ||
801 :     $compounds_borrowed_to_tokens{$prod}->{$token_id} > 0)
802 :     {
803 :     print STDERR "\t\tToken can visit $prod\n" unless !$debug;
804 :     $prods_are_ok = 1;
805 :     last;
806 :     }
807 :     }
808 :    
809 :     if ($prods_are_ok)
810 :     {
811 :     print STDERR "\t\tToken is OK\n" unless !$debug;
812 :     $ok_tokens{$token_id} = $compounds_to_tokens{$cpd}->{$token_id};
813 :     $num_ok_tokens += $compounds_to_tokens{$cpd}->{$token_id};
814 :     }
815 :     }
816 :    
817 :     map { $tokens_available{$_}->{$cpd} = $ok_tokens{$_} } keys %ok_tokens;
818 :     $tokens_needed{$cpd} = $stoich;
819 :    
820 :     if ($main && $num_ok_tokens >= 1)
821 :     {
822 :     if (! defined $path_outputs{$cpd} || defined $scenario_cycles{$cpd})
823 :     {
824 :     print STDERR "\t\tgot at least one token on main compound\n" unless !$debug;
825 :     $reaction_can_run = 1;
826 :     }
827 :     }
828 :     elsif (defined $path_inputs{$cpd})
829 :     {
830 :     print STDERR "\t\tno tokens available for path input: $cpd\n" unless !$debug;
831 :     next rxn;
832 :     }
833 :     elsif (defined $path_outputs{$cpd} && ! defined $scenario_cycles{$cpd})
834 :     {
835 :     print STDERR "\t\tno tokens available for path output: $cpd\n" unless !$debug;
836 :     next rxn;
837 :     }
838 :     }
839 :     }
840 :    
841 :     if ($reaction_can_run)
842 :     {
843 :     $reactions_to_tokens_available{$reaction} = \%tokens_available;
844 :     $reactions_to_tokens_needed{$reaction} = \%tokens_needed;
845 :     $reactions_to_main_substrates{$reaction} = \%main_substrates;
846 :     print STDERR "\tReaction $reaction can run\n" unless !$debug;
847 :     }
848 :     }
849 :    
850 :     # keep track of tokens used that will be used to run rxns. Clone tokens if necessary.
851 :     my %reactions_to_tokens_to_use;
852 :     my %tokens_to_use_to_reactions;
853 :     my %copy_of_compounds_to_tokens; # for determining which reaction uses which tokens
854 :    
855 :     foreach my $cpd (keys %compounds_to_tokens)
856 :     {
857 :     my %cpd_token_ids = %{$compounds_to_tokens{$cpd}};
858 :     my %new_cpd_token_ids;
859 :     map { $new_cpd_token_ids{$_} = $cpd_token_ids{$_} } keys %cpd_token_ids;
860 :     $copy_of_compounds_to_tokens{$cpd} = \%new_cpd_token_ids;
861 :     }
862 :    
863 :     foreach my $reaction (keys %reactions_to_tokens_available)
864 :     {
865 :     print STDERR "\tPreparing to run reaction $reaction\n" unless !$debug;
866 :    
867 :     # assemble tokens to run reaction, cloning ones that were used by
868 :     # other reactions during this cycle if necessary
869 :     my %tokens_available = %{$reactions_to_tokens_available{$reaction}};
870 :     my %tokens_needed = %{$reactions_to_tokens_needed{$reaction}};
871 :     my @final_tokens_to_use; # list of maps from substrates to tokens to use
872 :     my %clone_history; # map from token ids to ids of their new clones
873 :    
874 :     # check to see if any available tokens are already commited to the reverse reaction
875 :     my @token_ids = keys %tokens_available;
876 :    
877 :     foreach my $token_id (keys %tokens_available)
878 :     {
879 :     if (($reaction =~ /(.*)_R/ && defined $tokens_to_use_to_reactions{$token_id}->{$1."_L"}) ||
880 :     ($reaction =~ /(.*)_L/ && defined $tokens_to_use_to_reactions{$token_id}->{$1."_R"}))
881 :     {
882 :     # clone the token
883 :     my $new_token_id = $token_id_counter++;
884 :     &clone_token($token_id, $new_token_id, \%tokens, \%compounds_to_tokens,
885 :     \%compounds_borrowed_to_tokens);
886 :    
887 :     foreach my $icpd (keys %compounds_to_tokens)
888 :     {
889 :     if ($compounds_to_tokens{$icpd}->{$token_id} > 0)
890 :     {
891 :     $copy_of_compounds_to_tokens{$icpd}->{$new_token_id} =
892 :     $compounds_to_tokens{$icpd}->{$token_id};
893 :     }
894 :     }
895 :    
896 :     print STDERR "\t\tCloned token '$token_id', new token is '$new_token_id'\n" unless !$debug;
897 :    
898 :     $tokens_available{$new_token_id} = $tokens_available{$token_id};
899 :     delete $tokens_available{$token_id};
900 :     }
901 :     }
902 :    
903 :     # first, assemble tokens that have all the main compounds they need to run
904 :     @token_ids = keys %tokens_available;
905 :    
906 :     foreach my $token_id (keys %tokens_available)
907 :     {
908 :     my $has_all_main_cpds = 1;
909 :     my $need_to_clone = 0;
910 :    
911 :     foreach my $cpd (keys %tokens_needed)
912 :     {
913 :     if ($reactions_to_main_substrates{$reaction}->{$cpd})
914 :     {
915 :     if ($compounds_to_tokens{$cpd}->{$token_id} >= $tokens_needed{$cpd})
916 :     {
917 :     if ($copy_of_compounds_to_tokens{$cpd}->{$token_id} < $tokens_needed{$cpd})
918 :     {
919 :     $need_to_clone = 1;
920 :     }
921 :     }
922 :     else
923 :     {
924 :     $has_all_main_cpds = 0;
925 :     last;
926 :     }
927 :     }
928 :     }
929 :    
930 :     if ($has_all_main_cpds)
931 :     {
932 :     print STDERR "\t\ttoken '$token_id' has all main compounds\n" unless !$debug;
933 :    
934 :     delete $tokens_available{$token_id};
935 :    
936 :     if ($need_to_clone)
937 :     {
938 :     my $new_token_id = $token_id_counter++;
939 :     &clone_token($token_id, $new_token_id, \%tokens,
940 :     \%compounds_to_tokens,
941 :     \%compounds_borrowed_to_tokens);
942 :    
943 :     foreach my $icpd (keys %compounds_to_tokens)
944 :     {
945 :     if ($compounds_to_tokens{$icpd}->{$token_id} > 0)
946 :     {
947 :     $copy_of_compounds_to_tokens{$icpd}->{$new_token_id} =
948 :     $compounds_to_tokens{$icpd}->{$token_id};
949 :     }
950 :     }
951 :    
952 :     print STDERR "\t\tCloned token '$token_id', new token is '$new_token_id'\n" unless !$debug;
953 :    
954 :     $clone_history{$token_id} = $new_token_id;
955 :     $token_id = $new_token_id;
956 :     }
957 :    
958 :     # now assemble the compound to tokens map
959 :     my %cpd_to_tokens;
960 :    
961 :     foreach my $cpd (keys %tokens_needed)
962 :     {
963 :     if ($reactions_to_main_substrates{$reaction}->{$cpd})
964 :     {
965 :     for (my $i = 0; $i < $tokens_needed{$cpd}; $i++)
966 :     {
967 :     push @{$cpd_to_tokens{$cpd}}, $token_id;
968 :     $copy_of_compounds_to_tokens{$cpd}->{$token_id}--;
969 :     }
970 :     }
971 :     }
972 :    
973 :     push @final_tokens_to_use, \%cpd_to_tokens;
974 :     }
975 :     }
976 :    
977 :     if (scalar keys %tokens_available > 0)
978 :     {
979 :     # try to merge left over available tokens in all combinations that fulfill
980 :     # needed main substrates. Create map from main substrates to tokens.
981 :     my %tokens_available_for_cpds;
982 :     my %tokens_to_use_for_cpds;
983 :    
984 :     foreach my $token_id (keys %tokens_available)
985 :     {
986 :     foreach my $cpd (keys %{$tokens_available{$token_id}})
987 :     {
988 :     if ($reactions_to_main_substrates{$reaction}->{$cpd})
989 :     {
990 :     $tokens_available_for_cpds{$cpd}->{$token_id} =
991 :     $tokens_available{$token_id}->{$cpd};
992 :     }
993 :     }
994 :     }
995 :    
996 :     foreach my $cpd (keys %{$reactions_to_main_substrates{$reaction}})
997 :     {
998 :     my %available_token_ids;
999 :    
1000 :     foreach my $token_id (keys %{$tokens_available_for_cpds{$cpd}})
1001 :     {
1002 :     my $updated_token_id = $token_id;
1003 :    
1004 :     # in case the token has already been cloned for another cpd in this reaction
1005 :     while ($clone_history{$updated_token_id})
1006 :     {
1007 :     $updated_token_id = $clone_history{$updated_token_id};
1008 :     }
1009 :    
1010 :     $available_token_ids{$updated_token_id} = 1;
1011 :    
1012 :     if ($token_id != $updated_token_id)
1013 :     {
1014 :     $tokens_available_for_cpds{$cpd}->{$updated_token_id} =
1015 :     $tokens_available_for_cpds{$cpd}->{$token_id};
1016 :     delete $tokens_available_for_cpds{$cpd}->{$token_id};
1017 :     }
1018 :     }
1019 :    
1020 :     my $num_tokens_needed = $tokens_needed{$cpd};
1021 :     my $num_available_tokens = 0;
1022 :    
1023 :     foreach my $token_id (keys %available_token_ids)
1024 :     {
1025 :     $num_available_tokens += $tokens_available_for_cpds{$cpd}->{$token_id};
1026 :     }
1027 :    
1028 :     # if not enough tokens are available to fill out sets, create new ones
1029 :     my $num_short_of_full;
1030 :    
1031 :     if ($num_available_tokens == 0)
1032 :     {
1033 :     $num_short_of_full = $num_tokens_needed;
1034 :     }
1035 :     else
1036 :     {
1037 :     $num_short_of_full = ($num_tokens_needed - ($num_available_tokens % $num_tokens_needed)) % $num_tokens_needed;
1038 :     }
1039 :    
1040 :     print STDERR "\t\tFor $cpd, $num_tokens_needed tokens are needed, $num_available_tokens are available. Need to create $num_short_of_full to fill out sets\n" unless !$debug;
1041 :    
1042 :     for (my $i = 0; $i < $num_short_of_full; $i++)
1043 :     {
1044 :     my $new_token_id = $token_id_counter++;
1045 :     my %new_token;
1046 :     print STDERR "\t\tCreating new token '$new_token_id' for $cpd\n" unless !$debug;
1047 :    
1048 :     $new_token{visited_reactions} = {};
1049 :     $new_token{visited_compounds} = {};
1050 :     $compounds_to_tokens{$cpd}->{$new_token_id}++;
1051 :     $copy_of_compounds_to_tokens{$cpd}->{$new_token_id}++;
1052 :     $tokens{$new_token_id} = \%new_token;
1053 :     $available_token_ids{$new_token_id} = 1;
1054 :     $tokens_available_for_cpds{$cpd}->{$new_token_id}++;
1055 :    
1056 :     # if it's not a path input, remember that we've "borrowed" it and
1057 :     # will need to pay it back.
1058 :     if (! defined $path_inputs{$cpd})
1059 :     {
1060 :     $new_token{visited_compounds}->{$cpd} = $loop_counter;
1061 :     $compounds_borrowed_to_tokens{$cpd}->{$new_token_id}++;
1062 :     }
1063 :     else
1064 :     {
1065 :     $new_token{token_path_inputs} = { $cpd => 1 };
1066 :     $new_token{visited_compounds}->{$cpd} = 0;
1067 :     }
1068 :     }
1069 :    
1070 :     # for main compounds, there may be more tokens available than needed,
1071 :     # so we may assemble multiple token sets.
1072 :     my %tokens_not_yet_used;
1073 :    
1074 :     foreach my $token_id (keys %available_token_ids)
1075 :     {
1076 :     if ($copy_of_compounds_to_tokens{$cpd}->{$token_id} > 0)
1077 :     {
1078 :     $tokens_not_yet_used{$token_id} = $tokens_available_for_cpds{$cpd}->{$token_id};
1079 :     }
1080 :     }
1081 :    
1082 :     my @token_sets_for_cpd;
1083 :    
1084 :     print STDERR "\t\tNeed $num_tokens_needed tokens for $cpd, '@{[ sort { $a <=> $b } keys %available_token_ids ]}' are usable, '@{[ sort { $a <=> $b } keys %tokens_not_yet_used ]}' are not yet used\n" unless !$debug;
1085 :    
1086 :     my @token_set;
1087 :     my $num_tokens_still_needed = $num_tokens_needed;
1088 :    
1089 :     foreach my $token_id (sort { $a <=> $b } keys %available_token_ids)
1090 :     {
1091 :     # need to clone the token if it is all used up
1092 :     if ($copy_of_compounds_to_tokens{$cpd}->{$token_id} == 0)
1093 :     {
1094 :     my $new_token_id = $token_id_counter++;
1095 :     &clone_token($token_id, $new_token_id, \%tokens,
1096 :     \%compounds_to_tokens,
1097 :     \%compounds_borrowed_to_tokens);
1098 :    
1099 :     foreach my $icpd (keys %compounds_to_tokens)
1100 :     {
1101 :     if ($compounds_to_tokens{$icpd}->{$token_id} > 0)
1102 :     {
1103 :     $copy_of_compounds_to_tokens{$icpd}->{$new_token_id} =
1104 :     $compounds_to_tokens{$icpd}->{$token_id};
1105 :     }
1106 :     }
1107 :    
1108 :     print STDERR "\t\tCloned token '$token_id' for $cpd, new token is '$new_token_id'\n" unless !$debug;
1109 :    
1110 :     $clone_history{$token_id} = $new_token_id;
1111 :     $token_id = $new_token_id;
1112 :     }
1113 :    
1114 :     while ($copy_of_compounds_to_tokens{$cpd}->{$token_id} > 0)
1115 :     {
1116 :     push @token_set, $token_id;
1117 :     $copy_of_compounds_to_tokens{$cpd}->{$token_id}--;
1118 :     $num_tokens_still_needed--;
1119 :    
1120 :     if ($num_tokens_still_needed == 0)
1121 :     {
1122 :     print STDERR "\t\tPushing token set '@token_set' for $cpd\n" unless !$debug;
1123 :     my @copy_of_token_set = @token_set;
1124 :     push @token_sets_for_cpd, \@copy_of_token_set;
1125 :     @token_set = ();
1126 :     $num_tokens_still_needed = $num_tokens_needed;
1127 :     }
1128 :     }
1129 :     }
1130 :    
1131 :     $tokens_to_use_for_cpds{$cpd} = \@token_sets_for_cpd;
1132 :     }
1133 :    
1134 :     # in case a token had to be cloned for this reaction after another
1135 :     # compound already determined to use it, check history and use new token id
1136 :     foreach my $cpd (keys %tokens_to_use_for_cpds)
1137 :     {
1138 :     foreach my $token_set (@{$tokens_to_use_for_cpds{$cpd}})
1139 :     {
1140 :     for (my $i = 0; $i < scalar @$token_set; $i++)
1141 :     {
1142 :     my $updated_token_id = $token_set->[$i];
1143 :    
1144 :     while ($clone_history{$updated_token_id})
1145 :     {
1146 :     $updated_token_id = $clone_history{$updated_token_id};
1147 :     }
1148 :    
1149 :     if ($token_set->[$i] != $updated_token_id)
1150 :     {
1151 :     print STDERR "\t\tReplacing '$token_set->[$i]' with '$updated_token_id' for $cpd\n" unless !$debug;
1152 :     splice @$token_set, $i, 1, ($updated_token_id);
1153 :     }
1154 :     }
1155 :     }
1156 :     }
1157 :    
1158 :     # we may have multiple sets of tokens to use for some compounds, so
1159 :     # make sure we're prepared for each combination by cloning sets as necessary
1160 :     my %cpd_to_token_set_index;
1161 :     my $num_combinations = 1;
1162 :     my %used_token_sets_for_cpd;
1163 :    
1164 :     foreach my $cpd (keys %tokens_to_use_for_cpds)
1165 :     {
1166 :     $cpd_to_token_set_index{$cpd} = 0;
1167 :     $num_combinations *= scalar @{$tokens_to_use_for_cpds{$cpd}};
1168 :     }
1169 :    
1170 :     for (my $i = 0; $i < $num_combinations; $i++)
1171 :     {
1172 :     my %combination_cpds_to_tokens;
1173 :    
1174 :     print STDERR "\t\tPreparing combination ", $i+1, " out of $num_combinations, indices are @{[ map { $cpd_to_token_set_index{$_} } sort keys %cpd_to_token_set_index ]}\n" unless !$debug;
1175 :    
1176 :     # keep track of who was cloned to what for this combination in case
1177 :     # different compounds are using the same tokens
1178 :     my %clone_history_this_combination;
1179 :    
1180 :     foreach my $cpd (sort keys %cpd_to_token_set_index)
1181 :     {
1182 :     my $token_set_index = $cpd_to_token_set_index{$cpd};
1183 :     my @token_set = @{$tokens_to_use_for_cpds{$cpd}->[$token_set_index]};
1184 :     my @new_token_set;
1185 :    
1186 :     # don't clone if this is the first time using this token set
1187 :     # should be able to do this mathematically
1188 :     if (! defined $used_token_sets_for_cpd{$cpd}->{$token_set_index})
1189 :     {
1190 :     @new_token_set = @token_set;
1191 :     $used_token_sets_for_cpd{$cpd}->{$token_set_index} = 1;
1192 :     }
1193 :     else
1194 :     {
1195 :     foreach my $token_id (@token_set)
1196 :     {
1197 :     my $new_token_id;
1198 :    
1199 :     if ($clone_history_this_combination{$token_id})
1200 :     {
1201 :     $new_token_id = $clone_history_this_combination{$token_id};
1202 :     }
1203 :     else
1204 :     {
1205 :     $new_token_id = $token_id_counter++;
1206 :     &clone_token($token_id, $new_token_id, \%tokens,
1207 :     \%compounds_to_tokens,
1208 :     \%compounds_borrowed_to_tokens);
1209 :    
1210 :     print STDERR "\t\t\tCloned token '$token_id' for $cpd, new token is '$new_token_id'\n" unless !$debug;
1211 :    
1212 :     $clone_history_this_combination{$token_id} = $new_token_id;
1213 :     }
1214 :    
1215 :     push @new_token_set, $new_token_id;
1216 :     }
1217 :     }
1218 :    
1219 :     push @{$combination_cpds_to_tokens{$cpd}}, @new_token_set;
1220 :     }
1221 :    
1222 :     push @final_tokens_to_use, \%combination_cpds_to_tokens;
1223 :    
1224 :     # move the cpd to token set indices for the next combination
1225 :     foreach my $cpd (sort keys %cpd_to_token_set_index)
1226 :     {
1227 :     if ($cpd_to_token_set_index{$cpd} < scalar @{$tokens_to_use_for_cpds{$cpd}} - 1)
1228 :     {
1229 :     $cpd_to_token_set_index{$cpd}++;
1230 :     last;
1231 :     }
1232 :     else
1233 :     {
1234 :     $cpd_to_token_set_index{$cpd} = 0;
1235 :     }
1236 :     }
1237 :     }
1238 :     }
1239 :    
1240 :     # last step, need to create place-holder tokens for the non-main substrates
1241 :     # and record the tokens being used for this reaction
1242 :     foreach my $combination (@final_tokens_to_use)
1243 :     {
1244 :     foreach my $cpd (keys %tokens_needed)
1245 :     {
1246 :     if (! defined $reactions_to_main_substrates{$reaction}->{$cpd})
1247 :     {
1248 :     my $new_token_id = $token_id_counter++;
1249 :     my %new_token;
1250 :     $new_token{visited_reactions} = {};
1251 :     $new_token{visited_compounds} = {};
1252 :     $compounds_to_tokens{$cpd}->{$new_token_id} += $tokens_needed{$cpd};
1253 :     $copy_of_compounds_to_tokens{$cpd}->{$new_token_id} += $tokens_needed{$cpd};
1254 :     $new_token{token_path_inputs} = { $cpd => $tokens_needed{$cpd} };
1255 :     $tokens{$new_token_id} = \%new_token;
1256 :    
1257 :     for (my $i = 0; $i < $tokens_needed{$cpd}; $i++)
1258 :     {
1259 :     push @{$combination->{$cpd}}, $new_token_id;
1260 :     }
1261 :     }
1262 :     }
1263 :    
1264 :     foreach my $cpd (keys %$combination)
1265 :     {
1266 :     map { $tokens_to_use_to_reactions{$_}->{$reaction} = 1 } @{$combination->{$cpd}};
1267 :     }
1268 :     }
1269 :    
1270 :     $reactions_to_tokens_to_use{$reaction} = \@final_tokens_to_use;
1271 :     }
1272 :    
1273 :     # keep track of tokens merged during this round
1274 :     my %token_merge_history;
1275 :    
1276 :     foreach my $reaction (keys %reactions_to_tokens_to_use)
1277 :     {
1278 :     # we may have multiple sets of tokens to use
1279 :     my @final_tokens_to_use = @{$reactions_to_tokens_to_use{$reaction}};
1280 :     my $num_combinations = scalar @final_tokens_to_use;
1281 :    
1282 :     # since we may be running this reaction several times with some of the same
1283 :     # tokens, don't process token merge history until all sets have been run
1284 :     my %token_merge_history_this_reaction;
1285 :    
1286 :     for (my $i = 0; $i < $num_combinations; $i++)
1287 :     {
1288 :     print STDERR "\tRunning reaction $reaction (", $i+1, " out of $num_combinations)\n" unless !$debug;
1289 :    
1290 :     # assemble list of tokens to use for this combination
1291 :     my %tokens_to_use;
1292 :     # find the most recently visited compound's step
1293 :     my $most_recent_step = 0;
1294 :    
1295 :     # loop through substrates
1296 :     foreach my $cpd (keys %{$reactions_to_tokens_to_use{$reaction}})
1297 :     {
1298 :     # remove the token ids from the real compounds_to_tokens map
1299 :     my @cpd_token_ids = @{$final_tokens_to_use[$i]->{$cpd}};
1300 :     print STDERR "\t\tFound tokens '@cpd_token_ids' for $cpd\n" unless !$debug;
1301 :    
1302 :     foreach my $token_id (@cpd_token_ids)
1303 :     {
1304 :     # check if token_id has been merged into a new token by a previous reaction
1305 :     while (defined $token_merge_history{$token_id})
1306 :     {
1307 :     $token_id = $token_merge_history{$token_id};
1308 :     }
1309 :    
1310 :     $compounds_to_tokens{$cpd}->{$token_id}--;
1311 :    
1312 :     if ($compounds_to_tokens{$cpd}->{$token_id} == 0)
1313 :     {
1314 :     delete $compounds_to_tokens{$cpd}->{$token_id};
1315 :     }
1316 :    
1317 :     $tokens_to_use{$token_id} = 1;
1318 :    
1319 :     if ($tokens{$token_id}->{visited_compounds}->{$cpd} > $most_recent_step)
1320 :     {
1321 :     $most_recent_step = $tokens{$token_id}->{visited_compounds}->{$cpd};
1322 :     }
1323 :     }
1324 :     }
1325 :    
1326 :     # process list of unique token ids
1327 :     my @tokens_to_use = sort { $a <=> $b } keys %tokens_to_use;
1328 :     my $go_forward_token_id;
1329 :    
1330 :     if (scalar @tokens_to_use == 1)
1331 :     {
1332 :     $go_forward_token_id = shift @tokens_to_use;
1333 :     print STDERR "\t\tGoing forward with '$go_forward_token_id'\n" unless !$debug;
1334 :     }
1335 :     else
1336 :     {
1337 :     $go_forward_token_id = $token_id_counter++;
1338 :     $tokens{$go_forward_token_id} = {};
1339 :     print STDERR "\t\tRemember to merge tokens '@tokens_to_use' into '$go_forward_token_id'\n" unless !$debug;
1340 :    
1341 :     # record the need to merge - we'll do it after processing all sets of substrates
1342 :     foreach my $token_id (@tokens_to_use)
1343 :     {
1344 :     push @{$token_merge_history_this_reaction{$token_id}}, $go_forward_token_id;
1345 :     }
1346 :     }
1347 :    
1348 :     my $go_forward_token = $tokens{$go_forward_token_id};
1349 :     my @prodArr = @{$reactions_to_product_arrays{$reaction}};
1350 :    
1351 :     # add current reaction and products to accumulated token history.
1352 :     # reaction is mapped to loop counter to maintain history of order of execution
1353 :     $go_forward_token->{visited_reactions}->{$reaction} = $loop_counter;
1354 :    
1355 :     foreach my $prod (@prodArr)
1356 :     {
1357 :     my $cpd = @{$prod}[0];
1358 :     my $stoich = @{$prod}[1];
1359 :     my $main = @{$prod}[2];
1360 :    
1361 :     # keep track of path outputs we've seen
1362 :     if (defined $path_outputs{$cpd})
1363 :     {
1364 :     $path_outputs{$cpd} += $stoich;
1365 :     }
1366 :    
1367 :     if ($main)
1368 :     {
1369 :     $go_forward_token->{visited_compounds}->{$cpd} = $most_recent_step + 1;
1370 :     }
1371 :    
1372 :     # push tokens
1373 :     $compounds_to_tokens{$cpd}->{$go_forward_token_id} += $stoich;
1374 :     }
1375 :     }
1376 :    
1377 :     # now process this reaction's token merge history
1378 :     foreach my $token_id (sort { $a <=> $b } keys %token_merge_history_this_reaction)
1379 :     {
1380 :     # find the unique set of up to date merge ids
1381 :     my @merge_ids = @{$token_merge_history_this_reaction{$token_id}};
1382 :     my %updated_merge_ids;
1383 :    
1384 :     foreach my $token_id (@merge_ids)
1385 :     {
1386 :     while (defined $token_merge_history{$token_id})
1387 :     {
1388 :     $token_id = $token_merge_history{$token_id};
1389 :     }
1390 :    
1391 :     $updated_merge_ids{$token_id} = 1;
1392 :     }
1393 :    
1394 :     print STDERR "\t\tupdated merge id list for '$token_id': '@{[ keys %updated_merge_ids ]}'\n" unless !$debug;
1395 :    
1396 :     my $wrap_up_token_id = $token_id_counter++;
1397 :     $tokens{$wrap_up_token_id} = {};
1398 :     my $wrap_up_token = $tokens{$wrap_up_token_id};
1399 :     my @tokens_to_merge;
1400 :     push @tokens_to_merge, keys %updated_merge_ids, $token_id;
1401 :    
1402 :     # merge from oldest to youngest to update visited_compounds history
1403 :     foreach my $itoken_id (sort { $a <=> $b } @tokens_to_merge)
1404 :     {
1405 :     print STDERR "\t\t\tmerging '$itoken_id' into '$wrap_up_token_id'\n" unless !$debug;
1406 :    
1407 :     my $itoken = $tokens{$itoken_id};
1408 :     map { $wrap_up_token->{visited_reactions}->{$_} = $itoken->{visited_reactions}->{$_} } keys %{$itoken->{visited_reactions}};
1409 :     map { $wrap_up_token->{visited_compounds}->{$_} = $itoken->{visited_compounds}->{$_} } keys %{$itoken->{visited_compounds}};
1410 :     map { $wrap_up_token->{token_path_inputs}->{$_} += $itoken->{token_path_inputs}->{$_} } keys %{$itoken->{token_path_inputs}};
1411 :     $wrap_up_token->{initial_pass} |= $itoken->{initial_pass};
1412 :     $token_merge_history{$itoken_id} = $wrap_up_token_id;
1413 :    
1414 :     # tokens might be spread across multiple compounds; change them all
1415 :     # to new id
1416 :     foreach my $cpd (keys %compounds_to_tokens)
1417 :     {
1418 :     if ($compounds_to_tokens{$cpd}->{$itoken_id} > 0)
1419 :     {
1420 :     $compounds_to_tokens{$cpd}->{$wrap_up_token_id} +=
1421 :     $compounds_to_tokens{$cpd}->{$itoken_id};
1422 :     }
1423 :     }
1424 :    
1425 :     # tokens might have borrowed compounds; change them all to new id
1426 :     foreach my $cpd (keys %compounds_borrowed_to_tokens)
1427 :     {
1428 :     if ($compounds_borrowed_to_tokens{$cpd}->{$itoken_id} > 0)
1429 :     {
1430 :     $compounds_borrowed_to_tokens{$cpd}->{$wrap_up_token_id} +=
1431 :     $compounds_borrowed_to_tokens{$cpd}->{$itoken_id};
1432 :     }
1433 :     }
1434 :     }
1435 :     }
1436 :     }
1437 :    
1438 :     # now delete the tokens that were used merged in these reactions
1439 :     foreach my $token_id (keys %token_merge_history)
1440 :     {
1441 :     foreach my $icpd (keys %compounds_to_tokens)
1442 :     {
1443 :     delete $compounds_to_tokens{$icpd}->{$token_id};
1444 :     }
1445 :    
1446 :     foreach my $icpd (keys %compounds_borrowed_to_tokens)
1447 :     {
1448 :     delete $compounds_borrowed_to_tokens{$icpd}->{$token_id};
1449 :     }
1450 :    
1451 :     print STDERR "\t\tDeleting token '$token_id'\n" unless !$debug;
1452 :     delete $tokens{$token_id};
1453 :     }
1454 :    
1455 :     print STDERR "\nBalancing tokens\n" unless !$debug;
1456 :    
1457 :     foreach my $token_id (keys %tokens)
1458 :     {
1459 :     &balance_borrowing_and_giving($token_id, \%compounds_to_tokens,
1460 :     \%compounds_borrowed_to_tokens);
1461 :     }
1462 :    
1463 :     &print_token_status([sort { $a <=> $b } keys %tokens], \%tokens, \%compounds_to_tokens, \%compounds_borrowed_to_tokens, $fig);
1464 :    
1465 :     print STDERR "\nChecking for done\n" unless !$debug;
1466 :    
1467 :     print STDERR "\n\ntoken ids: @{[ sort { $a <=> $b } map { $_ if ! defined $tokens{$_}->{done} } keys %tokens ]}\n" unless !$debug;
1468 :    
1469 :     # we're done when all the main compounds in initial-pass tokens
1470 :     # have reached path outputs and repaid their borrowed tokens,
1471 :     # or have reached a dead end.
1472 :     # Check if we're done pushing and borrowing tokens first.
1473 :     my %not_done_tokens;
1474 :    
1475 :     foreach my $token_id (keys %tokens)
1476 :     {
1477 :     if (! defined $tokens{$token_id}->{done})
1478 :     {
1479 :     if (! &check_token_for_done($token_id, \%compounds_to_tokens,
1480 :     \%compounds_borrowed_to_tokens,
1481 :     \%all_compounds_to_main, \%path_outputs,
1482 :     \%scenario_cycles, \%tokens, $fig, $outputs_lists))
1483 :     {
1484 :     $not_done_tokens{$token_id} = 1;
1485 :     }
1486 :     }
1487 :     }
1488 :    
1489 :     print STDERR "\nChecking if we can pay back borrowed compounds from other tokens\n" unless !$debug;
1490 :    
1491 :     foreach my $bcpd (keys %compounds_borrowed_to_tokens)
1492 :     {
1493 :     next if $path_outputs{$bcpd}; # tokens must manage their own path outputs
1494 :    
1495 :     if (scalar keys %{$compounds_borrowed_to_tokens{$bcpd}} > 0 &&
1496 :     scalar keys %{$compounds_to_tokens{$bcpd}} > 0)
1497 :     {
1498 :     my %borrowers_to_givers;
1499 :    
1500 :     foreach my $borrower_id (keys %{$compounds_borrowed_to_tokens{$bcpd}})
1501 :     {
1502 :     next if defined $tokens{$borrower_id}->{done}; # don't repay deadenders
1503 :    
1504 :     my $num_needed = $compounds_borrowed_to_tokens{$bcpd}->{$borrower_id};
1505 :    
1506 :     giver: foreach my $giver_id (keys %{$compounds_to_tokens{$bcpd}})
1507 :     {
1508 :     next if defined $tokens{$giver_id}->{done};
1509 :    
1510 :     my $num_to_give = $compounds_to_tokens{$bcpd}->{$giver_id};
1511 :    
1512 :     print STDERR "\tToken '$giver_id' has $num_to_give $bcpd to give to '$borrower_id', which needs $num_needed\n" unless !$debug;
1513 :    
1514 :     my $borrower = $tokens{$borrower_id};
1515 :     my $giver = $tokens{$giver_id};
1516 :    
1517 :     # check whether the giver and borrower have conflicting histories
1518 :     foreach my $visited_reaction (keys %{$giver->{visited_reactions}})
1519 :     {
1520 :     if (($visited_reaction =~ /(.*)_R/ &&
1521 :     defined $borrower->{visited_reactions}->{$1."_L"}) ||
1522 :     ($visited_reaction =~ /(.*)_L/ &&
1523 :     defined $borrower->{visited_reactions}->{$1."_R"}))
1524 :     {
1525 :     print STDERR "\t\tConflict on $visited_reaction\n" unless !$debug;
1526 :     next giver;
1527 :     }
1528 :    
1529 :     }
1530 :    
1531 :     push @{$borrowers_to_givers{$borrower_id}}, $giver_id;
1532 :     }
1533 :     }
1534 :    
1535 :     # we have a list of potential givers for each borrower for this compound.
1536 :     # Now figure out who the lucky givers will be.
1537 :     my %givers_to_borrowers;
1538 :    
1539 :     foreach my $borrower_id (keys %borrowers_to_givers)
1540 :     {
1541 :     my %lucky_givers;
1542 :     my $num_needed = $compounds_borrowed_to_tokens{$bcpd}->{$borrower_id};
1543 :    
1544 :     print STDERR "\tCollecting lucky givers for '$borrower_id' for $bcpd, need $num_needed\n" unless !$debug;
1545 :    
1546 :     # check potential givers starting with those with the most
1547 :     my @potential_givers = reverse sort { $compounds_to_tokens{$bcpd}->{$a} <=> $compounds_to_tokens{$bcpd}->{$b} } @{$borrowers_to_givers{$borrower_id}};
1548 :    
1549 :     foreach my $giver_id (@potential_givers)
1550 :     {
1551 :     if (! defined $lucky_givers{$giver_id})
1552 :     {
1553 :     my $num_to_give = $compounds_to_tokens{$bcpd}->{$giver_id};
1554 :     print STDERR "\t\t'$giver_id' has $num_to_give to give\n" unless !$debug;
1555 :     $lucky_givers{$giver_id} = 1;
1556 :     $num_needed -= $num_to_give;
1557 :     last if $num_needed <= 0;
1558 :     }
1559 :     }
1560 :    
1561 :     foreach my $giver_id (keys %lucky_givers)
1562 :     {
1563 :     push @{$givers_to_borrowers{$giver_id}}, $borrower_id;
1564 :     }
1565 :     }
1566 :    
1567 :    
1568 :     foreach my $orig_giver_id (keys %givers_to_borrowers)
1569 :     {
1570 :     my @borrowers_list = @{$givers_to_borrowers{$orig_giver_id}};
1571 :     my @givers_list = ($orig_giver_id);
1572 :    
1573 :     # clone enough givers so that every borrower gets one. Last borrower
1574 :     # gets the orgiinal giver.
1575 :    
1576 :     while (scalar @borrowers_list > scalar @givers_list)
1577 :     {
1578 :     my $new_giver_id = $token_id_counter++;
1579 :     &clone_token($orig_giver_id, $new_giver_id, \%tokens,
1580 :     \%compounds_to_tokens, \%compounds_borrowed_to_tokens);
1581 :     push @givers_list, $new_giver_id;
1582 :     }
1583 :    
1584 :     for (my $k = 0; $k < scalar @givers_list; $k++)
1585 :     {
1586 :     my $giver_id = $givers_list[$k];
1587 :     my $giver = $tokens{$giver_id};
1588 :     my $borrower_id = $borrowers_list[$k];
1589 :     my $borrower = $tokens{$borrower_id};
1590 :    
1591 :     print STDERR "\n\tMerging '$giver_id' into '$borrower_id'\n" unless !$debug;
1592 :    
1593 :     # bump the borrower's visited reactions and compounds counters forward,
1594 :     # then merge the giver's visited reactions and compounds, unless
1595 :     # the borrower has already visited them
1596 :     map { $borrower->{visited_reactions}->{$_} += $loop_counter } keys %{$borrower->{visited_reactions}};
1597 :     map { $borrower->{visited_compounds}->{$_} += $loop_counter } keys %{$borrower->{visited_compounds}};
1598 :     map { $borrower->{visited_reactions}->{$_} = $giver->{visited_reactions}->{$_} unless defined $borrower->{visited_reactions}->{$_} } keys %{$giver->{visited_reactions}};
1599 :     map { $borrower->{visited_compounds}->{$_} = $giver->{visited_compounds}->{$_} unless defined $borrower->{visited_compounds}->{$_} } keys %{$giver->{visited_compounds}};
1600 :     map { $borrower->{token_path_inputs}->{$_} += $giver->{token_path_inputs}->{$_} } keys %{$giver->{token_path_inputs}};
1601 :     $borrower->{initial_pass} |= $giver->{initial_pass};
1602 :    
1603 :     foreach my $icpd (keys %compounds_to_tokens)
1604 :     {
1605 :     if ($compounds_to_tokens{$icpd}->{$giver_id} > 0)
1606 :     {
1607 :     $compounds_to_tokens{$icpd}->{$borrower_id} += $compounds_to_tokens{$icpd}->{$giver_id};
1608 :     delete $compounds_to_tokens{$icpd}->{$giver_id};
1609 :     }
1610 :     }
1611 :    
1612 :     foreach my $icpd (keys %compounds_borrowed_to_tokens)
1613 :     {
1614 :     if ($compounds_borrowed_to_tokens{$icpd}->{$giver_id} > 0)
1615 :     {
1616 :     $compounds_borrowed_to_tokens{$icpd}->{$borrower_id} += $compounds_borrowed_to_tokens{$icpd}->{$giver_id};
1617 :     delete $compounds_borrowed_to_tokens{$icpd}->{$giver_id};
1618 :     }
1619 :     }
1620 :    
1621 :     delete $tokens{$giver_id};
1622 :     delete $not_done_tokens{$giver_id};
1623 :    
1624 :     &balance_borrowing_and_giving($borrower_id, \%compounds_to_tokens,
1625 :     \%compounds_borrowed_to_tokens);
1626 :    
1627 :     &print_token_status([$borrower_id], \%tokens, \%compounds_to_tokens,
1628 :     \%compounds_borrowed_to_tokens, $fig);
1629 :    
1630 :     if (&check_token_for_done($borrower_id, \%compounds_to_tokens,
1631 :     \%compounds_borrowed_to_tokens,
1632 :     \%all_compounds_to_main, \%path_outputs,
1633 :     \%scenario_cycles,\%tokens, $fig, $outputs_lists))
1634 :     {
1635 :     delete $not_done_tokens{$borrower_id};
1636 :     }
1637 :     }
1638 :     }
1639 :     }
1640 :     }
1641 :    
1642 :     # Now check if we've reached a dead end, either a compound we can't push or a
1643 :     # borrowed compound we can't repay. Also determine whether there is a reaction
1644 :     # to run that can move an initial pass token forward.
1645 :    
1646 :     print STDERR "\nChecking for dead ends\n" unless !$debug;
1647 :    
1648 :     my $found_reaction_for_initial_pass_token = 0;
1649 :    
1650 :     check: foreach my $token_id (keys %not_done_tokens)
1651 :     {
1652 :     # determine which compounds the token is sitting on, and whether a
1653 :     # reaction can proceed from those compounds that isn't a loop
1654 :     # back to compounds already visited
1655 :     my %visited_compounds = %{$tokens{$token_id}->{visited_compounds}};
1656 :     my %visited_reactions = %{$tokens{$token_id}->{visited_reactions}};
1657 :     my $dead_end_cpd;
1658 :    
1659 :     print STDERR "\tChecking if '$token_id' can run\n" unless !$debug;
1660 :    
1661 :     my $found_reaction_for_token = 0;
1662 :    
1663 :     substrate: foreach my $cpd (keys %compounds_to_tokens)
1664 :     {
1665 :     if ($compounds_to_tokens{$cpd}->{$token_id} > 0 &&
1666 :     ! defined $path_outputs{$cpd})
1667 :     {
1668 :     print STDERR "\t\tChecking substrate $cpd (main: $all_compounds_to_main{$cpd})\n" unless !$debug;
1669 :    
1670 :     foreach my $reaction (@{$substrates_to_reactions->{$cpd}})
1671 :     {
1672 :     next if ($reaction =~ /(.*)_R/ && defined $visited_reactions{$1."_L"}) ||
1673 :     ($reaction =~ /(.*)_L/ && defined $visited_reactions{$1."_R"});
1674 :    
1675 :     print STDERR "\t\t\tChecking reaction $reaction\n" unless !$debug;
1676 :    
1677 :     my $prods_are_ok = 0;
1678 :    
1679 :     foreach my $prod (@{$reactions_to_products->{$reaction}})
1680 :     {
1681 :     print STDERR "\t\t\t\tChecking product $prod\n" unless !$debug;
1682 :    
1683 :     if (! defined $visited_compounds{$prod} ||
1684 :     $visited_compounds{$prod} >= $visited_compounds{$cpd} ||
1685 :     defined $path_outputs{$prod} ||
1686 :     $compounds_borrowed_to_tokens{$prod}->{$token_id} > 0)
1687 :     {
1688 :     $prods_are_ok = $prod;
1689 :     last;
1690 :     }
1691 :     }
1692 :    
1693 :     if ($prods_are_ok)
1694 :     {
1695 :     print STDERR "\tToken '$token_id' can run $reaction on $cpd to produce $prods_are_ok\n" unless !$debug;
1696 :     $found_reaction_for_token = 1;
1697 :    
1698 :     if ($tokens{$token_id}->{initial_pass})
1699 :     {
1700 :     $found_reaction_for_initial_pass_token = 1;
1701 :     }
1702 :    
1703 :     next substrate;
1704 :     }
1705 :     }
1706 :    
1707 :     # didn't find a reaction for this substrate
1708 :     $dead_end_cpd = $cpd if $all_compounds_to_main{$cpd} &&
1709 :     ! defined $path_outputs{$cpd};
1710 :     }
1711 :     }
1712 :    
1713 :     if ($dead_end_cpd)
1714 :     {
1715 :     print STDERR "\tToken '$token_id' has reached a dead end on $dead_end_cpd\n" unless !$debug;
1716 :     $tokens{$token_id}->{done} = "dead end on $dead_end_cpd";
1717 :     }
1718 :     elsif (! $found_reaction_for_token)
1719 :     {
1720 :     # didn't find any reaction to run.
1721 :     # check to see if there are borrowed compounds to repay
1722 :     product: foreach my $cpd (keys %compounds_borrowed_to_tokens)
1723 :     {
1724 :     if ($compounds_borrowed_to_tokens{$cpd}->{$token_id} > 0)
1725 :     {
1726 :     print STDERR "\t\tChecking product $cpd\n" unless !$debug;
1727 :    
1728 :     rxn: foreach my $reaction (@{$products_to_reactions->{$cpd}})
1729 :     {
1730 :     next if ($reaction =~ /(.*)_R/ && defined $visited_reactions{$1."_L"})
1731 :     || ($reaction =~ /(.*)_L/ && defined $visited_reactions{$1."_R"});
1732 :    
1733 :     print STDERR "\t\t\tChecking reaction $reaction\n" unless !$debug;
1734 :    
1735 :     my $substrates_are_ok = 0;
1736 :    
1737 :     foreach my $sub (@{$reactions_to_substrates->{$reaction}})
1738 :     {
1739 :     print STDERR "\t\t\t\tChecking substrate $sub\n" unless !$debug;
1740 :    
1741 :     if (defined $path_outputs{$sub})
1742 :     {
1743 :     # don't run reactions that use up outputs
1744 :     next rxn;
1745 :     }
1746 :    
1747 :     if (! defined $visited_compounds{$sub} ||
1748 :     $visited_compounds{$sub} <= $visited_compounds{$cpd} ||
1749 :     defined $scenario_cycles{$sub})
1750 :     {
1751 :     $substrates_are_ok = $sub;
1752 :     }
1753 :     }
1754 :    
1755 :     if ($substrates_are_ok)
1756 :     {
1757 :     print STDERR "\tToken '$token_id' can wait for $reaction on $substrates_are_ok to produce borrowed compound $cpd\n" unless !$debug;
1758 :     $found_reaction_for_token = 1;
1759 :    
1760 :     if ($tokens{$token_id}->{initial_pass})
1761 :     {
1762 :     $found_reaction_for_initial_pass_token = 1;
1763 :     }
1764 :    
1765 :     last product;
1766 :     }
1767 :     }
1768 :    
1769 :     # didn't find a reaction for this product
1770 :     $dead_end_cpd = $cpd if $all_compounds_to_main{$cpd} &&
1771 :     ! defined $path_outputs{$cpd};
1772 :     }
1773 :     }
1774 :    
1775 :     # didn't find any reaction to run.
1776 :     if (! $found_reaction_for_token)
1777 :     {
1778 :     if ($dead_end_cpd)
1779 :     {
1780 :     print STDERR "\tToken '$token_id' has reached a dead end on borrowed compound $dead_end_cpd\n" unless !$debug;
1781 :     $tokens{$token_id}->{done} = "dead end on borrowed compound $dead_end_cpd";
1782 :     }
1783 :     else
1784 :     {
1785 :     # nothing to push, borrow or do
1786 :     print STDERR "\tToken '$token_id' has reached a dead end\n" unless !$debug;
1787 :     $tokens{$token_id}->{done} = "dead end";
1788 :     }
1789 :     }
1790 :     }
1791 :     }
1792 :    
1793 :     if($find_first)
1794 :     {
1795 :     foreach my $token_id (sort { $tokens{$a}->{done} <=> $tokens{$b}->{done} } keys %tokens)
1796 :     {
1797 :     if($tokens{$token_id}->{done} == 1)
1798 :     {
1799 :     $done = 1;
1800 :     }
1801 :     }
1802 :     }
1803 :     # is there an initial pass token that can make progress?
1804 :     if ($found_reaction_for_initial_pass_token)
1805 :     {
1806 :     if (scalar keys %reactions_to_tokens_to_use == 0)
1807 :     {
1808 :     # Couldn't run any reactions this time around.
1809 :     # Push more tokens through from the beginning of the path to
1810 :     # supply more substrates.
1811 :     $add_path_inputs = 1;
1812 :     print STDERR "\nSupplying more path inputs to push stalled tokens\n" unless !$debug;
1813 :     }
1814 :     }
1815 :     else
1816 :     {
1817 :     $done = 1;
1818 :     }
1819 :    
1820 :     $loop_counter++;
1821 :    
1822 :     if ($loop_counter >= $infinite_loop_check)
1823 :     {
1824 :     $data_results{"infinite"} = 1;
1825 :     print STDERR "Encountered an infinite loop\n" unless !$debug;
1826 :     $done = 1;
1827 :     }
1828 :     }
1829 :    
1830 :     # reverse %compounds_to_tokens, since all tokens should be at path outputs now
1831 :     my %tokens_to_compounds;
1832 :    
1833 :     foreach my $cpd (keys %compounds_to_tokens)
1834 :     {
1835 :     foreach my $token_id (keys %{$compounds_to_tokens{$cpd}})
1836 :     {
1837 :     my $num_tokens = $compounds_to_tokens{$cpd}->{$token_id};
1838 :     $tokens_to_compounds{$token_id}->{$cpd} = $num_tokens if $num_tokens > 0;
1839 :     }
1840 :     }
1841 :    
1842 :     print STDERR "\n\ntoken ids: @{[ sort { $a <=> $b } map { $_ if ! defined $tokens{$_}->{done} } keys %tokens ]}\n" unless !$debug;
1843 :    
1844 :     my $path_counter = 1;
1845 :    
1846 :     foreach my $token_id (sort { $tokens{$a}->{done} <=> $tokens{$b}->{done} } keys %tokens)
1847 :     {
1848 :     my $token = $tokens{$token_id};
1849 :     my %visited_reactions = %{$token->{visited_reactions}};
1850 :     my @path = sort { $visited_reactions{$a} <=> $visited_reactions{$b} }
1851 :     keys %visited_reactions;
1852 :     my %visited_compounds = %{$token->{visited_compounds}};
1853 :     my @compounds = sort { $visited_compounds{$a} <=> $visited_compounds{$b} }
1854 :     keys %visited_compounds;
1855 :    
1856 :     print STDERR "Adding token id: $token_id\n" unless !$debug;
1857 :    
1858 :    
1859 :     #each key in data_results is a token which points to an array
1860 :     # [0]=initial pass [1]=0/1 if its done [2]=reaction path [3]=compounds
1861 :     # [4]= html string of inputs
1862 :     # [5]= html string of outputs
1863 :     # [6]= html string of borrowed compounds
1864 :     # [7] = array of path input compounds
1865 :     # [8] = array of path output compounds
1866 :    
1867 :    
1868 :     $data_results{$token_id} = [$token->{initial_pass},$token->{done},\@path,\@compounds,[],[],[],$token->{token_path_inputs},$tokens_to_compounds{$token_id}];
1869 :    
1870 :     foreach my $input (keys %{$token->{token_path_inputs}})
1871 :     {
1872 :     my $input_stoich = $token->{token_path_inputs}->{$input};
1873 :     my $output_stoich = $tokens_to_compounds{$token_id}->{$input};
1874 :    
1875 :     # don't balance scenario cycled compounds until final assembly
1876 :     if ($scenario_cycles{$input} && ! $create_assembly)
1877 :     {
1878 :     my @names = $fig->names_of_compound($input);
1879 :     push @{$data_results{$token_id}->[4]}, "\t\t$input_stoich\t$input $names[0]\n";
1880 :     next;
1881 :     }
1882 :    
1883 :     if ($input_stoich > $output_stoich)
1884 :     {
1885 :     delete $tokens_to_compounds{$token_id}->{$input};
1886 :     $input_stoich -= $output_stoich;
1887 :     $token->{token_path_inputs}->{$input} -= $output_stoich;
1888 :     my @names = $fig->names_of_compound($input);
1889 :     push @{$data_results{$token_id}->[4]}, "\t\t$input_stoich\t$input $names[0]\n";
1890 :     }
1891 :     elsif ($output_stoich > $input_stoich)
1892 :     {
1893 :     delete $token->{token_path_inputs}->{$input};
1894 :     $tokens_to_compounds{$token_id}->{$input} -= $input_stoich;
1895 :     }
1896 :     else
1897 :     {
1898 :     delete $token->{token_path_inputs}->{$input};
1899 :     delete $tokens_to_compounds{$token_id}->{$input};
1900 :     }
1901 :     }
1902 :    
1903 :     foreach my $output (keys %{$tokens_to_compounds{$token_id}})
1904 :     {
1905 :     my @names = $fig->names_of_compound($output);
1906 :     push @{$data_results{$token_id}->[5]},"\t\t$tokens_to_compounds{$token_id}->{$output}\t$output $names[0]\n";
1907 :     }
1908 :    
1909 :     if ($token->{done} != 1)
1910 :     {
1911 :     foreach my $cpd (sort keys %compounds_borrowed_to_tokens)
1912 :     {
1913 :     my $num = $compounds_borrowed_to_tokens{$cpd}->{$token_id};
1914 :     my @names = $fig->names_of_compound($cpd);
1915 :     push @{$data_results{$token_id}->[6]},"\t\t$num $cpd\t$names[0]\n" if ($num > 0);
1916 :     }
1917 :     }
1918 :     }
1919 :    
1920 :    
1921 :    
1922 :    
1923 :     return \%data_results;
1924 :     }
1925 :    
1926 :     sub balance_borrowing_and_giving
1927 :     {
1928 :     my ($token_id, $compounds_to_tokens, $compounds_borrowed_to_tokens) = @_;
1929 :    
1930 :     my %merged_compounds;
1931 :     map { $merged_compounds{$_} = 1 } keys %{$compounds_to_tokens};
1932 :     map { $merged_compounds{$_} = 1 } keys %{$compounds_borrowed_to_tokens};
1933 :    
1934 :     foreach my $icpd (keys %merged_compounds)
1935 :     {
1936 :     my $inum_to_give = $compounds_to_tokens->{$icpd}->{$token_id};
1937 :    
1938 :     if ($inum_to_give > 0)
1939 :     {
1940 :     my $inum_needed = $compounds_borrowed_to_tokens->{$icpd}->{$token_id};
1941 :    
1942 :     if ($inum_to_give == $inum_needed)
1943 :     {
1944 :     delete $compounds_borrowed_to_tokens->{$icpd}->{$token_id};
1945 :     delete $compounds_to_tokens->{$icpd}->{$token_id};
1946 :     }
1947 :     elsif ($inum_to_give > $inum_needed)
1948 :     {
1949 :     delete $compounds_borrowed_to_tokens->{$icpd}->{$token_id};
1950 :     $compounds_to_tokens->{$icpd}->{$token_id} -= $inum_needed;
1951 :     }
1952 :     else
1953 :     {
1954 :     $compounds_borrowed_to_tokens->{$icpd}->{$token_id} -= $inum_to_give;
1955 :     delete $compounds_to_tokens->{$icpd}->{$token_id};
1956 :     }
1957 :     }
1958 :     }
1959 :     }
1960 :    
1961 :     sub print_token_status
1962 :     {
1963 :     my ($token_id_list, $tokens, $compounds_to_tokens, $compounds_borrowed_to_tokens, $fig) = @_;
1964 :    
1965 :     print STDERR "\nToken status:\n" unless !$debug;
1966 :    
1967 :     foreach my $token_id (@$token_id_list)
1968 :     {
1969 :     next if defined $tokens->{$token_id}->{done};
1970 :    
1971 :     print STDERR "\n\ttoken: '$token_id', initial: $tokens->{$token_id}->{initial_pass}\n" unless !$debug;
1972 :    
1973 :     foreach my $cpd (sort keys %{$tokens->{$token_id}->{token_path_inputs}})
1974 :     {
1975 :     my $num = $tokens->{$token_id}->{token_path_inputs}->{$cpd};
1976 :     my @names = $fig->names_of_compound($cpd);
1977 :     print STDERR "\t\tInput: $num $cpd\t$names[0]\n" unless !$debug;
1978 :     }
1979 :    
1980 :     foreach my $cpd (sort keys %$compounds_to_tokens)
1981 :     {
1982 :     my $num = $compounds_to_tokens->{$cpd}->{$token_id};
1983 :     my @names = $fig->names_of_compound($cpd);
1984 :     print STDERR "\t\tStatus: $num $cpd\t$names[0]\n" if ($num > 0 && $debug);
1985 :     }
1986 :    
1987 :     foreach my $cpd (sort keys %$compounds_borrowed_to_tokens)
1988 :     {
1989 :     my $num = $compounds_borrowed_to_tokens->{$cpd}->{$token_id};
1990 :     my @names = $fig->names_of_compound($cpd);
1991 :     print STDERR "\t\tBorrowed: $num $cpd\t$names[0]\n" if ($num > 0 && $debug);
1992 :     }
1993 :    
1994 :     my %visited_compounds = %{$tokens->{$token_id}->{visited_compounds}};
1995 :     print STDERR "\t\tvisited_compounds: @{[map { ($_, $visited_compounds{$_} ) } sort { $visited_compounds{$a} <=> $visited_compounds{$b} } keys %visited_compounds]}\n" unless !$debug;
1996 :    
1997 :     my %visited_reactions = %{$tokens->{$token_id}->{visited_reactions}};
1998 :     print STDERR "\t\tvisited_reactions: @{[map { ($_, $visited_reactions{$_} ) } sort { $visited_reactions{$a} <=> $visited_reactions{$b} } keys %visited_reactions]}\n" unless !$debug;
1999 :    
2000 :     }
2001 :    
2002 :     print STDERR "\n" unless !$debug;
2003 :     }
2004 :    
2005 :     sub clone_token
2006 :     {
2007 :     my ($clone_id, $new_token_id, $tokens, $compounds_to_tokens, $compounds_borrowed_to_tokens) = @_;
2008 :     my (%new_token, %new_visited_reactions, %new_visited_compounds, %new_token_path_inputs);
2009 :    
2010 :     $tokens->{$new_token_id} = \%new_token;
2011 :    
2012 :     my $clone_token = $tokens->{$clone_id};
2013 :     map {$new_visited_reactions{$_} = $clone_token->{visited_reactions}->{$_}} keys %{$clone_token->{visited_reactions}};
2014 :     map {$new_visited_compounds{$_} = $clone_token->{visited_compounds}->{$_}} keys %{$clone_token->{visited_compounds}};
2015 :     map {$new_token_path_inputs{$_} = $clone_token->{token_path_inputs}->{$_}} keys %{$clone_token->{token_path_inputs}};
2016 :    
2017 :     $new_token{visited_reactions} = \%new_visited_reactions;
2018 :     $new_token{visited_compounds} = \%new_visited_compounds;
2019 :     $new_token{token_path_inputs} = \%new_token_path_inputs;
2020 :     $new_token{initial_pass} = $clone_token->{initial_pass};
2021 :    
2022 :     # tokens might be spread across multiple compounds
2023 :     foreach my $icpd (keys %$compounds_to_tokens)
2024 :     {
2025 :     if ($compounds_to_tokens->{$icpd}->{$clone_id} > 0)
2026 :     {
2027 :     $compounds_to_tokens->{$icpd}->{$new_token_id} =
2028 :     $compounds_to_tokens->{$icpd}->{$clone_id};
2029 :     }
2030 :     }
2031 :    
2032 :     # tokens might have borrowed compounds
2033 :     foreach my $icpd (keys %$compounds_borrowed_to_tokens)
2034 :     {
2035 :     if ($compounds_borrowed_to_tokens->{$icpd}->{$clone_id} > 0)
2036 :     {
2037 :     $compounds_borrowed_to_tokens->{$icpd}->{$new_token_id} =
2038 :     $compounds_borrowed_to_tokens->{$icpd}->{$clone_id};
2039 :     }
2040 :     }
2041 :    
2042 :     return \%new_token;
2043 :     }
2044 :    
2045 :     sub check_token_for_done
2046 :     {
2047 :     my ($token_id, $compounds_to_tokens, $compounds_borrowed_to_tokens, $all_compounds_to_main,
2048 :     $path_outputs, $scenario_cycles, $tokens, $fig, $outputs_lists) = @_;
2049 :    
2050 :     my $token_is_done_pushing = 1;
2051 :     my $token_is_done_borrowing = 1;
2052 :    
2053 :     # first determine if there is a main compound that isn't a path output
2054 :     foreach my $cpd (keys %$compounds_to_tokens)
2055 :     {
2056 :     if ($compounds_to_tokens->{$cpd}->{$token_id} > 0)
2057 :     {
2058 :     # also check if scenario cycle compounds need to be pushed
2059 :     if ($all_compounds_to_main->{$cpd} &&
2060 :     (! defined $path_outputs->{$cpd} ||
2061 :     ($scenario_cycles->{$cpd} &&
2062 :     $tokens->{$token_id}->{visited_compounds}->{$cpd} == 0)))
2063 :     {
2064 :     my @names = $fig->names_of_compound($cpd);
2065 :     print STDERR "\ttoken '$token_id' needs to push $cpd $names[0]\n" unless !$debug;
2066 :     $token_is_done_pushing = 0;
2067 :     }
2068 :     }
2069 :     }
2070 :    
2071 :     # now determine if one of the output lists has been satisfied
2072 :     if ($token_is_done_pushing)
2073 :     {
2074 :     my $found_a_list = 0;
2075 :    
2076 :     foreach my $cpd_list (@$outputs_lists)
2077 :     {
2078 :     print STDERR "\t\tchecking outputs_list @$cpd_list\n" unless !$debug;
2079 :    
2080 :     $found_a_list = 1;
2081 :    
2082 :     foreach my $cpd (@$cpd_list)
2083 :     {
2084 :     if (! $compounds_to_tokens->{$cpd}->{$token_id} > 0)
2085 :     {
2086 :     $found_a_list = 0;
2087 :     last;
2088 :     }
2089 :     }
2090 :    
2091 :     last if $found_a_list;
2092 :     }
2093 :    
2094 :     if (! $found_a_list)
2095 :     {
2096 :     print STDERR "\ttoken '$token_id' hasn't satisfied output compound list\n" unless !$debug;
2097 :     $token_is_done_pushing = 0;
2098 :     }
2099 :     }
2100 :    
2101 :     if ($token_is_done_pushing)
2102 :     {
2103 :     foreach my $cpd (keys %$compounds_borrowed_to_tokens)
2104 :     {
2105 :     if ($compounds_borrowed_to_tokens->{$cpd}->{$token_id} > 0)
2106 :     {
2107 : dejongh 1.17 # I don't know why I put in this "if" statement, so I'm going to do the same thing
2108 :     # in both cases, but print a distinguishing debug statement in case I need to find
2109 :     # the occurrences later
2110 : olson 1.1 if ($all_compounds_to_main{$cpd})
2111 :     {
2112 :     my @names = $fig->names_of_compound($cpd);
2113 : dejongh 1.17 print STDERR "\ttoken '$token_id' has borrowed $cpd $names[0] (main)\n" unless !$debug;
2114 :     $token_is_done_borrowing = 0;
2115 :     }
2116 :     else
2117 :     {
2118 :     my @names = $fig->names_of_compound($cpd);
2119 :     print STDERR "\ttoken '$token_id' has borrowed $cpd $names[0] (not main)\n" unless !$debug;
2120 : olson 1.1 $token_is_done_borrowing = 0;
2121 :     }
2122 :     }
2123 :     }
2124 :     }
2125 :    
2126 :     if (! $token_is_done_pushing || ! $token_is_done_borrowing)
2127 :     {
2128 :     return 0;
2129 :     }
2130 :     else
2131 :     {
2132 :     $tokens->{$token_id}->{done} = 1;
2133 :     print STDERR "\tToken '$token_id' is done\n" unless !$debug;
2134 :     return 1;
2135 :     }
2136 :     }
2137 :    
2138 :    
2139 :     sub write_fluxanalyzer_files
2140 :     {
2141 :     my ($dir, $path_inputs, $path_outputs, $path_array,
2142 :     $all_reactions,$reactions_to_substrate_arrays,$reactions_to_product_arrays,
2143 :     $cidToName) = @_;
2144 :    
2145 :     my $x_pos = 10;
2146 :     my $y_pos = 20;
2147 :    
2148 :     #Write the inputs/outputs to a seperate file, along with $stoich and if its main
2149 :     open(A_INPUT, ">$dir/inputs_main");
2150 :    
2151 :     print A_INPUT map {"$_\t$all_compounds_to_main{$_}\n"} keys %$path_inputs;
2152 :    
2153 :     close(A_INPUT);
2154 :    
2155 :     open(A_OUTPUT, ">$dir/outputs_main");
2156 :    
2157 :     print A_OUTPUT map {"$_\t$all_compounds_to_main{$_}\n"} keys %$path_outputs;
2158 :    
2159 :     close(A_OUTPUT);
2160 :    
2161 :    
2162 :    
2163 :     open(REACTIONS, ">$dir/reactions");
2164 :     open(INPUTS, ">$dir/inputs");
2165 :     open(OUTPUTS, ">$dir/outputs");
2166 :     open(PATH, ">$dir/path_info");
2167 :    
2168 :     my @inputs = keys %$path_inputs;
2169 :    
2170 :     foreach my $elem(@{$path_array}){
2171 :     print PATH $elem ."\n";
2172 :     }
2173 :     close(PATH);
2174 :    
2175 :     foreach my $cpd (@inputs)
2176 :     {
2177 :     my $toPrint = $cpd."up\t = 1 $cpd \t| \t";
2178 :     $toPrint .= $path_inputs->{$cpd};
2179 :     $toPrint .= " \t0 100 0 \t$x_pos $y_pos 1 1\t0.01\n";
2180 :     print REACTIONS $toPrint;
2181 :     $y_pos += 20;
2182 :    
2183 :     if ($y_pos == 300)
2184 :     {
2185 :     $x_pos += 60;
2186 :     $y_pos = 20;
2187 :     }
2188 :     }
2189 :    
2190 :     $x_pos += 60;
2191 :     $y_pos = 20;
2192 :    
2193 :     my (@display_array);
2194 :     foreach my $rxn (keys %$all_reactions)
2195 :     {
2196 :     my $direction = $all_reactions->{$rxn};
2197 :    
2198 :     my (@substrate_array, @product_array);
2199 :    
2200 :     push @display_array, $rxn;
2201 :    
2202 :     if ($direction eq "L")
2203 :     {
2204 :     @product_array = @{$reactions_to_substrate_arrays->{$rxn."_L"}};
2205 :     @substrate_array = @{$reactions_to_product_arrays->{$rxn."_L"}};
2206 :     }
2207 :     else
2208 :     {
2209 :     @substrate_array = @{$reactions_to_substrate_arrays->{$rxn."_R"}};
2210 :     @product_array = @{$reactions_to_product_arrays->{$rxn."_R"}};
2211 :     }
2212 :    
2213 :     foreach my $subTuple (@substrate_array)
2214 :     {
2215 :     my @temp = $fig->names_of_compound($subTuple->[0]);
2216 :     $cidToName->{$subTuple->[0]} = $temp [0] if ! defined $cidToName->{$subTuple->[0]};
2217 :     }
2218 :    
2219 :     foreach my $prodTuple(@product_array)
2220 :     {
2221 :     my @temp = $fig->names_of_compound($prodTuple->[0]);
2222 :     $cidToName->{$prodTuple->[0]} = $temp[0] if ! defined $cidToName->{$prodTuple->[0]};
2223 :     }
2224 :    
2225 :     if($direction eq "R" || $direction eq "B")
2226 :     {
2227 :     #write data in a strign for copying to file later
2228 :     my $toFile = '';
2229 :     $toFile.= $rxn."\t";
2230 :    
2231 :     #add all the substrates
2232 :     foreach my $curSub(@substrate_array){
2233 :     $toFile .= $curSub -> [1].' '. $curSub -> [0].' + ';
2234 :     }
2235 :    
2236 :     ##chop off the +
2237 :     chop($toFile);
2238 :     chop($toFile);
2239 :    
2240 :     $toFile.='= ';
2241 :    
2242 :     #add all the products
2243 :     foreach my $curProd(@product_array){
2244 :     $toFile .= $curProd -> [1].' '. $curProd -> [0].' + ';
2245 :     }
2246 :    
2247 :     #chop off the plus
2248 :     chop($toFile);
2249 :     chop($toFile);
2250 :    
2251 :     $toFile.="\t|\t#\t";
2252 :    
2253 :     if($direction eq "B"){
2254 :     $toFile.="-Inf";
2255 :     }
2256 :     else{
2257 :     $toFile.="0";
2258 :     }
2259 :    
2260 :     $toFile.=" Inf 0\t$x_pos $y_pos 1 1 \t0.01\n";
2261 :    
2262 :     print REACTIONS $toFile;
2263 :     }
2264 :     elsif($direction eq "L")
2265 :     {
2266 :     #write data in a strign for copying to file later
2267 :     my $toFile = '';
2268 :     $toFile.= $rxn."\t";
2269 :    
2270 :     #add all the substrates
2271 :     foreach my $curProd(@product_array){
2272 :     $toFile .= $curProd -> [1].' '. $curProd -> [0]." + ";
2273 :     }
2274 :    
2275 :     ##chop off the +
2276 :     chop($toFile);
2277 :     chop($toFile);
2278 :    
2279 :     $toFile.="= ";
2280 :    
2281 :     #add all the products
2282 :     foreach my $curSubstrate(@substrate_array){
2283 :     $toFile .= $curSubstrate -> [1].' '. $curSubstrate -> [0]." + ";
2284 :     }
2285 :    
2286 :     #chop off the plus
2287 :     chop($toFile);
2288 :     chop($toFile);
2289 :    
2290 :     $toFile.="\t|\t#\t";
2291 :    
2292 :     $toFile.="0";
2293 :    
2294 :     $toFile.=" Inf 0\t$x_pos $y_pos 1 1 \t0.01\n";
2295 :    
2296 :     print REACTIONS $toFile;
2297 :     }
2298 :    
2299 :     $y_pos += 20;
2300 :    
2301 :     if ($y_pos == 300)
2302 :     {
2303 :     $x_pos += 60;
2304 :     $y_pos = 20;
2305 :     }
2306 :     }
2307 :    
2308 :     my @outputs = keys %$path_outputs;
2309 :    
2310 :     $x_pos += 60;
2311 :     $y_pos = 20;
2312 :    
2313 :     foreach my $cpd (@inputs)
2314 :     {
2315 :     print INPUTS $cpd, "\t", $path_inputs->{$cpd}, "\t", $cidToName->{$cpd}, "\n";
2316 :     }
2317 :    
2318 :     foreach my $cpd (@outputs)
2319 :     {
2320 :     print OUTPUTS $cpd, "\t", $path_outputs->{$cpd}, "\t", $cidToName->{$cpd}, "\n";
2321 :     my $toPrint = $cpd."ex\t 1 $cpd = \t| \t# \t0 100 0 \t$x_pos $y_pos 1 1\t0.01\n";
2322 :     print REACTIONS $toPrint;
2323 :     $y_pos += 20;
2324 :    
2325 :     if ($y_pos == 300)
2326 :     {
2327 :     $x_pos += 60;
2328 :     $y_pos = 20;
2329 :     }
2330 :     }
2331 :    
2332 :     $x_pos += 60;
2333 :     $y_pos = 20;
2334 :    
2335 :     #print the macromolecule_synthesis and assembly file
2336 :     open(MACRO_SYTH,">$dir/macromolecule_synthesis");
2337 :     open(ASSEM,">$dir/assembly");
2338 :     my $toPrint = "M1 = ";
2339 :    
2340 :     foreach my $cpd (keys %$path_outputs)
2341 :     {
2342 :     $toPrint.="$path_outputs->{$cpd} $cpd + ";
2343 :     print ASSEM "$cpd\tM1\t-100 -100 1\n";
2344 :     $y_pos += 25;
2345 :     }
2346 :    
2347 :     chop $toPrint;
2348 :     chop $toPrint;
2349 :     chop $toPrint;
2350 :     print MACRO_SYTH $toPrint;
2351 :     close(MACRO_SYTH);
2352 :     close(ASSEM);
2353 :    
2354 :     #Print the metabolites for these subsystems.
2355 :     open(METABOLITES,">$dir/metabolites");
2356 :    
2357 :     foreach my $cid (keys %$cidToName)
2358 :     {
2359 :     my $name = $cidToName->{$cid};
2360 :     $name =~ s/\s/-/g;
2361 :     print METABOLITES $cid."\t".$name."\t0.001\t0\n";
2362 :     }
2363 :     close(METABOLITES);
2364 :    
2365 :     #Print the macromolucules file
2366 :     open(MACRO,">$dir/macromolecules");
2367 :     print MACRO "M1 \tM1 \t1 \t-100 -100 1 1\n";
2368 :     close(MACRO);
2369 :    
2370 :     $x_pos += 60;
2371 :     $y_pos = 20;
2372 :    
2373 :     print REACTIONS "mue\t\t\t|\t#\t0 100 0\t$x_pos $y_pos 1 1\t0.01\n";
2374 :    
2375 :     #close reaction equation file
2376 :     close(REACTIONS);
2377 :    
2378 :     # FluxAnalyzer requires this file
2379 :     open(APP, ">$dir/app_para.m");
2380 :     print APP "epsilon=1e-10;\nbasic_color=[0.7 0.7 0.7];\ncr_color=[0.5 0.5 1];\nbr_color=[1 0.2 0.2];\nnbr_color=[0.2 1 0.2];\ntext_color=[0 0 0];\nmacro_synth_color=[0 0 1];\nmacro_color=[0.6 0.6 1];\nbox_reaction_width=[0.12];\nbox_reaction_height=[0.06];\nbox_macro_width=[0.08];\nbox_macro_height=[0.06];\nfontsize_reaction=[11];\nfontsize_macro=[11];\nfluxmaps={'Fluxmap','dummy.pcx'};\n";
2381 :     close(APP);
2382 :     }
2383 :    
2384 :     sub write_final_fluxanalyzer_files
2385 :     {
2386 :     my ($dir, $path_inputs, $path_outputs, $all_reactions, $transport_reactions,
2387 :     $reactions_to_substrate_arrays,$reactions_to_product_arrays,
2388 :     $cidToName,$bioMass,$minSubstrates) = @_;
2389 :    
2390 :     #Write the inputs/outputs to a seperate file, along with $stoich and if its main
2391 :     open(A_INPUT, ">$dir/inputs_main");
2392 :    
2393 :     print A_INPUT map {"$_\t$all_compounds_to_main{$_}\n"} keys %$path_inputs;
2394 :    
2395 :     close(A_INPUT);
2396 :    
2397 :     open(A_OUTPUT, ">$dir/outputs_main");
2398 :    
2399 :     print A_OUTPUT map {"$_\t$all_compounds_to_main{$_}\n"} keys %$path_outputs;
2400 :    
2401 :     close(A_OUTPUT);
2402 :    
2403 :     my %open_transports;
2404 :    
2405 :     foreach my $cpd (keys %$minSubstrates)
2406 :     {
2407 :     map { $open_transports{$_} = 1 } @{$minSubstrates->{$cpd}};
2408 :     }
2409 :    
2410 :     open(REACTIONS, ">$dir/reactions");
2411 :     open(INPUTS, ">$dir/inputs");
2412 :     open(OUTPUTS, ">$dir/outputs");
2413 :    
2414 :     my @inputs = keys %$path_inputs;
2415 :    
2416 :     my $x_pos = 10;
2417 :     my $y_pos = 30;
2418 :    
2419 :     foreach my $cpd (@inputs)
2420 :     {
2421 :     my $toPrint = $cpd."up\t = 1 $cpd \t| \t";
2422 :    
2423 :     if (defined $minSubstrates->{$cpd})
2424 :     {
2425 :     $toPrint .= "#";
2426 :     }
2427 :     else
2428 :     {
2429 :     $toPrint .= "0";
2430 :     }
2431 :    
2432 :     $toPrint .= " \t0 Inf 0 \t$x_pos $y_pos 1 1\t0.01\n";
2433 :     print REACTIONS $toPrint;
2434 :     $y_pos += 30;
2435 :    
2436 :     if ($y_pos > 600)
2437 :     {
2438 :     $x_pos += 60;
2439 :     $y_pos = 30;
2440 :     }
2441 :     }
2442 :    
2443 :     $x_pos += 60;
2444 :     $y_pos = 30;
2445 :    
2446 :     foreach my $rxn (keys %$all_reactions)
2447 :     {
2448 :     my $direction = $all_reactions->{$rxn};
2449 :    
2450 :     my (@substrate_array, @product_array);
2451 :    
2452 :     if ($direction eq "L")
2453 :     {
2454 :     @product_array = @{$reactions_to_substrate_arrays->{$rxn."_L"}};
2455 :     @substrate_array = @{$reactions_to_product_arrays->{$rxn."_L"}};
2456 :     }
2457 :     else
2458 :     {
2459 :     @substrate_array = @{$reactions_to_substrate_arrays->{$rxn."_R"}};
2460 :     @product_array = @{$reactions_to_product_arrays->{$rxn."_R"}};
2461 :     }
2462 :    
2463 :     foreach my $subTuple (@substrate_array)
2464 :     {
2465 :     my @temp = $fig->names_of_compound($subTuple->[0]);
2466 :     $cidToName->{$subTuple->[0]} = $temp [0] if ! defined $cidToName->{$subTuple->[0]};
2467 :     }
2468 :    
2469 :     foreach my $prodTuple(@product_array)
2470 :     {
2471 :     my @temp = $fig->names_of_compound($prodTuple->[0]);
2472 :     $cidToName->{$prodTuple->[0]} = $temp[0] if ! defined $cidToName->{$prodTuple->[0]};
2473 :     }
2474 :    
2475 :     if($direction eq "R" || $direction eq "B")
2476 :     {
2477 :     #write data in a strign for copying to file later
2478 :     my $toFile = '';
2479 :     $toFile.= $rxn."\t";
2480 :    
2481 :     #add all the substrates
2482 :     foreach my $curSub(@substrate_array){
2483 :     $toFile .= $curSub -> [1].' '. $curSub -> [0].' + ';
2484 :     }
2485 :    
2486 :     ##chop off the +
2487 :     chop($toFile);
2488 :     chop($toFile);
2489 :    
2490 :     $toFile.='= ';
2491 :    
2492 :     #add all the products
2493 :     my $found_prod = 0;
2494 :     foreach my $curProd(@product_array){
2495 :     $toFile .= $curProd -> [1].' '. $curProd -> [0].' + ';
2496 :     $found_prod = 1;
2497 :     }
2498 :    
2499 :     if ($found_prod)
2500 :     {
2501 :     #chop off the plus
2502 :     chop($toFile);
2503 :     chop($toFile);
2504 :     }
2505 :    
2506 :     if (defined $transport_reactions->{$rxn})
2507 :     {
2508 :     if (defined $open_transports{$rxn})
2509 :     {
2510 :     $toFile.="\t|\t#\t";
2511 :     }
2512 :     else
2513 :     {
2514 :     $toFile.="\t|\t0\t";
2515 :     }
2516 :    
2517 :     if($direction eq "B"){
2518 :     $toFile.="-Inf";
2519 :     }
2520 :     else{
2521 :     $toFile.="0";
2522 :     }
2523 :    
2524 :     $toFile.=" Inf 0\t$x_pos $y_pos 1 1 \t0.01\n";
2525 :     $y_pos += 30;
2526 :     }
2527 :     elsif ($rxn =~ /sink/)
2528 :     {
2529 :     $toFile.="\t|\t0\t";
2530 :    
2531 :     if($direction eq "B"){
2532 :     $toFile.="-0.00001";
2533 :     }
2534 :     else{
2535 :     $toFile.="0";
2536 :     }
2537 :    
2538 :     $toFile.=" 0.00001 0\t$x_pos $y_pos 1 1 \t0.01\n";
2539 :     $y_pos += 30;
2540 :     }
2541 :     else
2542 :     {
2543 :     $toFile.="\t|\t#\t";
2544 :    
2545 :     if($direction eq "B"){
2546 :     $toFile.="-Inf";
2547 :     }
2548 :     else{
2549 :     $toFile.="0";
2550 :     }
2551 :    
2552 :     $toFile.=" Inf 0\t-10 -10 1 1 \t0.01\n";
2553 :     }
2554 :    
2555 :     print REACTIONS $toFile;
2556 :     }
2557 :     elsif($direction eq "L")
2558 :     {
2559 :     #write data in a strign for copying to file later
2560 :     my $toFile = '';
2561 :     $toFile.= $rxn."\t";
2562 :    
2563 :     #add all the substrates
2564 :     foreach my $curProd(@product_array){
2565 :     $toFile .= $curProd -> [1].' '. $curProd -> [0]." + ";
2566 :     }
2567 :    
2568 :     ##chop off the +
2569 :     chop($toFile);
2570 :     chop($toFile);
2571 :    
2572 :     $toFile.="= ";
2573 :    
2574 :     #add all the products
2575 :     foreach my $curSubstrate(@substrate_array){
2576 :     $toFile .= $curSubstrate -> [1].' '. $curSubstrate -> [0]." + ";
2577 :     }
2578 :    
2579 :     #chop off the plus
2580 :     chop($toFile);
2581 :     chop($toFile);
2582 :    
2583 :    
2584 :     $toFile.="\t|\t#\t";
2585 :     $toFile.="0";
2586 :     $toFile.=" Inf 0\t-10 -10 1 1 \t0.01\n";
2587 :    
2588 :     print REACTIONS $toFile;
2589 :     }
2590 :    
2591 :     if ($y_pos > 600)
2592 :     {
2593 :     $x_pos += 60;
2594 :     $y_pos = 30;
2595 :     }
2596 :     }
2597 :    
2598 :     my @outputs = keys %$path_outputs;
2599 :    
2600 :     $x_pos += 60;
2601 :     $y_pos = 30;
2602 :    
2603 :     foreach my $cpd (@inputs)
2604 :     {
2605 :     print INPUTS $cpd, "\t", $path_inputs->{$cpd}, "\t", $cidToName->{$cpd}, "\n";
2606 :     }
2607 :    
2608 :     foreach my $cpd (@outputs)
2609 :     {
2610 :     print OUTPUTS $cpd, "\t", $path_outputs->{$cpd}, "\t", $cidToName->{$cpd}, "\n";
2611 :     my $toPrint = $cpd."ex\t 1 $cpd = \t| \t# \t0 Inf 0 \t$x_pos $y_pos 1 1\t0.01\n";
2612 :     print REACTIONS $toPrint;
2613 :     $y_pos += 30;
2614 :    
2615 :     if ($y_pos > 600)
2616 :     {
2617 :     $x_pos += 60;
2618 :     $y_pos = 30;
2619 :     }
2620 :     }
2621 :    
2622 :     $x_pos += 60;
2623 :     $y_pos = 30;
2624 :    
2625 :     #print the macromolecule_synthesis and assembly file
2626 :     open(MACRO,">$dir/macromolecules");
2627 :     open(MACRO_SYTH,">$dir/macromolecule_synthesis");
2628 :     open(ASSEM,">$dir/assembly");
2629 :    
2630 :     my $toPrint = "M1 = ";
2631 :    
2632 :     foreach my $cpd (keys %$bioMass)
2633 :     {
2634 :     $toPrint .= "$bioMass->{$cpd} $cpd + ";
2635 :     print ASSEM "$cpd\tM1\t-100 -100 1\n";
2636 :     }
2637 :    
2638 :     chop $toPrint;
2639 :     chop $toPrint;
2640 :     chop $toPrint;
2641 :    
2642 :     print MACRO "M1 \tM1 \t1 \t-100 -100 1 1\n";
2643 :     print MACRO_SYTH $toPrint, "\n";
2644 :    
2645 :     close(MACRO);
2646 :     close(MACRO_SYTH);
2647 :     close(ASSEM);
2648 :    
2649 :     #Print the metabolites for these subsystems.
2650 :     open(METABOLITES,">$dir/metabolites");
2651 :    
2652 :     foreach my $cid (keys %$cidToName)
2653 :     {
2654 :     my $name = $cidToName->{$cid};
2655 :     $name =~ s/\s/-/g;
2656 :     print METABOLITES $cid."\t".$name."\t0.001\t0\n";
2657 :     }
2658 :     close(METABOLITES);
2659 :    
2660 :     $x_pos += 60;
2661 :     $y_pos = 30;
2662 :    
2663 :     print REACTIONS "mue\t\t\t|\t#\t0 100 0\t$x_pos $y_pos 1 1\t0.01\n";
2664 :    
2665 :     #close reaction equation file
2666 :     close(REACTIONS);
2667 :    
2668 :     # FluxAnalyzer requires this file
2669 :     open(APP, ">$dir/app_para.m");
2670 :     print APP "epsilon=1e-10;\nbasic_color=[0.7 0.7 0.7];\ncr_color=[0.5 0.5 1];\nbr_color=[1 0.2 0.2];\nnbr_color=[0.2 1 0.2];\ntext_color=[0 0 0];\nmacro_synth_color=[0 0 1];\nmacro_color=[0.6 0.6 1];\nbox_reaction_width=[0.06];\nbox_reaction_height=[0.03];\nbox_macro_width=[0.06];\nbox_macro_height=[0.03];\nfontsize_reaction=[11];\nfontsize_macro=[11];\nfluxmaps={'Fluxmap','dummy_medium.pcx'};\n";
2671 :     close(APP);
2672 :     }
2673 :    
2674 :    
2675 :     sub clear_arrays
2676 :     {
2677 :     undef %reactions_to_substrate_arrays;
2678 :     undef %reactions_to_product_arrays;
2679 :     undef %all_compounds_to_main;
2680 :     undef %all_reactions;
2681 :     undef %scenario_cycles;
2682 :     undef @all_outputs_lists;
2683 :     undef %all_inputs;
2684 :     undef %all_outputs;
2685 :    
2686 :     %reactions_to_substrate_arrays = ();
2687 :     %reactions_to_product_arrays = ();
2688 :     %all_compounds_to_main = ();
2689 :     %all_reactions = ();
2690 :     %scenario_cycles = ();
2691 :     @all_outputs_lists = ();
2692 :     %all_inputs = ();
2693 :     %all_outputs = ();
2694 :     }
2695 :    
2696 :     sub load_superset_file
2697 :     {
2698 :     %superset_to_ss = ();
2699 :     %ss_to_superset = ();
2700 :    
2701 : dejongh 1.20 open(FILE,"<$FIG_Config::global/Models/hope_supersets.txt") or die("Failed to open hope_supersets.txt");
2702 : olson 1.1 while(<FILE>)
2703 :     {
2704 :     my @line = split(/\t/,$_);
2705 :     map { s/ /_/g } @line;
2706 :     map { s/\"//g } @line;
2707 :     map { chomp } @line;
2708 :     $superset_to_ss{$line[0]} = [] if !defined $superset_to_ss{$line[0]};
2709 :     $ss_to_superset{$line[1]} = $line[0];
2710 :     push(@{$superset_to_ss{$line[0]}},$line[1]);
2711 :     }
2712 :     close(FILE);
2713 :    
2714 : dejongh 1.24 return \%ss_to_superset
2715 : olson 1.1 }
2716 :    
2717 :     # This function runs a given scenario that is defined the specified subsystem
2718 :     # It returns the data as it is found by process_paths in a hash reference
2719 :     # This is used internally by model.pm, and shouldn't be called externally
2720 :    
2721 :     sub load_scenario
2722 :     {
2723 : dejongh 1.19 my ($genome,$ssa,$scenario) = @_;
2724 : olson 1.1
2725 :     #load up the arrays with the info we need
2726 : dejongh 1.19 process_init($ssa,$scenario,$genome,0);
2727 : olson 1.1
2728 :     #assume all path inputs and outputs are main
2729 :     map { $all_compounds_to_main{$_} = 1 } keys %all_inputs;
2730 :     map { $all_compounds_to_main{$_} = 1 } keys %all_outputs;
2731 :    
2732 :     }
2733 :    
2734 :     sub internal_scenario
2735 :     {
2736 : dejongh 1.19 my ($genome,$ssa,$scenario,$find_first) = @_;
2737 : olson 1.1
2738 : olson 1.2 print STDERR "\nIn internal_scenario with '$genome', '$ssa', '$scenario', '$find_first'\n" if $debug;
2739 : olson 1.1
2740 : dejongh 1.19 load_scenario($genome,$ssa,$scenario);
2741 : olson 1.1
2742 :     return execute_paths([],$find_first,[],[]);
2743 :     }
2744 :    
2745 :     sub run_scenario
2746 :     {
2747 : dejongh 1.19 my($genome,$superset,$subsystem,$scenario,$find_first) = @_;
2748 : dejongh 1.23 my $scenario_dir = get_scenario_directory($genome) . "/PathInfo/$superset/$subsystem/$scenario";
2749 : olson 1.2
2750 :     system("rm", "-rf", $scenario_dir);
2751 : dejongh 1.18 &FIG::verify_dir($scenario_dir);
2752 : olson 1.2
2753 : olson 1.1 #make sure the arrays are empty to start out
2754 :     &clear_arrays;
2755 : dejongh 1.19 write_scenario(internal_scenario($genome,$subsystem,$scenario,$find_first),$scenario_dir);
2756 : olson 1.1 }
2757 :    
2758 : olson 1.2 sub compare_scenario
2759 :     {
2760 : dejongh 1.17 my($genome,$superset,$ss_name,$scenario_name,$dont_copy) = @_;
2761 : olson 1.2 my @genome_paths;
2762 : dejongh 1.23 my $scenario_dir_all = get_scenario_directory('All') . "/PathInfo/$superset/$ss_name/$scenario_name";
2763 : olson 1.2 my $subsystem = $fig->get_subsystem($ss_name);
2764 : dejongh 1.17 my @additional_reactions = $subsystem->get_hope_additional_reactions($scenario_name);
2765 : olson 1.2 my %additional_reactions;
2766 :     map { $additional_reactions{$_} = 1 } @additional_reactions;
2767 :    
2768 : dejongh 1.4 my %ss_reactions;
2769 :    
2770 :     if ($genome eq "All")
2771 : olson 1.2 {
2772 : dejongh 1.17 my %all_reactions = $subsystem->get_hope_reactions;
2773 : dejongh 1.4 foreach my $role (keys %all_reactions)
2774 : olson 1.2 {
2775 : dejongh 1.4 map { $ss_reactions{$_} = 1 } @{$all_reactions{$role}};
2776 : olson 1.2 }
2777 :     }
2778 : dejongh 1.4 else
2779 : olson 1.2 {
2780 : dejongh 1.24 my %reactions_for_genome = get_reactions_for_genome_in_subsystem($genome,$ss_name);
2781 : parrello 1.32 map { $ss_reactions{$_} = 1 } keys %reactions_for_genome if keys %reactions_for_genome;
2782 : olson 1.2 }
2783 :    
2784 :     # first find paths in the All directory that should be valid for the genome
2785 :     # based on the reactions associated with it in the subsystems
2786 : dejongh 1.17 opendir (DIR_ALL,$scenario_dir_all) or return [];
2787 : dejongh 1.18 my @sub_dirs = readdir DIR_ALL;
2788 : olson 1.2 close DIR_ALL;
2789 :    
2790 :     my %paths_all;
2791 :    
2792 :     for my $path (@sub_dirs)
2793 :     {
2794 : dejongh 1.18 next if $path !~ /path/; # skip . and ..
2795 : olson 1.2 my $match = 1;
2796 :     open (PATH, "$scenario_dir_all/$path/path_info");
2797 :     my @reactions = <PATH>;
2798 :     close PATH;
2799 :    
2800 :     my $reaction_string = "";
2801 :    
2802 :     foreach my $reaction (sort @reactions)
2803 :     {
2804 :     if ($reaction =~ /(R\d\d\d\d\d)/)
2805 :     {
2806 :     $reaction_string .= $reaction;
2807 :    
2808 :     if (! exists($ss_reactions{$1}) && ! exists($additional_reactions{$1}))
2809 :     {
2810 :     $match = 0;
2811 :     }
2812 :     }
2813 :     }
2814 :    
2815 :     $paths_all{$reaction_string} = $match;
2816 :     push @genome_paths, $path if $match;
2817 :     }
2818 :    
2819 :     # now check all paths found for this particular organism, and make sure they
2820 :     # are in the appropriate All subdirectory
2821 : dejongh 1.23 my $scenario_dir_genome = get_scenario_directory($genome). "/PathInfo/$superset/$ss_name/$scenario_name";
2822 : olson 1.2
2823 : dejongh 1.17 opendir (DIR_GENOME,$scenario_dir_genome) or return \@genome_paths;
2824 : dejongh 1.18 @sub_dirs = readdir DIR_GENOME;
2825 : olson 1.2 close DIR_GENOME;
2826 :    
2827 :     my $path_counter = scalar keys %paths_all;
2828 :    
2829 :     for my $path (@sub_dirs)
2830 :     {
2831 : dejongh 1.18 next if $path !~ /path/; # skip . and ..
2832 : olson 1.2 my $match = 1;
2833 :     open (PATH, "$scenario_dir_genome/$path/path_info");
2834 :     my @reactions = <PATH>;
2835 :     close PATH;
2836 :    
2837 :     my $reaction_string = "";
2838 :    
2839 :     foreach my $reaction (sort @reactions)
2840 :     {
2841 :     if ($reaction =~ /(R\d\d\d\d\d)/)
2842 :     {
2843 :     $reaction_string .= $reaction;
2844 :     }
2845 :     }
2846 :    
2847 :     if (! exists($paths_all{$reaction_string}))
2848 :     {
2849 : dejongh 1.17 if ($dont_copy)
2850 :     {
2851 :     print STDERR "$scenario_dir_genome/$path not found in All\n";
2852 :     }
2853 :     else
2854 :     {
2855 :     $path_counter++;
2856 :     my $new_path_name = "path_".$path_counter;
2857 :     my $temp_sdg = $scenario_dir_genome;
2858 :     $temp_sdg =~ s/\(/\\\(/g;
2859 :     $temp_sdg =~ s/\)/\\\)/g;
2860 :     my $temp_sda = $scenario_dir_all;
2861 :     $temp_sda =~ s/\(/\\\(/g;
2862 :     $temp_sda =~ s/\)/\\\)/g;
2863 :     `cp -R $temp_sdg/$path $temp_sda/$new_path_name`;
2864 :     push @genome_paths, $new_path_name;
2865 :     print STDERR "Copied $temp_sdg/$path to $temp_sda/$new_path_name\n";
2866 :     }
2867 :     }
2868 :    
2869 :     unless ($dont_copy)
2870 :     {
2871 :     # remove genome-specific paths
2872 :     rmtree("$scenario_dir_genome/$path");
2873 :     }
2874 :     }
2875 :    
2876 :     unless ($dont_copy)
2877 :     {
2878 : dejongh 1.26 # copy each genome-specific path, with the same name as the path in the "All" directory
2879 : dejongh 1.17 foreach my $path (@genome_paths)
2880 :     {
2881 : olson 1.2 my $temp_sdg = $scenario_dir_genome;
2882 :     $temp_sdg =~ s/\(/\\\(/g;
2883 :     $temp_sdg =~ s/\)/\\\)/g;
2884 :     my $temp_sda = $scenario_dir_all;
2885 :     $temp_sda =~ s/\(/\\\(/g;
2886 :     $temp_sda =~ s/\)/\\\)/g;
2887 : dejongh 1.26 `cp -r $temp_sda/$path $temp_sdg`;
2888 : dejongh 1.17 }
2889 :     }
2890 : olson 1.2
2891 :     return \@genome_paths;
2892 :     }
2893 :    
2894 : olson 1.1 sub write_scenario
2895 :     {
2896 :     my($scenario_data,$scenario_dir) = @_;
2897 :     delete $scenario_data->{"infinite"};
2898 :     my $path_count = 1;
2899 :     my @list_of_done_tokens=();
2900 :    
2901 :     print STDERR "Paths: ", keys %{$scenario_data}, "\n" if $debug;
2902 :    
2903 :     foreach my $try_path (keys %{$scenario_data})
2904 :     {
2905 :     print STDERR "\t Checking $try_path\n" if $debug;
2906 :     if($scenario_data->{$try_path}->[1] != 1)
2907 :     {
2908 :     print STDERR "\t Token $try_path is not complete\n" if $debug;
2909 :     next;
2910 :     }
2911 :     print STDERR "These are the contents of the list: " , @list_of_done_tokens, "\n" if $debug;
2912 :     if(scalar @list_of_done_tokens == 0){
2913 :     push @list_of_done_tokens, $try_path;
2914 :     next;
2915 :     }
2916 :    
2917 :     # Check this token's ($try_path) values against the values of the keys($elem) stored in
2918 :     # list_of_done_tokens. If they match, don't add it to the finished token list, if it doesn't
2919 :     # match, add it.
2920 :    
2921 :     print STDERR "These are the contents of the list: " , @list_of_done_tokens, "\n" if $debug;
2922 :    
2923 :     my $found_match = 0;
2924 :    
2925 :     foreach my $elem(@list_of_done_tokens){
2926 :     my @done_reactions = @{$scenario_data->{$elem}->[2]};
2927 :    
2928 :     print STDERR "This is the path we're trying: " , $try_path, "\t","This is the path already in the array: " , $elem , "\nThis is the size of the array: " . @list_of_done_tokens ."\n" if $debug;
2929 :    
2930 :     # if the list of reactions match, they represent the same path
2931 :    
2932 :     my @path_reactions = @{$scenario_data->{$try_path}->[2]};
2933 :     my (%diff_reactions_1, %diff_reactions_2);
2934 :    
2935 :     map {$diff_reactions_1{$_} = 1} @path_reactions;
2936 :     map {delete $diff_reactions_1{$_}} @done_reactions;
2937 :     map {$diff_reactions_2{$_} = 1} @done_reactions;
2938 :     map {delete $diff_reactions_2{$_}} @path_reactions;
2939 :    
2940 :     if (scalar keys %diff_reactions_1 == 0 && scalar keys %diff_reactions_2 == 0)
2941 :     {
2942 :     print STDERR "They match.\n" if $debug;
2943 :     $found_match = 1;
2944 :     last;
2945 :     }
2946 :     }
2947 :    
2948 :     if (! $found_match)
2949 :     {
2950 :     push @list_of_done_tokens, $try_path;
2951 :     print STDERR $try_path, " Added to the array\n" if $debug;
2952 :    
2953 :     }
2954 :     }
2955 :    
2956 :     foreach my $path (@list_of_done_tokens){
2957 :     if($scenario_data->{$path}->[1] != 1)
2958 :     {
2959 :     next;
2960 :     }
2961 :     if(@list_of_done_tokens == 0){
2962 :     push @list_of_done_tokens, $path;
2963 :     }
2964 :    
2965 :     #create input/output info
2966 :     my $input_hash = $scenario_data->{$path}->[7];
2967 :     my $output_hash = $scenario_data->{$path}->[8];
2968 :     my $reaction_path = $scenario_data->{$path}->[2];
2969 :     my @reaction_array = @$reaction_path;
2970 :    
2971 :     print STDERR "\nInputs:\n" if $debug;
2972 :     print STDERR map{"$_ => $input_hash->{$_}" } keys %$input_hash, "\n" if $debug;
2973 :     print STDERR "\nOutputs:\n" if $debug;
2974 :     print STDERR map{"$_ => $output_hash->{$_}"} keys %$output_hash, "\n" if $debug;
2975 :    
2976 :     # divide stoichiometry by greatest common denominator
2977 :     my ($min_stoich, @all_stoichs);
2978 :    
2979 :     map { push @all_stoichs, $input_hash->{$_}; $min_stoich = $input_hash->{$_} if $input_hash->{$_} < $min_stoich || $min_stoich == 0 } keys %{$input_hash};
2980 :     map { push @all_stoichs, $output_hash->{$_}; $min_stoich = $output_hash->{$_} if $output_hash->{$_} < $min_stoich || $min_stoich == 0 } keys %{$output_hash};
2981 :    
2982 :     my ($gcd, @gcd_candidates);
2983 :    
2984 :     outer: for ($gcd = $min_stoich; $gcd > 1; $gcd--)
2985 :     {
2986 :     foreach my $stoich (@all_stoichs)
2987 :     {
2988 :     next outer if $stoich % $gcd != 0;
2989 :     }
2990 :    
2991 :     last; # found a gcd
2992 :     }
2993 :    
2994 :     map { $input_hash->{$_} /= $gcd } keys %{$input_hash};
2995 :     map { $output_hash->{$_} /= $gcd } keys %{$output_hash};
2996 :    
2997 :     mkdir "$scenario_dir/path_$path_count";
2998 :    
2999 :     open(FILE, ">$scenario_dir/path_$path_count/path_info");
3000 :     foreach my $elem(@reaction_array){
3001 :     print FILE scalar @reaction_array, "\t", $elem, "\n";
3002 :     }
3003 :     close(FILE);
3004 :    
3005 :     &write_fluxanalyzer_files("$scenario_dir/path_$path_count",$input_hash,
3006 :     $output_hash, \@reaction_array,\%all_reactions,
3007 :     \%reactions_to_substrate_arrays,
3008 :     \%reactions_to_product_arrays,
3009 :     {});
3010 :     $path_count++;
3011 :     }
3012 : dejongh 1.25 print STDERR "no paths for $scenario_dir\n" if $path_count == 1 && $debug;
3013 : olson 1.1 @list_of_done_tokens=();
3014 :     undef %{$scenario_data};
3015 :    
3016 :     }
3017 :    
3018 :    
3019 :     sub load_subsystem
3020 :     {
3021 :     my ($genome,$ss_name) = @_;
3022 :     my $subsystem = $fig->get_subsystem($ss_name);
3023 :     my @ss_scenarios = $subsystem->get_hope_scenario_names;
3024 :     foreach my $name (@ss_scenarios)
3025 :     {
3026 :     load_scenario($genome,$ss_name,$name);
3027 :     }
3028 :     }
3029 :    
3030 :     sub internal_subsystem
3031 :     {
3032 :     my ($genome,$ss_name,$find_first) = @_;
3033 :     my %scenario_to_paths;
3034 :    
3035 :     my $subsystem = $fig->get_subsystem($ss_name);
3036 :     my @ss_scenarios = $subsystem->get_hope_scenario_names;
3037 :    
3038 :     foreach my $name (@ss_scenarios)
3039 :     {
3040 :     $scenario_to_paths{$name} = internal_scenario($genome,$ss_name,$name,$find_first);
3041 :     }
3042 :    
3043 :     return \%scenario_to_paths;
3044 :     }
3045 :    
3046 :     sub run_subsystem
3047 :     {
3048 :     my ($genome,$superset,$subsystem,$find_first) = @_;
3049 :    
3050 :     my $subsystem_obj = $fig->get_subsystem($subsystem);
3051 : olson 1.10
3052 :     if (!$subsystem_obj)
3053 :     {
3054 :     warn "Cannot open subsystem $subsystem\n";
3055 :     return;
3056 :     }
3057 :    
3058 : olson 1.1 my @ss_scenarios = $subsystem_obj->get_hope_scenario_names;
3059 : olson 1.10
3060 : dejongh 1.23 my $dir = get_scenario_directory($genome) . "/PathInfo/$superset/$subsystem";
3061 : olson 1.10 system("rm", "-rf", $dir);
3062 :     &FIG::verify_dir($dir);
3063 : olson 1.1
3064 :     foreach my $name (@ss_scenarios)
3065 :     {
3066 :     run_scenario($genome,$superset,$subsystem,$name,$find_first);
3067 :     }
3068 :     }
3069 :    
3070 : olson 1.2 sub compare_subsystem
3071 :     {
3072 : dejongh 1.17 my ($genome,$superset,$subsystem,$dont_copy) = @_;
3073 : olson 1.2 my %genome_scenarios;
3074 :    
3075 :     my $subsystem_obj = $fig->get_subsystem($subsystem);
3076 : olson 1.11
3077 :     if (!$subsystem_obj)
3078 :     {
3079 :     warn "Cannot open subsystem $subsystem\n";
3080 :     return;
3081 :     }
3082 :    
3083 : olson 1.2 my @ss_scenarios = $subsystem_obj->get_hope_scenario_names;
3084 :    
3085 :     foreach my $name (@ss_scenarios)
3086 :     {
3087 : dejongh 1.17 $genome_scenarios{$name} = compare_scenario($genome,$superset,$subsystem,$name,$dont_copy);
3088 : olson 1.2 }
3089 :    
3090 :     return \%genome_scenarios;
3091 :     }
3092 :    
3093 : olson 1.1 sub load_superset
3094 :     {
3095 :     my($genome, $superset_name) = @_;
3096 :    
3097 :     my @subsystems = @{$superset_to_ss{$superset_name}};
3098 :     foreach my $ss_name (@subsystems)
3099 :     {
3100 :     load_subsystem($genome,$ss_name);
3101 :     }
3102 :     }
3103 :    
3104 :     sub internal_superset
3105 :     {
3106 :     my($genome, $superset_name,$find_first) = @_;
3107 :    
3108 :     my @subsystems = @{$superset_to_ss{$superset_name}};
3109 :    
3110 :     my %supersets_data;
3111 :    
3112 :     foreach my $ss_name (@subsystems)
3113 :     {
3114 :     $supersets_data{$ss_name} = internal_subsystem($genome,$ss_name,$find_first);
3115 :     }
3116 :    
3117 :     return \%supersets_data;
3118 :     }
3119 :    
3120 :     sub run_superset
3121 :     {
3122 :     my($genome, $superset_name,$find_first) = @_;
3123 :    
3124 :     my @subsystems = @{$superset_to_ss{$superset_name}};
3125 :    
3126 : olson 1.10
3127 : dejongh 1.23 my $dir = get_scenario_directory($genome) . "/PathInfo/$superset_name";
3128 : olson 1.10 system("rm", "-rf", $dir);
3129 :     &FIG::verify_dir($dir);
3130 : olson 1.1
3131 :     foreach my $ss_name (@subsystems)
3132 :     {
3133 : dejongh 1.30 print STDERR "Running Scenarios for $genome for subsystem $ss_name\n";
3134 : olson 1.1 run_subsystem($genome,$superset_name,$ss_name,$find_first);
3135 :     }
3136 :     }
3137 :    
3138 : olson 1.2 sub compare_superset
3139 :     {
3140 : dejongh 1.17 my($genome, $superset_name, $dont_copy) = @_;
3141 : olson 1.2
3142 :     my @subsystems = @{$superset_to_ss{$superset_name}};
3143 :     my %genome_subsystems;
3144 :    
3145 :     foreach my $ss_name (@subsystems)
3146 :     {
3147 :     print STDERR "Comparing Scenarios for $genome in subsystem $ss_name\n";
3148 : dejongh 1.17 $genome_subsystems{$ss_name} = compare_subsystem($genome,$superset_name,$ss_name,$dont_copy);
3149 : olson 1.2 }
3150 :    
3151 :     return \%genome_subsystems;
3152 :     }
3153 :    
3154 : olson 1.1
3155 :     sub load_supersets
3156 :     {
3157 :     my($genome) = @_;
3158 : dejongh 1.20 &load_superset_file;
3159 :    
3160 : olson 1.1 foreach my $superset (keys %superset_to_ss)
3161 :     {
3162 :     load_superset($genome,$superset);
3163 :     }
3164 :    
3165 :     return (\%all_reactions,\%reactions_to_substrate_arrays,\%reactions_to_product_arrays);
3166 :    
3167 :     }
3168 :    
3169 :     sub run_supersets
3170 :     {
3171 :     my($genome,$find_first) = @_;
3172 : dejongh 1.20 &load_superset_file;
3173 :    
3174 : dejongh 1.23 my $dir = get_scenario_directory($genome) . "/PathInfo";
3175 : dejongh 1.20 system("rm", "-rf", $dir);
3176 :     &FIG::verify_dir($dir);
3177 :    
3178 : olson 1.1 foreach my $superset (keys %superset_to_ss)
3179 :     {
3180 :     run_superset($genome,$superset,$find_first);
3181 :     }
3182 :    
3183 :     }
3184 :    
3185 : olson 1.2 sub compare_supersets
3186 :     {
3187 : dejongh 1.17 my($genome, $dont_copy) = @_;
3188 : olson 1.2 my %genome_supersets;
3189 :    
3190 : dejongh 1.20 &load_superset_file;
3191 :    
3192 : olson 1.2 foreach my $superset (keys %superset_to_ss)
3193 :     {
3194 : dejongh 1.17 $genome_supersets{$superset} = compare_superset($genome,$superset,$dont_copy);
3195 : olson 1.2 }
3196 :    
3197 :     return \%genome_supersets;
3198 :     }
3199 :    
3200 : olson 1.1 sub run_genome_report
3201 :     {
3202 :     my ($genome) = @_;
3203 :     my @string_out;
3204 :     push @string_out,"Genome $genome\n";
3205 :     #get all the subsystems this genome is involved in
3206 :     foreach my $superset (keys %superset_to_ss)
3207 :     {
3208 :     foreach my $name (@{$superset_to_ss{$superset}})
3209 :     {
3210 :     push @string_out, @{print_ss_report($name,internal_subsystem($genome,$name,0))};
3211 :     }
3212 :     }
3213 :     return \@string_out;
3214 :     }
3215 :    
3216 :     sub print_ss_report
3217 :     {
3218 :     my ($ss_name,$scenario_to_paths) = @_;
3219 :     my @output = ();
3220 :     my %scenario_path_count;
3221 :    
3222 :     push(@output,"\tSubsystem $ss_name\n");
3223 :    
3224 :     foreach my $scenario (keys %$scenario_to_paths)
3225 :     {
3226 :     if($scenario_to_paths->{$scenario}->{"infinite"})
3227 :     {
3228 :     push(@output,"\tWarning: Possible Infinite loop\n");
3229 :     }
3230 :     delete $scenario_to_paths->{$scenario}->{"infinite"};
3231 :     foreach my $token (keys %{$scenario_to_paths->{$scenario}})
3232 :     {
3233 :     $scenario_path_count{$scenario}++ if ($scenario_to_paths->{$scenario}->{$token}->[1]);
3234 :     }
3235 :     }
3236 :    
3237 :     push @output , map { "\t\t$_ has $scenario_path_count{$_} path(s).\n" } keys %scenario_path_count;
3238 :    
3239 :     return \@output;
3240 :     }
3241 :    
3242 :    
3243 :     sub internal_assembly
3244 :     {
3245 :     #lets get the genome, and a array reference to the paths we want to build togather
3246 :     my ($paths,$input_path,$output_path,$one_path) = @_;
3247 :    
3248 :    
3249 : olson 1.2 print STDERR $paths."\n" if $debug;
3250 : olson 1.1 clear_arrays();
3251 :    
3252 :     #This gets us an array of arrays, each subarray holds
3253 : dejongh 1.12 # [0] = genome [1] = Scenarios [2] = superset [3] = subsystem [4] = scenario [5] = path
3254 : olson 1.1 # OR [0] = genome [1] = assembly [2] = path_name
3255 :     my @assembly_scenarios = @{parse_assembly_scenarios($paths)};
3256 :    
3257 :     #split and the input/output path for later as well
3258 :     #these should only return one path array...so just grab that one
3259 :     my @input_arr;
3260 :     my @output_arr;
3261 :    
3262 :     if($input_path != undef && $output_path != undef )
3263 :     {
3264 :     @input_arr = @{parse_assembly_scenarios($input_path)};
3265 :     @output_arr = @{parse_assembly_scenarios($output_path)};
3266 :     }
3267 :    
3268 :     #load all the kegg information for each 'scenario' from the paths we have selected
3269 :     foreach my $scenario (@assembly_scenarios)
3270 :     {
3271 : dejongh 1.12 if(scalar @$scenario > 5) #this is a normal scenario path
3272 : olson 1.1 {
3273 : dejongh 1.12 print STDERR "Checking $scenario->[3] $scenario->[4] $scenario->[5] \n" if $debug;
3274 :     process_init($scenario->[3],$scenario->[4],$scenario->[0],1);
3275 : olson 1.1 }
3276 :     else #This is a assembly path
3277 :     {
3278 :     #read in the input/output compounds and mark main's correctly
3279 : dejongh 1.20 my $genome = shift @$scenario;
3280 : dejongh 1.23 my $path = get_scenario_directory($genome) . "/" . join "/" , @$scenario;
3281 : olson 1.1
3282 :     open(M_IN,"$path/inputs_main");
3283 :     while(<M_IN>)
3284 :     {
3285 :     my @line = split(/\t/,$_);
3286 :     $all_compounds_to_main{$line[0]} = $line[1];
3287 :     }
3288 :     close(M_IN);
3289 :    
3290 :     open(M_OUT,"$path/outputs_main");
3291 :     while(<M_OUT>)
3292 :     {
3293 :     my @line = split(/\t/,$_);
3294 :     $all_compounds_to_main{$line[0]} = $line[1];
3295 :     }
3296 :     close(M_OUT);
3297 :     }
3298 :     }
3299 :     #assume all path inputs and outputs are main
3300 :     map { $all_compounds_to_main{$_} = 1 } keys %all_inputs;
3301 :     map { $all_compounds_to_main{$_} = 1 } keys %all_outputs;
3302 :    
3303 :     print STDERR "Inputs: " if $debug;
3304 : olson 1.2 print STDERR map { $_."\n" } keys %all_inputs if $debug;
3305 : olson 1.1 print STDERR "Outputs: " if $debug;
3306 : olson 1.2 print STDERR map { $_."\n" } keys %all_outputs if $debug;
3307 : olson 1.1
3308 :     #run process paths
3309 :     return execute_paths(\@assembly_scenarios,$one_path,\@input_arr,\@output_arr);
3310 :     }
3311 :    
3312 :     sub run_assembly
3313 :     {
3314 :     my ($paths,$genome,$write_name,$one_path) = @_;
3315 :     print STDERR "\nThis is the passed information.\n";
3316 :     print STDERR $paths , "\n" , @$paths, "\n";
3317 :     print STDERR $genome . "\n";
3318 :     print STDERR $write_name . "\n";
3319 : dejongh 1.23 my $dir = get_scenario_directory($genome) . "/Assemblies/$write_name";
3320 : olson 1.10 system("rm", "-rf", $dir);
3321 :     &FIG::verify_dir($dir);
3322 :     write_scenario(internal_assembly($paths,[],[],$one_path),$dir);
3323 : olson 1.1 }
3324 :    
3325 :     sub expand_paths
3326 :     {
3327 :     my ($paths) = @_;
3328 :    
3329 :     my @final_paths;
3330 :    
3331 :     foreach my $path (@$paths)
3332 :     {
3333 :     if($path eq "" || $path eq "//")
3334 :     {
3335 :     next;
3336 :     }
3337 : dejongh 1.12 my $length = 6;
3338 :     if($path =~ /Assemblies/)
3339 : olson 1.1 {
3340 :     $length = 4;
3341 :     }
3342 :     print STDERR "Expanding $path \n" if $debug;
3343 :     my @parts = split "/", $path;
3344 :     shift @parts; # get ride of the first blank entry from /$genome
3345 : parrello 1.32 $length =$length - scalar @parts;
3346 : dejongh 1.23 my $genome = shift @parts;
3347 : olson 1.1 $path = join "/" , @parts;
3348 : dejongh 1.23 print STDERR "Length : $length Path: $path\n" if $debug;
3349 :     my @temp = @{expand_recursive($genome, $path, $length)};
3350 : olson 1.1 push @final_paths , @temp if scalar @temp > 0;
3351 :     }
3352 :    
3353 :    
3354 :     return \@final_paths;
3355 :     }
3356 :    
3357 :     sub expand_recursive
3358 :     {
3359 : dejongh 1.23 my ($genome,$path,$count) = @_;
3360 : olson 1.1 my @sub_dirs;
3361 : dejongh 1.23 if($path =~ m\Assemblies$\ || $path =~ m\Analysis$\ || $path =~ m\Curation$\)
3362 : olson 1.1 {
3363 :     return [];
3364 :     }
3365 :     if($count !=0)
3366 :     {
3367 : dejongh 1.20 #read this path, and pull out all the sub-directories.
3368 : dejongh 1.23 my $model_dir = get_scenario_directory($genome) . "/$path/";
3369 : olson 1.10
3370 :     print STDERR "reading directory $model_dir\n" if $debug;
3371 : formsma 1.16 opendir (DIR, $model_dir) or die("$model_dir");
3372 : olson 1.10 print STDERR "reading directory $model_dir\n" if $debug;
3373 : dejongh 1.18 @sub_dirs = readdir DIR;
3374 : olson 1.1 close DIR;
3375 :     print STDERR "Found: " if $debug;
3376 :     print STDERR @sub_dirs , "," if $debug;
3377 :     }
3378 :     else
3379 :     {
3380 :    
3381 :     return [$path];
3382 :     }
3383 :     $count--;
3384 :     my @to_return;
3385 :     foreach my $sub_path (@sub_dirs)
3386 :     {
3387 : dejongh 1.18 next if $sub_path =~ /^\.$/ || $sub_path =~ /^\.\.$/;
3388 : olson 1.1 print STDERR "Calling on $sub_path , $count \n" if $debug;
3389 : dejongh 1.23 push @to_return , @{expand_recursive($genome, "$path/$sub_path",$count)};
3390 : olson 1.1 }
3391 :     print STDERR "Returning" ,@to_return if $debug;
3392 :     return \@to_return;
3393 :    
3394 :     }
3395 :    
3396 :     sub parse_assembly_scenarios
3397 :     {
3398 :     my ($paths) = @_;
3399 :    
3400 :     $paths = expand_paths($paths);
3401 :     my @array_of_path_arrays;
3402 :    
3403 :     foreach my $path (@$paths)
3404 :     {
3405 :     my @parts = split "/", $path;
3406 :     #shift @parts;
3407 :     push(@array_of_path_arrays, \@parts);
3408 :     }
3409 :    
3410 :     return \@array_of_path_arrays;
3411 :     }
3412 :    
3413 :    
3414 :     sub write_selected_scenarios
3415 :     {
3416 :     my($checked,$genome,$ssa,$sc_name) = @_;
3417 :     my (@tempArray);
3418 :    
3419 :     #Load this scenario again with all of its rxns and cpds
3420 :    
3421 :     model::clear_arrays();
3422 :    
3423 :     model::process_init($ssa,$sc_name,$genome,0);
3424 :    
3425 :     #assume all path inputs and outputs are main
3426 :     map { $all_compounds_to_main{$_} = 1 } keys %all_inputs;
3427 :     map { $all_compounds_to_main{$_} = 1 } keys %all_outputs;
3428 :    
3429 :     model::create_reactions({},{},{},{});
3430 :    
3431 :     ##End of scenario loading
3432 :    
3433 :    
3434 :     #setup the filesystem to store the scenario/paths
3435 :     my $superset = $ss_to_superset{$ssa};
3436 : dejongh 1.23 my $base_dir = get_scenario_directory($genome) . "/PathInfo/$superset/$ssa/$sc_name/";
3437 : olson 1.2 system("rm", "-rf", $base_dir);
3438 : olson 1.10 &FIG::verify_dir($base_dir);
3439 : olson 1.1
3440 :     #for the selected paths, lets gather their cpd from the checkbox and write
3441 :     #the path that we need
3442 :     foreach my $path (@$checked)
3443 :     {
3444 :     my (%input_hash, %output_hash, $path_name);
3445 :     #process the strings to get the information from the parameters
3446 :     my @items = split(";", $path);
3447 :     $path_name = $items[0];
3448 :     #next we have the input compounds ids/stoich/main
3449 :     map { if ($_ =~ /(.*):(.*):(.*)/)
3450 :     { $input_hash{$1}+= $2 } } split ",", $items[1];
3451 :     #the third part has the output compounds ids/stoich/main
3452 :     map { if ($_ =~ /(.*):(.*):(.*)/)
3453 :     { $output_hash{$1} += $2 } } split ",",$items[2];
3454 :     #the fourth part has the strings of reactions visited
3455 :     @tempArray = split("#" , $items[3]);
3456 :    
3457 :    
3458 :     print STDERR "\nInputs:\n" if $debug;
3459 :     print STDERR map{"$_ => $input_hash{$_}" } keys %input_hash, "\n" if $debug;
3460 :     print STDERR "\nOutputs:\n" if $debug;
3461 :     print STDERR map{"$_ => $output_hash{$_}"} keys %output_hash, "\n" if $debug;
3462 :    
3463 :    
3464 :     # divide stoichiometry by greatest common denominator
3465 :     my ($min_stoich, @all_stoichs);
3466 :    
3467 :     map { push @all_stoichs, $input_hash{$_}; $min_stoich = $input_hash{$_} if $input_hash{$_} < $min_stoich || $min_stoich == 0 } keys %input_hash;
3468 :     map { push @all_stoichs, $output_hash{$_}; $min_stoich = $output_hash{$_} if $output_hash{$_} < $min_stoich || $min_stoich == 0 } keys %output_hash;
3469 :    
3470 :     my ($gcd, @gcd_candidates);
3471 :    
3472 :     outer: for ($gcd = $min_stoich; $gcd > 1; $gcd--)
3473 :     {
3474 :     foreach my $stoich (@all_stoichs)
3475 :     {
3476 :     next outer if $stoich % $gcd != 0;
3477 :     }
3478 :    
3479 :     last; # found a gcd
3480 :     }
3481 :    
3482 :     map { $input_hash{$_} /= $gcd } keys %input_hash;
3483 :     map { $output_hash{$_} /= $gcd } keys %output_hash;
3484 :    
3485 :     mkdir "$base_dir/$path_name";
3486 :     &write_fluxanalyzer_files("$base_dir/$path_name",\%input_hash,
3487 :     \%output_hash,\@tempArray,\%all_reactions,
3488 :     \%reactions_to_substrate_arrays,
3489 :     \%reactions_to_product_arrays,
3490 :     {});
3491 :     }
3492 :     return $base_dir;
3493 :     }
3494 :    
3495 :     #This write function assumes that we have just run a assembly (and we haven't cleared the arrays)
3496 :     # becuase the write_fluxanalyzer_files function is depended on those global arrays for the reactions.
3497 :    
3498 :     sub write_assembly
3499 :     {
3500 :     my($input,$genome,$name) = @_;
3501 :    
3502 :     my $paths = $input->[1];
3503 :     my $file_paths = $input->[0];
3504 :     print STDERR "Paths: @$paths \n File_Dirs: @$file_paths \n" if $debug;
3505 :     my @tempArray;
3506 :    
3507 :     chomp $genome;
3508 :     chomp $name;
3509 :     #setup the filesystem to store the assembly
3510 : dejongh 1.23 my $base_dir = get_scenario_directory($genome) . "/Assemblies/$name";
3511 : olson 1.10 system("rm", "-rf", $base_dir);
3512 :     &FIG::verify_dir($base_dir);
3513 :    
3514 : olson 1.1
3515 :     ##Here we want to reload the cpd and rxn info so we can write it later
3516 :     clear_arrays();
3517 :    
3518 :     #This gets us an array of arrays, each subarray holds
3519 : dejongh 1.12 # [0] = genome [1] = Scenarios [2] = superset [3] = subsystem [4] = scenario [5] = path
3520 : olson 1.1 my @assembly_scenarios = @{parse_assembly_scenarios($file_paths)};
3521 :    
3522 :     #load all the kegg information for each 'scenario' from the paths we have selected
3523 :     foreach my $scenario (@assembly_scenarios)
3524 :     {
3525 : dejongh 1.12 print STDERR "Checking $scenario->[3] $scenario->[4] $scenario->[5] \n" if $debug;
3526 :     process_init($scenario->[3],$scenario->[4],$scenario->[0],1);
3527 : olson 1.1 }
3528 :     #assume all path inputs and outputs are main
3529 :     map { $all_compounds_to_main{$_} = 1 } keys %all_inputs;
3530 :     map { $all_compounds_to_main{$_} = 1 } keys %all_outputs;
3531 :    
3532 :     create_assembly_reactions({},{},{},{});
3533 :    
3534 :     ##End of rxn,cpd loading
3535 :    
3536 :     #for the selected paths, lets gather their cpd from the checkbox and write
3537 :     #the path that we need
3538 :     foreach my $path (@$paths)
3539 :     {
3540 :     my (%input_hash, %output_hash, $path_name);
3541 :     #process the strings to get the information from the parameters
3542 :     my @items = split(";", $path);
3543 :     $path_name = $items[0];
3544 :     #next we have the input compounds ids/stoich/main
3545 :     map { if ($_ =~ /(.*):(.*):(.*)/)
3546 :     { $input_hash{$1}+= $2 } } split ",", $items[1];
3547 :     #the third part has the output compounds ids/stoich/main
3548 :     map { if ($_ =~ /(.*):(.*):(.*)/)
3549 :     { $output_hash{$1} += $2 } } split ",",$items[2];
3550 :     #the fourth part has the strings of reactions visited
3551 :     @tempArray = split("#" , $items[3]);
3552 :    
3553 :     print STDERR "\nInputs:\n" if $debug;
3554 :     print STDERR map{"$_ => $input_hash{$_}" } keys %input_hash, "\n" if $debug;
3555 :     print STDERR "\nOutputs:\n" if $debug;
3556 :     print STDERR map{"$_ => $output_hash{$_}"} keys %output_hash, "\n" if $debug;
3557 :    
3558 :     # divide stoichiometry by greatest common denominator
3559 :     my ($min_stoich, @all_stoichs);
3560 :    
3561 :     map { push @all_stoichs, $input_hash{$_}; $min_stoich = $input_hash{$_} if $input_hash{$_} < $min_stoich || $min_stoich == 0 } keys %input_hash;
3562 :     map { push @all_stoichs, $output_hash{$_}; $min_stoich = $output_hash{$_} if $output_hash{$_} < $min_stoich || $min_stoich == 0 } keys %output_hash;
3563 :    
3564 :     my ($gcd, @gcd_candidates);
3565 :    
3566 :     outer: for ($gcd = $min_stoich; $gcd > 1; $gcd--)
3567 :     {
3568 :     foreach my $stoich (@all_stoichs)
3569 :     {
3570 :     next outer if $stoich % $gcd != 0;
3571 :     }
3572 :    
3573 :     last; # found a gcd
3574 :     }
3575 :    
3576 :     map { $input_hash{$_} /= $gcd } keys %input_hash;
3577 :     map { $output_hash{$_} /= $gcd } keys %output_hash;
3578 :    
3579 :     print STDERR "Making directory $base_dir" if $debug;
3580 :    
3581 :    
3582 :     mkdir "$base_dir/$path_name";
3583 :     &write_fluxanalyzer_files("$base_dir/$path_name",\%input_hash,\%output_hash,
3584 :     \@tempArray,\%all_reactions,
3585 :     \%reactions_to_substrate_arrays,
3586 :     \%reactions_to_product_arrays,
3587 :     {});
3588 :     }
3589 :     }
3590 :    
3591 :     sub show_path_results
3592 :     {
3593 :     my ($data_results,$html,$cgi) = @_;
3594 :    
3595 :     print STDERR "Starting Results Display\n";
3596 :    
3597 :    
3598 :     #Display infinite loop warning if the indicator is on
3599 :     if($data_results->{"infinite"})
3600 :     {
3601 :     push(@$html, "<h3>Warning: Looks like an infinite loop</h3>");
3602 :     }
3603 :    
3604 :     #Delete the infinite loop indicator, so we don't need to have a if statment to check for it
3605 :     delete $data_results->{"infinite"};
3606 :    
3607 :     my $path_counter = 1;
3608 :     my $reactionPath;
3609 :     foreach my $token_id (sort { $data_results->{$a}->[1] <=> $data_results->{$b}->[1] }keys %$data_results)
3610 :     {
3611 :     print STDERR "Token id : $token_id\n";
3612 :    
3613 :     if(!($token_id =~ /^\d/))
3614 :     {
3615 :     next;
3616 :     }
3617 :    
3618 :     my @path = @{$data_results->{$token_id}->[2]};
3619 :     my @compounds = @{$data_results->{$token_id}->[3]};
3620 :    
3621 :     push(@$html, "<pre>Token: $token_id\tInitial Pass: $data_results->{$token_id}->[0]\tDone:$data_results->{$token_id}->[1]\n\tReactions: @path\n\tVisted Compounds: @compounds\n\tPath Inputs\n@{$data_results->{$token_id}->[4]}\n\tPath Outputs\n@{$data_results->{$token_id}->[5]}\n\tBorrowed\n@{$data_results->{$token_id}->[6]}\n</pre>");
3622 :    
3623 :    
3624 :     if ($data_results->{$token_id}->[1] == 1)
3625 :     {
3626 :     my $path_name = "path_".$path_counter++;
3627 :     my @tempArray;
3628 :     foreach my $elem(@{$data_results->{$token_id}->[2]}){
3629 :     push @tempArray, $elem;
3630 :     }
3631 :     $reactionPath = join "#", @tempArray;
3632 :     my $checkbox=$cgi->checkbox(-name=>"$path_name", -label=>'',
3633 :     -value=>"$path_name;@{[ join ',', map { $_ . ':' . $data_results->{$token_id}->[7]->{$_} . ':' . $all_compounds_to_main{$_} } keys %{$data_results->{$token_id}->[7]} ]};@{[ join ',', map { $_ . ':' . $data_results->{$token_id}->[8]->{$_} . ':' . $all_compounds_to_main{$_} } keys %{$data_results->{$token_id}->[8]} ]};$reactionPath");
3634 :     push @$html, $checkbox, "&nbsp;$path_name", $cgi->br;
3635 :     }
3636 :    
3637 :     push @$html, "<hr>";
3638 :    
3639 :    
3640 :     }
3641 :    
3642 :    
3643 :     push @$html, $cgi->hidden(-name=>'reaction_info',
3644 :     -value=>$reactionPath);
3645 :    
3646 :     }
3647 :    
3648 :     sub show_path_results_two
3649 :     {
3650 :     my ($data_results,$html,$cgi) = @_;
3651 :    
3652 :     print STDERR "Starting Results Display\n";
3653 :    
3654 :    
3655 :     #Display infinite loop warning if the indicator is on
3656 :     if($data_results->{"infinite"})
3657 :     {
3658 :     push(@$html, "<h3>Warning: Looks like an infinite loop</h3>");
3659 :     }
3660 :    
3661 :     #Delete the infinite loop indicator, so we don't need to have a if statment to check for it
3662 :     delete $data_results->{"infinite"};
3663 :    
3664 :     my $path_counter = 1;
3665 :     my $reactionPath;
3666 :     foreach my $token_id (sort { $data_results->{$a}->[1] <=> $data_results->{$b}->[1] }keys %$data_results)
3667 :     {
3668 :     print STDERR "Token id : $token_id\n";
3669 :    
3670 :     if(!($token_id =~ /^\d/))
3671 :     {
3672 :     next;
3673 :     }
3674 :    
3675 :     my @path = @{$data_results->{$token_id}->[2]};
3676 :     my @compounds = @{$data_results->{$token_id}->[3]};
3677 :    
3678 :     push(@$html, "<pre>Token: $token_id\tInitial Pass: $data_results->{$token_id}->[0]\tDone:$data_results->{$token_id}->[1]\n\tReactions: @path\n\tVisted Compounds: @compounds\n\tPath Inputs\n@{$data_results->{$token_id}->[4]}\n\tPath Outputs\n@{$data_results->{$token_id}->[5]}\n\tBorrowed\n@{$data_results->{$token_id}->[6]}\n</pre>");
3679 :    
3680 :    
3681 :     if ($data_results->{$token_id}->[1] == 1)
3682 :     {
3683 :     my $path_name = "path_".$path_counter++;
3684 :     my @tempArray;
3685 :     foreach my $elem(@{$data_results->{$token_id}->[2]}){
3686 :     push @tempArray, $elem;
3687 :     }
3688 :     $reactionPath = join "#", @tempArray;
3689 :     my $checkbox=$cgi->checkbox(-name=>"checked", -label=>'',
3690 :     -value=>"$path_name;@{[ join ',', map { $_ . ':' . $data_results->{$token_id}->[7]->{$_} . ':' . $all_compounds_to_main{$_} } keys %{$data_results->{$token_id}->[7]} ]};@{[ join ',', map { $_ . ':' . $data_results->{$token_id}->[8]->{$_} . ':' . $all_compounds_to_main{$_} } keys %{$data_results->{$token_id}->[8]} ]};$reactionPath");
3691 :     push @$html, $checkbox, "&nbsp;$path_name", $cgi->br;
3692 :     }
3693 :    
3694 :     push @$html, "<hr>";
3695 :     }
3696 :    
3697 :     }
3698 :    
3699 :    
3700 :     sub set_loop_max
3701 :     {
3702 :     my ($number) = @_;
3703 :     $loop_max = $number;
3704 :     }
3705 :    
3706 :     sub set_loop_max_assembly
3707 :     {
3708 :     my ($number) = @_;
3709 :     $loop_max_assembly = $number;
3710 :     }
3711 :    
3712 : dejongh 1.8 sub analyze_scenario_connections
3713 :     {
3714 :     my ($genome_id) = @_;
3715 : dejongh 1.23 my $scenario_dir = get_scenario_directory($genome_id);
3716 :     my $pathinfo_dir = $scenario_dir . "/PathInfo";
3717 : dejongh 1.26 my @paths = `find $pathinfo_dir -type d -name "path_*"`;
3718 : dejongh 1.8 my %inputs;
3719 :     my %outputs;
3720 : olson 1.1
3721 : dejongh 1.12 foreach my $dir (@paths)
3722 : dejongh 1.8 {
3723 : dejongh 1.12 chomp $dir;
3724 : dejongh 1.8 my ($cat, $subsys, $scenario);
3725 : olson 1.1
3726 : dejongh 1.23 if ($dir =~ (/\/PathInfo\/(.*)\/(.*)\/(.*)\//)){
3727 : dejongh 1.8 $cat = $1;
3728 :     $subsys = $2;
3729 :     $scenario = $3;
3730 :     }
3731 :     else
3732 :     {
3733 :     next;
3734 :     }
3735 : olson 1.1
3736 : dejongh 1.8 open(M_INPUTS,$dir."/inputs_main") or die("Failed to open $dir/inputs_main");
3737 :    
3738 :     while (<M_INPUTS>)
3739 :     {
3740 :     chomp;
3741 :     my ($cpd, $main) = split "\t" , $_;
3742 :     my $info = join "\t", $subsys, $scenario;
3743 :     $inputs{$cpd}->{$info} = 1 if $main eq "1";
3744 :     }
3745 :     close M_INPUTS;
3746 :    
3747 :     open(M_OUTPUTS, $dir."/outputs_main");
3748 :    
3749 :     while (<M_OUTPUTS>)
3750 :     {
3751 :     chomp;
3752 :     my ($cpd, $main) = split "\t" , $_;
3753 :     my $info = join "\t", $subsys, $scenario;
3754 :     $outputs{$cpd}->{$info} = 1 if $main eq "1";
3755 :     }
3756 :     close M_OUTPUTS;
3757 : dejongh 1.12 }
3758 : dejongh 1.8
3759 : dejongh 1.23 my $analysis_dir = $scenario_dir . "/Analysis";
3760 : dejongh 1.12 mkdir $analysis_dir;
3761 : dejongh 1.8
3762 : dejongh 1.12 open (IN_CONN, ">$analysis_dir/inputs_to_scenarios");
3763 :     foreach my $cpd (sort keys %inputs)
3764 :     {
3765 :     map { print IN_CONN "$cpd\t$_\n"; } keys %{$inputs{$cpd}};
3766 :     }
3767 :     close IN_CONN;
3768 : dejongh 1.8
3769 : dejongh 1.12 open (OUT_CONN, ">$analysis_dir/outputs_to_scenarios");
3770 :     foreach my $cpd (sort keys %outputs)
3771 :     {
3772 :     map { print OUT_CONN "$cpd\t$_\n"; } keys %{$outputs{$cpd}};
3773 : dejongh 1.8 }
3774 : dejongh 1.12 close OUT_CONN;
3775 : olson 1.1 }
3776 :    
3777 : formsma 1.14 sub predict_pegs_used
3778 :     {
3779 :     my ($genome_id) = @_;
3780 :     my @scenarios = @{Scenario->get_genome_scenarios("All",1)};
3781 :     unless(scalar(@scenarios))
3782 :     {
3783 :     return undef;
3784 :     }
3785 :     my %reaction_to_pegs;
3786 :    
3787 :     my @ss_names;
3788 : dejongh 1.24 &load_superset_file;
3789 : formsma 1.14 foreach (keys %superset_to_ss)
3790 :     {
3791 :     foreach my $subsys (@{$superset_to_ss{$_}})
3792 :     {
3793 :     push @ss_names, $subsys;
3794 :     }
3795 :    
3796 :     }
3797 :     foreach my $subsystem_name(@ss_names)
3798 :     {
3799 :     my $subsystem = $fig->get_subsystem($subsystem_name);
3800 :     next if(!defined $subsystem);
3801 : dejongh 1.24 my %reactions_for_ss = get_reactions_for_genome_in_subsystem($genome_id,$subsystem_name);
3802 : parrello 1.32 next if(! keys %reactions_for_ss);
3803 : dejongh 1.17 foreach my $reaction (keys %reactions_for_ss)
3804 : formsma 1.14 {
3805 :     if(defined $reaction_to_pegs{$reaction})
3806 :     {
3807 : dejongh 1.17 push @{$reaction_to_pegs{$reaction}} , @{$reactions_for_ss{$reaction}};
3808 : formsma 1.14 }
3809 :     else
3810 :     {
3811 : dejongh 1.17 $reaction_to_pegs{$reaction} = $reactions_for_ss{$reaction};
3812 : formsma 1.14 }
3813 :     }
3814 :     }
3815 :    
3816 :     my %peg_to_scenario;
3817 :     foreach my $scenario (@scenarios)
3818 :     {
3819 :     my @scenario_reactions = @{$scenario->get_reaction_ids};
3820 :     my $path_valid = 1;
3821 :     my %pegs;
3822 :     foreach my $reaction (@scenario_reactions)
3823 :     {
3824 :     if(!defined $reaction_to_pegs{$reaction})
3825 :     {
3826 :     $path_valid = 0;
3827 :     last;
3828 :     }
3829 :     else
3830 :     {
3831 :     map {$pegs{$_} = 1} @{$reaction_to_pegs{$reaction}};
3832 :     }
3833 :     }
3834 :     if($path_valid)
3835 :     {
3836 :     foreach my $peg (keys %pegs)
3837 :     {
3838 :     if(!defined $peg_to_scenario{$peg})
3839 :     {
3840 :     $peg_to_scenario{$peg} = [$scenario->get_id()];
3841 :     }
3842 :     else
3843 :     {
3844 :     push @{$peg_to_scenario{$peg}} , $scenario->get_id();
3845 :     }
3846 :     }
3847 :     }
3848 :     }
3849 :    
3850 :     return \%peg_to_scenario;
3851 :     }
3852 : dejongh 1.19
3853 : dejongh 1.21 sub analyze_network_for_genome
3854 :     {
3855 :     my ($genome_id, $transports, $biomass) = @_;
3856 :    
3857 :     my ($node_type, $node_connections);
3858 :     my $debug = 0;
3859 :     my $unreachable = 0;
3860 :     chomp $genome_id;
3861 :    
3862 :     my @transports = split ("_", $transports);
3863 :     my @biomass = split ("_", $biomass);
3864 :    
3865 :     my %compound_id_to_name;
3866 :     open(IN,"$FIG_Config::global/Models/All/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt");
3867 :     while(<IN>)
3868 :     {
3869 :     chomp($_);
3870 :     my($id,$name) = split("\t");
3871 :     $compound_id_to_name{$id} = $name;
3872 :     }
3873 :     close(IN);
3874 :    
3875 :     Scenario::set_fig($fig,$genome_id);
3876 :     print STDERR "Running check_model_hope on genome $genome_id\n" if($debug);
3877 :     print STDERR "\tStage 1 - Loading Scenarios, transports and biomass\n" if($debug);
3878 :     my @scenarios = @{Scenario->get_genome_scenarios($genome_id,0)};
3879 :    
3880 :     my %scenarios_used;
3881 :     my %scenarios_id;
3882 :    
3883 :     foreach (@scenarios)
3884 :     {
3885 :     $scenarios_used{$_->get_scenario_name()} = 0;
3886 :     $scenarios_id{$_->get_id()} = 0;
3887 :     }
3888 :    
3889 :     my %biomass_lib_cpds;
3890 :     foreach my $biomass_cpd (@biomass)
3891 :     {
3892 :     $biomass_lib_cpds{$biomass_cpd} = 1;
3893 :     }
3894 :    
3895 :     my (%transportable_cpds, %transportable_in_biomass);
3896 :     foreach my $transport_cpd (@transports)
3897 :     {
3898 :     if (defined $biomass_lib_cpds{$transport_cpd})
3899 :     {
3900 :     $transportable_in_biomass{$transport_cpd} = 1;
3901 :     }
3902 :     else
3903 :     {
3904 :     $transportable_cpds{$transport_cpd} = 1;
3905 :     }
3906 :     }
3907 :    
3908 :     print STDERR "\tStage 1 - Complete\n" if($debug);
3909 :     print STDERR "\tStage 2 - Run Analysis Tools\n" if($debug);
3910 :    
3911 :     # Total list of what compounds we reach.
3912 :     my %list_generated;
3913 :    
3914 :     # Here we are running all the scenarios we have inputs for
3915 :     # to analyze what biomass we can reach. We are also identifying
3916 :     # any Scenarios that are not being utilized.
3917 :    
3918 :     my %current_inputs;
3919 :     foreach (keys %transportable_cpds)
3920 :     {
3921 :     # keep track of the paths that lead to each compound
3922 :     $current_inputs{$_}->{"[]transport of ".$_."/".$compound_id_to_name{$_}."[$_]"} = 1;
3923 :     }
3924 :    
3925 :     while(1)
3926 :     {
3927 :     my ($temp_in,$break) = generate_new_inputs(\%current_inputs);
3928 :    
3929 :     foreach my $cpd (keys %$temp_in)
3930 :     {
3931 :     map { $current_inputs{$cpd}->{$_} = 1 } keys %{$temp_in->{$cpd}};
3932 :     }
3933 :    
3934 :     map { $list_generated{$_} = 1 } keys %$temp_in;
3935 :    
3936 :     print STDERR "Rechecking with inputs\n" if($debug);
3937 :     if($break)
3938 :     {
3939 :     last;
3940 :     }
3941 :     }
3942 :     print STDERR "Finished checking organism scenarios\n" if($debug);
3943 :    
3944 :     my %biomass_reached;
3945 :    
3946 :     foreach my $found (keys %list_generated)
3947 :     {
3948 :     if(defined $biomass_lib_cpds{$found})
3949 :     {
3950 :     $biomass_reached{$found} = $current_inputs{$found};
3951 :     }
3952 :     }
3953 :    
3954 :     foreach (keys %transportable_cpds)
3955 :     {
3956 :     push (@{$node_type->{'transport'}}, $_."/".$compound_id_to_name{$_});
3957 :     }
3958 :    
3959 :     print STDERR "\tStage 2 - Complete!\n" if($debug);
3960 :     print STDERR "\tStage 3 - Loading All Scenarios\n" if($debug);
3961 :    
3962 :     @scenarios = @{Scenario->get_genome_scenarios("All",0)};
3963 :    
3964 :     my %old_scenarios_used = %scenarios_used;
3965 :     my %old_scenarios_id = %scenarios_id;
3966 :    
3967 :     my %old_biomass_reached;
3968 :     foreach my $cid (keys %biomass_reached)
3969 :     {
3970 :     my %inner_hash_copy = %{$biomass_reached{$cid}};
3971 :     $old_biomass_reached{$cid} = \%inner_hash_copy
3972 :     }
3973 :    
3974 :     my %old_current_inputs;
3975 :     foreach my $cid (keys %current_inputs)
3976 :     {
3977 :     my %inner_hash_copy = %{$current_inputs{$cid}};
3978 :     $old_current_inputs{$cid} = \%inner_hash_copy;
3979 :     }
3980 :    
3981 :    
3982 : parrello 1.32 %list_generated = ();
3983 : dejongh 1.21
3984 :     foreach (@scenarios)
3985 :     {
3986 :     $scenarios_used{$_->get_scenario_name()} = 0 if ! exists $scenarios_used{$_->get_scenario_name()};
3987 :     $scenarios_id{$_->get_id()} = 0 if ! exists $scenarios_id{$_->get_id()};
3988 :     }
3989 :    
3990 :     while(1)
3991 :     {
3992 :     my ($temp_in,$break) = generate_new_inputs(\%current_inputs);
3993 :    
3994 :     foreach my $cpd (keys %$temp_in)
3995 :     {
3996 :     map { $current_inputs{$cpd}->{$_} = 1 } keys %{$temp_in->{$cpd}};
3997 :     }
3998 :    
3999 :     map { $list_generated{$_} = 1 } keys %$temp_in;
4000 :    
4001 :     print STDERR "Rechecking with inputs\n" if($debug);
4002 :     if($break)
4003 :     {
4004 :     last;
4005 :     }
4006 :     }
4007 :     print STDERR "Finished checking all scenarios\n" if($debug);
4008 :    
4009 : parrello 1.32 %biomass_reached = ();
4010 : dejongh 1.21
4011 :     foreach my $found (keys %list_generated)
4012 :     {
4013 :     if(defined $biomass_lib_cpds{$found})
4014 :     {
4015 :     $biomass_reached{$found} = $current_inputs{$found};
4016 :     }
4017 :     }
4018 :    
4019 :     my (%cpds_to_print);
4020 :     my %org_cpds_to_check;
4021 :     map { $org_cpds_to_check{$_} = 1 } keys %old_biomass_reached;
4022 :     foreach my $cpd (sort keys %biomass_reached)
4023 :     {
4024 :     next if exists $old_biomass_reached{$cpd};
4025 :    
4026 :     $cpds_to_print{$cpd} = 1;
4027 :     my $done = 0;
4028 :    
4029 :     my @cpd_list = ($cpd);
4030 :     while (! $done)
4031 :     {
4032 :     my @new_substrates = collect_substrates(\%current_inputs, @cpd_list);
4033 :     @cpd_list = ();
4034 :     foreach my $substrate (@new_substrates)
4035 :     {
4036 :     if (exists $old_current_inputs{$substrate})
4037 :     {
4038 :     $org_cpds_to_check{$substrate} = 1;
4039 :     next;
4040 :     }
4041 :     push (@cpd_list, $substrate) if ! exists $cpds_to_print{$substrate};
4042 :     $cpds_to_print{$substrate} = 1;
4043 :     }
4044 :     $done = 1 if scalar @new_substrates == 0;
4045 :     }
4046 :     }
4047 :    
4048 :     foreach my $cpd (keys %cpds_to_print)
4049 :     {
4050 :     foreach my $path (keys %{$current_inputs{$cpd}})
4051 :     {
4052 :     my ($substrates,$scenario,$products) = $path =~ /\[(.*)\](.*)\[(.*)\]/;
4053 :    
4054 :     foreach my $substrate (split ",", $substrates)
4055 :     {
4056 :     foreach my $product (split ",", $products)
4057 :     {
4058 :     $node_connections->{$substrate."/".$compound_id_to_name{$substrate}}->$scenario = 1;
4059 :     $node_connections->{$scenario}->$product."/".$compound_id_to_name{$product} = 1;
4060 :     }
4061 :     }
4062 :     }
4063 :     }
4064 :    
4065 :     %cpds_to_print = ();
4066 :    
4067 :     foreach my $cpd (sort keys %org_cpds_to_check)
4068 :     {
4069 :     $cpds_to_print{$cpd} = 1;
4070 :     my $done = 0;
4071 :    
4072 :     my @cpd_list = ($cpd);
4073 :     while (! $done)
4074 :     {
4075 :     my @new_substrates = collect_substrates(\%old_current_inputs, @cpd_list);
4076 :     @cpd_list = ();
4077 :     foreach my $substrate (@new_substrates)
4078 :     {
4079 :     push (@cpd_list, $substrate) if ! exists $cpds_to_print{$substrate};
4080 :     $cpds_to_print{$substrate} = 1;
4081 :     }
4082 :     $done = 1 if scalar @new_substrates == 0;
4083 :     }
4084 :     }
4085 :    
4086 :     foreach my $cpd (keys %cpds_to_print)
4087 :     {
4088 :     foreach my $path (keys %{$old_current_inputs{$cpd}})
4089 :     {
4090 :     my ($substrates,$scenario,$products) = $path =~ /\[(.*)\](.*)\[(.*)\]/;
4091 :    
4092 :     foreach my $substrate (split ",", $substrates)
4093 :     {
4094 :     foreach my $product (split ",", $products)
4095 :     {
4096 :     $node_connections->{$substrate."/".$compound_id_to_name{$substrate}}->$scenario = 1;
4097 :     $node_connections->{$scenario}->$product."/".$compound_id_to_name{$product} = 1;
4098 :     }
4099 :     }
4100 :     }
4101 :     }
4102 :    
4103 :     foreach (sort keys %biomass_lib_cpds)
4104 :     {
4105 :     if(! exists $biomass_reached{$_} && ! exists $old_biomass_reached{$_})
4106 :     {
4107 :     print STDERR "Biomass $_/".$compound_id_to_name{$_}." not reachable using current transports" if ($debug);
4108 :     $unreachable++;
4109 :     if (exists $transportable_in_biomass{$_})
4110 :     {
4111 :     print STDERR " (transportable)" if ($debug);
4112 :     push (@{$node_type->{'biomass not reached (transportable)'}}, $_."/".$compound_id_to_name{$_});
4113 :     }
4114 :     else
4115 :     {
4116 :     push (@{$node_type->{'biomass not reached'}}, $_."/".$compound_id_to_name{$_});
4117 :     }
4118 :     print STDERR "\n" if ($debug);
4119 :     }
4120 :     elsif ($old_biomass_reached{$_})
4121 :     {
4122 :     push (@{$node_type->{'biomass reached'}}, $_."/".$compound_id_to_name{$_});
4123 :     }
4124 :     else
4125 :     {
4126 :     # do the paths all include a non-organism scenario as the last step?
4127 :     my $last_not_org = 1;
4128 :    
4129 :     foreach my $path_set (keys %{$biomass_reached{$_}})
4130 :     {
4131 :     my @paths = split ";", $path_set;
4132 :     my $last_path = $paths[-1];
4133 :     my ($substrates,$scenario,$products) = $last_path =~ /\[(.*)\](.*)\[(.*)\]/;
4134 :    
4135 :     if (exists $old_scenarios_used{$scenario})
4136 :     {
4137 :     $last_not_org = 0;
4138 :     last;
4139 :     }
4140 :     }
4141 :    
4142 :     if ($last_not_org == 1)
4143 :     {
4144 :     push (@{$node_type->{'biomass reached all last_not_org'}}, $_."/".$compound_id_to_name{$_});
4145 :     }
4146 :     else
4147 :     {
4148 :     push (@{$node_type->{'biomass reached all'}}, $_."/".$compound_id_to_name{$_});
4149 :     }
4150 :     }
4151 :     }
4152 :    
4153 :     foreach my $scen_name (sort keys %scenarios_used)
4154 :     {
4155 :     if($scenarios_used{$scen_name} == 1)
4156 :     {
4157 :     if (exists ($old_scenarios_used{$scen_name}))
4158 :     {
4159 :     push (@{$node_type->{'scenario present'}}, $scen_name);
4160 :     }
4161 :     else
4162 :     {
4163 :     push (@{$node_type->{'scenario present all'}}, $scen_name);
4164 :     }
4165 :     }
4166 :     }
4167 :    
4168 :     foreach my $cpd (keys %current_inputs)
4169 :     {
4170 :     # compounds present with All scenarios that aren't biomass
4171 :     if (! exists($old_current_inputs{$cpd}) && ! exists ($biomass_reached{$cpd}))
4172 :     {
4173 :     push (@{$node_type->{'compound present all'}}, $cpd."/".$compound_id_to_name{$cpd});
4174 :     }
4175 :     }
4176 :    
4177 :     print STDERR "\tStage 3 - Complete!\n" if($debug);
4178 :    
4179 :     return ($node_type, $node_connections);
4180 :    
4181 :     sub generate_new_inputs
4182 :     {
4183 :     my ($inputs) = @_;
4184 :     my %new_inputs;
4185 :     my $break = 1;
4186 :    
4187 :     my %used_scenario_paths_this_time;
4188 :    
4189 :     foreach my $scenario (@scenarios)
4190 :     {
4191 :     if($scenarios_id{$scenario->get_id()} == 1) # || $scenarios_used{$scenario->get_scenario_name()} != 0)
4192 :     {
4193 :     next;
4194 :     }
4195 :    
4196 :     my $sids = join ",", sort @{$scenario->get_substrates};
4197 :     my $pids = join ",", sort @{$scenario->get_products};
4198 :     my $cid_scenario_name = "[".$sids."]".$scenario->get_scenario_name()."[".$pids."]";
4199 :    
4200 :     my @matched;
4201 :    
4202 :     foreach my $substrate (@{$scenario->get_substrates()})
4203 :     {
4204 :     if(exists $inputs->{$substrate})
4205 :     {
4206 :     push @matched, $substrate;
4207 :     }
4208 :     }
4209 :    
4210 :     if (scalar @matched != scalar(@{$scenario->get_substrates()}))
4211 :     {
4212 :     next;
4213 :     }
4214 :     else
4215 :     {
4216 :     print STDERR "Scenario ".$cid_scenario_name." ran.\n" if($debug);
4217 :     $break = 0;
4218 :     $scenarios_id{$scenario->get_id()} = 1;
4219 :     $scenarios_used{$scenario->get_scenario_name()} = 1;
4220 :    
4221 :     # local determination that we've used an equivalent path already in this
4222 :     # time through the subroutine, so no need to process further
4223 :     next if (exists $used_scenario_paths_this_time{$cid_scenario_name});
4224 :     $used_scenario_paths_this_time{$cid_scenario_name} = 1;
4225 :    
4226 :     foreach my $product (@{$scenario->get_products()})
4227 :     {
4228 :     $new_inputs{$product}->{$cid_scenario_name} = 1;
4229 :     }
4230 :     }
4231 :     }
4232 :     return (\%new_inputs,$break);
4233 :     }
4234 :    
4235 :     sub collect_substrates
4236 :     {
4237 :     my $input_hash = shift @_;
4238 :     my @cpds = @_;
4239 :     my %substrate_hash;
4240 :    
4241 :     foreach my $cpd (@cpds)
4242 :     {
4243 :     foreach my $path (keys %{$input_hash->{$cpd}})
4244 :     {
4245 :     my ($substrates,$scenario,$products) = $path =~ /\[(.*)\](.*)\[(.*)\]/;
4246 :    
4247 :     foreach my $substrate (split ",", $substrates)
4248 :     {
4249 :     $substrate_hash{$substrate} = 1;
4250 :     }
4251 :     }
4252 :     }
4253 :     return keys %substrate_hash;
4254 :    
4255 :    
4256 :     }
4257 :     }
4258 :    
4259 : dejongh 1.24 sub get_reactions_for_genome_in_subsystem
4260 : dejongh 1.19 {
4261 :     my ($genome_id, $ss_name) = @_;
4262 : dejongh 1.24 my $subsystem = $fig->get_subsystem($ss_name);
4263 :     my %genome_reactions = ();
4264 :     my %hope_reactions = $subsystem->get_hope_reactions;
4265 :    
4266 : dejongh 1.31 foreach my $role (keys %hope_reactions)
4267 : dejongh 1.24 {
4268 : dejongh 1.31 my @peg_list = $fig->seqs_with_role($role, "", $genome_id);
4269 :    
4270 :     if (scalar @peg_list > 0)
4271 : dejongh 1.24 {
4272 : dejongh 1.31 foreach my $reaction (@{$hope_reactions{$role}})
4273 : dejongh 1.24 {
4274 : dejongh 1.31 push @{$genome_reactions{$reaction}}, @peg_list;
4275 : dejongh 1.24 }
4276 : dejongh 1.31 }
4277 : dejongh 1.24 }
4278 : dejongh 1.19
4279 : dejongh 1.20 # assume that the hope_supersets.txt has been loaded
4280 :     my $superset = $ss_to_superset{$ss_name};
4281 : dejongh 1.23 my $pathdir_org = get_scenario_directory($genome_id);
4282 :     my $file_path = "$pathdir_org/Curation/PathInfo/$superset/$ss_name/included_reactions_for_org";
4283 : dejongh 1.19 if (open (INCLUDEDRXNS, $file_path))
4284 :     {
4285 :     my @included_reactions = (<INCLUDEDRXNS>);
4286 :     close INCLUDEDRXNS;
4287 :     chomp (@included_reactions);
4288 : dejongh 1.24 map { $genome_reactions{$_} = 1 } @included_reactions;
4289 : dejongh 1.19 }
4290 : dejongh 1.31
4291 : dejongh 1.24 return %genome_reactions;
4292 : dejongh 1.19 }
4293 : dejongh 1.20
4294 : dejongh 1.23 sub get_scenario_directory
4295 : dejongh 1.20 {
4296 :     my ($genome) = @_;
4297 :     # want to match "All" and "All.bk", etc.
4298 : dejongh 1.23 return $genome =~ "^All" ? $FIG_Config::global."/Models/$genome/Scenarios" : $fig->scenario_directory($genome);
4299 : dejongh 1.20
4300 :     }
4301 :    
4302 : dejongh 1.28 # return three values:
4303 :     # 1) reference to a hash of invalid roles, i.e., roles no longer in the subsystem.
4304 :     # 2) reference to a hash of invalid reactions, i.e., reactions no longer in KEGG.
4305 :     # 3) flag set to main classification value if the subsystem classification has changed
4306 :     # to Experimental or none.
4307 :     sub check_hope_info
4308 :     {
4309 :     my ($ss) = @_;
4310 :     my (%invalid_roles, %invalid_reactions, $invalid_classification);
4311 :     my $sub = $fig->get_subsystem($ss);
4312 :    
4313 :     if (! defined $sub)
4314 :     {
4315 :     return ( {}, {}, "Subsystem no longer exists" );
4316 :     }
4317 :    
4318 :     my $classification = $sub->get_classification->[0];
4319 :     $invalid_classification = $classification if ($classification eq '' || $classification =~ /Experimental/);
4320 :     my %hope_reactions = $sub->get_hope_reactions;
4321 :     my @roles = $sub->get_roles;
4322 :     my %roles;
4323 :     my %reactions;
4324 :    
4325 :     map { $roles{$_} = 1 } @roles;
4326 :    
4327 :     foreach my $role (keys %hope_reactions)
4328 :     {
4329 :     if (! exists($roles{$role}))
4330 :     {
4331 :     $invalid_roles{$role} = $hope_reactions{$role};
4332 :     }
4333 :    
4334 :     map {push @{$reactions{$_}}, $role} @{$hope_reactions{$role}};
4335 :     }
4336 :    
4337 :     foreach my $name ($sub->get_hope_scenario_names)
4338 :     {
4339 :     my @additional_reactions = $sub->get_hope_additional_reactions($name);
4340 :     map { push @{$reactions{$_}}, "ADDITIONAL_REACTIONS" } @additional_reactions;
4341 :     }
4342 :    
4343 :    
4344 :     foreach my $reaction (keys %reactions)
4345 :     {
4346 : dejongh 1.29 unless ($reaction eq "not_in_KEGG" || $fig->valid_reaction_id($reaction))
4347 : dejongh 1.28 {
4348 :     $invalid_reactions{$reaction} = $reactions{$reaction};
4349 :     }
4350 :     }
4351 :    
4352 :     return (\%invalid_roles, \%invalid_reactions, $invalid_classification);
4353 :    
4354 :     }
4355 :    
4356 : olson 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3