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

Annotation of /FigKernelPackages/FIGMODELmodel.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : chenry 1.18 use strict;
2 :     package FIGMODELmodel;
3 : chenry 1.12
4 : chenry 1.18 =head1 FIGMODELmodel object
5 :     =head2 Introduction
6 :     Module for manipulating model objects.
7 :     =head2 Core Object Methods
8 :    
9 :     =head3 new
10 :     Definition:
11 :     FIGMODELmodel = FIGMODELmodel->new();
12 :     Description:
13 :     This is the constructor for the FIGMODELmodel object.
14 :     =cut
15 :     sub new {
16 : chenry 1.21 my ($class,$figmodel,$id,$metagenome) = @_;
17 : chenry 1.18 #Error checking first
18 :     if (!defined($figmodel)) {
19 :     print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n";
20 :     return undef;
21 :     }
22 :     my $self = {_figmodel => $figmodel};
23 :     bless $self;
24 :     if (!defined($id)) {
25 :     $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,undef):id must be defined to create a model object");
26 :     return undef;
27 :     }
28 :     #Checking that the id exists
29 : chenry 1.21 if (!defined($metagenome) || $metagenome != 1) {
30 :     my $modelHandle = $self->figmodel()->database()->get_object_manager("model");
31 :     if (defined($modelHandle)) {
32 :     #If the id is a number, we get the model row by index
33 :     if ($id =~ m/^\d+$/) {
34 :     my $objects = $modelHandle->get_objects();
35 :     if (defined($objects->[$id])) {
36 :     $self->{_data} = $objects->[$id];
37 :     $self->figmodel()->{_models}->{$id} = $self;
38 :     }
39 :     } else {
40 :     my $objects = $modelHandle->get_objects({id => $id});
41 :     if (defined($objects->[0])) {
42 :     $self->{_data} = $objects->[0];
43 :     } elsif (!defined($objects->[0]) && $id =~ m/(.+)(V[^V]+)/) {
44 :     $objects = $modelHandle->get_objects({id => $1});
45 :     if (defined($objects->[0])) {
46 :     $self->{_selectedversion} = $2;
47 :     $self->{_data} = $objects->[0];
48 :     }
49 :     }
50 :     }
51 :     }
52 :     if (defined($self->{_data})) {
53 :     $self->{_modeltype} = "genome";
54 :     }
55 : chenry 1.18 }
56 : chenry 1.21 if (!defined($self->{_data})) {
57 :     my $modelHandle = $self->figmodel()->database()->get_object_manager("mgmodel");
58 :     if (defined($modelHandle)) {
59 :     #If the id is a number, we get the model row by index
60 :     if ($id =~ m/^\d+$/) {
61 :     my $objects = $modelHandle->get_objects();
62 :     if (defined($objects->[$id])) {
63 :     $self->{_data} = $objects->[$id];
64 :     $self->figmodel()->{_models}->{"MG".$id} = $self;
65 :     }
66 :     } else {
67 :     my $objects = $modelHandle->get_objects({id => $id});
68 :     if (defined($objects->[0])) {
69 :     $self->{_data} = $objects->[0];
70 :     } elsif (!defined($objects->[0]) && $id =~ m/(.+)(V[^V]+)/) {
71 :     $objects = $modelHandle->get_objects({id => $1});
72 :     if (defined($objects->[0])) {
73 :     $self->{_selectedversion} = $2;
74 :     $self->{_data} = $objects->[0];
75 :     }
76 :     }
77 : chenry 1.18 }
78 :     }
79 : chenry 1.21 if (defined($self->{_data})) {
80 :     $self->{_modeltype} = "metagenome";
81 :     }
82 : chenry 1.18 }
83 :     if (!defined($self->{_data})) {
84 : chenry 1.21 $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!");
85 : chenry 1.18 return undef;
86 :     }
87 :     $self->figmodel()->{_models}->{$self->id()} = $self;
88 :     return $self;
89 :     }
90 :    
91 :     =head3 config
92 :     Definition:
93 :     ref::key value = FIGMODELmodel->config(string::key);
94 :     Description:
95 :     Trying to avoid using calls that assume configuration data is stored in a particular manner.
96 :     Call this function to get file paths etc.
97 :     =cut
98 :     sub config {
99 :     my ($self,$key) = @_;
100 :     return $self->figmodel()->config($key);
101 :     }
102 :    
103 : chenry 1.27 =head3 error_message
104 :     Definition:
105 :     string:message text = FIGMODELmodel->error_message(string::message);
106 :     Description:
107 :     =cut
108 :     sub error_message {
109 :     my ($self,$message) = @_;
110 :     return $self->figmodel()->error_message("FIGMODELmodel:".$self->id().":".$message);
111 :     }
112 :    
113 : chenry 1.18 =head3 fail
114 :     Definition:
115 :     -1 = FIGMODELmodel->fail();
116 :     Description:
117 :     Standard return for failed functions.
118 :     =cut
119 :     sub fail {
120 :     my ($self) = @_;
121 :     return $self->figmodel()->fail();
122 :     }
123 :    
124 :     =head3 success
125 :     Definition:
126 :     1 = FIGMODELmodel->success();
127 :     Description:
128 :     Standard return for successful functions.
129 :     =cut
130 :     sub success {
131 :     my ($self) = @_;
132 :     return $self->figmodel()->success();
133 :     }
134 :    
135 :     =head3 figmodel
136 :     Definition:
137 :     FIGMODEL = FIGMODELmodel->figmodel();
138 :     Description:
139 :     Returns a FIGMODEL object
140 :     =cut
141 :     sub figmodel {
142 :     my ($self) = @_;
143 :     return $self->{"_figmodel"};
144 :     }
145 :    
146 :     =head3 fig
147 :     Definition:
148 :     FIGMODEL = FIGMODELmodel->fig();
149 :     Description:
150 :     Returns a FIG object
151 :     =cut
152 :     sub fig {
153 :     my ($self) = @_;
154 :     if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {
155 :     if ($self->source() =~ /^RAST/) {
156 :     $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());
157 :     } else {
158 :     $self->{"_fig"} = $self->figmodel()->fig();
159 :     }
160 :     }
161 :     return $self->{"_fig"};
162 :     }
163 :    
164 :     =head3 aquireModelLock
165 :    
166 :     Definition:
167 :    
168 :     FIGMODELmodel->aquireModelLock();
169 :    
170 :     Description:
171 :    
172 :     Locks the database for alterations relating to the current model object
173 :    
174 :     =cut
175 :     sub aquireModelLock {
176 :     my ($self) = @_;
177 :     $self->figmodel()->database()->genericLock($self->id());
178 :     }
179 :    
180 :     =head3 releaseModelLock
181 :    
182 :     Definition:
183 :    
184 :     FIGMODELmodel->releaseModelLock();
185 :    
186 :     Description:
187 :    
188 :     Unlocks the database for alterations relating to the current model object
189 :    
190 :     =cut
191 :     sub releaseModelLock {
192 :     my ($self) = @_;
193 :     $self->figmodel()->database()->genericUnlock($self->id());
194 :     }
195 :    
196 :     =head3 mgdata
197 :     Definition:
198 :     FIGMODEL = FIGMODELmodel->mgdata();
199 :     Description:
200 :     Returns a mgrast database object
201 :     =cut
202 :     sub mgdata {
203 :     my ($self) = @_;
204 : chenry 1.22 if (!defined($self->{_mgdata})) {
205 :     my $mgrastdbhandle = $self->figmodel()->database()->get_object_manager("mgjob");
206 :     my $objects = $mgrastdbhandle->get_objects( { 'genome_id' => $self->genome() } );
207 :     $self->{_mgdata} = $objects->[0];
208 : chenry 1.18 }
209 :     return $self->{_mgdata};
210 :     }
211 :    
212 :     =head3 id
213 :     Definition:
214 :     string = FIGMODELmodel->id();
215 :     Description:
216 :     Returns model id
217 :     =cut
218 :     sub id {
219 :     my ($self) = @_;
220 :     return $self->{_data}->id();
221 :     }
222 :    
223 : chenry 1.21 =head3 get_model_type
224 :     Definition:
225 :     string = FIGMODELmodel->get_model_type();
226 :     Description:
227 :     Returns the type of the model
228 :     =cut
229 :     sub get_model_type {
230 :     my ($self) = @_;
231 :     return $self->{_modeltype};
232 :     }
233 :    
234 : chenry 1.18 =head3 owner
235 :     Definition:
236 :     string = FIGMODELmodel->owner();
237 :     Description:
238 :     Returns the username for the model owner
239 :     =cut
240 :     sub owner {
241 :     my ($self) = @_;
242 :     return $self->{_data}->owner();
243 :     }
244 :    
245 : chenry 1.27 =head3 users
246 :     Definition:
247 :     string = FIGMODELmodel->users();
248 :     Description:
249 :     =cut
250 :     sub users {
251 :     my ($self) = @_;
252 :     return $self->{_data}->users();
253 :     }
254 :    
255 : chenry 1.18 =head3 name
256 :     Definition:
257 :     string = FIGMODELmodel->name();
258 :     Description:
259 :     Returns the name of the organism or metagenome sample being modeled
260 :     =cut
261 :     sub name {
262 :     my ($self) = @_;
263 :    
264 :     if (!defined($self->{_name})) {
265 :     my $source = $self->source();
266 :     if ($source =~ /^MGRAST/) {
267 :     $self->{_name} = "NA";
268 :     if (defined($self->mgdata())) {
269 :     $self->{_name} = $self->mgdata()->genome_name;
270 :     }
271 :     } elsif ($source !~ /^RAST/) {
272 :     $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());
273 :     } else {
274 :     $self->{_name} = $self->figmodel()->get_genome_stats($self->genome())->{NAME}->[0];
275 :     }
276 :     }
277 :    
278 :     return $self->{_name};
279 :     }
280 :    
281 :     =head3 get_reaction_class
282 :     Definition:
283 :     string = FIGMODELmodel->get_reaction_class(string::reaction ID);
284 :     Description:
285 :     Returns reaction class
286 :     =cut
287 :     sub get_reaction_class {
288 :     my ($self,$reaction,$nohtml,$brief_flux) = @_;
289 :    
290 :     if (!-e $self->directory()."ReactionClassification-".$self->id().".tbl") {
291 :     if (!defined($self->{_reaction_classes})) {
292 :     $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id()."-Complete.tbl",";","|",0,["REACTION"]);
293 :     if (!defined($self->{_reaction_classes})) {
294 :     return undef;
295 :     }
296 :     }
297 :    
298 :     my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
299 :     if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
300 :     my $class;
301 :     my $min = $ClassRow->{MIN}->[0];
302 :     my $max = $ClassRow->{MAX}->[0];
303 :     if ($ClassRow->{CLASS}->[0] eq "Positive") {
304 :     $class = "Essential =>";
305 :     $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
306 :     } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {
307 :     $class = "Essential <=";
308 : dejongh 1.25 $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $class.="<br>[Flux: ".$max." to ".$min."]<br>";
309 : chenry 1.18 } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {
310 :     $class = "Active =>";
311 :     $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
312 :     } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {
313 :     $class = "Active <=";
314 : dejongh 1.25 $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $class.="<br>[Flux: ".$max." to ".$min."]<br>";
315 : chenry 1.18 } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {
316 :     $class = "Active <=>";
317 :     $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
318 :     } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {
319 :     $class = "Inactive";
320 :     } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {
321 :     $class = "Disconnected";
322 :     }
323 :    
324 :     if (!defined($nohtml) || $nohtml ne "1") {
325 :     $class = "<span title=\"Flux:".$min." to ".$max."\">".$class."</span>";
326 :     }
327 :    
328 :     return $class;
329 :     }
330 :     return undef;
331 :     }
332 :    
333 :     if (!defined($self->{_reaction_classes})) {
334 :     $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().".tbl",";","|",0,["REACTION"]);
335 :     if (!defined($self->{_reaction_classes})) {
336 :     return undef;
337 :     }
338 :     }
339 :    
340 :     my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
341 :     my $classstring = "";
342 :     if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
343 :     for (my $i=0; $i < @{$ClassRow->{CLASS}};$i++) {
344 :     if (length($classstring) > 0) {
345 :     $classstring .= "<br>";
346 :     }
347 :     my $NewClass;
348 :     my $min = $ClassRow->{MIN}->[$i];
349 :     my $max = $ClassRow->{MAX}->[$i];
350 :     if ($ClassRow->{CLASS}->[$i] eq "Positive") {
351 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";
352 :     $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
353 :     } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {
354 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";
355 : dejongh 1.25 $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $NewClass.="<br>[Flux: ".$max." to ".$min."]<br>";
356 : chenry 1.18 } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {
357 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";
358 :     $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
359 :     } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {
360 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";
361 : dejongh 1.25 $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $NewClass.="<br>[Flux: ".$max." to ".$min."]<br>";
362 : chenry 1.18 } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {
363 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";
364 :     $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
365 :     } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {
366 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";
367 :     } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {
368 :     $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";
369 :     }
370 :    
371 :     if (!defined($nohtml) || $nohtml ne "1") {
372 :     $NewClass = "<span title=\"Flux:".$min." to ".$max."\">".$NewClass."</span>";
373 :     }
374 :     $classstring .= $NewClass;
375 :     }
376 :     }
377 :     return $classstring;
378 :     }
379 :    
380 :     =head3 get_biomass
381 :     Definition:
382 :     string = FIGMODELmodel->get_biomass();
383 :     Description:
384 :     Returns data for the biomass reaction
385 :     =cut
386 :     sub get_biomass {
387 :     my ($self) = @_;
388 :    
389 :     if (!defined($self->{_biomass})) {
390 :     my $rxntbl = $self->reaction_table();
391 :     if (defined($rxntbl)) {
392 :     for (my $i=0; $i < $rxntbl->size(); $i++) {
393 :     if ($rxntbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {
394 :     $self->{_biomass} = $rxntbl->get_row($i)->{"LOAD"}->[0];
395 :     last;
396 :     }
397 :     }
398 :     }
399 :     }
400 :    
401 :     return $self->get_reaction_data($self->{_biomass});
402 :     }
403 :    
404 :     =head3 get_reaction_data
405 :     Definition:
406 :     string = FIGMODELmodel->get_reaction_data(string::reaction ID);
407 :     Description:
408 :     Returns model reaction data
409 :     =cut
410 :     sub get_reaction_data {
411 :     my ($self,$reaction) = @_;
412 :     if (!defined($self->reaction_table())) {
413 :     return undef;
414 :     }
415 :     if ($reaction =~ m/^\d+$/) {
416 : chenry 1.21 return $self->reaction_table()->get_row($reaction);
417 : chenry 1.18 }
418 :     return $self->reaction_table()->get_row_by_key($reaction,"LOAD");
419 :     }
420 :    
421 :     =head3 get_reaction_id
422 :     Definition:
423 :     string = FIGMODELmodel->get_reaction_id(string::reaction ID);
424 :     Description:
425 :     Returns model reaction id
426 :     =cut
427 :     sub get_reaction_id {
428 :     my ($self,$reaction) = @_;
429 :    
430 :     my $Data = $self->get_reaction_data($reaction);
431 :     if (!defined($Data) || !defined($Data->{LOAD}->[0])) {
432 :     return undef;
433 :     }
434 :     return $Data->{LOAD}->[0];
435 :     }
436 :    
437 :     =head3 load_model_table
438 :    
439 :     Definition: FIGMODELTable = FIGMODELmodel->load_model_table(string:table name,0/1:refresh the table));
440 :    
441 :     Description: Returns the table specified by the input filename. Table will be stored in a file in the model directory.
442 :    
443 :     =cut
444 :     sub load_model_table {
445 :     my ($self,$name,$refresh) = @_;
446 : chenry 1.19 if (defined($refresh) && $refresh == 1) {
447 : chenry 1.18 delete $self->{"_".$name};
448 :     }
449 :     if (!defined($self->{"_".$name})) {
450 : chenry 1.19 my $tbldef = $self->figmodel()->config($name);
451 : chenry 1.18 if (!defined($tbldef)) {
452 :     return undef;
453 :     }
454 : chenry 1.19 my $itemDelim = "|";
455 :     if (defined($tbldef->{itemdelimiter}->[0])) {
456 :     $itemDelim = $tbldef->{itemdelimiter}->[0];
457 :     if ($itemDelim eq "SC") {
458 :     $itemDelim = ";";
459 :     }
460 :     }
461 :     my $columnDelim = "\t";
462 :     if (defined($tbldef->{columndelimiter}->[0])) {
463 :     $columnDelim = $tbldef->{columndelimiter}->[0];
464 :     if ($columnDelim eq "SC") {
465 :     $columnDelim = ";";
466 :     }
467 :     }
468 :     my $suffix = ".tbl";
469 :     if (defined($tbldef->{filename_suffix}->[0])) {
470 :     $suffix = $tbldef->{filename_suffix}->[0];
471 :     }
472 :     my $filename = $self->directory().$name."-".$self->id().$self->selected_version().$suffix;
473 :     if (defined($tbldef->{filename_prefix}->[0])) {
474 :     if ($tbldef->{filename_prefix}->[0] eq "NONE") {
475 :     $filename = $self->directory().$self->id().$self->selected_version().$suffix;
476 :     } else {
477 :     $filename = $self->directory().$tbldef->{filename_prefix}->[0]."-".$self->id().$self->selected_version().$suffix;
478 :     }
479 :     }
480 :     if (-e $filename) {
481 :     $self->{"_".$name} = $self->figmodel()->database()->load_table($filename,$columnDelim,$itemDelim,$tbldef->{headingline}->[0],$tbldef->{hashcolumns});
482 :     } else {
483 : chenry 1.18 if (defined($tbldef->{prefix})) {
484 : chenry 1.19 $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim,join(@{$tbldef->{prefix}},"\n"));
485 : chenry 1.18 } else {
486 : chenry 1.19 $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim);
487 : chenry 1.18 }
488 :     }
489 :     }
490 :     return $self->{"_".$name};
491 :     }
492 :    
493 :     =head3 create_table_prototype
494 :    
495 :     Definition:
496 :     FIGMODELTable::table = FIGMODELmodel->create_table_prototype(string::table);
497 :     Description:
498 :     Returns a empty FIGMODELTable with all the metadata associated with the input table name
499 :    
500 :     =cut
501 :     sub create_table_prototype {
502 :     my ($self,$TableName) = @_;
503 :     #Checking if the table definition exists in the FIGMODELconfig file
504 : chenry 1.21 my $tbldef = $self->figmodel()->config($TableName);
505 :     if (!defined($tbldef)) {
506 : chenry 1.18 $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:Definition not found for ".$TableName);
507 :     return undef;
508 :     }
509 :     #Checking that this is a database table
510 : chenry 1.21 if (!defined($tbldef->{tabletype}) || $tbldef->{tabletype}->[0] ne "ModelTable") {
511 : chenry 1.18 $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:".$TableName." is not a model table!");
512 :     return undef;
513 :     }
514 : chenry 1.21 #Setting default values for table parameters
515 :     my $prefix;
516 :     if (defined($tbldef->{prefix})) {
517 :     $prefix = join("\n",@{$self->config($TableName)->{prefix}})."\n";
518 : chenry 1.18 }
519 : chenry 1.21 my $itemDelim = "|";
520 :     if (defined($tbldef->{itemdelimiter}->[0])) {
521 :     $itemDelim = $tbldef->{itemdelimiter}->[0];
522 :     if ($itemDelim eq "SC") {
523 :     $itemDelim = ";";
524 :     }
525 :     }
526 :     my $columnDelim = "\t";
527 :     if (defined($tbldef->{columndelimiter}->[0])) {
528 :     $columnDelim = $tbldef->{columndelimiter}->[0];
529 :     if ($columnDelim eq "SC") {
530 :     $columnDelim = ";";
531 :     }
532 :     }
533 :     my $suffix = ".tbl";
534 :     if (defined($tbldef->{filename_suffix}->[0])) {
535 :     $suffix = $tbldef->{filename_suffix}->[0];
536 :     }
537 :     my $filename = $self->directory().$TableName."-".$self->id().$self->selected_version().$suffix;
538 :     if (defined($tbldef->{filename_prefix}->[0])) {
539 :     if ($tbldef->{filename_prefix}->[0] eq "NONE") {
540 :     $filename = $self->directory().$self->id().$self->selected_version().$suffix;
541 :     } else {
542 :     $filename = $self->directory().$tbldef->{filename_prefix}->[0]."-".$self->id().$self->selected_version().$suffix;
543 :     }
544 : chenry 1.18 }
545 : chenry 1.21 #Creating the table prototype
546 :     my $tbl = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim,$prefix);
547 : chenry 1.18 return $tbl;
548 :     }
549 :    
550 :     =head3 get_reaction_number
551 :     Definition:
552 :     int = FIGMODELmodel->get_reaction_number();
553 :     Description:
554 :     Returns the number of reactions in the model
555 :     =cut
556 :     sub get_reaction_number {
557 :     my ($self) = @_;
558 :     if (!defined($self->reaction_table())) {
559 :     return 0;
560 :     }
561 :     return $self->reaction_table()->size();
562 :     }
563 :    
564 :     =head3 reaction_table
565 :     Definition:
566 :     FIGMODELTable = FIGMODELmodel->reaction_table();
567 :     Description:
568 :     Returns FIGMODELTable with the reaction list for the model
569 :     =cut
570 :     sub reaction_table {
571 : chenry 1.19 my ($self,$clear) = @_;
572 : chenry 1.22 if (defined($clear) && $clear == 1) {
573 :     delete $self->{_reaction_table};
574 : chenry 1.21 }
575 : chenry 1.22 if (!defined($self->{_reaction_table})) {
576 :     $self->{_reaction_table} = $self->load_model_table("ModelReactions",$clear);
577 :     my $classTbl = $self->reaction_class_table();
578 :     if (defined($classTbl)) {
579 :     for (my $i=0; $i < $classTbl->size(); $i++) {
580 :     my $row = $classTbl->get_row($i);
581 :     if (defined($row->{REACTION})) {
582 :     my $rxnRow = $self->{_reaction_table}->get_row_by_key($row->{"REACTION"}->[0],"LOAD");
583 :     if (defined($row->{MEDIA})) {
584 :     for (my $j=0; $j < @{$row->{MEDIA}};$j++) {
585 :     my $class = "Active <=>";
586 :     if ($row->{CLASS}->[$j] eq "Positive") {
587 :     $class = "Essential =>";
588 :     } elsif ($row->{CLASS}->[$j] eq "Negative") {
589 :     $class = "Essential <=";
590 :     } elsif ($row->{CLASS}->[$j] eq "Blocked") {
591 :     $class = "Inactive";
592 :     } elsif ($row->{CLASS}->[$j] eq "Positive variable") {
593 :     $class = "Active =>";
594 :     } elsif ($row->{CLASS}->[$j] eq "Negative variable") {
595 :     $class = "Active <=";
596 :     } elsif ($row->{CLASS}->[$j] eq "Variable") {
597 :     $class = "Active <=>";
598 :     } elsif ($row->{CLASS}->[$j] eq "Dead") {
599 :     $class = "Dead";
600 :     }
601 :     push(@{$rxnRow->{PREDICTIONS}},$row->{MEDIA}->[$j].":".$class);
602 : chenry 1.19 }
603 :     }
604 : chenry 1.18 }
605 :     }
606 :     }
607 :     }
608 : chenry 1.21 return $self->{_reaction_table};
609 :     }
610 :    
611 :     =head3 essentials_table
612 :     Definition:
613 :     FIGMODELTable = FIGMODELmodel->essentials_table();
614 :     Description:
615 :     Returns FIGMODELTable with the essential genes for the model
616 :     =cut
617 :     sub essentials_table {
618 :     my ($self,$clear) = @_;
619 :     my $tbl = $self->load_model_table("ModelEssentialGenes",$clear);
620 : chenry 1.19 return $tbl;
621 :     }
622 : chenry 1.18
623 : chenry 1.19 =head3 model_history
624 :     Definition:
625 :     FIGMODELTable = FIGMODELmodel->model_history();
626 :     Description:
627 :     Returns FIGMODELTable with the history of model changes
628 :     =cut
629 :     sub model_history {
630 :     my ($self,$clear) = @_;
631 :     return $self->load_model_table("ModelHistory",$clear);
632 : chenry 1.18 }
633 :    
634 :     =head3 feature_table
635 :     Definition:
636 :     FIGMODELTable = FIGMODELmodel->feature_table();
637 :     Description:
638 :     Returns FIGMODELTable with the feature list for the model
639 :     =cut
640 :     sub feature_table {
641 :     my ($self) = @_;
642 :    
643 :     if (!defined($self->{_feature_data})) {
644 :     #Getting the genome feature list
645 :     my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
646 :     if (!defined($FeatureTable)) {
647 :     print STDERR "FIGMODELmodel:feature_table:Could not get features for genome ".$self->genome()." in database!";
648 :     return undef;
649 :     }
650 :     #Getting the reaction table for the model
651 :     my $rxnTable = $self->reaction_table();
652 :     if (!defined($rxnTable)) {
653 : chenry 1.20 print STDERR "FIGMODELmodel:feature_table:Could not get reaction table for model ".$self->id()." in database!";
654 : chenry 1.18 return undef;
655 : chenry 1.20 }
656 : chenry 1.18 #Cloning the feature table
657 :     $self->{_feature_data} = $FeatureTable->clone_table_def();
658 :     $self->{_feature_data}->add_headings(($self->id()."REACTIONS",$self->id()."PREDICTIONS"));
659 : chenry 1.20 for (my $i=0; $i < $rxnTable->size(); $i++) {
660 :     my $Row = $rxnTable->get_row($i);
661 :     if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
662 :     foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
663 :     my $temp = $GeneSet;
664 :     $temp =~ s/\+/|/g;
665 :     $temp =~ s/\sAND\s/|/gi;
666 :     $temp =~ s/\sOR\s/|/gi;
667 :     $temp =~ s/[\(\)\s]//g;
668 :     my @GeneList = split(/\|/,$temp);
669 :     foreach my $Gene (@GeneList) {
670 :     my $FeatureRow = $self->{_feature_data}->get_row_by_key("fig|".$self->genome().".".$Gene,"ID");
671 : chenry 1.18 if (!defined($FeatureRow)) {
672 :     $FeatureRow = $FeatureTable->get_row_by_key("fig|".$self->genome().".".$Gene,"ID");
673 :     if (defined($FeatureRow)) {
674 :     $self->{_feature_data}->add_row($FeatureRow);
675 :     }
676 :     }
677 :     if (defined($FeatureRow)) {
678 :     $self->{_feature_data}->add_data($FeatureRow,$self->id()."REACTIONS",$Row->{"LOAD"}->[0],1);
679 :     }
680 : chenry 1.20 }
681 :     }
682 :     }
683 :     }
684 :     #Loading predictions
685 : chenry 1.21 my $esstbl = $self->essentials_table();
686 :     for (my $i=0; $i < $self->{_feature_data}->size(); $i++) {
687 :     my $Row = $self->{_feature_data}->get_row($i);
688 :     if ($Row->{ID}->[0] =~ m/(peg\.\d+)/) {
689 :     my $gene = $1;
690 :     my @rows = $esstbl->get_rows_by_key($gene,"ESSENTIAL GENES");
691 :     my $mediahash;
692 :     for (my $j=0; $j < $esstbl->size(); $j++) {
693 :     $mediahash->{$esstbl->get_row($j)->{MEDIA}->[0]} = 0;
694 :     }
695 :     for (my $j=0; $j < @rows; $j++) {
696 :     $mediahash->{$rows[$j]->{MEDIA}->[0]} = 1;
697 :     }
698 :     my @mediaList = keys(%{$mediahash});
699 :     for (my $j=0; $j < @mediaList; $j++) {
700 :     if ($mediahash->{$mediaList[$j]} == 1) {
701 :     push(@{$Row->{$self->id()."PREDICTIONS"}},$mediaList[$j].":essential");
702 :     } else {
703 :     push(@{$Row->{$self->id()."PREDICTIONS"}},$mediaList[$j].":nonessential");
704 : chenry 1.20 }
705 :     }
706 : chenry 1.21 }
707 : chenry 1.20 }
708 : chenry 1.18 }
709 :     return $self->{_feature_data};
710 :     }
711 :    
712 :     =head3 reaction_class_table
713 :     Definition:
714 :     FIGMODELTable = FIGMODELmodel->reaction_class_table();
715 :     Description:
716 :     Returns FIGMODELTable with the reaction class data, and creates the table file if it does not exist
717 :     =cut
718 :     sub reaction_class_table {
719 : chenry 1.19 my ($self,$clear) = @_;
720 :     return $self->load_model_table("ModelReactionClasses",$clear);
721 : chenry 1.18 }
722 :    
723 :     =head3 compound_class_table
724 :     Definition:
725 :     FIGMODELTable = FIGMODELmodel->compound_class_table();
726 :     Description:
727 :     Returns FIGMODELTable with the compound class data, and creates the table file if it does not exist
728 :     =cut
729 :     sub compound_class_table {
730 : chenry 1.19 my ($self,$clear) = @_;
731 :     return $self->load_model_table("ModelCompoundClasses",$clear);
732 : chenry 1.18 }
733 :    
734 :     =head3 get_essential_genes
735 :     Definition:
736 :     [string::peg ID] = FIGMODELmodel->get_essential_genes(string::media condition);
737 :     Description:
738 :     Returns an reference to an array of the predicted essential genes during growth in the input media condition
739 :     =cut
740 :     sub get_essential_genes {
741 :     my ($self,$media) = @_;
742 : chenry 1.21 my $tbl = $self->essentials_table();
743 :     my $row = $tbl->get_row_by_key($media,"MEDIA");
744 :     if (defined($row)) {
745 :     return $row->{"ESSENTIAL GENES"};
746 : chenry 1.18 }
747 : chenry 1.21 return undef;
748 : chenry 1.18 }
749 :    
750 :     =head3 compound_table
751 :     Definition:
752 :     FIGMODELTable = FIGMODELmodel->compound_table();
753 :     Description:
754 :     Returns FIGMODELTable with the compound list for the model
755 :     =cut
756 :     sub compound_table {
757 :     my ($self) = @_;
758 :    
759 :     if (!defined($self->{_compound_table})) {
760 : chenry 1.21 $self->{_compound_table} = $self->create_table_prototype("ModelCompounds");
761 :     #Loading the reactions
762 :     my $ReactionTable = $self->figmodel()->database()->get_table("REACTIONS");
763 :     my $BiomassTable = $self->figmodel()->database()->get_table("BIOMASS");
764 :     #Loading the model
765 :     my $ModelTable = $self->reaction_table();
766 :     #Checking that the tables were loaded
767 :     if (!defined($ModelTable) || !defined($ReactionTable)) {
768 :     return undef;
769 :     }
770 :     #Finding the biomass reaction
771 :     for (my $i=0; $i < $ModelTable->size(); $i++) {
772 :     my $ID = $ModelTable->get_row($i)->{"LOAD"}->[0];
773 :     my $Row = $ReactionTable->get_row_by_key($ID,"DATABASE");
774 :     my $IsBiomass = 0;
775 :     if (!defined($Row)) {
776 :     $Row = $BiomassTable->get_row_by_key($ID,"DATABASE");
777 :     $IsBiomass = 1;
778 :     }
779 :     if (defined($Row->{"EQUATION"}->[0])) {
780 :     $_ = $Row->{"EQUATION"}->[0];
781 :     my @OriginalArray = /(cpd\d\d\d\d\d[\[\w]*)/g;
782 :     foreach my $Compound (@OriginalArray) {
783 :     my $ID = substr($Compound,0,8);
784 :     my $NewRow = $self->{_compound_table}->get_row_by_key($ID,"DATABASE",1);
785 :     if ($IsBiomass == 1) {
786 :     $self->{_compound_table}->add_data($NewRow,"BIOMASS",$Row->{"DATABASE"}->[0],1);
787 :     }
788 :     if (length($Compound) > 8) {
789 :     #print $Compound."\t".$Row->{"EQUATION"}->[0]."\t".$Row->{"DATABASE"}->[0]."\n";
790 :     my $Compartment = substr($Compound,8,1);
791 :     $self->{_compound_table}->add_data($NewRow,"COMPARTMENTS",$Compartment,1);
792 :     $self->{_compound_table}->add_data($NewRow,"TRANSPORTERS",$Row->{"DATABASE"}->[0],1);
793 :     }
794 :     }
795 :     }
796 :     }
797 : chenry 1.18 }
798 : chenry 1.21
799 : chenry 1.18 return $self->{_compound_table};
800 :     }
801 :    
802 :     =head3 get_compound_data
803 :     Definition:
804 :     {string:key=>[string]:values} = FIGMODELmodel->get_compound_data(string::compound ID);
805 :     Description:
806 :     Returns model compound data
807 :     =cut
808 :     sub get_compound_data {
809 :     my ($self,$compound) = @_;
810 :    
811 :     if (!defined($self->compound_table())) {
812 :     return undef;
813 :     }
814 :     if ($compound =~ m/^\d+$/) {
815 :     return $self->compound_table()->get_row($compound);
816 :     }
817 :     return $self->compound_table()->get_row_by_key($compound,"DATABASE");
818 :     }
819 :    
820 :     =head3 get_feature_data
821 :     Definition:
822 :     {string:key=>[string]:values} = FIGMODELmodel->get_feature_data(string::feature ID);
823 :     Description:
824 :     Returns model feature data
825 :     =cut
826 :     sub get_feature_data {
827 :     my ($self,$feature) = @_;
828 :     if (!defined($self->feature_table())) {
829 :     return undef;
830 :     }
831 :     if ($feature =~ m/^\d+$/) {
832 :     return $self->feature_table()->get_row($feature);
833 :     }
834 :     if ($feature =~ m/(peg\.\d+)/) {
835 :     $feature = $1;
836 :     }
837 :     return $self->feature_table()->get_row_by_key("fig|".$self->genome().".".$feature,"ID");
838 :     }
839 :    
840 :     =head3 genome
841 :     Definition:
842 :     string = FIGMODELmodel->genome();
843 :     Description:
844 :     Returns model genome
845 :     =cut
846 :     sub genome {
847 :     my ($self) = @_;
848 :     return $self->{_data}->genome();
849 :     }
850 :    
851 :     =head3 source
852 :     Definition:
853 :     string = FIGMODELmodel->source();
854 :     Description:
855 :     Returns model source
856 :     =cut
857 :     sub source {
858 :     my ($self) = @_;
859 :     return $self->{_data}->source();
860 :     }
861 :    
862 :     =head3 rights
863 :     Definition:
864 :     1/0 = FIGMODELmodel->rights(string::username);
865 :     Description:
866 :     Returns 1 if the input user can view the model, and zero otherwise
867 :     =cut
868 :     sub rights {
869 :     my ($self,$username) = @_;
870 :     if ($self->public()) {
871 :     return 1;
872 :     }
873 :     if (!defined($username) || $username eq "NONE") {
874 :     return 0;
875 :     }
876 :     if (!defined($self->{_userrights}->{$username})) {
877 :     $self->{_userrights}->{$self->{_data}->owner()} = 1;
878 :     my @users = split(/\|/,$self->{_data}->users());
879 :     for (my $i=0; $i < @users;$i++) {
880 :     $self->{_userrights}->{$users[$i]} = 1;
881 :     }
882 :     }
883 :     return $self->{_userrights}->{$username};
884 :     }
885 :    
886 :     =head3 public
887 :     Definition:
888 :     1/0 = FIGMODELmodel->public();
889 :     Description:
890 :     Returns 1 if the model is public, and zero otherwise
891 :     =cut
892 :     sub public {
893 :     my ($self) = @_;
894 :     if ($self->{_data}->users() eq "all") {
895 :     return 1;
896 :     }
897 :     return 0;
898 :     }
899 :    
900 :     =head3 directory
901 :     Definition:
902 :     string = FIGMODELmodel->directory();
903 :     Description:
904 :     Returns model directory
905 :     =cut
906 :     sub directory {
907 :     my ($self) = @_;
908 :    
909 :     if (!defined($self->{_directory})) {
910 :     my $userdirectory = "";
911 :     if ($self->owner() ne "master") {
912 :     $userdirectory = $self->owner()."/";
913 :     }
914 :     my $source = $self->source();
915 :     if ($source =~ /^MGRAST/) {
916 : chenry 1.21 $self->{_directory} = $self->figmodel()->config("mgrast model directory")->[0].$userdirectory.$self->genome()."/";
917 : chenry 1.18 } elsif ($source =~ /^RAST/) {
918 :     $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
919 :     } elsif ($source =~ /^SEED/) {
920 :     $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
921 :     } elsif ($source =~ /^PM/) {
922 :     if (length($userdirectory) == 0) {
923 :     $self->{_directory} = $self->figmodel()->config("imported model directory")->[0].$self->id()."/";
924 :     } else {
925 :     $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->id()."/";
926 :     }
927 :     }
928 :     }
929 :    
930 :     return $self->{_directory};
931 :     }
932 :    
933 :     =head3 filename
934 :     Definition:
935 :     string = FIGMODELmodel->filename();
936 :     Description:
937 :     Returns model filename
938 :     =cut
939 :     sub filename {
940 :     my ($self) = @_;
941 :    
942 :     return $self->directory().$self->id().$self->selected_version().".txt";
943 :     }
944 :    
945 :     =head3 version
946 :     Definition:
947 :     string = FIGMODELmodel->version();
948 :     Description:
949 :     Returns the version of the model
950 :     =cut
951 :     sub version {
952 :     my ($self) = @_;
953 :    
954 :     if (!defined($self->{_version})) {
955 :     if (!defined($self->{_selectedversion})) {
956 :     $self->{_version} = "V".$self->{_data}->version().".".$self->{_data}->autocompleteVersion();
957 :     } else {
958 :     $self->{_version} = $self->{_selectedversion};
959 :     }
960 :     }
961 :     return $self->{_version};
962 :     }
963 :    
964 :     =head3 selected_version
965 :     Definition:
966 :     string = FIGMODELmodel->selected_version();
967 :     Description:
968 :     Returns the selected version of the model
969 :     =cut
970 :     sub selected_version {
971 :     my ($self) = @_;
972 :    
973 :     if (!defined($self->{_selectedversion})) {
974 :     return "";
975 :     }
976 :     return $self->{_selectedversion};
977 :     }
978 :    
979 :     =head3 modification_time
980 :     Definition:
981 :     string = FIGMODELmodel->modification_time();
982 :     Description:
983 :     Returns the selected version of the model
984 :     =cut
985 :     sub modification_time {
986 :     my ($self) = @_;
987 :     return $self->{_data}->modificationDate();
988 :     }
989 :    
990 :     =head3 gene_reactions
991 :     Definition:
992 :     string = FIGMODELmodel->gene_reactions();
993 :     Description:
994 :     Returns the number of reactions added by the gap filling
995 :     =cut
996 :     sub gene_reactions {
997 :     my ($self) = @_;
998 :     return ($self->{_data}->reactions() - $self->{_data}->autoCompleteReactions() - $self->{_data}->spontaneousReactions() - $self->{_data}->gapFillReactions());
999 :     }
1000 :    
1001 :     =head3 total_compounds
1002 :     Definition:
1003 :     string = FIGMODELmodel->total_compounds();
1004 :     Description:
1005 :     Returns the number of compounds in the model
1006 :     =cut
1007 :     sub total_compounds {
1008 :     my ($self) = @_;
1009 :     return $self->{_data}->compounds();
1010 :     }
1011 :    
1012 :     =head3 gapfilling_reactions
1013 :     Definition:
1014 :     string = FIGMODELmodel->gapfilling_reactions();
1015 :     Description:
1016 :     Returns the number of reactions added by the gap filling
1017 :     =cut
1018 :     sub gapfilling_reactions {
1019 :     my ($self) = @_;
1020 :     return ($self->{_data}->autoCompleteReactions()+$self->{_data}->gapFillReactions());
1021 :     }
1022 :    
1023 :     =head3 total_reactions
1024 :     Definition:
1025 :     string = FIGMODELmodel->total_reactions();
1026 :     Description:
1027 :     Returns the total number of reactions in the model
1028 :     =cut
1029 :     sub total_reactions {
1030 :     my ($self) = @_;
1031 :     return $self->{_data}->reactions();
1032 :     }
1033 :    
1034 :     =head3 model_genes
1035 :     Definition:
1036 :     string = FIGMODELmodel->model_genes();
1037 :     Description:
1038 :     Returns the number of genes mapped to one or more reactions in the model
1039 :     =cut
1040 :     sub model_genes {
1041 :     my ($self) = @_;
1042 :     return $self->{_data}->associatedGenes();
1043 :     }
1044 :    
1045 :     =head3 class
1046 :     Definition:
1047 :     string = FIGMODELmodel->class();
1048 :     Description:
1049 :     Returns the class of the model: gram positive, gram negative, other
1050 :     =cut
1051 :     sub class {
1052 :     my ($self) = @_;
1053 :     return $self->{_data}->cellwalltype();
1054 :     }
1055 :    
1056 : chenry 1.19 sub autocompleteMedia {
1057 : chenry 1.22 my ($self,$newMedia) = @_;
1058 : chenry 1.24 if (defined($newMedia)) {
1059 :     return $self->{_data}->autoCompleteMedia($newMedia);
1060 :     }
1061 :     return $self->{_data}->autoCompleteMedia();
1062 : chenry 1.22 }
1063 :    
1064 :     sub biomassReaction {
1065 :     my ($self,$newBiomass) = @_;
1066 :     if (!defined($newBiomass)) {
1067 :     return $self->{_data}->biomassReaction();
1068 :     } else {
1069 :     #Figuring out what the old biomass is
1070 :     my $oldBiomass = $self->{_data}->biomassReaction();
1071 : chenry 1.27 $self->{_data}->biomassReaction($newBiomass);
1072 :     #Changing the biomass reaction in the model file
1073 :     my $rxnTbl = $self->reaction_table();
1074 :     for (my $i=0; $i < $rxnTbl->size(); $i++) {
1075 :     my $row = $rxnTbl->get_row($i);
1076 :     if ($row->{LOAD}->[0] =~ m/^bio/) {
1077 :     $row->{LOAD}->[0] = $newBiomass;
1078 :     }
1079 :     }
1080 :     $rxnTbl->save();
1081 : chenry 1.22 if ($newBiomass ne $oldBiomass) {
1082 :     #Figuring out if the new biomass exists
1083 : chenry 1.27 my $handle = $self->figmodel()->database()->get_object_manager("bof");
1084 :     my $objects = $handle->get_objects({id=>$newBiomass});
1085 : chenry 1.22 if (!defined($objects) || !defined($objects->[0])) {
1086 :     print STDERR "Could not find new biomass reaction ".$newBiomass."\n";
1087 :     return $oldBiomass;
1088 :     }
1089 :     }
1090 : chenry 1.27 return $newBiomass;
1091 : chenry 1.22 }
1092 : chenry 1.19 }
1093 :    
1094 : chenry 1.29 =head3 growth
1095 :     Definition:
1096 :     double = FIGMODELmodel->growth();
1097 :     Description:
1098 :     =cut
1099 :     sub growth {
1100 :     my ($self,$inGrowth) = @_;
1101 :     if (!defined($inGrowth)) {
1102 :     return $self->{_data}->growth();
1103 :     } else {
1104 :     return $self->{_data}->growth($inGrowth);
1105 :     }
1106 :     }
1107 :    
1108 :     =head3 noGrowthCompounds
1109 :     Definition:
1110 :     string = FIGMODELmodel->noGrowthCompounds();
1111 :     Description:
1112 :     =cut
1113 :     sub noGrowthCompounds {
1114 :     my ($self,$inCompounds) = @_;
1115 :     if (!defined($inCompounds)) {
1116 :     return $self->{_data}->noGrowthCompounds();
1117 :     } else {
1118 :     return $self->{_data}->noGrowthCompounds($inCompounds);
1119 :     }
1120 :     }
1121 :    
1122 : chenry 1.18 =head3 taxonomy
1123 :     Definition:
1124 :     string = FIGMODELmodel->taxonomy();
1125 :     Description:
1126 :     Returns model taxonomy or biome if this is an metagenome model
1127 :     =cut
1128 :     sub taxonomy {
1129 :     my ($self) = @_;
1130 :    
1131 :     if (!defined($self->{_taxonomy})) {
1132 :     my $source = $self->source();
1133 :     if ($source =~ /^MGRAST/) {
1134 :     $self->{_taxonomy} = "NA";
1135 :     my $mgmodeldata = $self->figmodel()->database()->mg_model_data($self->genome());
1136 :     if (defined($mgmodeldata)) {
1137 :     $self->{_taxonomy} = $mgmodeldata->biome();
1138 :     }
1139 :     } else {
1140 :     $self->{_taxonomy} = $self->fig()->taxonomy_of($self->genome());
1141 :     }
1142 :     }
1143 :    
1144 :     return $self->{_taxonomy};
1145 :     }
1146 :    
1147 :     =head3 genome_size
1148 :     Definition:
1149 :     string = FIGMODELmodel->genome_size();
1150 :     Description:
1151 :     Returns size of the modeled genome in KB
1152 :     =cut
1153 :     sub genome_size {
1154 :     my ($self) = @_;
1155 :    
1156 :     if (!defined($self->{_genome_size})) {
1157 :     my $source = $self->source();
1158 :     if ($source =~ /^MGRAST/) {
1159 :     $self->{_genome_size} = "NA";
1160 :     if (defined($self->mgdata())) {
1161 :     $self->{_genome_size} = $self->mgdata()->size;
1162 :     }
1163 :     } else {
1164 :     $self->{_genome_size} = $self->fig()->genome_szdna($self->genome());
1165 :     }
1166 :     }
1167 :    
1168 :     return $self->{_genome_size};
1169 :     }
1170 :    
1171 :     =head3 genome_genes
1172 :     Definition:
1173 :     string = FIGMODELmodel->genome_genes();
1174 :     Description:
1175 :     Returns the number of genes in the modeled genome
1176 :     =cut
1177 :     sub genome_genes {
1178 :     my ($self) = @_;
1179 :    
1180 :     if (!defined($self->{_genome_genes})) {
1181 :     my $source = $self->source();
1182 :     if ($source =~ /^MGRAST/) {
1183 :     $self->{_genome_genes} = "NA";
1184 :     if (defined($self->mgdata())) {
1185 :     $self->{_genome_genes} = $self->mgdata()->genome_contig_count;
1186 :     }
1187 :     } else {
1188 :     $self->{_genome_genes} = $self->figmodel()->get_genome_stats($self->genome())->{"TOTAL GENES"}->[0];
1189 :     }
1190 :     }
1191 :    
1192 :     return $self->{_genome_genes};
1193 :     }
1194 :    
1195 :     =head3 run_default_model_predictions
1196 :     Definition:
1197 :     0/1::status = FIGMODELmodel->run_default_model_predictions(string::media ID);
1198 :     Description:
1199 :     =cut
1200 :     sub run_default_model_predictions {
1201 :     my ($self,$Media) = @_;
1202 :    
1203 :     #Assuming complete media if none is provided
1204 :     if (!defined($Media)) {
1205 : chenry 1.21 $Media = $self->autocompleteMedia();
1206 : chenry 1.18 }
1207 :    
1208 :     #Predicting essentiality
1209 :     my $result = $self->figmodel()->RunFBASimulation($self->id(),"SINGLEKO",undef,undef,[$self->id()],[$Media]);
1210 :     #Checking that the table is defined and the output file exists
1211 :     if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {
1212 : chenry 1.21 my $tbl = $self->essentials_table();
1213 :     my $row = $tbl->get_row_by_key($Media,"MEDIA",1);
1214 :     $row->{"ESSENTIAL GENES"} = $result->get_row(0)->{"ESSENTIALGENES"};
1215 :     $tbl->save();
1216 : chenry 1.18 } else {
1217 :     $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not identify essential reactions for model ".$self->id().$self->selected_version().".");
1218 :     return $self->figmodel()->fail();
1219 :     }
1220 :    
1221 :     #Classifying reactions and compounds
1222 :     my $tbl = $self->classify_model_reactions($Media);
1223 :     if (!defined($tbl)) {
1224 :     $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not classify reactions for model ".$self->id().$self->selected_version().".");
1225 :     return $self->figmodel()->fail();
1226 :     }
1227 :     $tbl->save();
1228 :    
1229 :     return $self->figmodel()->success();
1230 :     }
1231 :    
1232 :     =head3 update_stats_for_gap_filling
1233 :     Definition:
1234 :     {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);
1235 :     Description:
1236 :     =cut
1237 :     sub update_stats_for_gap_filling {
1238 :     my ($self,$gapfilltime) = @_;
1239 :     $self->{_data}->autoCompleteTime($gapfilltime);
1240 :     $self->{_data}->autocompleteDate(time());
1241 :     $self->{_data}->modificationDate(time());
1242 :     my $version = $self->{_data}->autocompleteVersion();
1243 :     $self->{_data}->autocompleteVersion($version+1);
1244 : chenry 1.22 $self->update_model_stats();
1245 : chenry 1.18 }
1246 :    
1247 :     =head3 update_stats_for_build
1248 :     Definition:
1249 :     {string => [string]} = FIGMODELmodel->update_stats_for_build();
1250 :     Description:
1251 :     =cut
1252 :     sub update_stats_for_build {
1253 :     my ($self) = @_;
1254 :     $self->{_data}->builtDate(time());
1255 :     $self->{_data}->modificationDate(time());
1256 :     my $version = $self->{_data}->version();
1257 :     $self->{_data}->version($version+1);
1258 : chenry 1.28 $self->update_model_stats();
1259 : chenry 1.18 }
1260 :    
1261 :     =head3 update_model_stats
1262 :     Definition:
1263 :     FIGMODELmodel->update_model_stats();
1264 :     Description:
1265 :     =cut
1266 :     sub update_model_stats {
1267 :     my ($self) = @_;
1268 :    
1269 :     #Getting reaction table
1270 :     my $rxntbl = $self->reaction_table();
1271 :     if (!defined($rxntbl)) {
1272 : chenry 1.27 die $self->error_message("update_model_stats:Could not load reaction list!");
1273 : chenry 1.18 }
1274 :     my $cpdtbl = $self->compound_table();
1275 :    
1276 :     #Calculating all necessary stats
1277 :     my %GeneHash;
1278 :     my %NonpegHash;
1279 :     my %CompoundHash;
1280 :     my $spontaneousReactions = 0;
1281 :     my $gapFillReactions = 0;
1282 :     my $biologReactions = 0;
1283 :     my $transporters = 0;
1284 :     my $autoCompleteReactions = 0;
1285 :     my $associatedSubsystemGenes = 0;
1286 :     for (my $i=0; $i < $rxntbl->size(); $i++) {
1287 :     my $Row = $rxntbl->get_row($i);
1288 :     if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
1289 :     my $ReactionRow = $self->figmodel()->get_reaction($Row->{"LOAD"}->[0]);
1290 :     if (defined($ReactionRow->{"EQUATION"}->[0])) {
1291 :     #Checking for extracellular metabolites which indicate that this is a transporter
1292 :     if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {
1293 :     $transporters++;
1294 :     }
1295 :     }
1296 :     #Identifying spontaneous/biolog/gapfilling/gene associated reactions
1297 :     if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {
1298 :     $biologReactions++;
1299 :     } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GROW/i) {
1300 :     $gapFillReactions++;
1301 :     } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {
1302 :     $spontaneousReactions++;
1303 :     } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {
1304 :     $autoCompleteReactions++;
1305 :     } else {
1306 :     foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
1307 :     $_ = $GeneSet;
1308 :     my @GeneList = /(peg\.\d+)/g;
1309 :     foreach my $Gene (@GeneList) {
1310 :     if ($Gene =~ m/(peg\.\d+)/) {
1311 :     $GeneHash{$1} = 1;
1312 :     } else {
1313 :     $NonpegHash{$Gene} = 1;
1314 :     }
1315 :     }
1316 :     }
1317 :     }
1318 :     }
1319 :     }
1320 :     my @genes = keys(%GeneHash);
1321 :     my @othergenes = keys(%NonpegHash);
1322 :    
1323 :     #Setting the reaction count
1324 :     $self->{_data}->reactions($rxntbl->size());
1325 :     #Setting the metabolite count
1326 : chenry 1.21 $self->{_data}->compounds($cpdtbl->size());
1327 : chenry 1.18 #Setting the gene count
1328 :     my $geneCount = @genes + @othergenes;
1329 :     $self->{_data}->associatedGenes($geneCount);
1330 :     #Setting remaining stats
1331 :     $self->{_data}->spontaneousReactions($spontaneousReactions);
1332 :     $self->{_data}->gapFillReactions($gapFillReactions);
1333 :     $self->{_data}->biologReactions($biologReactions);
1334 :     $self->{_data}->transporters($transporters);
1335 :     $self->{_data}->autoCompleteReactions($autoCompleteReactions);
1336 :     $self->{_data}->associatedSubsystemGenes($associatedSubsystemGenes);
1337 :     #Setting the model class
1338 :     my $class = "";
1339 :     for (my $i=0; $i < @{$self->figmodel()->config("class list")}; $i++) {
1340 :     if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i]))) {
1341 :     if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->id()})) {
1342 :     $class = $self->figmodel()->config("class list")->[$i];
1343 :     last;
1344 :     }
1345 :     if ($class eq "" && defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->genome()})) {
1346 :     $class = $self->figmodel()->config("class list")->[$i];
1347 :     }
1348 :     }
1349 :     }
1350 :     if ($class eq "") {
1351 :     $class = $self->figmodel()->get_genome_stats($self->genome())->{CLASS}->[0];
1352 :     }
1353 :     if ($class eq "") {
1354 :     $class = "unknown";
1355 :     }
1356 :     $self->{_data}->cellwalltype($class);
1357 :     }
1358 :    
1359 :     =head3 GapFillModel
1360 :     Definition:
1361 :     (success/fail) = FIGMODELmodel->GapFillModel();
1362 :     Description:
1363 :     This function performs an optimization to identify the minimal set of reactions that must be added to a model in order for biomass to be produced by the biomass reaction in the model.
1364 :     Before running the gap filling, the existing model is backup in the same directory with the current version numbers appended.
1365 :     If the model has been gap filled previously, the previous gap filling reactions are removed prior to running the gap filling again.
1366 :     =cut
1367 :     sub GapFillModel {
1368 :     my ($self,$donotclear) = @_;
1369 :    
1370 :     #Setting status of model to gap filling
1371 :     my $OrganismID = $self->genome();
1372 :     $self->set_status(1,"Auto completion running");
1373 :     my $UniqueFilename = $self->figmodel()->filename();
1374 :     my $StartTime = time();
1375 : chenry 1.19
1376 :     #Reading original reaction table
1377 :     my $OriginalRxn = $self->reaction_table();
1378 :     #Clearing the table
1379 :     $self->reaction_table(1);
1380 : chenry 1.18
1381 :     #Removing any gapfilling reactions that may be currently present in the model
1382 :     if (!defined($donotclear) || $donotclear != 1) {
1383 :     my $ModelTable = $self->reaction_table();
1384 :     for (my $i=0; $i < $ModelTable->size(); $i++) {
1385 : chenry 1.19 if (!defined($ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0]) || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "AUTOCOMPLETION" || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "GAP FILLING" || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/ || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/GROWMATCH/) {
1386 : chenry 1.18 $ModelTable->delete_row($i);
1387 :     $i--;
1388 :     }
1389 :     }
1390 :     $ModelTable->save();
1391 :     }
1392 :    
1393 :     #Calling the MFAToolkit to run the gap filling optimization
1394 : chenry 1.22 my $Media = $self->autocompleteMedia();
1395 :     if ($Media eq "Complete") {
1396 :     system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1397 :     } else {
1398 : chenry 1.18 #Loading media, changing bounds, saving media as a test media
1399 :     my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);
1400 :     for (my $i=0; $i < $MediaTable->size(); $i++) {
1401 :     if ($MediaTable->get_row($i)->{"Min"}->[0] < 0) {
1402 :     $MediaTable->get_row($i)->{"Min"}->[0] = -10000;
1403 :     }
1404 :     if ($MediaTable->get_row($i)->{"Max"}->[0] > 0) {
1405 :     $MediaTable->get_row($i)->{"Max"}->[0] = 10000;
1406 :     }
1407 :     }
1408 :     $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1409 : chenry 1.20 #print $self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0,"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef)."\n";
1410 : chenry 1.18 system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0,"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1411 :     unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1412 :     }
1413 :    
1414 :     #Parse the solutions file for the model and get the reaction list from it
1415 :     my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);
1416 : chenry 1.19
1417 : chenry 1.18 #Clearing the mfatoolkit output and log file
1418 :     $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");
1419 :     if (!defined($SolutionData) || $SolutionData->size() == 0) {
1420 :     $self->set_status(1,"Auto completion failed: auto completion solution not found");
1421 :     $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");
1422 :     return $self->fail();
1423 :     }
1424 :    
1425 :     #Looking for the last printed recursive MILP solution
1426 :     for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
1427 :     if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {
1428 :     my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);
1429 :     my @TempThree = split(/\|/,$AllSolutions);
1430 :     if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {
1431 :     my @TempFour = split(/,/,$1);
1432 :     my $DirectionList;
1433 :     my $ReactionList;
1434 :     for (my $j=0; $j < @TempFour; $j++) {
1435 :     my $ID = "";
1436 :     my $Sign = "=>";
1437 :     if ($TempFour[$j] =~ m/\-(rxn\d\d\d\d\d)/) {
1438 :     $ID = $1;
1439 :     $Sign = "<=";
1440 :     } elsif ($TempFour[$j] =~ m/\+(rxn\d\d\d\d\d)/) {
1441 :     $ID = $1;
1442 :     }
1443 :     if ($ID ne "") {
1444 :     if ($self->figmodel()->reversibility_of_reaction($ID) ne $Sign) {
1445 :     $Sign = "<=>";
1446 :     }
1447 :     push(@{$DirectionList},$Sign);
1448 :     push(@{$ReactionList},$ID);
1449 :     }
1450 :     }
1451 : chenry 1.19 $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"AUTOCOMPLETION",0,1);
1452 : chenry 1.18 }
1453 :     last;
1454 :     }
1455 :     }
1456 : chenry 1.19
1457 : chenry 1.18 #Updating model stats with gap filling results
1458 :     my $ElapsedTime = time() - $StartTime;
1459 : chenry 1.19 $self->reaction_table(1);
1460 :     $self->calculate_model_changes($OriginalRxn,"AUTOCOMPLETION");
1461 :    
1462 : chenry 1.18 #Determining why each gap filling reaction was added
1463 :     $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);
1464 :     if (!defined($donotclear) || $donotclear != 1) {
1465 :     if ($self->id() !~ m/MGRast/) {
1466 :     $self->update_stats_for_gap_filling($ElapsedTime);
1467 :     }
1468 :     } else {
1469 :     $self->update_model_stats();
1470 :     }
1471 :     #Printing the updated SBML file
1472 :     $self->PrintSBMLFile();
1473 :     $self->PrintModelLPFile($self->id());
1474 :     $self->set_status(1,"Auto completion successfully finished");
1475 :     $self->run_default_model_predictions();
1476 :    
1477 :     return $self->success();
1478 :     }
1479 :    
1480 : chenry 1.19 =head3 calculate_model_changes
1481 :     Definition:
1482 :     FIGMODELmodel->calculate_model_changes(FIGMODELTable:original reaction table,string:modification cause);
1483 :     Description:
1484 :    
1485 :     =cut
1486 :    
1487 :     sub calculate_model_changes {
1488 : chenry 1.20 my ($self,$originalReactions,$cause,$tbl,$version) = @_;
1489 : chenry 1.19 my $modTime = time();
1490 : chenry 1.20 if (!defined($version)) {
1491 :     $version = $self->selected_version();
1492 :     }
1493 : chenry 1.19 my $user = $self->figmodel()->user();
1494 :     #Getting the history table
1495 :     my $histTbl = $self->model_history();
1496 :     #Checking for differences
1497 : chenry 1.20 if (!defined($tbl)) {
1498 :     $tbl = $self->reaction_table();
1499 :     }
1500 : chenry 1.19 for (my $i=0; $i < $tbl->size(); $i++) {
1501 :     my $row = $tbl->get_row($i);
1502 :     my $orgRow = $originalReactions->get_row_by_key($row->{LOAD}->[0],"LOAD");
1503 :     if (!defined($orgRow)) {
1504 :     $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => $row->{DIRECTIONALITY}, GeneChange => $row->{"ASSOCIATED PEG"}, Action => ["ADDED"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1505 :     } else {
1506 :     my $geneChanges;
1507 :     my $directionChange;
1508 :     if ($orgRow->{"DIRECTIONALITY"}->[0] ne $row->{"DIRECTIONALITY"}->[0]) {
1509 :     $directionChange = $orgRow->{"DIRECTIONALITY"}->[0]." to ".$row->{"DIRECTIONALITY"}->[0];
1510 :     }
1511 :     for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
1512 :     my $match = 0;
1513 : chenry 1.20 if (defined($orgRow->{"ASSOCIATED PEG"})) {
1514 :     for (my $k=0; $k < @{$orgRow->{"ASSOCIATED PEG"}}; $k++) {
1515 :     if ($row->{"ASSOCIATED PEG"}->[$j] eq $orgRow->{"ASSOCIATED PEG"}->[$k]) {
1516 :     $match = 1;
1517 :     }
1518 : chenry 1.19 }
1519 :     }
1520 :     if ($match == 0) {
1521 :     push(@{$geneChanges},"Added ".$row->{"ASSOCIATED PEG"}->[$j]);
1522 :     }
1523 :     }
1524 : chenry 1.20 if (defined($orgRow->{"ASSOCIATED PEG"})) {
1525 :     for (my $k=0; $k < @{$orgRow->{"ASSOCIATED PEG"}}; $k++) {
1526 :     my $match = 0;
1527 :     if (defined($row->{"ASSOCIATED PEG"})) {
1528 :     for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
1529 :     if ($row->{"ASSOCIATED PEG"}->[$j] eq $orgRow->{"ASSOCIATED PEG"}->[$k]) {
1530 :     $match = 1;
1531 :     }
1532 :     }
1533 :     }
1534 :     if ($match == 0) {
1535 :     push(@{$geneChanges},"Removed ".$orgRow->{"ASSOCIATED PEG"}->[$k]);
1536 : chenry 1.19 }
1537 :     }
1538 :     }
1539 :     if ((defined($directionChange) && length($directionChange) > 0) || defined($geneChanges) && @{$geneChanges} > 0) {
1540 :     $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => [$directionChange], GeneChange => $geneChanges, Action => ["CHANGE"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1541 :     }
1542 :     }
1543 :     }
1544 :     #Looking for removed reactions
1545 :     for (my $i=0; $i < $originalReactions->size(); $i++) {
1546 :     my $row = $originalReactions->get_row($i);
1547 :     my $orgRow = $tbl->get_row_by_key($row->{LOAD}->[0],"LOAD");
1548 :     if (!defined($orgRow)) {
1549 :     $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => $row->{DIRECTIONALITY}, GeneChange => $row->{"ASSOCIATED PEG"}, Action => ["REMOVED"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1550 :     }
1551 :     }
1552 :     $histTbl->save();
1553 :     }
1554 :    
1555 : chenry 1.18 =head3 GapGenModel
1556 :     Definition:
1557 :     FIGMODELmodel->GapGenModel();
1558 :     Description:
1559 :     Runs the gap generation algorithm to correct a single false positive prediction. Results are loaded into a table.
1560 :     =cut
1561 :    
1562 :     sub GapGenModel {
1563 :     my ($self,$Media,$KOList,$NoKOList,$Experiment,$SolutionLimit) = @_;
1564 :    
1565 :     #Enforcing nonoptional arguments
1566 :     if (!defined($Media)) {
1567 :     return undef;
1568 :     }
1569 :     if (!defined($KOList)) {
1570 :     $KOList->[0] = "none";
1571 :     }
1572 :     if (!defined($NoKOList)) {
1573 :     $NoKOList->[0] = "none";
1574 :     }
1575 :     if (!defined($Experiment)) {
1576 :     $Experiment= "ReactionKO";
1577 :     }
1578 :     if (!defined($SolutionLimit)) {
1579 :     $SolutionLimit = "10";
1580 :     }
1581 :    
1582 :     #Translating the KO lists into arrays
1583 :     if (ref($KOList) ne "ARRAY") {
1584 :     my $temp = $KOList;
1585 :     $KOList = ();
1586 :     push(@{$KOList},split(/[,;]/,$temp));
1587 :     }
1588 :     my $noKOHash;
1589 :     if (defined($NoKOList) && ref($NoKOList) ne "ARRAY") {
1590 :     my $temp = $NoKOList;
1591 :     $NoKOList = ();
1592 :     push(@{$NoKOList},split(/[,;]/,$temp));
1593 :     foreach my $rxn (@{$NoKOList}) {
1594 :     $noKOHash->{$rxn} = 1;
1595 :     }
1596 :     }
1597 :    
1598 :     #Checking if solutions exist for the input parameters
1599 :     $self->aquireModelLock();
1600 :     my $tbl = $self->load_model_table("GapGenSolutions");
1601 :     my $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1602 :     my $solutions;
1603 :     if (defined($solutionRow)) {
1604 :     #Checking if any solutions conform to the no KO list
1605 :     foreach my $solution (@{$solutionRow->{Solutions}}) {
1606 :     my @reactions = split(/,/,$solution);
1607 :     my $include = 1;
1608 :     foreach my $rxn (@reactions) {
1609 :     if ($rxn =~ m/(rxn\d\d\d\d\d)/) {
1610 :     if (defined($noKOHash->{$1})) {
1611 :     $include = 0;
1612 :     }
1613 :     }
1614 :     }
1615 :     if ($include == 1) {
1616 :     push(@{$solutions},$solution);
1617 :     }
1618 :     }
1619 :     } else {
1620 :     $solutionRow = {Media => [$Media],Experiment => [$Experiment],KOlist => [join(",",@{$KOList})]};
1621 :     $tbl->add_row($solutionRow);
1622 :     $self->figmodel()->database()->save_table($tbl);
1623 :     }
1624 :     $self->releaseModelLock();
1625 :    
1626 :     #Returning solution list of solutions were found
1627 :     if (defined($solutions) && @{$solutions} > 0) {
1628 :     return $solutions;
1629 :     }
1630 :    
1631 :     #Getting unique filename
1632 :     my $Filename = $self->figmodel()->filename();
1633 :    
1634 :     #Running the gap generation
1635 :     system($self->figmodel()->GenerateMFAToolkitCommandLineCall($Filename,$self->id().$self->selected_version(),$Media,["GapGeneration"],{"Recursive MILP solution limit" => $SolutionLimit ,"Reactions that should always be active" => join(";",@{$NoKOList}),"Reactions to knockout" => join(";",@{$KOList}),"Reactions that are always blocked" => "none"},"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log",undef,undef));
1636 :     my $ProblemReport = $self->figmodel()->LoadProblemReport($Filename);
1637 :     if (!defined($ProblemReport)) {
1638 :     $self->figmodel()->error_message("FIGMODEL:GapGenerationAlgorithm;No problem report;".$Filename.";".$self->id().$self->selected_version().";".$Media.";".$KOList.";".$NoKOList);
1639 :     return undef;
1640 :     }
1641 :    
1642 :     #Clearing the output folder and log file
1643 :     $self->figmodel()->clearing_output($Filename,"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log");
1644 :    
1645 :     #Saving the solution
1646 :     $self->aquireModelLock();
1647 :     $tbl = $self->load_model_table("GapGenSolutions");
1648 :     $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1649 :     for (my $j=0; $j < $ProblemReport->size(); $j++) {
1650 :     if ($ProblemReport->get_row($j)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^)]+)/) {
1651 :     my @SolutionList = split(/\|/,$1);
1652 :     for (my $k=0; $k < @SolutionList; $k++) {
1653 :     if ($SolutionList[$k] =~ m/(\d+):(.+)/) {
1654 :     push(@{$solutionRow->{Solutions}},$2);
1655 :     push(@{$solutions},$2);
1656 :     }
1657 :     }
1658 :     }
1659 :     }
1660 :     $self->figmodel()->database()->save_table($tbl);
1661 :     $self->releaseModelLock();
1662 :    
1663 :     return $solutions;
1664 :     }
1665 :    
1666 :     =head3 datagapfill
1667 :     Definition:
1668 :     success()/fail() = FIGMODELmodel->datagapfill();
1669 :     Description:
1670 :     Run gapfilling on the input run specifications
1671 :     =cut
1672 :     sub datagapfill {
1673 :     my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;
1674 :     my $UniqueFilename = $self->figmodel()->filename();
1675 :     if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {
1676 :     system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id().$self->selected_version(),"NoBounds",["DataGapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0],"Gap filling runs" => join(";",@{$GapFillingRunSpecs})},"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,undef));
1677 :     #Checking that the solution exists
1678 :     if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {
1679 :     $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");
1680 :     $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
1681 :     return undef;
1682 :     }
1683 :     my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);
1684 :     if (defined($TansferFileSuffix)) {
1685 :     system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
1686 :     }
1687 :     #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted
1688 :     $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
1689 :     return $GapFillResultTable;
1690 :     }
1691 :     if (defined($TansferFileSuffix)) {
1692 :     $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
1693 :     }
1694 :     return undef;
1695 :     }
1696 :    
1697 :     =head3 TestSolutions
1698 :     Definition:
1699 :     $model->TestSolutions($ModelID,$NumProcessors,$ProcessorIndex,$GapFill);
1700 :     Description:
1701 :     Example:
1702 :     =cut
1703 :    
1704 :     sub TestSolutions {
1705 :     my ($self,$OriginalErrorFilename,$GapFillResultTable) = @_;
1706 :     #Getting the filename
1707 :     my $UniqueFilename = $self->figmodel()->filename();
1708 :     #Reading in the original error matrix which has the headings for the original model simulation
1709 :     my $OriginalErrorData;
1710 :     if (!defined($OriginalErrorFilename) || !-e $self->directory().$OriginalErrorFilename) {
1711 :     my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
1712 :     $OriginalErrorData = [$HeadingVector,$Errorvector];
1713 :     } else {
1714 :     $OriginalErrorData = $self->figmodel()->database()->load_single_column_file($self->directory().$OriginalErrorFilename,"");
1715 :     }
1716 :     my $HeadingHash;
1717 :     my @HeadingArray = split(/;/,$OriginalErrorData->[0]);
1718 :     my @OrigErrorArray = split(/;/,$OriginalErrorData->[1]);
1719 :     for (my $i=0; $i < @HeadingArray; $i++) {
1720 :     my @SubArray = split(/:/,$HeadingArray[$i]);
1721 :     $HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]} = $i;
1722 :     }
1723 :     #Scanning through the gap filling solutions
1724 :     my $TempVersion = "V".$UniqueFilename;
1725 :     my $ErrorMatrixLines;
1726 :     for (my $i=0; $i < $GapFillResultTable->size(); $i++) {
1727 :     print "Starting problem solving ".$i."\n";
1728 :     my $ErrorLine = $GapFillResultTable->get_row($i)->{"Experiment"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution index"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution cost"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0];
1729 :     #Integrating solution into test model
1730 :     my $ReactionArray;
1731 :     my $DirectionArray;
1732 :     my @ReactionList = split(/,/,$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0]);
1733 :     my %SolutionHash;
1734 :     for (my $k=0; $k < @ReactionList; $k++) {
1735 :     if ($ReactionList[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) {
1736 :     my $Reaction = $2;
1737 :     my $Sign = $1;
1738 :     if (defined($SolutionHash{$Reaction})) {
1739 :     $SolutionHash{$Reaction} = "<=>";
1740 :     } elsif ($Sign eq "-") {
1741 :     $SolutionHash{$Reaction} = "<=";
1742 :     } elsif ($Sign eq "+") {
1743 :     $SolutionHash{$Reaction} = "=>";
1744 :     } else {
1745 :     $SolutionHash{$Reaction} = $Sign;
1746 :     }
1747 :     }
1748 :     }
1749 :     @ReactionList = keys(%SolutionHash);
1750 :     for (my $k=0; $k < @ReactionList; $k++) {
1751 :     push(@{$ReactionArray},$ReactionList[$k]);
1752 :     push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});
1753 :     }
1754 :     print "Integrating solution!\n";
1755 :     $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1);
1756 :     $self->PrintModelLPFile();
1757 :     #Running the model against all available experimental data
1758 :     print "Running test model!\n";
1759 :     my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
1760 :    
1761 :     @HeadingArray = split(/;/,$HeadingVector);
1762 :     my @ErrorArray = @OrigErrorArray;
1763 :     my @TempArray = split(/;/,$Errorvector);
1764 :     for (my $j=0; $j < @HeadingArray; $j++) {
1765 :     my @SubArray = split(/:/,$HeadingArray[$j]);
1766 :     $ErrorArray[$HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]}] = $TempArray[$j];
1767 :     }
1768 :     $ErrorLine .= ";".$FalsePostives."/".$FalseNegatives.";".join(";",@ErrorArray);
1769 :     push(@{$ErrorMatrixLines},$ErrorLine);
1770 :     print "Finishing problem solving ".$i."\n";
1771 :     }
1772 :     #Clearing out the test model
1773 :     if (-e $self->directory().$self->id().$TempVersion.".txt") {
1774 :     unlink($self->directory().$self->id().$TempVersion.".txt");
1775 :     unlink($self->directory()."SimulationOutput".$self->id().$TempVersion.".txt");
1776 :     }
1777 :     return $ErrorMatrixLines;
1778 :     }
1779 :    
1780 :     =head3 status
1781 :     Definition:
1782 :     int::model status = FIGMODELmodel->status();
1783 :     Description:
1784 :     Returns the current status of the SEED model associated with the input genome ID.
1785 :     model status = 1: model exists
1786 :     model status = 0: model is being built
1787 :     model status = -1: model does not exist
1788 :     model status = -2: model build failed
1789 :     =cut
1790 :     sub status {
1791 :     my ($self) = @_;
1792 : chenry 1.19 return $self->{_data}->status();
1793 : chenry 1.18 }
1794 :    
1795 :     =head3 message
1796 :     Definition:
1797 :     string::model message = FIGMODELmodel->message();
1798 :     Description:
1799 :     Returns a message associated with the models current status
1800 :     =cut
1801 :     sub message {
1802 :     my ($self) = @_;
1803 : chenry 1.19 return $self->{_data}->message();
1804 : chenry 1.18 }
1805 :    
1806 :     =head3 set_status
1807 :     Definition:
1808 :     (success/fail) = FIGMODELmodel->set_status(int::new status,string::status message);
1809 :     Description:
1810 :     Changes the current status of the SEED model
1811 :     new status = 1: model exists
1812 :     new status = 0: model is being built
1813 :     new status = -1: model does not exist
1814 :     new status = -2: model build failed
1815 :     =cut
1816 :     sub set_status {
1817 :     my ($self,$NewStatus,$Message) = @_;
1818 : chenry 1.19 $self->{_data}->status($NewStatus);
1819 :     $self->{_data}->message($Message);
1820 : chenry 1.18 return $self->config("SUCCESS")->[0];
1821 :     }
1822 :    
1823 : chenry 1.29 =head3 generate_gpr_hash
1824 : chenry 1.18 Definition:
1825 : chenry 1.29 FIGMODELmodel->generate_gpr_hash();
1826 : chenry 1.18 Description:
1827 :     =cut
1828 :    
1829 : chenry 1.29 sub generate_gpr_hash {
1830 :     my ($self) = @_;
1831 : chenry 1.18 #Getting genome stats
1832 :     my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1833 :     my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
1834 :     if (!defined($FeatureTable)) {
1835 : chenry 1.29 $self->error_message("generate_gpr_hash:genome features could not be accessed!");
1836 :     return undef;
1837 : chenry 1.18 }
1838 :     #Checking that the number of genes exceeds the minimum size
1839 :     if ($FeatureTable->size() < $self->config("minimum genome size for modeling")->[0]) {
1840 : chenry 1.29 $self->error_message("generate_gpr_hash:genome rejected as too small for modeling!");
1841 :     return undef;
1842 : chenry 1.18 }
1843 :     #Sorting GenomeData by gene location on the chromosome
1844 : chenry 1.20 my $ftrTbl = $self->figmodel()->database()->get_table("ROLERXNMAPPING");
1845 : chenry 1.18 $FeatureTable->sort_rows("MIN LOCATION");
1846 :     my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);
1847 :     my %LastGenePosition;
1848 :     my $GeneRoles;
1849 :     for (my $j=0; $j < $FeatureTable->size(); $j++) {
1850 :     my $CurrentRow = $FeatureTable->get_row($j);
1851 :     #"ID","ALIASES","MIN LOCATION","MAX LOCATION","ROLES","SUBSYSTEMS","SUBSYSTEM CLASS"
1852 :     if (defined($CurrentRow)) {
1853 :     my $GeneID = $CurrentRow->{"ID"}->[0];
1854 :     if ($GeneID =~ m/(peg\.\d+)/) {
1855 :     $GeneID = $1;
1856 :     }
1857 :     foreach my $Role (@{$CurrentRow->{"ROLES"}}) {
1858 :     if ($self->figmodel()->role_is_valid($Role) != 0) {
1859 :     push(@{$GeneRoles->{$GeneID}},$Role);
1860 :     my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);
1861 :     if ($ReactionHashArrayRef != 0) {
1862 :     foreach my $Mapping (@{$ReactionHashArrayRef}) {
1863 :     if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
1864 :     if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) {
1865 :     if ($Mapping->{"MASTER"}->[0] eq 1) {
1866 :     #Creating a complex if consecutive genes have been assigned to the same reaction
1867 :     $ComplexHash->{$Mapping->{"REACTION"}->[0]}->{$Mapping->{"COMPLEX"}->[0]}->{$Role}->{$GeneID} = 1;
1868 :     if (!defined($LastGenePosition{$Mapping->{"REACTION"}->[0]})) {
1869 :     $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j;
1870 :     push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID);
1871 :     } elsif (($j-$LastGenePosition{$Mapping->{"REACTION"}->[0]}) < 3 && $LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) {
1872 :     my $CurrentComplex = pop(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}});
1873 :     push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$CurrentComplex."+".$GeneID);
1874 :     } elsif ($LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) {
1875 :     push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID);
1876 :     }
1877 :     $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j;
1878 :     #Adding a subsystem for the reaction
1879 :     if ($self->figmodel()->subsystem_is_valid($Mapping->{"SUBSYSTEM"}->[0]) == 1) {
1880 :     ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},$Mapping->{"SUBSYSTEM"}->[0]);
1881 :     if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 1) {
1882 :     if ($Mapping->{"SOURCE"}->[0] =~ m/Hope\sFiles/) {
1883 :     $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 1;
1884 :     } elsif ($Mapping->{"SOURCE"}->[0] =~ m/SEED/) {
1885 :     $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 2;
1886 :     } elsif (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) {
1887 :     $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3;
1888 :     }
1889 :     }
1890 :     }
1891 :     #Handling confidence
1892 :     if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) {
1893 :     if ($Mapping->{"SOURCE"}->[0] =~ m/MATT/) {
1894 :     $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3;
1895 :     } elsif ($Mapping->{"SOURCE"}->[0] =~ m/CHRIS/) {
1896 :     $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 4;
1897 :     } else {
1898 :     $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 5;
1899 :     }
1900 :     }
1901 :     #Parsing sources
1902 :     ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},split(/\|/,$Mapping->{"SOURCE"}->[0]));
1903 :     } else {
1904 :     push(@{$SuggestedMappings},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role);
1905 :     }
1906 :     } else {
1907 :     push(@{$UnrecognizedReactions},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role);
1908 :     }
1909 :     }
1910 :     }
1911 :     }
1912 :     }
1913 :     }
1914 :     }
1915 :     }
1916 :     #Creating nonadjacent complexes
1917 :     my @ReactionList = keys(%{$ReactionHash});
1918 :     foreach my $Reaction (@ReactionList) {
1919 :     #If multiple genes are assigned to the reaction, we check if they should should be in a complex
1920 :     if (@{$ReactionHash->{$Reaction}->{"GENES"}} > 0 && defined($ComplexHash->{$Reaction})) {
1921 :     my $GeneArray;
1922 :     foreach my $Complex (keys(%{$ComplexHash->{$Reaction}})) {
1923 :     my %ComplexComponents;
1924 :     foreach my $CurrentGeneSet (@{$ReactionHash->{$Reaction}->{"GENES"}}) {
1925 :     my @GeneList = split(/\+/,$CurrentGeneSet);
1926 :     my %RoleHash;
1927 :     foreach my $Gene (@GeneList) {
1928 :     foreach my $Role (@{$GeneRoles->{$Gene}}) {
1929 :     if (defined($ComplexHash->{$Reaction}->{$Complex}->{$Role})) {
1930 :     $RoleHash{$Role} = 1;
1931 :     }
1932 :     }
1933 :     }
1934 :     if (keys(%RoleHash) > 0) {
1935 :     if (!defined($ComplexComponents{join("|",sort(keys(%RoleHash)))})) {
1936 :     my @RoleList = keys(%RoleHash);
1937 :     my @ComplexList = keys(%ComplexComponents);
1938 :     foreach my $ComplexSet (@ComplexList) {
1939 :     my @RoleList = split(/\|/,$ComplexSet);
1940 :     my $Match = 0;
1941 :     foreach my $SingleRole (@RoleList) {
1942 :     if (defined($RoleHash{$SingleRole})) {
1943 :     $Match = 1;
1944 :     last;
1945 :     }
1946 :     }
1947 :     if ($Match == 1) {
1948 :     foreach my $SingleRole (@RoleList) {
1949 :     $RoleHash{$SingleRole} = 1
1950 :     }
1951 :     push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},@{$ComplexComponents{$ComplexSet}});
1952 :     delete $ComplexComponents{$ComplexSet};
1953 :     }
1954 :     }
1955 :     }
1956 :     push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},$CurrentGeneSet);
1957 :     }
1958 :     }
1959 :     my @Position;
1960 :     my @Options;
1961 :     my $Count = 0;
1962 :     foreach my $RoleSet (keys(%ComplexComponents)) {
1963 :     push(@Position,0);
1964 :     push(@{$Options[$Count]},@{$ComplexComponents{$RoleSet}});
1965 :     $Count++;
1966 :     }
1967 :     my $Done = 0;
1968 :     $Count = 0;
1969 :     my $NewRelationship;
1970 :     while($Done == 0) {
1971 :     #Creating complex with current indecies
1972 :     $NewRelationship->[$Count] = $Options[0]->[$Position[0]];
1973 :     for (my $i=1; $i < @Position; $i++) {
1974 :     $NewRelationship->[$Count] .= "+".$Options[$i]->[$Position[$i]];
1975 :     }
1976 :     $NewRelationship->[$Count] = join("+",$self->figmodel()->remove_duplicates(split(/\+/,$NewRelationship->[$Count])));
1977 :     $Count++;
1978 :     #Iterating indecies
1979 :     my $CurrentIndex = 0;
1980 :     while($CurrentIndex >= 0) {
1981 :     if ($CurrentIndex >= @Position) {
1982 :     $CurrentIndex = -1000;
1983 :     } elsif ($Position[$CurrentIndex]+1 == @{$Options[$CurrentIndex]}) {
1984 :     $Position[$CurrentIndex] = -1;
1985 :     $CurrentIndex++;
1986 :     } else {
1987 :     $Position[$CurrentIndex]++;
1988 :     $CurrentIndex--;
1989 :     }
1990 :     }
1991 :     if ($CurrentIndex == -1000) {
1992 :     $Done = 1;
1993 :     }
1994 :     }
1995 :     push(@{$GeneArray},@{$NewRelationship});
1996 :     }
1997 :     @{$ReactionHash->{$Reaction}->{"GENES"}} = $self->figmodel()->remove_duplicates(@{$GeneArray});
1998 :     }
1999 :     }
2000 : chenry 1.29 return $ReactionHash;
2001 :     }
2002 : chenry 1.18
2003 : chenry 1.29 =head3 CreateSingleGenomeReactionList
2004 :     Definition:
2005 :     FIGMODELmodel->CreateSingleGenomeReactionList();
2006 :     Description:
2007 :     This function uses fig calls to obtain a list of genes and functions for a genome, and it uses a file mapping reactions and functional roles to produce a reaction list.
2008 :     Example:
2009 :     =cut
2010 :    
2011 :     sub CreateSingleGenomeReactionList {
2012 :     my ($self,$RunGapFilling) = @_;
2013 :     #Creating directory
2014 :     if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/") {
2015 :     system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/");
2016 :     } elsif ($self->owner() eq "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->genome()."/") {
2017 :     system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->genome()."/");
2018 :     }
2019 :     if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/") {
2020 :     system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/");
2021 :     }
2022 :     #Checking that the number of genes exceeds the minimum size
2023 :     my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
2024 :     if ($FeatureTable->size() < $self->config("minimum genome size for modeling")->[0]) {
2025 :     $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome rejected as too small for modeling!");
2026 :     return $self->fail();
2027 :     }
2028 :     #Setting up needed variables
2029 :     my $OriginalModelTable = undef;
2030 :     if ($self->status() == 0) {
2031 :     $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList:Model is already being built. Canceling current build.");
2032 :     return $self->fail();
2033 :     }elsif ($self->status() == 1) {
2034 :     $OriginalModelTable = $self->reaction_table();
2035 :     $self->set_status(0,"Rebuilding preliminary reconstruction");
2036 :     } else {
2037 :     $self->set_status(0,"Preliminary reconstruction");
2038 :     }
2039 :     #Getting reaction hash
2040 :     my $ReactionHash = $self->generate_gpr_hash();
2041 :     if (!defined($ReactionHash)) {
2042 :     $self->figmodel()->error_message("CreateSingleGenomeReactionList:unable to create reaction hash!");
2043 :     return $self->fail();
2044 :     }
2045 : chenry 1.18 #Getting the reaction table
2046 :     my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");
2047 :     #Creating the model reaction table
2048 :     my $NewModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n");
2049 : chenry 1.29 my @ReactionList = keys(%{$ReactionHash});
2050 : chenry 1.18 foreach my $ReactionID (@ReactionList) {
2051 :     #Getting the thermodynamic reversibility from the database
2052 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
2053 :     my $Subsystem = "NONE";
2054 :     if (defined($ReactionHash->{$ReactionID}->{"SUBSYSTEMS"})) {
2055 :     $Subsystem = join("|",@{$ReactionHash->{$ReactionID}->{"SUBSYSTEMS"}});
2056 :     }
2057 :     my $Source = "NONE";
2058 :     if (defined($ReactionHash->{$ReactionID}->{"SOURCE"})) {
2059 :     $Source = join("|",@{$ReactionHash->{$ReactionID}->{"SOURCE"}});
2060 :     }
2061 :     $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => [join("|",@{$ReactionHash->{$ReactionID}->{"GENES"}})],"SUBSYSTEM" => [$Subsystem],"CONFIDENCE" => [$ReactionHash->{$ReactionID}->{"CONFIDENCE"}->[0]],"REFERENCE" => [$Source],"NOTES" => ["NONE"]});
2062 :     }
2063 : chenry 1.20 #Getting feature rows for features that are lumped
2064 : chenry 1.29 my $ftrTbl = $self->figmodel()->database()->get_table("ROLERXNMAPPING");
2065 : chenry 1.20 my @rows = $ftrTbl->get_rows_by_key("2","MASTER");
2066 :     for (my $i=0; $i < @rows; $i++) {
2067 :     my $rxn = $rows[$i]->{REACTION}->[0];
2068 :     my $role = $rows[$i]->{ROLE}->[0];
2069 :     my @orgrows = $FeatureTable->get_rows_by_key($role,"ROLES");
2070 :     my $genes;
2071 :     for (my $j=0; $j < @orgrows; $j++) {
2072 :     if ($orgrows[$j]->{ID}->[0] =~ m/(peg\.\d+)/) {
2073 :     push(@{$genes},$1);
2074 :     }
2075 :     }
2076 :     if (defined($genes) && @{$genes} > 0) {
2077 :     my $row = $NewModelTable->get_row_by_key($rxn,"LOAD");
2078 :     my $newGeneAssociations;
2079 :     for (my $k=0; $k < @{$genes}; $k++) {
2080 :     for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
2081 :     my $peg = $genes->[$k];
2082 :     if ($row->{"ASSOCIATED PEG"}->[$j] !~ m/$peg/) {
2083 :     push(@{$newGeneAssociations},$row->{"ASSOCIATED PEG"}->[$j]."+".$peg);
2084 :     }
2085 :     }
2086 :     }
2087 :     $row->{"ASSOCIATED PEG"} = $newGeneAssociations;
2088 :     }
2089 :     }
2090 : chenry 1.18 #Adding the spontaneous and universal reactions
2091 :     foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
2092 :     #Getting the thermodynamic reversibility from the database
2093 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
2094 :     #Checking if the reaction is already in the model
2095 :     if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) {
2096 :     $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]});
2097 :     }
2098 :     }
2099 :     foreach my $ReactionID (@{$self->config("universal reactions")}) {
2100 :     #Getting the thermodynamic reversibility from the database
2101 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
2102 :     #Checking if the reaction is already in the model
2103 :     if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) {
2104 :     $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]});
2105 :     }
2106 :     }
2107 : chenry 1.27 #Creating biomass reaction for model
2108 :     my $biomassID = $self->BuildSpecificBiomassReaction();
2109 :     if ($biomassID !~ m/bio\d\d\d\d\d/) {
2110 : chenry 1.18 $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
2111 : chenry 1.27 die $self->error_message("CreateSingleGenomeReactionList: Could not generate biomass reaction!");
2112 :     }
2113 :     #Getting the biomass reaction PPO object
2114 :     my $bioRxn = $self->figmodel()->database()->get_object("bof",{id=>$biomassID});
2115 :     if (!defined($bioRxn)) {
2116 :     die $self->error_message("CreateSingleGenomeReactionList: Could not find biomass reaction ".$biomassID."!");
2117 :     }
2118 :     #Getting the list of essential reactions for biomass reaction
2119 :     my $ReactionList;
2120 :     my $essentialReactions = $bioRxn->essentialRxn();
2121 :     if (defined($essentialReactions) && $essentialReactions =~ m/rxn\d\d\d\d\d/) {
2122 :     push(@{$ReactionList},split(/\|/,$essentialReactions));
2123 :     if ($essentialReactions !~ m/$biomassID/) {
2124 :     push(@{$ReactionList},$biomassID);
2125 :     }
2126 : chenry 1.30 } else {
2127 :     push(@{$ReactionList},$biomassID);
2128 : chenry 1.18 }
2129 :     #Adding biomass reactions to the model table
2130 :     foreach my $BOFReaction (@{$ReactionList}) {
2131 :     #Getting the thermodynamic reversibility from the database
2132 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($BOFReaction);
2133 :     #Checking if the reaction is already in the model
2134 :     if (!defined($NewModelTable->get_row_by_key($BOFReaction,"LOAD"))) {
2135 :     if ($BOFReaction =~ m/bio/) {
2136 :     $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["BOF"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [1],"REFERENCE" => ["Biomass objective function"],"NOTES" => ["NONE"]});
2137 :     } else {
2138 :     $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["INITIAL GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Initial gap filling"],"NOTES" => ["NONE"]});
2139 :     }
2140 :     }
2141 :     }
2142 :     #Completing any incomplete reactions sets
2143 :     my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS");
2144 :     for (my $i=0; $i < $ReactionSetTable->size(); $i++) {
2145 :     if (defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
2146 :     foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) {
2147 :     if (!defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
2148 :     #Getting the thermodynamic reversibility from the database
2149 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction);
2150 :     $NewModelTable->add_row({"LOAD" => [$Reaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["REACTION SET GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Added due to presence of ".$ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0]],"NOTES" => ["NONE"]});
2151 :     }
2152 :     }
2153 :     }
2154 :     }
2155 : chenry 1.20 #If an original model exists, we copy the gap filling solution from that model
2156 :     if (defined($OriginalModelTable)) {
2157 :     for (my $i=0; $i < $OriginalModelTable->size(); $i++) {
2158 :     if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ && $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] ne "INITIAL GAP FILLING") {
2159 :     my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD");
2160 :     if (!defined($Row)) {
2161 :     $NewModelTable->add_row($OriginalModelTable->get_row($i));
2162 :     }
2163 :     }
2164 :     }
2165 :     }
2166 : chenry 1.18 #Now we compare the model to the previous model to determine if any differences exist that aren't gap filling reactions
2167 :     if (defined($OriginalModelTable)) {
2168 :     my $PerfectMatch = 1;
2169 :     my $ReactionCount = 0;
2170 :     for (my $i=0; $i < $OriginalModelTable->size(); $i++) {
2171 :     #We only check that nongapfilling reactions exist in the new model
2172 :     if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] !~ m/GAP/ || $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "INITIAL GAP FILLING") {
2173 :     $ReactionCount++;
2174 :     my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD");
2175 :     if (defined($Row)) {
2176 :     #We check that the reaction directionality is identical
2177 :     if ($Row->{"DIRECTIONALITY"}->[0] ne $OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]) {
2178 :     if (defined($OriginalModelTable->get_row($i)->{"NOTES"}->[0]) && $OriginalModelTable->get_row($i)->{"NOTES"}->[0] =~ m/Directionality\sswitched\sfrom\s([^\s])/) {
2179 :     if ($1 ne $Row->{"DIRECTIONALITY"}->[0]) {
2180 :     print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$1." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n";
2181 :     $PerfectMatch = 0;
2182 :     last;
2183 :     }
2184 :     } else {
2185 :     print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n";
2186 :     $PerfectMatch = 0;
2187 :     last;
2188 :     }
2189 :     }
2190 :     #We check that the genes assigned to the reaction are identical
2191 :     if ($PerfectMatch == 1 && @{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}} != @{$Row->{"ASSOCIATED PEG"}}) {
2192 :     print "Gene associatation mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}}." vs ".@{$Row->{"ASSOCIATED PEG"}}."\n";
2193 :     $PerfectMatch = 0;
2194 :     last;
2195 :     }
2196 :     if ($PerfectMatch == 1) {
2197 :     my @GeneSetOne = sort(@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}});
2198 :     my @GeneSetTwo = sort(@{$Row->{"ASSOCIATED PEG"}});
2199 :     for (my $j=0; $j < @GeneSetOne; $j++) {
2200 :     if ($GeneSetOne[$j] ne $GeneSetTwo[$j]) {
2201 :     print "Gene mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$GeneSetOne[$j]." vs ".$GeneSetTwo[$j]."\n";
2202 :     $PerfectMatch = 0;
2203 :     $i = $OriginalModelTable->size();
2204 :     last;
2205 :     }
2206 :     }
2207 :     }
2208 :     } else {
2209 :     print "Original model contains an extra reaction:".$OriginalModelTable->get_row($i)->{"LOAD"}->[0]."\n";
2210 :     $PerfectMatch = 0;
2211 :     last;
2212 :     }
2213 :     }
2214 :     }
2215 :     if ($PerfectMatch == 1 && $ReactionCount == $NewModelTable->size()) {
2216 :     #Bailing out of function as the model has not changed
2217 :     $self->set_status(1,"rebuild canceled because model has not changed");
2218 :     return $self->success();
2219 :     }
2220 :     }
2221 :     #Saving the new model to file
2222 :     $self->set_status(1,"Preliminary reconstruction complete");
2223 :     $self->figmodel()->database()->save_table($NewModelTable);
2224 : chenry 1.19 $self->reaction_table(1);
2225 : chenry 1.18 #Updating the model stats table
2226 :     $self->update_stats_for_build();
2227 :     $self->PrintSBMLFile();
2228 : chenry 1.28 $self->PrintModelLPFile();
2229 :     $self->PrintModelSimpleReactionTable();
2230 : chenry 1.19 if (defined($OriginalModelTable)) {
2231 :     $self->calculate_model_changes($OriginalModelTable,"REBUILD");
2232 :     }
2233 : chenry 1.18 #Adding model to gapfilling queue
2234 :     if (defined($RunGapFilling) && $RunGapFilling == 1) {
2235 :     $self->set_status(1,"Autocompletion queued");
2236 :     $self->figmodel()->add_job_to_queue("gapfillmodel?".$self->id(),"QSUB","cplex",$self->owner(),"BACK");
2237 :     }
2238 :     return $self->success();
2239 :     }
2240 :    
2241 :     =head3 CreateMetaGenomeReactionList
2242 :     Definition:
2243 :     (success/fail) = FIGMODELmodel->CreateMetaGenomeReactionList();
2244 :     Description:
2245 :     This is the code called to create or update the reaction list for a metgenome model
2246 :     =cut
2247 :    
2248 :     sub CreateMetaGenomeReactionList {
2249 :     my ($self) = @_;
2250 :     #Checking if the metagenome file exists
2251 :     if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {
2252 :     $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
2253 : chenry 1.21 return $self->fail();
2254 : chenry 1.18 }
2255 :     #Loading metagenome data
2256 :     my $MGRASTData = $self->figmodel()->database()->load_multiple_column_file($self->config("raw MGRAST directory")->[0].$self->genome().".summary","\t");
2257 :     if (!defined($MGRASTData)) {
2258 :     $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
2259 : chenry 1.21 return $self->fail();
2260 : chenry 1.18 }
2261 :     #Setting up needed variables
2262 :     my $OriginalModelTable = undef;
2263 :     #Checking status
2264 :     if ($self->status() < 0) {
2265 :     $self->set_status(0,"Preliminary reconstruction");
2266 :     } elsif ($self->status() == 0) {
2267 :     $self->error_message("FIGMODEL->CreateModelReactionList:Model is already being built. Canceling current build.");
2268 :     return $self->fail();
2269 :     } else {
2270 :     $OriginalModelTable = $self->reaction_table();
2271 :     $self->ArchiveModel();
2272 :     $self->set_status(0,"Rebuilding preliminary reconstruction");
2273 :     }
2274 : chenry 1.21 #Creating a hash of escores and pegs associated with each role
2275 :     my $rolePegHash;
2276 :     my $roleEscores;
2277 :     for (my $i=0; $i < @{$MGRASTData};$i++) {
2278 :     #MD5,PEG,number of sims,role,sim e-scores,max escore,min escore,ave escore,stdev escore,ave exponent,stddev exponent
2279 :     $rolePegHash->{$MGRASTData->[$i]->[3]}->{substr($MGRASTData->[$i]->[1],4)} = 1;
2280 :     push(@{$roleEscores->{$MGRASTData->[$i]->[3]}},split(/;/,$MGRASTData->[$i]->[4]));
2281 :     }
2282 : chenry 1.18 #Getting the reaction table
2283 : chenry 1.21 my $ReactionTable = $self->figmodel()->database()->get_table("REACTIONS");
2284 : chenry 1.18 #Creating model table
2285 : chenry 1.21 my $ModelTable = $self->create_table_prototype("ModelReactions");
2286 :     print $ModelTable->filename();
2287 :     my @roles = keys(%{$rolePegHash});
2288 :     for (my $i=0; $i < @roles; $i++) {
2289 :     my $min = -1;
2290 :     my $max = -1;
2291 :     my $count = @{$roleEscores->{$roles[$i]}};
2292 :     my $ave = 0;
2293 :     my $stdev = 0;
2294 :     my $aveexp = 0;
2295 :     my $stdevexp = 0;
2296 :     for (my $j=0; $j < @{$roleEscores->{$roles[$i]}}; $j++) {
2297 :     if ($roleEscores->{$roles[$i]} < $min || $min == -1) {
2298 :     $min = $roleEscores->{$roles[$i]};
2299 :     }
2300 :     if ($roleEscores->{$roles[$i]} > $max || $max == -1) {
2301 :     $max = $roleEscores->{$roles[$i]};
2302 :     }
2303 :     $ave += $roleEscores->{$roles[$i]}->[$j];
2304 :     if ($roleEscores->{$roles[$i]}->[$j] =~ m/e(-\d+$)/) {
2305 :     $aveexp += $1;
2306 :     }
2307 :     }
2308 :     $ave = $ave/$count;
2309 :     $aveexp = $aveexp/$count;
2310 :     for (my $j=0; $j < @{$roleEscores->{$roles[$i]}}; $j++) {
2311 :     $stdev += ($roleEscores->{$roles[$i]}->[$j]-$ave)*($roleEscores->{$roles[$i]}->[$j]-$ave);
2312 :     if ($roleEscores->{$roles[$i]}->[$j] =~ m/e(-\d+$)/) {
2313 :     $stdevexp += ($1-$aveexp)*($1-$aveexp);
2314 :     }
2315 :     }
2316 :     $stdev = sqrt($stdev/$count);
2317 :     $stdevexp = sqrt($stdevexp/$count);
2318 : chenry 1.18 #Checking for subsystems
2319 : chenry 1.21 my $GeneSubsystems = $self->figmodel()->subsystems_of_role($roles[$i]);
2320 : chenry 1.18 #Checking if there are reactions associated with this role
2321 : chenry 1.21 my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($roles[$i]);
2322 : chenry 1.18 if ($ReactionHashArrayRef != 0) {
2323 :     foreach my $Mapping (@{$ReactionHashArrayRef}) {
2324 :     if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
2325 :     if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) {
2326 :     if ($Mapping->{"MASTER"}->[0] eq 1) {
2327 :     #Checking if the reaction is already in the model
2328 :     my $ReactionRow = $ModelTable->get_row_by_key($Mapping->{"REACTION"}->[0],"LOAD");
2329 :     if (!defined($ReactionRow)) {
2330 :     $ReactionRow = {"LOAD" => [$Mapping->{"REACTION"}->[0]],"DIRECTIONALITY" => [$self->figmodel()->reversibility_of_reaction($Mapping->{"REACTION"}->[0])],"COMPARTMENT" => ["c"]};
2331 :     $ModelTable->add_row($ReactionRow);
2332 :     }
2333 : chenry 1.21 my %pegHash = %{$rolePegHash->{$roles[$i]}};
2334 :     if (defined($ReactionRow->{"ASSOCIATED PEG"})) {
2335 :     for (my $j=0; $j < @{$ReactionRow->{"ASSOCIATED PEG"}}; $j++) {
2336 :     $pegHash{$ReactionRow->{"ASSOCIATED PEG"}->[$j]} = 1;
2337 :     }
2338 :     }
2339 :     delete $ReactionRow->{"ASSOCIATED PEG"};
2340 :     push(@{$ReactionRow->{"ASSOCIATED PEG"}},keys(%pegHash));
2341 :     push(@{$ReactionRow->{"REFERENCE"}},$count.":".$ave.":".$stdev.":".$aveexp.":".$stdevexp.":".$min.":".$max);
2342 : chenry 1.18 if (defined($GeneSubsystems)) {
2343 :     push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});
2344 :     }
2345 :     }
2346 :     }
2347 :     }
2348 :     }
2349 :     }
2350 :     }
2351 :    
2352 :     #Adding the spontaneous and universal reactions
2353 :     foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
2354 :     #Getting the thermodynamic reversibility from the database
2355 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
2356 :     #Checking if the reaction is already in the model
2357 :     if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) {
2358 :     $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]});
2359 :     }
2360 :     }
2361 :     foreach my $ReactionID (@{$self->config("universal reactions")}) {
2362 :     #Getting the thermodynamic reversibility from the database
2363 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
2364 :     #Checking if the reaction is already in the model
2365 :     if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) {
2366 :     $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]});
2367 :     }
2368 :     }
2369 :    
2370 :     #Completing any incomplete reactions sets
2371 :     my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS");
2372 :     for (my $i=0; $i < $ReactionSetTable->size(); $i++) {
2373 :     if (defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
2374 :     foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) {
2375 :     if (!defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
2376 :     #Getting the thermodynamic reversibility from the database
2377 :     my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction);
2378 :     $ModelTable->add_row({"LOAD" => [$Reaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["REACTION SET GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Added due to presence of ".$ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0]],"NOTES" => ["NONE"]});
2379 :     }
2380 :     }
2381 :     }
2382 :     }
2383 :    
2384 :     #Clearing the previous model from the cache
2385 :     $self->figmodel()->database()->ClearDBModel($self->id(),1);
2386 :     $ModelTable->save();
2387 :    
2388 :     return $self->success();
2389 :     }
2390 :    
2391 :     =head3 ArchiveModel
2392 :     Definition:
2393 :     (success/fail) = FIGMODELmodel->ArchiveModel();
2394 :     Description:
2395 :     This function archives the specified model in the model directory with the current version numbers appended.
2396 :     This function is used to preserve old versions of models prior to overwriting so new versions may be compared with old versions.
2397 :     =cut
2398 :     sub ArchiveModel {
2399 :     my ($self) = @_;
2400 :    
2401 :     #Checking that the model file exists
2402 :     if (!(-e $self->filename())) {
2403 :     $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");
2404 :     return $self->fail();
2405 :     }
2406 :    
2407 :     #Copying the model file
2408 :     system("cp ".$self->filename()." ".$self->directory().$self->id().$self->version().".txt");
2409 :     }
2410 :    
2411 :     =head3 PrintModelDataToFile
2412 :     Definition:
2413 :     (success/fail) = FIGMODELmodel->PrintModelDataToFile();
2414 :     Description:
2415 :     This function uses the MFAToolkit to print out all of the compound and reaction data for the input model.
2416 :     Some of the data printed by the toolkit is calculated internally in the toolkit and not stored in any files, so this data can only be retrieved through this
2417 :     function. The LoadModel function for example would not give you this data.
2418 :     =cut
2419 :     sub PrintModelDataToFile {
2420 :     my($self) = @_;
2421 :    
2422 :     #Running the MFAToolkit on the model file
2423 :     my $OutputIndex = $self->figmodel()->filename();
2424 :     my $Command = $self->config("MFAToolkit executable")->[0]." parameterfile ../Parameters/Printing.txt resetparameter output_folder ".$OutputIndex.'/ LoadCentralSystem "'.$self->filename().'"';
2425 :     system($Command);
2426 :    
2427 :     #Copying the model file printed by the toolkit out of the output directory and into the model directory
2428 :     if (!-e $self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().".txt") {
2429 :     $self->figmodel()->error_message("New model file not created due to an error. Check that the input modelfile exists.");
2430 :     $self->figmodel()->cleardirectory($OutputIndex);
2431 :     return $self->fail();
2432 :     }
2433 :    
2434 :     $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().'.txt" "'.$self->directory().$self->id().$self->selected_version().'Data.txt"';
2435 :     system($Command);
2436 :     $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex.'/ErrorLog0.txt" "'.$self->directory().'ModelErrors.txt"';
2437 :     system($Command);
2438 :     $self->figmodel()->cleardirectory($OutputIndex);
2439 :     return $self->success();
2440 :     }
2441 :    
2442 :     =head2 Analysis Functions
2443 :    
2444 :     =head3 run_microarray_analysis
2445 :     Definition:
2446 :     int::status = FIGMODEL->run_microarray_analysis(string::media,string::job id,string::gene calls);
2447 :     Description:
2448 :     Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2449 :     =cut
2450 :     sub run_microarray_analysis {
2451 :     my ($self,$media,$label,$index,$genecall) = @_;
2452 :     $genecall =~ s/_/:/g;
2453 :     $genecall =~ s/\//;/g;
2454 :     my $uniqueFilename = $self->figmodel()->filename();
2455 :     my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($uniqueFilename,$self->id(),$media,["ProductionMFA","ShewenellaExperiment"],{"Microarray assertions" => $label.";".$index.";".$genecall,"MFASolver" => "CPLEX","Network output location" => "/scratch/"},"MicroarrayAnalysis-".$uniqueFilename.".txt",undef,$self->selected_version());
2456 :     system($command);
2457 :     my $filename = $self->figmodel()->config("MFAToolkit output directory")->[0].$uniqueFilename."/MicroarrayOutput-".$index.".txt";
2458 : chenry 1.20 if (-e $filename) {
2459 :     my $output = $self->figmodel()->database()->load_single_column_file($filename);
2460 :     if (defined($output->[0])) {
2461 :     my @array = split(/;/,$output->[0]);
2462 :     $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2463 :     return ($array[0],$array[1],$array[8].":".$array[2],$array[9].":".$array[3],$array[10].":".$array[4],$array[11].":".$array[5],$array[12].":".$array[6],$array[13].":".$array[7]);
2464 :     }
2465 :     print STDERR $filename." is empty!";
2466 :     }
2467 :     print STDERR $filename." not found!";
2468 :     $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2469 :    
2470 : chenry 1.18 return undef;
2471 :     }
2472 :    
2473 :     =head3 find_minimal_pathways
2474 :     Definition:
2475 :     int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
2476 :     Description:
2477 :     Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2478 :     =cut
2479 :     sub find_minimal_pathways {
2480 :     my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2481 :    
2482 :     #Setting default media
2483 :     if (!defined($media)) {
2484 :     $media = "Complete";
2485 :     }
2486 :    
2487 :     #Setting default solution number
2488 :     if (!defined($solutionnum)) {
2489 :     $solutionnum = "5";
2490 :     }
2491 :    
2492 :     #Setting additional exchange fluxes
2493 :     if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2494 :     if ($self->id() eq "iAF1260") {
2495 :     $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2496 :     } else {
2497 :     $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2498 :     }
2499 :     }
2500 :    
2501 :     #Translating objective
2502 :     my $objectivestring;
2503 :     if ($objective eq "ALL") {
2504 :     #Getting the list of universal building blocks
2505 :     my $buildingblocks = $self->config("universal building blocks");
2506 :     my @objectives = keys(%{$buildingblocks});
2507 :     #Getting the nonuniversal building blocks
2508 :     my $otherbuildingblocks = $self->config("nonuniversal building blocks");
2509 :     my @array = keys(%{$otherbuildingblocks});
2510 :     if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
2511 :     my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
2512 :     if (defined($equation)) {
2513 :     for (my $i=0; $i < @array; $i++) {
2514 :     if (CORE::index($equation,$array[$i]) > 0) {
2515 :     push(@objectives,$array[$i]);
2516 :     }
2517 :     }
2518 :     }
2519 :     }
2520 :     for (my $i=0; $i < @objectives; $i++) {
2521 :     $self->find_minimal_pathways($media,$objectives[$i]);
2522 :     }
2523 :     return;
2524 :     } elsif ($objective eq "ENERGY") {
2525 :     $objectivestring = "MAX;FLUX;rxn00062;c;1";
2526 :     } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2527 :     if ($objective =~ m/\[(\w)\]/) {
2528 :     $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2529 :     $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2530 :     } else {
2531 :     $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2532 :     $additionalexchange .= ";".$objective."[c]:-100:0";
2533 :     }
2534 :     } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2535 :     my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2536 :     for (my $i=0; $i < @{$Products};$i++) {
2537 :     my $temp = $Products->[$i]->{"DATABASE"}->[0];
2538 :     if ($additionalexchange !~ m/$temp/) {
2539 :     #$additionalexchange .= ";".$temp."[c]:-100:0";
2540 :     }
2541 :     }
2542 :     for (my $i=0; $i < @{$Reactants};$i++) {
2543 :     print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2544 :     $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2545 :     print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2546 :     }
2547 :     return;
2548 :     }
2549 :    
2550 :     #Adding additional drains
2551 :     if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
2552 :     $additionalexchange .= ";cpd15666[c]:0:100";
2553 :     } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2554 :     $additionalexchange .= ";cpd12370[c]:0:100";
2555 :     } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2556 :     $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2557 :     }
2558 :    
2559 :     #Running MFAToolkit
2560 :     my $filename = $self->figmodel()->filename();
2561 :     my $command;
2562 :     if (defined($AllReversible) && $AllReversible == 1) {
2563 :     $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"Make all reactions reversible in MFA"=>1, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "CPLEX","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2564 :     } else {
2565 :     $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"Make all reactions reversible in MFA"=>0, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "CPLEX","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2566 :     }
2567 :     system($command);
2568 :    
2569 :     #Loading problem report
2570 :     my $results = $self->figmodel()->LoadProblemReport($filename);
2571 :     #Clearing output
2572 :     $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2573 :     if (!defined($results)) {
2574 :     print STDERR $objective." pathway results not found!\n";
2575 :     return;
2576 :     }
2577 :    
2578 :     #Parsing output
2579 :     my @Array;
2580 :     my $row = $results->get_row(1);
2581 :     if (defined($row->{"Notes"}->[0])) {
2582 :     $_ = $row->{"Notes"}->[0];
2583 :     @Array = /\d+:([^\|]+)\|/g;
2584 :     }
2585 :    
2586 :     #Writing output to file
2587 :     $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
2588 :     }
2589 :    
2590 :     =head3 find_minimal_pathways
2591 :     Definition:
2592 :     int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
2593 :     Description:
2594 :     Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2595 :     =cut
2596 :     sub find_minimal_pathways_two {
2597 :     my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2598 :    
2599 :     #Setting default media
2600 :     if (!defined($media)) {
2601 :     $media = "Complete";
2602 :     }
2603 :    
2604 :     #Setting default solution number
2605 :     if (!defined($solutionnum)) {
2606 :     $solutionnum = "5";
2607 :     }
2608 :    
2609 :     #Setting additional exchange fluxes
2610 :     if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2611 :     if ($self->id() eq "iAF1260") {
2612 :     $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2613 :     } else {
2614 :     $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2615 :     }
2616 :     }
2617 :    
2618 :     #Translating objective
2619 :     my $objectivestring;
2620 :     if ($objective eq "ALL") {
2621 :     #Getting the list of universal building blocks
2622 :     my $buildingblocks = $self->config("universal building blocks");
2623 :     my @objectives = keys(%{$buildingblocks});
2624 :     #Getting the nonuniversal building blocks
2625 :     my $otherbuildingblocks = $self->config("nonuniversal building blocks");
2626 :     my @array = keys(%{$otherbuildingblocks});
2627 :     if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
2628 :     my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
2629 :     if (defined($equation)) {
2630 :     for (my $i=0; $i < @array; $i++) {
2631 :     if (CORE::index($equation,$array[$i]) > 0) {
2632 :     push(@objectives,$array[$i]);
2633 :     }
2634 :     }
2635 :     }
2636 :     }
2637 :     for (my $i=0; $i < @objectives; $i++) {
2638 :     $self->find_minimal_pathways($media,$objectives[$i]);
2639 :     }
2640 :     return;
2641 :     } elsif ($objective eq "ENERGY") {
2642 :     $objectivestring = "MAX;FLUX;rxn00062;c;1";
2643 :     } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2644 :     if ($objective =~ m/\[(\w)\]/) {
2645 :     $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2646 :     $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2647 :     } else {
2648 :     $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2649 :     $additionalexchange .= ";".$objective."[c]:-100:0";
2650 :     }
2651 :     } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2652 :     my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2653 :     for (my $i=0; $i < @{$Products};$i++) {
2654 :     my $temp = $Products->[$i]->{"DATABASE"}->[0];
2655 :     if ($additionalexchange !~ m/$temp/) {
2656 :     #$additionalexchange .= ";".$temp."[c]:-100:0";
2657 :     }
2658 :     }
2659 :     for (my $i=0; $i < @{$Reactants};$i++) {
2660 :     print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2661 :     $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2662 :     print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2663 :     }
2664 :     return;
2665 :     }
2666 :    
2667 :     #Adding additional drains
2668 :     if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
2669 :     $additionalexchange .= ";cpd15666[c]:0:100";
2670 :     } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2671 :     $additionalexchange .= ";cpd12370[c]:0:100";
2672 :     } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2673 :     $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2674 :     }
2675 :    
2676 :     #Running MFAToolkit
2677 :     my $filename = $self->figmodel()->filename();
2678 :     my $command;
2679 :     if (defined($AllReversible) && $AllReversible == 1) {
2680 :     $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"use simple variable and constraint names"=>1,"Make all reactions reversible in MFA"=>1, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "SCIP","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2681 :     } else {
2682 :     $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"use simple variable and constraint names"=>1,"Make all reactions reversible in MFA"=>0, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "SCIP","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());
2683 :     }
2684 :     print $command."\n";
2685 :     system($command);
2686 :    
2687 :     #Loading problem report
2688 :     my $results = $self->figmodel()->LoadProblemReport($filename);
2689 :     #Clearing output
2690 :     $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2691 :     if (!defined($results)) {
2692 :     print STDERR $objective." pathway results not found!\n";
2693 :     return;
2694 :     }
2695 :    
2696 :     #Parsing output
2697 :     my @Array;
2698 :     my $row = $results->get_row(1);
2699 :     if (defined($row->{"Notes"}->[0])) {
2700 :     $_ = $row->{"Notes"}->[0];
2701 :     @Array = /\d+:([^\|]+)\|/g;
2702 :     }
2703 :    
2704 :     #Writing output to file
2705 :     $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
2706 :     }
2707 :    
2708 :     sub combine_minimal_pathways {
2709 :     my ($self) = @_;
2710 :    
2711 :     my $tbl;
2712 :     if (-e $self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl") {
2713 :     $tbl = FIGMODELTable::load_table($self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",";","|",0,["Objective","Media","Reversible"]);
2714 :     } else {
2715 :     $tbl = FIGMODELTable->new(["Objective","Media","Reactions","Reversible","Shortest path","Number of essentials","Essentials","Length"],$self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",["Objective","Media","Reversible"],";","|");
2716 :     }
2717 :     my @files = glob($self->directory()."MinimalPathways-*");
2718 :     for (my $i=0; $i < @files;$i++) {
2719 :     if ($files[$i] =~ m/MinimalPathways\-(\S+)\-(cpd\d\d\d\d\d)\-(\w+)\-(\d)\-/ || $files[$i] =~ m/MinimalPathways\-(\S+)\-(ENERGY)\-(\w+)\-(\d)\-/) {
2720 :     my $reactions = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2721 :     if (defined($reactions) && @{$reactions} > 0 && length($reactions->[0]) > 0) {
2722 :     my $newrow = {"Objective"=>[$2],"Media"=>[$1],"Reversible"=>[$4]};
2723 :     my $row = $tbl->get_table_by_key($newrow->{"Objective"}->[0],"Objective")->get_table_by_key($newrow->{"Media"}->[0],"Media")->get_row_by_key($newrow->{"Reversible"}->[0],"Reversible");
2724 :     if (!defined($row)) {
2725 :     $row = $tbl->add_row($newrow);
2726 :     }
2727 :     $row->{Reactions} = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2728 :     delete($row->{"Shortest path"});
2729 :     delete($row->{"Number of essentials"});
2730 :     delete($row->{"Essentials"});
2731 :     delete($row->{"Length"});
2732 :     for (my $j=0; $j < @{$row->{Reactions}}; $j++) {
2733 :     my @array = split(/,/,$row->{Reactions}->[$j]);
2734 :     $row->{"Length"}->[$j] = @array;
2735 :     if (!defined($row->{"Shortest path"}->[0]) || $row->{"Length"}->[$j] < $row->{"Shortest path"}->[0]) {
2736 :     $row->{"Shortest path"}->[0] = $row->{"Length"}->[$j];
2737 :     }
2738 :     $row->{"Number of essentials"}->[0] = 0;
2739 :     for (my $k=0; $k < @array;$k++) {
2740 :     if ($array[$k] =~ m/(rxn\d\d\d\d\d)/) {
2741 :     my $class = $self->get_reaction_class($1,1);
2742 :     my $temp = $row->{Media}->[0].":Essential";
2743 :     if ($class =~ m/$temp/) {
2744 :     $row->{"Number of essentials"}->[$j]++;
2745 :     if (!defined($row->{"Essentials"}->[$j]) && length($row->{"Essentials"}->[$j]) > 0) {
2746 :     $row->{"Essentials"}->[$j] = $array[$k];
2747 :     } else {
2748 :     $row->{"Essentials"}->[$j] .= ",".$array[$k];
2749 :     }
2750 :     }
2751 :     }
2752 :     }
2753 :     }
2754 :     }
2755 :     }
2756 :     }
2757 :     $tbl->save();
2758 :     }
2759 :    
2760 :     =head3 calculate_growth
2761 :     Definition:
2762 :     string::growth = FIGMODELmodel->calculate_growth(string:media);
2763 :     Description:
2764 :     Calculating growth in the input media
2765 :     =cut
2766 :     sub calculate_growth {
2767 : chenry 1.19 my ($self,$Media,$outputDirectory,$InParameters,$saveLPFile) = @_;
2768 :     #Setting the Media
2769 :     if (!defined($Media) || length($Media) == 0) {
2770 :     $Media = $self->autocompleteMedia();
2771 :     }
2772 :     #Setting parameters for the run
2773 :     my $DefaultParameters = $self->figmodel()->defaultParameters();
2774 :     if (defined($InParameters)) {
2775 :     my @parameters = keys(%{$InParameters});
2776 :     for (my $i=0; $i < @parameters; $i++) {
2777 :     $DefaultParameters->{$parameters[$i]} = $InParameters->{$parameters[$i]};
2778 :     }
2779 :     }
2780 :     $DefaultParameters->{"optimize metabolite production if objective is zero"} = 1;
2781 :     #Setting filenames
2782 : chenry 1.18 my $UniqueFilename = $self->figmodel()->filename();
2783 : chenry 1.19 if (!defined($outputDirectory)) {
2784 :     $outputDirectory = $self->config("database message file directory")->[0];
2785 :     }
2786 :     my $fluxFilename = $outputDirectory."Fluxes-".$self->id()."-".$Media.".txt";
2787 :     my $cpdFluxFilename = $outputDirectory."CompoundFluxes-".$self->id()."-".$Media.".txt";
2788 :     #Running FBA
2789 : chenry 1.20 #print $self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],$DefaultParameters,$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version())."\n";
2790 : chenry 1.19 system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],$DefaultParameters,$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version()));
2791 :     #Saving LP file if requested
2792 :     if (defined($saveLPFile) && $saveLPFile == 1 && -e $self->figmodel()->{"MFAToolkit output directory"}->[0].$UniqueFilename."/CurrentProblem.lp") {
2793 : chenry 1.20 system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/CurrentProblem.lp ".$self->directory().$self->id().".lp");
2794 : chenry 1.19 }
2795 : chenry 1.18 my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);
2796 :     my $Result;
2797 :     if (defined($ProblemReport)) {
2798 :     my $Row = $ProblemReport->get_row(0);
2799 :     if (defined($Row) && defined($Row->{"Objective"}->[0])) {
2800 : chenry 1.19 if ($Row->{"Objective"}->[0] < 0.00000001 || $Row->{"Objective"}->[0] == 1e7) {
2801 :     $Result = "NOGROWTH";
2802 :     if (defined($Row->{"Individual metabolites with zero production"}->[0]) && $Row->{"Individual metabolites with zero production"}->[0] =~ m/cpd\d\d\d\d\d/) {
2803 :     $Result .= ":".$Row->{"Individual metabolites with zero production"}->[0];
2804 :     }
2805 : chenry 1.29 $self->growth(0);
2806 :     $self->noGrowthCompounds($Row->{"Individual metabolites with zero production"}->[0]);
2807 : chenry 1.18 } else {
2808 : chenry 1.19 if (-e $self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionReactionData0.txt") {
2809 :     system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionReactionData0.txt ".$fluxFilename);
2810 :     system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionCompoundData0.txt ".$cpdFluxFilename);
2811 :     }
2812 : chenry 1.18 $Result = $Row->{"Objective"}->[0];
2813 : chenry 1.29 $self->growth($Result);
2814 :     $self->noGrowthCompounds("NONE");
2815 : chenry 1.18 }
2816 :     }
2817 :     }
2818 : chenry 1.19 #Deleting files if necessary
2819 :     if ($self->figmodel()->config("preserve all log files")->[0] ne "yes") {
2820 :     $self->figmodel()->cleardirectory($UniqueFilename);
2821 :     unlink($self->figmodel()->config("database message file directory")->[0].$self->id()."-".$Media."-GrowthTest.txt");
2822 :     }
2823 :     #Returning result
2824 : chenry 1.18 return $Result;
2825 :     }
2826 :    
2827 :     =head3 classify_model_reactions
2828 :     Definition:
2829 :     (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);
2830 :     Description:
2831 :     This function uses the MFAToolkit to minimize and maximize the flux through every reaction in the input model during minimal growth on the input media.
2832 :     The results are returned in a hash of strings where the keys are the reaction IDs and the strings are structured as follows: "Class;Min flux;Max flux".
2833 :     Possible values for "Class" include:
2834 :     1.) Positive: these reactions are essential in the forward direction.
2835 :     2.) Negative: these reactions are essential in the reverse direction.
2836 :     3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.
2837 :     4.) Negative variable: these reactions are nonessential, but they only ever proceed in the reverse direction.
2838 :     5.) Variable: these reactions are nonessential and proceed in the forward or reverse direction.
2839 :     6.) Blocked: these reactions never carry any flux at all in the media condition tested.
2840 :     7.) Dead: these reactions are disconnected from the network.
2841 :     =cut
2842 :     sub classify_model_reactions {
2843 :     my ($self,$Media,$SaveChanges) = @_;
2844 :    
2845 :     #Getting unique file for printing model output
2846 :     my $UniqueFilename = $self->figmodel()->filename();
2847 :     #Running the MFAToolkit
2848 :     system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],{"identify dead ends" => 1,"find tight bounds" => 1,"MFASolver" => "GLPK"},"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
2849 :     #Reading in the output bounds file
2850 :     my ($ReactionTB,$CompoundTB,$DeadCompounds,$DeadEndCompounds,$DeadReactions);
2851 :     if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") {
2852 :     $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]);
2853 :     }
2854 :     if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt") {
2855 :     $CompoundTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt",";","|",1,["DATABASE ID"]);
2856 :     }
2857 :     if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt") {
2858 :     $DeadReactions = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt",""));
2859 :     }
2860 :     if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt") {
2861 :     $DeadCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt",""));
2862 :     }
2863 :     if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt") {
2864 :     $DeadEndCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt",""));
2865 :     }
2866 :     if (!defined($ReactionTB) && !defined($CompoundTB)) {
2867 :     print STDERR "FIGMODEL:ClassifyModelReactions: Classification file not found when classifying reactions in ".$self->id().$self->selected_version()." with ".$Media." media. Most likely the model did not grow.\n";
2868 :     return (undef,undef);
2869 :     }
2870 :    
2871 :     #Clearing output
2872 :     $self->figmodel()->clearing_output($UniqueFilename,"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
2873 :     #Creating the table objects that will hold the results of the reaction classification
2874 :     my $rxnclasstable = $self->reaction_class_table();
2875 :     my $cpdclasstable = $self->compound_class_table();
2876 :     #Loading the compound table
2877 :     if (defined($CompoundTB)) {
2878 :     for (my $i=0; $i < $CompoundTB->size(); $i++) {
2879 :     my $Row = $CompoundTB->get_row($i);
2880 :     if (defined($Row->{"DATABASE ID"})) {
2881 :     #Getting the compound row
2882 :     my $CpdRow = $cpdclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0].$Row->{COMPARTMENT}->[0],"COMPOUND",1);
2883 :     #Setting row values
2884 :     my $Max = 0;
2885 :     my $Min = 0;
2886 :     my $Class = "Unknown";
2887 :     if (defined($DeadCompounds) && defined($DeadCompounds->{$Row->{"DATABASE ID"}->[0]})) {
2888 :     $Class = "Dead";
2889 :     } elsif (defined($DeadEndCompounds) && defined($DeadEndCompounds->{$Row->{"DATABASE ID"}->[0]})) {
2890 :     $Class = "Deadend";
2891 :     } elsif (defined($Row->{"Min DRAIN_FLUX"}) && defined($Row->{"Max DRAIN_FLUX"}) && $Row->{"Min DRAIN_FLUX"}->[0] ne "1e+07") {
2892 :     $Max = $Row->{"Max DRAIN_FLUX"}->[0];
2893 :     $Min = $Row->{"Min DRAIN_FLUX"}->[0];
2894 :     if ($Row->{"Min DRAIN_FLUX"}->[0] > 0.00000001) {
2895 :     $Class = "Positive";
2896 :     } elsif ($Row->{"Max DRAIN_FLUX"}->[0] < -0.00000001) {
2897 :     $Class = "Negative";
2898 :     } elsif ($Row->{"Min DRAIN_FLUX"}->[0] < -0.00000001) {
2899 :     if ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) {
2900 :     $Class = "Variable";
2901 :     } else {
2902 :     $Max = 0;
2903 :     $Class = "Negative variable";
2904 :     }
2905 :     } elsif ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) {
2906 :     $Min = 0;
2907 :     $Class = "Positive variable";
2908 :     } else {
2909 :     $Min = 0;
2910 :     $Max = 0;
2911 :     $Class = "Blocked";
2912 :     }
2913 :     }
2914 :     my $index = 0;
2915 :     if (defined($CpdRow->{MEDIA})) {
2916 :     for (my $i=0; $i < @{$CpdRow->{MEDIA}};$i++) {
2917 :     $index++;
2918 :     if ($CpdRow->{MEDIA}->[$i] eq $Media) {
2919 :     $index = $i;
2920 :     last;
2921 :     }
2922 :     }
2923 :     }
2924 :     $CpdRow->{MIN}->[$index] = $Min;
2925 :     $CpdRow->{MAX}->[$index] = $Max;
2926 :     $CpdRow->{CLASS}->[$index] = $Class;
2927 :     $CpdRow->{MEDIA}->[$index] = $Media;
2928 :     }
2929 :     }
2930 :     if (!defined($SaveChanges) || $SaveChanges == 1) {
2931 :     $cpdclasstable->save();
2932 :     }
2933 :     }
2934 :     if (defined($ReactionTB)) {
2935 :     for (my $i=0; $i < $ReactionTB->size(); $i++) {
2936 :     my $Row = $ReactionTB->get_row($i);
2937 :     if (defined($Row->{"DATABASE ID"})) {
2938 :     #Getting the compound row
2939 :     my $Compartment = "c";
2940 :     if (defined($Row->{COMPARTMENT}->[0])) {
2941 :     $Compartment = $Row->{COMPARTMENT}->[0];
2942 :     }
2943 :     my $RxnRow = $rxnclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0],"REACTION",1);
2944 :     my $Max = 0;
2945 :     my $Min = 0;
2946 :     my $Class = "Unknown";
2947 :     if (defined($DeadReactions) && defined($DeadReactions->{$Row->{"DATABASE ID"}->[0]})) {
2948 :     $Class = "Dead";
2949 :     } elsif (defined($Row->{"Min FLUX"}) && defined($Row->{"Max FLUX"})) {
2950 :     $Max = $Row->{"Max FLUX"}->[0];
2951 :     $Min = $Row->{"Min FLUX"}->[0];
2952 :     if ($Row->{"Min FLUX"}->[0] > 0.00000001) {
2953 :     $Class = "Positive";
2954 :     } elsif ($Row->{"Max FLUX"}->[0] < -0.00000001) {
2955 :     $Class = "Negative";
2956 :     } elsif ($Row->{"Min FLUX"}->[0] < -0.00000001) {
2957 :     if ($Row->{"Max FLUX"}->[0] > 0.00000001) {
2958 :     $Class = "Variable";
2959 :     } else {
2960 :     $Max = 0;
2961 :     $Class = "Negative variable";
2962 :     }
2963 :     } elsif ($Row->{"Max FLUX"}->[0] > 0.00000001) {
2964 :     $Min = 0;
2965 :     $Class = "Positive variable";
2966 :     } else {
2967 :     $Min = 0;
2968 :     $Max = 0;
2969 :     $Class = "Blocked";
2970 :     }
2971 :     }
2972 :     my $index = 0;
2973 :     if (defined($RxnRow->{MEDIA})) {
2974 :     for (my $i=0; $i < @{$RxnRow->{MEDIA}};$i++) {
2975 :     $index++;
2976 :     if ($RxnRow->{MEDIA}->[$i] eq $Media) {
2977 :     $index = $i;
2978 :     last;
2979 :     }
2980 :     }
2981 :     }
2982 :     $RxnRow->{MIN}->[$index] = $Min;
2983 :     $RxnRow->{MAX}->[$index] = $Max;
2984 :     $RxnRow->{CLASS}->[$index] = $Class;
2985 :     $RxnRow->{MEDIA}->[$index] = $Media;
2986 :     }
2987 :     }
2988 :     if (!defined($SaveChanges) || $SaveChanges == 1) {
2989 :     $rxnclasstable->save();
2990 :     }
2991 :     }
2992 :     return ($rxnclasstable,$cpdclasstable);
2993 :     }
2994 :    
2995 :     =head3 RunAllStudiesWithDataFast
2996 :     Definition:
2997 :     (integer::false positives,integer::false negatives,integer::correct negatives,integer::correct positives,string::error vector,string heading vector) = FIGMODELmodel->RunAllStudiesWithDataFast(string::experiment,0/1::print result);
2998 :     Description:
2999 :     Simulates every experimental condition currently available for the model.
3000 :     =cut
3001 :    
3002 :     sub RunAllStudiesWithDataFast {
3003 :     my ($self,$Experiment,$PrintResults) = @_;
3004 :    
3005 :     #Printing lp and key file for model
3006 :     if (!-e $self->directory()."FBA-".$self->id().$self->selected_version().".lp") {
3007 :     $self->PrintModelLPFile();
3008 :     }
3009 :     my $UniqueFilename = $self->figmodel()->filename();
3010 :    
3011 :     #Determing the simulations that need to be run
3012 :     my $ExperimentalDataTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment);
3013 :     #Creating the table of jobs to submit
3014 :     my $JobArray = $self->GetSimulationJobTable($ExperimentalDataTable,$Experiment,$UniqueFilename);
3015 :     #Printing the job file
3016 :     if (!-d $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/") {
3017 :     system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/");
3018 :     }
3019 :     $JobArray->save();
3020 :    
3021 :     #Running simulations
3022 :     system($self->config("mfalite executable")->[0]." ".$self->config("Reaction database directory")->[0]."masterfiles/MediaTable.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Jobfile.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
3023 :     #Parsing the results
3024 :     my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
3025 :     if (!defined($Results)) {
3026 :     $self->figmodel()->error_message("FIGMODELmodel:RunAllStudiesWithDataFast:Could not find simulation results: ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
3027 :     return undef;
3028 :     }
3029 :     my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector,$SimulationResults) = $self->EvaluateSimulationResults($Results,$ExperimentalDataTable);
3030 :     #Printing results to file
3031 :     $self->figmodel()->database()->save_table($SimulationResults,undef,undef,undef,"False negatives\tFalse positives\tCorrect negatives\tCorrect positives\n".$FalseNegatives."\t".$FalsePostives."\t".$CorrectNegatives."\t".$CorrectPositives."\n");
3032 :     $self->figmodel()->clearing_output($UniqueFilename);
3033 :    
3034 :     return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector);
3035 :     }
3036 :    
3037 :     =head3 GetSimulationJobTable
3038 :     Definition:
3039 :     my $JobTable = $model->GetSimulationJobTable($Experiment,$PrintResults,$Version);
3040 :     Description:
3041 :     =cut
3042 :    
3043 :     sub GetSimulationJobTable {
3044 :     my ($self,$SimulationTable,$Experiment,$Folder) = @_;
3045 :    
3046 :     #Determing the simulations that need to be run
3047 :     if (!defined($SimulationTable)) {
3048 :     $SimulationTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment);
3049 :     if (!defined($SimulationTable)) {
3050 :     return undef;
3051 :     }
3052 :     }
3053 :    
3054 :     #Creating the job table
3055 :     my $JobTable = $self->figmodel()->CreateJobTable($Folder);
3056 :     for (my $i=0; $i < $SimulationTable->size(); $i++) {
3057 :     if ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Gene\sKO/) {
3058 :     my $Row = $JobTable->get_row_by_key("Gene KO","LABEL",1);
3059 :     $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
3060 :     } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Media\sgrowth/) {
3061 :     my $Row = $JobTable->get_row_by_key("Growth phenotype","LABEL",1);
3062 :     $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
3063 :     } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Interval\sKO/) {
3064 :     my $Row = $JobTable->get_row_by_key($SimulationTable->get_row($i)->{"Heading"}->[0],"LABEL",1);
3065 :     $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
3066 :     $JobTable->add_data($Row,"GENE KO",$SimulationTable->get_row($i)->{"Experiment type"}->[0],1);
3067 :     }
3068 :     }
3069 :    
3070 :     #Filling in model specific elements of the job table
3071 :     for (my $i=0; $i < $JobTable->size(); $i++) {
3072 :     if ($JobTable->get_row($i)->{"LABEL"}->[0] =~ m/Gene\sKO/) {
3073 :     $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "SINGLEKO";
3074 :     $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 1;
3075 :     } else {
3076 :     $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "GROWTH";
3077 :     $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 0;
3078 :     }
3079 :     $JobTable->get_row($i)->{"LP FILE"}->[0] = $self->directory()."FBA-".$self->id().$self->selected_version();
3080 :     $JobTable->get_row($i)->{"MODEL"}->[0] = $self->directory().$self->id().$self->selected_version().".txt";
3081 :     $JobTable->get_row($i)->{"SAVE FLUXES"}->[0] = 0;
3082 :     }
3083 :    
3084 :     return $JobTable;
3085 :     }
3086 :    
3087 :     =head3 EvaluateSimulationResults
3088 :     Definition:
3089 :     (integer::false positives,integer::false negatives,integer::correct negatives,integer::correct positives,string::error vector,string heading vector,FIGMODELtable::simulation results) = FIGMODELmodel->EvaluateSimulationResults(FIGMODELtable::raw simulation results,FIGMODELtable::experimental data);
3090 :     Description:
3091 :     Compares simulation results with experimental data to produce a table indicating where predictions are incorrect.
3092 :     =cut
3093 :    
3094 :     sub EvaluateSimulationResults {
3095 :     my ($self,$Results,$ExperimentalDataTable) = @_;
3096 :    
3097 :     #Comparing experimental results with simulation results
3098 :     my $SimulationResults = FIGMODELTable->new(["Run result","Experiment type","Media","Experiment ID","Reactions knocked out"],$self->directory()."SimulationOutput".$self->id().$self->selected_version().".txt",["Experiment ID","Media"],"\t",",",undef);
3099 :     my $FalsePostives = 0;
3100 :     my $FalseNegatives = 0;
3101 :     my $CorrectNegatives = 0;
3102 :     my $CorrectPositives = 0;
3103 :     my @Errorvector;
3104 :     my @HeadingVector;
3105 :     my $ReactionKOWithGeneHash;
3106 :     for (my $i=0; $i < $Results->size(); $i++) {
3107 :     if ($Results->get_row($i)->{"LABEL"}->[0] eq "Gene KO") {
3108 :     if (defined($Results->get_row($i)->{"REACTION KO WITH GENES"})) {
3109 :     for (my $j=0; $j < @{$Results->get_row($i)->{"REACTION KO WITH GENES"}}; $j++) {
3110 :     my @Temp = split(/:/,$Results->get_row($i)->{"REACTION KO WITH GENES"}->[$j]);
3111 :     if (defined($Temp[1]) && length($Temp[1]) > 0) {
3112 :     $ReactionKOWithGeneHash->{$Temp[0]} = $Temp[1];
3113 :     }
3114 :     }
3115 :     }
3116 :     if ($Results->get_row($i)->{"OBJECTIVE"}->[0] == 0) {
3117 :     for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) {
3118 :     my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading");
3119 :     if (defined($Row)) {
3120 :     my $KOReactions = "none";
3121 :     if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) {
3122 :     $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]};
3123 :     }
3124 :     push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
3125 :     my $Status = "Unknown";
3126 :     if ($Row->{"Growth"}->[0] > 0) {
3127 :     $Status = "False negative";
3128 :     $FalseNegatives++;
3129 :     push(@Errorvector,3);
3130 :     } else {
3131 :     $Status = "False positive";
3132 :     $FalsePostives++;
3133 :     push(@Errorvector,2);
3134 :     }
3135 :     $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
3136 :     }
3137 :     }
3138 :     } else {
3139 :     for (my $j=0; $j < @{$Results->get_row($i)->{"ESSENTIALGENES"}}; $j++) {
3140 :     #print $j."\t".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]."\n";
3141 :     my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j],"Heading");
3142 :     if (defined($Row)) {
3143 :     my $KOReactions = "none";
3144 :     if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]})) {
3145 :     $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]};
3146 :     }
3147 :     push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
3148 :     my $Status = "Unknown";
3149 :     if ($Row->{"Growth"}->[0] > 0) {
3150 :     $Status = "False negative";
3151 :     $FalseNegatives++;
3152 :     push(@Errorvector,3);
3153 :     } else {
3154 :     $Status = "Correct negative";
3155 :     $CorrectNegatives++;
3156 :     push(@Errorvector,1);
3157 :     }
3158 :     $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
3159 :     }
3160 :     }
3161 :     for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) {
3162 :     my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading");
3163 :     if (defined($Row)) {
3164 :     my $KOReactions = "none";
3165 :     if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) {
3166 :     $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]};
3167 :     }
3168 :     push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
3169 :     my $Status = "Unknown";
3170 :     if ($Row->{"Growth"}->[0] > 0) {
3171 :     $Status = "Correct positive";
3172 :     $CorrectPositives++;
3173 :     push(@Errorvector,0);
3174 :     } else {
3175 :     $Status = "False positive";
3176 :     $FalsePostives++;
3177 :     push(@Errorvector,2);
3178 :     }
3179 :     $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
3180 :     }
3181 :     }
3182 :     }
3183 :     } elsif ($Results->get_row($i)->{"LABEL"}->[0] eq "Growth phenotype") {
3184 :     my $Row = $ExperimentalDataTable->get_row_by_key("Media growth:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"MEDIA"}->[0],"Heading");
3185 :     if (defined($Row)) {
3186 :     push(@HeadingVector,$Row->{"Heading"}->[0].":none");
3187 :     my $Status = "Unknown";
3188 :     if ($Row->{"Growth"}->[0] > 0) {
3189 :     if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
3190 :     $Status = "Correct positive";
3191 :     $CorrectPositives++;
3192 :     push(@Errorvector,0);
3193 :     } else {
3194 :     $Status = "False negative";
3195 :     $FalseNegatives++;
3196 :     push(@Errorvector,3);
3197 :     }
3198 :     } else {
3199 :     if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
3200 :     $Status = "False positive";
3201 :     $FalsePostives++;
3202 :     push(@Errorvector,2);
3203 :     } else {
3204 :     $Status = "Correct negative";
3205 :     $CorrectNegatives++;
3206 :     push(@Errorvector,1);
3207 :     }
3208 :     }
3209 :     $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Media growth"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Media"}->[0]],"Reactions knocked out" => ["none"]});
3210 :     }
3211 :     } elsif ($Results->get_row($i)->{"LABEL"}->[0] =~ m/Interval\sKO/ && defined($Results->get_row($i)->{"KOGENES"}->[0])) {
3212 :     my $Row = $ExperimentalDataTable->get_row_by_key($Results->get_row($i)->{"LABEL"}->[0],"Heading");
3213 :     if (defined($Row)) {
3214 :     my $Status = "Unknown";
3215 :     if ($Row->{"Growth"}->[0] > 0) {
3216 :     if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
3217 :     $Status = "Correct positive";
3218 :     $CorrectPositives++;
3219 :     push(@Errorvector,0);
3220 :     } else {
3221 :     $Status = "False negative";
3222 :     $FalseNegatives++;
3223 :     push(@Errorvector,3);
3224 :     }
3225 :     } else {
3226 :     if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
3227 :     $Status = "False positive";
3228 :     $FalsePostives++;
3229 :     push(@Errorvector,2);
3230 :     } else {
3231 :     $Status = "Correct negative";
3232 :     $CorrectNegatives++;
3233 :     push(@Errorvector,1);
3234 :     }
3235 :     }
3236 :     $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Interval KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => ["none"]});
3237 :     }
3238 :     }
3239 :     }
3240 :    
3241 :     return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,join(";",@Errorvector),join(";",@HeadingVector),$SimulationResults);
3242 :     }
3243 :    
3244 :     =head3 InspectSolution
3245 :     Definition:
3246 :     $model->InspectSolution(string::gene knocked out,string::media condition,[string]::list of reactions);
3247 :     Description:
3248 :     =cut
3249 :    
3250 :     sub InspectSolution {
3251 :     my ($self,$GeneKO,$Media,$ReactionList) = @_;
3252 :    
3253 :     #Getting a directory for the results
3254 :     my $UniqueFilename = $self->figmodel()->filename();
3255 :     system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/");
3256 :     my $TempVersion = "V".$UniqueFilename;
3257 :    
3258 :     #Setting gene ko to none if no genes are to be knocked out
3259 :     if ($GeneKO !~ m/^peg\./) {
3260 :     $GeneKO = "none";
3261 :     }
3262 :    
3263 :     #Implementing the input solution in the test model
3264 :     my $ReactionArray;
3265 :     my $DirectionArray;
3266 :     my %SolutionHash;
3267 :     for (my $k=0; $k < @{$ReactionList}; $k++) {
3268 :     if ($ReactionList->[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) {
3269 :     my $Reaction = $2;
3270 :     my $Sign = $1;
3271 :     if (defined($SolutionHash{$Reaction})) {
3272 :     $SolutionHash{$Reaction} = "<=>";
3273 :     } elsif ($Sign eq "-") {
3274 :     $SolutionHash{$Reaction} = "<=";
3275 :     } elsif ($Sign eq "+") {
3276 :     $SolutionHash{$Reaction} = "=>";
3277 :     } else {
3278 :     $SolutionHash{$Reaction} = $Sign;
3279 :     }
3280 :     }
3281 :     }
3282 :     my @TempList = keys(%SolutionHash);
3283 :     for (my $k=0; $k < @TempList; $k++) {
3284 :     push(@{$ReactionArray},$TempList[$k]);
3285 :     push(@{$DirectionArray},$SolutionHash{$TempList[$k]});
3286 :     }
3287 :    
3288 :     print "Integrating solution!\n";
3289 :     $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"SolutionInspection",1,1);
3290 :    
3291 :     #Printing lp and key file for model
3292 :     $self->PrintModelLPFile();
3293 :    
3294 :     #Running FBA on the test model
3295 :     my $JobTable = $self->figmodel()->CreateJobTable($UniqueFilename);
3296 :     $JobTable->add_row({"LABEL" => ["TEST"],"RUNTYPE" => ["GROWTH"],"LP FILE" => [$self->directory()."FBA-".$self->id().$TempVersion],"MODEL" => [$self->directory().$self->id().$TempVersion.".txt"],"MEDIA" => [$Media],"REACTION KO" => ["none|".join("|",@{$ReactionList})],"GENE KO" => [$GeneKO],"SAVE FLUXES" => [0],"SAVE NONESSENTIALS" => [0]});
3297 :     $JobTable->save();
3298 :    
3299 :     #Running simulations
3300 :     system($self->config("mfalite executable")->[0]." ".$self->config("Reaction database directory")->[0]."masterfiles/MediaTable.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Jobfile.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
3301 :    
3302 :     #Parsing the results
3303 :     my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
3304 :     if (!defined($Results)) {
3305 :     $self->figmodel()->error_message("FIGMODELmodel:InspectSolution:Could not load problem report ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
3306 :     return undef;
3307 :     }
3308 :    
3309 :     #Making sure that the model grew with all reactions present
3310 :     my $Found = 0;
3311 :     for (my $i=0; $i < $Results->size(); $i++) {
3312 :     if (defined($Results->get_row($i)->{"KOGENES"}->[0]) && defined($Results->get_row($i)->{"KOREACTIONS"}->[0]) && $Results->get_row($i)->{"KOREACTIONS"}->[0] eq "none" && $Results->get_row($i)->{"KOGENES"}->[0] eq $GeneKO && $Results->get_row($i)->{"OBJECTIVE"}->[0] > 0.00001) {
3313 :     $Found = 1;
3314 :     }
3315 :     }
3316 :     if ($Found == 0) {
3317 :     print "Solution no longer valid\n";
3318 :     return undef;
3319 :     }
3320 :    
3321 :     #Making sure all of the reactions added are still necessary
3322 :     my $FinalReactionList;
3323 :     for (my $k=0; $k < $Results->size(); $k++) {
3324 :     if (defined($Results->get_row($k)->{"KOGENES"}->[0]) && $Results->get_row($k)->{"KOGENES"}->[0] eq $GeneKO) {
3325 :     if (defined($Results->get_row($k)->{"KOREACTIONS"}->[0]) && $Results->get_row($k)->{"KOREACTIONS"}->[0] =~ m/rxn\d\d\d\d\d/ && $Results->get_row($k)->{"OBJECTIVE"}->[0] < 0.000001) {
3326 :     push(@{$FinalReactionList},$Results->get_row($k)->{"KOREACTIONS"}->[0]);
3327 :     }
3328 :     }
3329 :     }
3330 :    
3331 :     #Deleting extra files created
3332 :     unlink($self->directory()."FBA-".$self->id().$TempVersion.".lp");
3333 :     unlink($self->directory()."FBA-".$self->id().$TempVersion.".key");
3334 :     unlink($self->directory().$self->id().$TempVersion.".txt");
3335 :    
3336 :     #Deleting the test model and the MFA folder
3337 :     $self->figmodel()->clearing_output($UniqueFilename);
3338 :    
3339 :     return $FinalReactionList;
3340 :     }
3341 :    
3342 :     =head3 GapFillingAlgorithm
3343 :    
3344 :     Definition:
3345 :     FIGMODELmodel->GapFillingAlgorithm();
3346 :    
3347 :     Description:
3348 :     This is a wrapper for running the gap filling algorithm on any model in the database.
3349 :     The algorithm performs a gap filling for any false negative prediction of the avialable experimental data.
3350 :     This function is threaded to improve efficiency: one thread does nothing but using the MFAToolkit to fill gaps for every false negative prediction.
3351 :     The other thread reads in the gap filling solutions, builds a test model for each solution, and runs the test model against all available experimental data.
3352 :     This function prints two important output files in the Model directory:
3353 :     1.) GapFillingOutput.txt: this is a summary of the results of the gap filling analysis
3354 :     2.) GapFillingErrorMatrix.txt: this lists the correct and incorrect predictions for each gapfilling solution implemented in a test model.
3355 :     =cut
3356 :    
3357 :     sub GapFillingAlgorithm {
3358 :     my ($self) = @_;
3359 :    
3360 :     #First the input model version and model filename should be simulated and the false negatives identified
3361 :     my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
3362 :    
3363 :     #Getting the filename
3364 :     my $UniqueFilename = $self->figmodel()->filename();
3365 :    
3366 :     #Printing the original performance vector
3367 :     $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-OPEM".".txt",[$HeadingVector,$Errorvector]);
3368 :    
3369 :     my $PreviousGapFilling;
3370 :     if (-e $self->directory().$self->id().$self->selected_version()."-GFS.txt") {
3371 :     #Backing up the old solution file
3372 :     system("cp ".$self->directory().$self->id().$self->selected_version()."-GFS.txt ".$self->directory().$self->id().$self->selected_version()."-OldGFS.txt");
3373 :     unlink($self->directory().$self->id().$self->selected_version()."-GFS.txt");
3374 :     }
3375 :     if (-e $self->directory().$self->id().$self->selected_version()."-OldGFS.txt") {
3376 :     #Reading in the solution file from the previous gap filling if it exists
3377 :     $PreviousGapFilling = $self->figmodel()->database()->load_table($self->directory().$self->id().$self->selected_version()."-OldGFS.txt",";",",",0,["Experiment"]);
3378 :     }
3379 :    
3380 :     #Now we use the simulation output to make the gap filling run data
3381 :     my @Errors = split(/;/,$Errorvector);
3382 :     my @Headings = split(/;/,$HeadingVector);
3383 :     my $GapFillingRunSpecs = "";
3384 :     my $Count = 0;
3385 :     my $RescuedPreviousResults;
3386 :     my $RunCount = 0;
3387 :     my $SolutionExistedCount = 0;
3388 :     my $AcceptedSolutions = 0;
3389 :     my $RejectedSolutions = 0;
3390 :     my $NoExistingSolutions = 0;
3391 :     for (my $i=0; $i < @Errors; $i++) {
3392 :     if ($Errors[$i] == 3) {
3393 :     my @HeadingDataArray = split(/:/,$Headings[$i]);
3394 :     if ($HeadingDataArray[2] !~ m/^peg\./ || $HeadingDataArray[3] ne "none") {
3395 :     my $SolutionFound = 0;
3396 :     if (defined($PreviousGapFilling) && defined($PreviousGapFilling->get_row_by_key($HeadingDataArray[2],"Experiment"))) {
3397 :     my @Rows = $PreviousGapFilling->get_rows_by_key($HeadingDataArray[2],"Experiment");
3398 :     for (my $j=0; $j < @Rows; $j++) {
3399 :     if ($HeadingDataArray[2] =~ m/^peg\./) {
3400 :     my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"});
3401 :     if (defined($ReactionList)) {
3402 :     print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n";
3403 :     $SolutionFound++;
3404 :     push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList}));
3405 :     $AcceptedSolutions++;
3406 :     } else {
3407 :     $RejectedSolutions++;
3408 :     }
3409 :     } else {
3410 :     my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"});
3411 :     if (defined($ReactionList)) {
3412 :     print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n";
3413 :     $SolutionFound++;
3414 :     push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList}));
3415 :     $AcceptedSolutions++;
3416 :     } else {
3417 :     $RejectedSolutions++;
3418 :     }
3419 :     }
3420 :     }
3421 :     } else {
3422 :     $NoExistingSolutions++;
3423 :     }
3424 :     if ($SolutionFound == 0) {
3425 :     $RunCount++;
3426 :     if (length($GapFillingRunSpecs) > 0) {
3427 :     $GapFillingRunSpecs .= ";";
3428 :     }
3429 :     $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].":".$HeadingDataArray[3];
3430 :     } else {
3431 :     $SolutionExistedCount++;
3432 :     }
3433 :     }
3434 :     $Count++;
3435 :     }
3436 :     }
3437 :    
3438 :     #Updating the growmatch progress table
3439 :     my $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1);
3440 :     $Row->{"INITIAL FP"}->[0] = $FalsePostives;
3441 :     $Row->{"INITIAL FN"}->[0] = $FalseNegatives;
3442 :     $Row->{"GF TIMING"}->[0] = time()."-";
3443 :     $Row->{"FN WITH SOL"}->[0] = $FalseNegatives-$NoExistingSolutions;
3444 :     $Row->{"FN WITH ACCEPTED SOL"}->[0] = $SolutionExistedCount;
3445 :     $Row->{"TOTAL ACCEPTED GF SOL"}->[0] = $AcceptedSolutions;
3446 :     $Row->{"TOTAL REJECTED GF SOL"}->[0] = $RejectedSolutions;
3447 :     $Row->{"FN WITH NO SOL"}->[0] = $NoExistingSolutions+$RejectedSolutions;
3448 :     $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM");
3449 :    
3450 :     #Running the gap filling once to correct all false negative errors
3451 :     my $SolutionsFound = 0;
3452 :     my $GapFillingArray;
3453 :     push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));
3454 :     my $GapFillingResults = $self->datagapfill($GapFillingArray,"GFS");
3455 :     if (defined($GapFillingResults)) {
3456 :     $SolutionsFound = 1;
3457 :     }
3458 :    
3459 :     if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {
3460 :     #Printing previous solutions to GFS file
3461 :     $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",$RescuedPreviousResults,1);
3462 :     $SolutionsFound = 1;
3463 :     }
3464 :    
3465 :     #Recording the finishing of the gapfilling
3466 :     $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1);
3467 :     $Row->{"GF TIMING"}->[0] .= time();
3468 :     $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM");
3469 :    
3470 :     if ($SolutionsFound == 1) {
3471 :     #Scheduling solution testing
3472 :     $self->figmodel()->add_job_to_queue("testsolutions?".$self->id().$self->selected_version()."?-1?GF","QSUB","fast","master","BACK");
3473 :     } else {
3474 :     $self->figmodel()->error_message("No false negative predictions found. Data gap filling not necessary!");
3475 :     }
3476 :    
3477 :     return $self->success();
3478 :     }
3479 :    
3480 :     =head3 SolutionReconciliation
3481 :     Definition:
3482 :     FIGMODELmodel->SolutionReconciliation();
3483 :     Description:
3484 :     This is a wrapper for running the solution reconciliation algorithm on any model in the database.
3485 :     The algorithm performs a reconciliation of any gap filling solutions to identify the combination of solutions that results in the optimal model.
3486 :     This function prints out one output file in the Model directory: ReconciliationOutput.txt: this is a summary of the results of the reconciliation analysis
3487 :     =cut
3488 :    
3489 :     sub SolutionReconciliation {
3490 :     my ($self,$GapFill,$Stage) = @_;
3491 :    
3492 :     #Setting the output filenames
3493 :     my $OutputFilename;
3494 :     my $OutputFilenameTwo;
3495 :     if ($GapFill == 1) {
3496 :     $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GFReconciliation.txt";
3497 :     $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GFSRS.txt";
3498 :     } else {
3499 :     $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GGReconciliation.txt";
3500 :     $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GGSRS.txt";
3501 :     }
3502 :    
3503 :     #In stage one, we run the reconciliation and create a test file to check combined solution performance
3504 :     if (!defined($Stage) || $Stage == 1) {
3505 :     my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3506 :     my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3507 :     $Row->{"GF RECONCILATION TIMING"}->[0] = time()."-";
3508 :     $GrowMatchTable->save();
3509 :     $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3510 :    
3511 :     #Getting a unique filename
3512 :     my $UniqueFilename = $self->figmodel()->filename();
3513 :    
3514 :     #Copying over the necessary files
3515 :     if ($GapFill == 1) {
3516 :     if (!-e $self->directory().$self->id().$self->selected_version()."-GFEM.txt") {
3517 :     print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GFEM.txt file not found. Could not reconcile!";
3518 :     return 0;
3519 :     }
3520 :     if (!-e $self->directory().$self->id().$self->selected_version()."-OPEM.txt") {
3521 :     print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-OPEM.txt file not found. Could not reconcile!";
3522 :     return 0;
3523 :     }
3524 :     system("cp ".$self->directory().$self->id().$self->selected_version()."-GFEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GFEM.txt");
3525 :     system("cp ".$self->directory().$self->id().$self->selected_version()."-OPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3526 :     #Backing up and deleting the existing reconciliation file
3527 :     if (-e $OutputFilename) {
3528 :     system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGFReconciliation.txt");
3529 :     unlink($OutputFilename);
3530 :     }
3531 :     } else {
3532 :     if (!-e $self->directory().$self->id().$self->selected_version()."-GGEM.txt") {
3533 :     print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGEM.txt file not found. Could not reconcile!";
3534 :     return 0;
3535 :     }
3536 :     if (!-e $self->directory().$self->id().$self->selected_version()."-GGOPEM.txt") {
3537 :     print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt file not found. Could not reconcile!";
3538 :     return 0;
3539 :     }
3540 :     system("cp ".$self->directory().$self->id().$self->selected_version()."-GGEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GGEM.txt");
3541 :     system("cp ".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3542 :     #Backing up and deleting the existing reconciliation file
3543 :     if (-e $OutputFilename) {
3544 :     system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGGReconciliation.txt");
3545 :     unlink($OutputFilename);
3546 :     }
3547 :     }
3548 :    
3549 :     #Running the reconciliation
3550 :     system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),"NONE",["SolutionReconciliation"],{"Solution data for model optimization" => $UniqueFilename},"Reconciliation".$UniqueFilename.".log",undef,$self->selected_version()));
3551 :     $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3552 :     $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3553 :     $Row->{"GF RECONCILATION TIMING"}->[0] .= time();
3554 :     $GrowMatchTable->save();
3555 :     $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3556 :    
3557 :     #Loading the problem report from the reconciliation run
3558 :     my $ReconciliatonOutput = $self->figmodel()->LoadProblemReport($UniqueFilename);
3559 :     print $UniqueFilename."\n";
3560 :     #Clearing output files
3561 :     $self->figmodel()->clearing_output($UniqueFilename,"Reconciliation".$UniqueFilename.".log");
3562 :     $ReconciliatonOutput->save("/home/chenry/Test.txt");
3563 :    
3564 :     #Checking the a problem report was found and was loaded
3565 :     if (!defined($ReconciliatonOutput) || $ReconciliatonOutput->size() < 1 || !defined($ReconciliatonOutput->get_row(0)->{"Notes"}->[0])) {
3566 :     print STDERR "FIGMODEL:SolutionReconciliation: MFAToolkit output from SolutionReconciliation of ".$self->id()." not found!\n\n";
3567 :     return 0;
3568 :     }
3569 :    
3570 :     #Processing the solutions
3571 :     my $SolutionCount = 0;
3572 :     my $ReactionSetHash;
3573 :     my $SingleReactionHash;
3574 :     my $ReactionDataHash;
3575 :     for (my $n=0; $n < $ReconciliatonOutput->size(); $n++) {
3576 :     if (defined($ReconciliatonOutput->get_row($n)->{"Notes"}->[0]) && $ReconciliatonOutput->get_row($n)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^;]+)/) {
3577 :     #Breaking up the solution into reaction sets
3578 :     my @ReactionSets = split(/\|/,$1);
3579 :     #Creating reaction lists for each set
3580 :     my $SolutionHash;
3581 :     for (my $i=0; $i < @ReactionSets; $i++) {
3582 :     if (length($ReactionSets[$i]) > 0) {
3583 :     my @Alternatives = split(/:/,$ReactionSets[$i]);
3584 :     for (my $j=1; $j < @Alternatives; $j++) {
3585 :     if (length($Alternatives[$j]) > 0) {
3586 :     push(@{$SolutionHash->{$Alternatives[$j]}},$Alternatives[0]);
3587 :     }
3588 :     }
3589 :     if (@Alternatives == 1) {
3590 :     $SingleReactionHash->{$Alternatives[0]}->{$SolutionCount} = 1;
3591 :     if (!defined($SingleReactionHash->{$Alternatives[0]}->{"COUNT"})) {
3592 :     $SingleReactionHash->{$Alternatives[0]}->{"COUNT"} = 0;
3593 :     }
3594 :     $SingleReactionHash->{$Alternatives[0]}->{"COUNT"}++;
3595 :     }
3596 :     }
3597 :     }
3598 :     #Identifying reactions sets and storing the sets in the reactions set hash
3599 :     foreach my $Solution (keys(%{$SolutionHash})) {
3600 :     my $SetKey = join(",",sort(@{$SolutionHash->{$Solution}}));
3601 :     if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount})) {
3602 :     $ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount} = 1;
3603 :     if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"})) {
3604 :     $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"} = 0;
3605 :     }
3606 :     $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"}++;
3607 :     }
3608 :     $ReactionSetHash->{$SetKey}->{$Solution}->{$SolutionCount} = 1;
3609 :     if (!defined($ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"})) {
3610 :     $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"} = 0;
3611 :     }
3612 :     $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"}++;
3613 :     }
3614 :     $SolutionCount++;
3615 :     }
3616 :     }
3617 :    
3618 :     #Handling the scenario where no solutions were found
3619 :     if ($SolutionCount == 0) {
3620 :     print STDERR "FIGMODEL:SolutionReconciliation: Reconciliation unsuccessful. No solution found.\n\n";
3621 :     return 0;
3622 :     }
3623 :    
3624 :     #Printing results without solution performance figures. Also printing solution test file
3625 :     open (RECONCILIATION, ">$OutputFilename");
3626 :     #Printing the file heading
3627 :     print RECONCILIATION "DATABASE;DEFINITION;REVERSIBLITY;DELTAG;DIRECTION;NUMBER OF SOLUTIONS";
3628 :     for (my $i=0; $i < $SolutionCount; $i++) {
3629 :     print RECONCILIATION ";Solution ".$i;
3630 :     }
3631 :     print RECONCILIATION "\n";
3632 :     #Printing the singlet reactions first
3633 :     my $Solutions;
3634 :     print RECONCILIATION "SINGLET REACTIONS\n";
3635 :     my @SingletReactions = keys(%{$SingleReactionHash});
3636 :     for (my $j=0; $j < $SolutionCount; $j++) {
3637 :     $Solutions->[$j]->{"BASE"} = $j;
3638 :     }
3639 :     for (my $i=0; $i < @SingletReactions; $i++) {
3640 :     my $ReactionData;
3641 :     if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3642 :     $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3643 :     } else {
3644 :     my $Direction = substr($SingletReactions[$i],0,1);
3645 :     if ($Direction eq "+") {
3646 :     $Direction = "=>";
3647 :     } else {
3648 :     $Direction = "<=";
3649 :     }
3650 :     my $Reaction = substr($SingletReactions[$i],1);
3651 :     $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3652 :     $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3653 :     $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3654 :     if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3655 :     $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3656 :     }
3657 :     if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3658 :     $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3659 :     }
3660 :     if (!defined($ReactionData->{"DELTAG"}->[0])) {
3661 :     $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3662 :     }
3663 :     $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3664 :     }
3665 :     print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$SingleReactionHash->{$SingletReactions[$i]}->{"COUNT"};
3666 :     for (my $j=0; $j < $SolutionCount; $j++) {
3667 :     print RECONCILIATION ";";
3668 :     if (defined($SingleReactionHash->{$SingletReactions[$i]}->{$j})) {
3669 :     $Solutions->[$j]->{$SingletReactions[$i]} = 1;
3670 :     $Solutions->[$j]->{"BASE"} = $j;
3671 :     print RECONCILIATION "|".$j."|";
3672 :     }
3673 :     }
3674 :     print RECONCILIATION "\n";
3675 :     }
3676 :     #Printing the reaction sets with alternatives
3677 :     print RECONCILIATION "Reaction sets with alternatives\n";
3678 :     my @ReactionSets = keys(%{$ReactionSetHash});
3679 :     foreach my $ReactionSet (@ReactionSets) {
3680 :     my $NewSolutions;
3681 :     my $BaseReactions;
3682 :     my $AltList = [$ReactionSet];
3683 :     push(@{$AltList},keys(%{$ReactionSetHash->{$ReactionSet}}));
3684 :     for (my $j=0; $j < @{$AltList}; $j++) {
3685 :     my $CurrentNewSolutions;
3686 :     my $Index;
3687 :     if ($j == 0) {
3688 :     print RECONCILIATION "NEW SET\n";
3689 :     } elsif ($AltList->[$j] ne $ReactionSet) {
3690 :     print RECONCILIATION "ALTERNATIVE SET\n";
3691 :     #For each base solution in which this set is represented, we copy the base solution to the new solution
3692 :     my $NewSolutionCount = 0;
3693 :     for (my $k=0; $k < $SolutionCount; $k++) {
3694 :     if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3695 :     if (defined($Solutions)) {
3696 :     $Index->{$k} = @{$Solutions} + $NewSolutionCount;
3697 :     } else {
3698 :     $Index->{$k} = $NewSolutionCount;
3699 :     }
3700 :     if (defined($NewSolutions) && @{$NewSolutions} > 0) {
3701 :     $Index->{$k} += @{$NewSolutions};
3702 :     }
3703 :     $CurrentNewSolutions->[$NewSolutionCount] = {};
3704 :     foreach my $Reaction (keys(%{$Solutions->[$k]})) {
3705 :     $CurrentNewSolutions->[$NewSolutionCount]->{$Reaction} = $Solutions->[$k]->{$Reaction};
3706 :     }
3707 :     $NewSolutionCount++;
3708 :     }
3709 :     }
3710 :     }
3711 :     if ($j == 0 || $AltList->[$j] ne $ReactionSet) {
3712 :     my @SingletReactions = split(/,/,$AltList->[$j]);
3713 :     for (my $i=0; $i < @SingletReactions; $i++) {
3714 :     #Adding base reactions to base solutions and set reactions the new solutions
3715 :     if ($j == 0) {
3716 :     push(@{$BaseReactions},$SingletReactions[$i]);
3717 :     } else {
3718 :     for (my $k=0; $k < @{$CurrentNewSolutions}; $k++) {
3719 :     $CurrentNewSolutions->[$k]->{$SingletReactions[$i]} = 1;
3720 :     }
3721 :     }
3722 :     #Getting reaction data and printing reaction in output file
3723 :     my $ReactionData;
3724 :     if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3725 :     $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3726 :     } else {
3727 :     my $Direction = substr($SingletReactions[$i],0,1);
3728 :     if ($Direction eq "+") {
3729 :     $Direction = "=>";
3730 :     } else {
3731 :     $Direction = "<=";
3732 :     }
3733 :     my $Reaction = substr($SingletReactions[$i],1);
3734 :     $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3735 :     $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3736 :     $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3737 :     if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3738 :     $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3739 :     }
3740 :     if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3741 :     $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3742 :     }
3743 :     if (!defined($ReactionData->{"DELTAG"}->[0])) {
3744 :     $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3745 :     }
3746 :     $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3747 :     }
3748 :     print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{"COUNT"};
3749 :     for (my $k=0; $k < $SolutionCount; $k++) {
3750 :     print RECONCILIATION ";";
3751 :     if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3752 :     if ($j == 0) {
3753 :     print RECONCILIATION "|".$k."|";
3754 :     } else {
3755 :     print RECONCILIATION "|".$Index->{$k}."|";
3756 :     }
3757 :     }
3758 :     }
3759 :     print RECONCILIATION "\n";
3760 :     }
3761 :     #Adding the current new solutions to the new solutions array
3762 :     if (defined($CurrentNewSolutions) && @{$CurrentNewSolutions} > 0) {
3763 :     push(@{$NewSolutions},@{$CurrentNewSolutions});
3764 :     }
3765 :     }
3766 :     }
3767 :     #Adding the base reactions to all existing solutions
3768 :     for (my $j=0; $j < @{$Solutions}; $j++) {
3769 :     if (defined($ReactionSetHash->{$ReactionSet}->{$ReactionSet}->{$Solutions->[$j]->{"BASE"}})) {
3770 :     foreach my $SingleReaction (@{$BaseReactions}) {
3771 :     $Solutions->[$j]->{$SingleReaction} = 1;
3772 :     }
3773 :     }
3774 :     }
3775 :     #Adding the new solutions to the set of existing solutions
3776 :     push(@{$Solutions},@{$NewSolutions});
3777 :     }
3778 :     close(RECONCILIATION);
3779 :     #Now printing a file that defines all of the solutions in a format the testsolutions function understands
3780 :     open (RECONCILIATION, ">$OutputFilenameTwo");
3781 :     print RECONCILIATION "Experiment;Solution index;Solution cost;Solution reactions\n";
3782 :     for (my $i=0; $i < @{$Solutions}; $i++) {
3783 :     delete($Solutions->[$i]->{"BASE"});
3784 :     print RECONCILIATION "SR".$i.";".$i.";10;".join(",",keys(%{$Solutions->[$i]}))."\n";
3785 :     }
3786 :     close(RECONCILIATION);
3787 :    
3788 :     $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3789 :     $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3790 :     $Row->{"GF RECON TESTING TIMING"}->[0] = time()."-";
3791 :     $Row->{"GF RECON SOLUTIONS"}->[0] = @{$Solutions};
3792 :     $GrowMatchTable->save();
3793 :     $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3794 :    
3795 :     #Scheduling the solution testing
3796 :     if ($GapFill == 1) {
3797 :     system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GFSR:BACK:fast:QSUB\"");
3798 :     } else {
3799 :     system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GGSR:BACK:fast:QSUB\"");
3800 :     }
3801 :     } else {
3802 :     #Reading in the solution testing results
3803 :     my $Data;
3804 :     if ($GapFill == 1) {
3805 :     $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GFSREM.txt","");
3806 :     } else {
3807 :     $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GGSREM.txt","");
3808 :     }
3809 :    
3810 :     #Reading in the preliminate reconciliation report
3811 :     my $OutputData = $self->figmodel()->database()->load_single_column_file($OutputFilename,"");
3812 :     #Replacing the file tags with actual performance data
3813 :     my $Count = 0;
3814 :     for (my $i=0; $i < @{$Data}; $i++) {
3815 :     if ($Data->[$i] =~ m/^SR(\d+);.+;(\d+\/\d+);/) {
3816 :     my $Index = $1;
3817 :     my $Performance = $Index."/".$2;
3818 :     for (my $j=0; $j < @{$OutputData}; $j++) {
3819 :     $OutputData->[$j] =~ s/\|$Index\|/$Performance/g;
3820 :     }
3821 :     }
3822 :     }
3823 :     $self->figmodel()->database()->print_array_to_file($OutputFilename,$OutputData);
3824 :    
3825 :     my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3826 :     my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3827 :     $Row->{"GF RECON TESTING TIMING"}->[0] .= time();
3828 :     $GrowMatchTable->save();
3829 :     $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3830 :     }
3831 :    
3832 :     return 1;
3833 :     }
3834 :    
3835 : chenry 1.27 =head3 DetermineCofactorLipidCellWallComponents
3836 : chenry 1.18 Definition:
3837 : chenry 1.27 {cofactor=>{string:compound id=>float:coefficient},lipid=>...cellWall=>} = FIGMODELmodel->DetermineCofactorLipidCellWallComponents();
3838 : chenry 1.18 Description:
3839 :     =cut
3840 : chenry 1.27 sub DetermineCofactorLipidCellWallComponents {
3841 : chenry 1.18 my ($self) = @_;
3842 : chenry 1.27 my $templateResults;
3843 :     my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
3844 :     my $Class = $self->{_data}->cellwalltype();
3845 :     my $Name = $self->name();
3846 :     my $translation = {COFACTOR=>"cofactor",LIPIDS=>"lipid","CELL WALL"=>"cellWall"};
3847 :     #Checking for phoenix variants
3848 :     my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");
3849 :     my $Phoenix = 0;
3850 :     my @Rows = $PhoenixVariantTable->get_rows_by_key($self->genome(),"GENOME");
3851 :     my $VariantHash;
3852 :     for (my $i=0; $i < @Rows; $i++) {
3853 :     $Phoenix = 1;
3854 :     if (defined($Rows[$i]->{"SUBSYSTEM"}) && defined($Rows[$i]->{"VARIANT"})) {
3855 :     $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];
3856 : chenry 1.18 }
3857 : chenry 1.27 }
3858 :     #Collecting genome data
3859 :     my $RoleHash;
3860 :     my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
3861 :     for (my $i=0; $i < $FeatureTable->size(); $i++) {
3862 :     for (my $j=0; $j < @{$FeatureTable->get_row($i)->{"ROLES"}}; $j++) {
3863 :     $RoleHash->{$FeatureTable->get_row($i)->{"ROLES"}->[$j]} = 1;
3864 :     my $Subsystems = $self->figmodel()->subsystems_of_role($FeatureTable->get_row($i)->{"ROLES"}->[$j]);
3865 :     if (defined($Subsystems)) {
3866 :     for (my $k=0; $k < @{$Subsystems}; $k++) {
3867 :     if ($Phoenix == 0) {
3868 :     $VariantHash->{$Subsystems->[$k]} = 1;
3869 : chenry 1.18 }
3870 :     }
3871 :     }
3872 :     }
3873 : chenry 1.27 }
3874 :     #Scanning through the template item by item and determinine which biomass components should be added
3875 :     my $includedHash;
3876 :     my $BiomassReactionTemplateTable = $self->figmodel()->database()->get_table("BIOMASSTEMPLATE");
3877 :     for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {
3878 :     my $Row = $BiomassReactionTemplateTable->get_row($i);
3879 :     if (defined($translation->{$Row->{CLASS}->[0]})) {
3880 :     my $coef = -1;
3881 :     if ($Row->{"REACTANT"}->[0] eq "NO") {
3882 :     $coef = 1;
3883 :     if ($Row->{"COEFFICIENT"}->[0] =~ m/cpd/) {
3884 :     $coef = $Row->{"COEFFICIENT"}->[0];
3885 :     }
3886 :     }
3887 : chenry 1.18 if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {
3888 : chenry 1.27 $includedHash->{$Row->{"ID"}->[0]} = 1;
3889 :     $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
3890 :     } elsif (defined($Row->{"INCLUSION CRITERIA"}->[0])) {
3891 : chenry 1.18 my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];
3892 :     my $End = 0;
3893 :     while ($End == 0) {
3894 :     if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {
3895 : chenry 1.27 print $Criteria." : ";
3896 : chenry 1.18 my $Start = "";
3897 :     my $End = "";
3898 :     my $Condition = $1;
3899 :     my $Data = $2;
3900 :     if ($1 ne "AND" && $1 ne "OR") {
3901 :     $Start = $1;
3902 :     $End = $4;
3903 :     $Condition = $2;
3904 :     $Data = $3;
3905 :     }
3906 :     my $Result = "YES";
3907 :     if ($Condition eq "OR") {
3908 :     $Result = "NO";
3909 :     }
3910 :     my @Array = split(/\|/,$Data);
3911 :     for (my $j=0; $j < @Array; $j++) {
3912 :     if ($Array[$j] eq "YES" && $Condition eq "OR") {
3913 :     $Result = "YES";
3914 :     last;
3915 :     } elsif ($Array[$j] eq "NO" && $Condition eq "AND") {
3916 :     $Result = "NO";
3917 :     last;
3918 :     } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {
3919 : chenry 1.27 if (defined($includedHash->{$1}) && $Condition eq "OR") {
3920 : chenry 1.18 $Result = "YES";
3921 :     last;
3922 : chenry 1.27 } elsif (!defined($includedHash->{$1}) && $Condition eq "AND") {
3923 : chenry 1.18 $Result = "NO";
3924 :     last;
3925 :     }
3926 :     } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {
3927 : chenry 1.27 if (!defined($includedHash->{$1}) && $Condition eq "OR") {
3928 : chenry 1.18 $Result = "YES";
3929 :     last;
3930 : chenry 1.27 } elsif (defined($includedHash->{$1}) && $Condition eq "AND") {
3931 : chenry 1.18 $Result = "NO";
3932 :     last;
3933 :     }
3934 :     } elsif ($Array[$j] =~ m/^NAME:(.+)/) {
3935 :     my $Comparison = $1;
3936 :     if ((!defined($Comparison) || !defined($Name) || $Name =~ m/$Comparison/) && $Condition eq "OR") {
3937 :     $Result = "YES";
3938 :     last;
3939 :     } elsif (defined($Comparison) && defined($Name) && $Name !~ m/$Comparison/ && $Condition eq "AND") {
3940 :     $Result = "NO";
3941 :     last;
3942 :     }
3943 :     } elsif ($Array[$j] =~ m/^!NAME:(.+)/) {
3944 :     my $Comparison = $1;
3945 :     if ((!defined($Comparison) || !defined($Name) || $Name !~ m/$Comparison/) && $Condition eq "OR") {
3946 :     $Result = "YES";
3947 :     last;
3948 :     } elsif (defined($Comparison) && defined($Name) && $Name =~ m/$Comparison/ && $Condition eq "AND") {
3949 :     $Result = "NO";
3950 :     last;
3951 :     }
3952 :     } elsif ($Array[$j] =~ m/^SUBSYSTEM:(.+)/) {
3953 :     my @SubsystemArray = split(/`/,$1);
3954 :     if (@SubsystemArray == 1) {
3955 :     if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "OR") {
3956 :     $Result = "YES";
3957 :     last;
3958 :     } elsif ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "AND") {
3959 :     $Result = "NO";
3960 :     last;
3961 :     }
3962 :     } else {
3963 :     my $Match = 0;
3964 :     for (my $k=1; $k < @SubsystemArray; $k++) {
3965 :     if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) {
3966 :     $Match = 1;
3967 :     last;
3968 :     }
3969 :     }
3970 :     if ($Match == 1 && $Condition eq "OR") {
3971 :     $Result = "YES";
3972 :     last;
3973 :     } elsif ($Match != 1 && $Condition eq "AND") {
3974 :     $Result = "NO";
3975 :     last;
3976 :     }
3977 :     }
3978 :     } elsif ($Array[$j] =~ m/^!SUBSYSTEM:(.+)/) {
3979 :     my @SubsystemArray = split(/`/,$1);
3980 :     if (@SubsystemArray == 1) {
3981 :     if ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "OR") {
3982 :     $Result = "YES";
3983 :     last;
3984 :     } elsif (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "AND") {
3985 :     $Result = "NO";
3986 :     last;
3987 :     }
3988 :     } else {
3989 :     my $Match = 0;
3990 :     for (my $k=1; $k < @SubsystemArray; $k++) {
3991 :     if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) {
3992 :     $Match = 1;
3993 :     last;
3994 :     }
3995 :     }
3996 :     if ($Match != 1 && $Condition eq "OR") {
3997 :     $Result = "YES";
3998 :     last;
3999 :     } elsif ($Match == 1 && $Condition eq "AND") {
4000 :     $Result = "NO";
4001 :     last;
4002 :     }
4003 :     }
4004 :     } elsif ($Array[$j] =~ m/^ROLE:(.+)/) {
4005 :     if (defined($RoleHash->{$1}) && $Condition eq "OR") {
4006 :     $Result = "YES";
4007 :     last;
4008 :     } elsif (!defined($RoleHash->{$1}) && $Condition eq "AND") {
4009 :     $Result = "NO";
4010 :     last;
4011 :     }
4012 :     } elsif ($Array[$j] =~ m/^!ROLE:(.+)/) {
4013 :     if (!defined($RoleHash->{$1}) && $Condition eq "OR") {
4014 :     $Result = "YES";
4015 :     last;
4016 :     } elsif (defined($RoleHash->{$1}) && $Condition eq "AND") {
4017 :     $Result = "NO";
4018 :     last;
4019 :     }
4020 :     } elsif ($Array[$j] =~ m/^CLASS:(.+)/) {
4021 :     if ($Class eq $1 && $Condition eq "OR") {
4022 :     $Result = "YES";
4023 :     last;
4024 :     } elsif ($Class ne $1 && $Condition eq "AND") {
4025 :     $Result = "NO";
4026 :     last;
4027 :     }
4028 :     } elsif ($Array[$j] =~ m/^!CLASS:(.+)/) {
4029 :     if ($Class ne $1 && $Condition eq "OR") {
4030 :     $Result = "YES";
4031 :     last;
4032 :     } elsif ($Class eq $1 && $Condition eq "AND") {
4033 :     $Result = "NO";
4034 :     last;
4035 :     }
4036 :     }
4037 :     }
4038 :     $Criteria = $Start.$Result.$End;
4039 :     print $Criteria."\n";
4040 :     } else {
4041 :     $End = 1;
4042 :     last;
4043 :     }
4044 :     }
4045 :     if ($Criteria eq "YES") {
4046 : chenry 1.27 $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
4047 :     $includedHash->{$Row->{"ID"}->[0]} = 1;
4048 : chenry 1.18 }
4049 :     }
4050 :     }
4051 : chenry 1.27 }
4052 :     my $types = ["cofactor","lipid","cellWall"];
4053 :     my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4054 :     for (my $i=0; $i < @{$types}; $i++) {
4055 :     my @list = keys(%{$templateResults->{$types->[$i]}});
4056 :     my $entries = 0;
4057 :     for (my $j=0; $j < @list; $j++) {
4058 :     if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4059 :     my $objs = $cpdMgr->get_objects({id=>$list[$j]});
4060 :     if (!defined($objs->[0]) || $objs->[0]->mass() == 0) {
4061 :     $templateResults->{$types->[$i]}->{$list[$j]} = -1e-5;
4062 : chenry 1.18 } else {
4063 : chenry 1.27 $entries++;
4064 : chenry 1.18 }
4065 : chenry 1.27 }
4066 :     }
4067 :     for (my $j=0; $j < @list; $j++) {
4068 :     if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4069 :     $templateResults->{$types->[$i]}->{$list[$j]} = -1/$entries;
4070 :     } elsif ($templateResults->{$types->[$i]}->{$list[$j]} =~ m/cpd/) {
4071 :     my $netCoef = 0;
4072 :     my @allcpd = split(/,/,$templateResults->{$types->[$i]}->{$list[$j]});
4073 :     for (my $k=0; $k < @allcpd; $k++) {
4074 :     if (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} ne "-1e-5") {
4075 :     $netCoef += (1/$entries);
4076 :     } elsif (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} eq "-1e-5") {
4077 :     $netCoef += 1e-5;
4078 :     }
4079 : chenry 1.18 }
4080 : chenry 1.27 $templateResults->{$types->[$i]}->{$list[$j]} = $netCoef;
4081 : chenry 1.18 }
4082 :     }
4083 : chenry 1.27 }
4084 :     return $templateResults;
4085 :     }
4086 :    
4087 :     =head3 BuildSpecificBiomassReaction
4088 :     Definition:
4089 :     FIGMODELmodel->BuildSpecificBiomassReaction();
4090 :     Description:
4091 :     =cut
4092 :     sub BuildSpecificBiomassReaction {
4093 :     my ($self) = @_;
4094 :     #Getting the database handle for biomass reactions
4095 :     my $bioMgr = $self->figmodel()->database()->get_object_manager("bof");
4096 :     #Checking if the current biomass reaction appears in more than on model, if not, this biomass reaction is conserved for this model
4097 :     my $biomassID = $self->biomassReaction();
4098 :     if ($biomassID =~ m/bio\d\d\d\d\d/) {
4099 :     my $mdlMgr = $self->figmodel()->database()->get_object_manager("model");
4100 :     my $mdlObs = $mdlMgr->get_objects({biomassReaction=>$biomassID});
4101 :     if (defined($mdlObs->[1])) {
4102 :     $biomassID = "NONE";
4103 :     }
4104 :     }
4105 :     #If the biomass ID is "NONE", then we create a new biomass reaction for the model
4106 :     my $bioObj;
4107 :     my $originalPackages = "";
4108 :     my $originalEssReactions = "";
4109 : chenry 1.28 if ($biomassID !~ m/bio\d\d\d\d\d/) {
4110 : chenry 1.27 #Getting the current largest ID
4111 :     $biomassID = $self->figmodel()->database()->check_out_new_id("bof");
4112 :     $bioObj = $bioMgr->create({
4113 :     id=>$biomassID,owner=>$self->owner(),users=>$self->users(),name=>"Biomass",equation=>"NONE",protein=>"0",
4114 :     energy=>"0",DNA=>"0",RNA=>"0",lipid=>"0",cellWall=>"0",cofactor=>"0",
4115 :     modificationDate=>time(),creationDate=>time(),
4116 :     cofactorPackage=>"NONE",lipidPackage=>"NONE",cellWallPackage=>"NONE",
4117 :     DNACoef=>"NONE",RNACoef=>"NONE",proteinCoef=>"NONE",lipidCoef=>"NONE",
4118 :     cellWallCoef=>"NONE",cofactorCoef=>"NONE",essentialRxn=>"NONE"});
4119 :     if (!defined($bioObj)) {
4120 :     die $self->error_message("BuildSpecificBiomassReaction():Could not create new biomass reaction ".$biomassID."!");
4121 :     }
4122 :     } else {
4123 :     #Getting the biomass DB handler from the database
4124 :     my $objs = $bioMgr->get_objects({id=>$biomassID});
4125 :     if (!defined($objs->[0])) {
4126 :     die $self->error_message("BuildSpecificBiomassReaction():Could not find biomass reaction ".$biomassID." in database!");
4127 :     }
4128 :     $bioObj = $objs->[0];
4129 :     $bioObj->owner($self->owner());
4130 :     $bioObj->users($self->users());
4131 :     if (defined($bioObj->essentialRxn())) {
4132 :     $originalEssReactions = $bioObj->essentialRxn();
4133 :     $originalPackages = $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage();
4134 :     }
4135 :     }
4136 :     #Getting genome stats
4137 :     my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
4138 :     my $Class = $self->{_data}->cellwalltype();
4139 :     #Setting global coefficients based on cell wall type
4140 :     my $biomassCompounds;
4141 :     my $compounds;
4142 :     if ($Class eq "Gram positive") {
4143 :     $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.323,cpd00052=>-0.199,cpd00062=>-0.215};
4144 :     $compounds->{protein} = {cpd00001=>1,cpd00023=>-0.0637,cpd00033=>-0.0999,cpd00035=>-0.0653,cpd00039=>-0.0790,cpd00041=>-0.0362,cpd00051=>-0.0472,cpd00053=>-0.0637,cpd00054=>-0.0529,cpd00060=>-0.0277,cpd00065=>-0.0133,cpd00066=>-0.0430,cpd00069=>-0.0271,cpd00084=>-0.0139,cpd00107=>-0.0848,cpd00119=>-0.0200,cpd00129=>-0.0393,cpd00132=>-0.0362,cpd00156=>-0.0751,cpd00161=>-0.0456,cpd00322=>-0.0660};
4145 :     $bioObj->protein("0.5284");
4146 :     $bioObj->DNA("0.026");
4147 :     $bioObj->RNA("0.0655");
4148 :     $bioObj->lipid("0.075");
4149 :     $bioObj->cellWall("0.25");
4150 :     $bioObj->cofactor("0.10");
4151 :     } else {
4152 :     $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.322,cpd00052=>-0.2,cpd00062=>-0.216};
4153 :     $compounds->{protein} = {cpd00001=>1,cpd00023=>-0.0492,cpd00033=>-0.1145,cpd00035=>-0.0961,cpd00039=>-0.0641,cpd00041=>-0.0451,cpd00051=>-0.0554,cpd00053=>-0.0492,cpd00054=>-0.0403,cpd00060=>-0.0287,cpd00065=>-0.0106,cpd00066=>-0.0347,cpd00069=>-0.0258,cpd00084=>-0.0171,cpd00107=>-0.0843,cpd00119=>-0.0178,cpd00129=>-0.0414,cpd00132=>-0.0451,cpd00156=>-0.0791,cpd00161=>-0.0474,cpd00322=>-0.0543};
4154 :     $bioObj->protein("0.563");
4155 :     $bioObj->DNA("0.031");
4156 :     $bioObj->RNA("0.21");
4157 :     $bioObj->lipid("0.093");
4158 :     $bioObj->cellWall("0.177");
4159 :     $bioObj->cofactor("0.039");
4160 :     }
4161 :     #Setting energy coefficient for all reactions
4162 :     $bioObj->energy("40");
4163 :     $compounds->{energy} = {cpd00002=>-1,cpd00001=>-1,cpd00008=>1,cpd00009=>1,cpd00067=>1};
4164 :     #Setting DNA coefficients based on GC content
4165 :     my $gc = $self->figmodel()->get_genome_gc_content($self->genome());
4166 :     $compounds->{DNA} = {cpd00012=>1,cpd00115=>0.5*(1-$gc),cpd00241=>0.5*$gc,cpd00356=>0.5*$gc,cpd00357=>0.5*(1-$gc)};
4167 :     #Setting Lipid,cell wall,and cofactor coefficients based on biomass template
4168 :     my $templateResults = $self->DetermineCofactorLipidCellWallComponents();
4169 :     $compounds->{cofactor} = $templateResults->{cofactor};
4170 :     $compounds->{lipid} = $templateResults->{lipid};
4171 :     $compounds->{cellWall} = $templateResults->{cellWall};
4172 :     #Getting package number for cofactor, lipid, and cell wall
4173 :     my $packages;
4174 :     my $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4175 :     my $packageTypes = ["Cofactor","Lipid","CellWall"];
4176 :     my $translation = {"Cofactor"=>"cofactor","Lipid"=>"lipid","CellWall"=>"cellWall"};
4177 :     for (my $i=0; $i < @{$packageTypes}; $i++) {
4178 :     my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4179 :     my $function = $translation->{$packageTypes->[$i]}."Package";
4180 :     if (@cpdList == 0) {
4181 :     $bioObj->$function("NONE");
4182 :     } else {
4183 :     my $cpdgrpObs = $cpdgrpMgr->get_objects({type=>$packageTypes->[$i]."Package"});
4184 :     for (my $j=0; $j < @{$cpdgrpObs}; $j++) {
4185 :     $packages->{$packageTypes->[$i]}->{$cpdgrpObs->[$j]->grouping()}->{$cpdgrpObs->[$j]->COMPOUND()} = 1;
4186 :     }
4187 :     my @packageList = keys(%{$packages->{$packageTypes->[$i]}});
4188 :     my $packageHash;
4189 :     for (my $j=0; $j < @packageList; $j++) {
4190 :     $packageHash->{join("|",sort(keys(%{$packages->{$packageTypes->[$i]}->{$packageList[$j]}})))} = $packageList[$j];
4191 :     }
4192 :     if (defined($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))})) {
4193 :     $bioObj->$function($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))});
4194 :     } else {
4195 :     my $newPackageID = $self->figmodel()->database()->check_out_new_id($packageTypes->[$i]."Pkg");
4196 :     $bioObj->$function($newPackageID);
4197 :     my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4198 :     for (my $j=0; $j < @cpdList; $j++) {
4199 :     $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4200 :     $cpdgrpMgr->create({COMPOUND=>$cpdList[$j],grouping=>$newPackageID,type=>$packageTypes->[$i]."Package"});
4201 :     }
4202 :     }
4203 :     }
4204 :     }
4205 :     #Filling in coefficient terms in database and calculating global reaction coefficients based on classification abundancies
4206 :     my $equationCompounds;
4207 :     my $types = ["RNA","DNA","protein","lipid","cellWall","cofactor","energy"];
4208 :     my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4209 :     for (my $i=0; $i < @{$types}; $i++) {
4210 :     my $coefString = "";
4211 :     my @compounds = sort(keys(%{$compounds->{$types->[$i]}}));
4212 :     #Building coefficient strings and determining net mass for component types
4213 :     my $netMass = 0;
4214 :     for (my $j=0; $j < @compounds; $j++) {
4215 :     my $objs = $cpdMgr->get_objects({id=>$compounds[$j]});
4216 :     my $mass = 0;
4217 :     if (defined($objs->[0]) && $objs->[0]->mass() != 0) {
4218 :     $mass = $objs->[0]->mass();
4219 :     $netMass += -$compounds->{$types->[$i]}->{$compounds[$j]}*$objs->[0]->mass();
4220 :     }
4221 :     if (!defined($equationCompounds->{$compounds[$j]})) {
4222 :     $equationCompounds->{$compounds[$j]}->{"coef"} = 0;
4223 :     $equationCompounds->{$compounds[$j]}->{"type"} = $types->[$i];
4224 :     $equationCompounds->{$compounds[$j]}->{"mass"} = $mass;
4225 :     }
4226 :     $coefString .= $compounds->{$types->[$i]}->{$compounds[$j]}."|";
4227 :     }
4228 :     $netMass = 0.001*$netMass;
4229 :     #Calculating coefficients for all component compounds
4230 :     for (my $j=0; $j < @compounds; $j++) {
4231 :     #Normalizing certain type coefficients by mass
4232 :     my $function = $types->[$i];
4233 :     my $fraction = $bioObj->$function();
4234 :     if ($types->[$i] ne "energy") {
4235 :     $fraction = $fraction/$netMass;
4236 :     }
4237 :     if ($compounds->{$types->[$i]}->{$compounds[$j]} eq 1e-5) {
4238 :     $fraction = 1;
4239 :     }
4240 :     $equationCompounds->{$compounds[$j]}->{"coef"} += $fraction*$compounds->{$types->[$i]}->{$compounds[$j]};
4241 :     }
4242 :     chop($coefString);
4243 :     if (length($coefString) == 0) {
4244 :     $coefString = "NONE";
4245 :     }
4246 :     my $function = $types->[$i]."Coef";
4247 :     if ($types->[$i] ne "energy") {
4248 :     $bioObj->$function($coefString);
4249 :     }
4250 :     }
4251 :     #Adding biomass to compound list
4252 : chenry 1.29 $equationCompounds->{cpd17041}->{coef} = -1;
4253 :     $equationCompounds->{cpd17041}->{type} = "macromolecule";
4254 :     $equationCompounds->{cpd17042}->{coef} = -1;
4255 :     $equationCompounds->{cpd17042}->{type} = "macromolecule";
4256 :     $equationCompounds->{cpd17043}->{coef} = -1;
4257 :     $equationCompounds->{cpd17043}->{type} = "macromolecule";
4258 : chenry 1.27 $equationCompounds->{cpd11416}->{coef} = 1;
4259 :     $equationCompounds->{cpd11416}->{type} = "macromolecule";
4260 :     #Building equation from hash and populating compound biomass table
4261 :     my @compoundList = keys(%{$equationCompounds});
4262 :     my ($reactants,$products);
4263 :     #Deleting existing Biomass Compound info
4264 :     my $cpdbofMgr = $self->figmodel()->database()->get_object_manager("cpdbof");
4265 :     my $matchingObjs = $cpdbofMgr->get_objects({BIOMASS=>$biomassID});
4266 :     for (my $i=0; $i < @{$matchingObjs}; $i++) {
4267 :     $matchingObjs->[$i]->delete();
4268 :     }
4269 :     my $typeCategories = {"macromolecule"=>"M","RNA"=>"R","DNA"=>"D","protein"=>"P","lipid"=>"L","cellWall"=>"W","cofactor"=>"C","energy"=>"E"};
4270 :     my $productmass = 0;
4271 :     my $reactantmass = 0;
4272 :     my $totalmass = 0;
4273 :     foreach my $compound (@compoundList) {
4274 :     if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4275 :     $totalmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4276 :     }
4277 :     if ($equationCompounds->{$compound}->{coef} < 0) {
4278 :     if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4279 :     $reactantmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4280 :     }
4281 :     $reactants->{$compound} = $self->figmodel()->format_coefficient(-1*$equationCompounds->{$compound}->{coef});
4282 :     } else {
4283 :     if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4284 :     $productmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4285 : chenry 1.18 }
4286 : chenry 1.27 $products->{$compound} = $self->figmodel()->format_coefficient($equationCompounds->{$compound}->{coef});
4287 : chenry 1.18 }
4288 : chenry 1.27 #Adding biomass reaction compounds to the biomass compound table
4289 :     $cpdbofMgr = $self->figmodel()->database()->get_object_manager("cpdbof");
4290 :     $cpdbofMgr->create({COMPOUND=>$compound,BIOMASS=>$biomassID,coefficient=>$equationCompounds->{$compound}->{coef},compartment=>"c",category=>$typeCategories->{$equationCompounds->{$compound}->{type}}});
4291 :     }
4292 :     print "Total mass = ".$totalmass.", Reactant mass = ".$reactantmass.", Product mass = ".$productmass."\n";
4293 :     my $Equation = "";
4294 :     my @ReactantList = sort(keys(%{$reactants}));
4295 :     for (my $i=0; $i < @ReactantList; $i++) {
4296 :     if (length($Equation) > 0) {
4297 :     $Equation .= " + ";
4298 :     }
4299 :     $Equation .= "(".$reactants->{$ReactantList[$i]}.") ".$ReactantList[$i];
4300 :     }
4301 :     $Equation .= " => ";
4302 :     my $First = 1;
4303 :     @ReactantList = sort(keys(%{$products}));
4304 :     for (my $i=0; $i < @ReactantList; $i++) {
4305 :     if ($First == 0) {
4306 :     $Equation .= " + ";
4307 :     }
4308 :     $First = 0;
4309 :     $Equation .= "(".$products->{$ReactantList[$i]}.") ".$ReactantList[$i];
4310 :     }
4311 :     $bioObj->equation($Equation);
4312 :     #Setting the biomass reaction of this model
4313 :     $self->biomassReaction($biomassID);
4314 :     $self->figmodel()->print_biomass_reaction_file($biomassID);
4315 :     #Checking if the biomass reaction remained unchanged
4316 :     if ($originalPackages ne "" && $originalPackages eq $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage()) {
4317 :     print "UNCHANGED!\n";
4318 :     $bioObj->essentialRxn($originalEssReactions);
4319 :     } else {
4320 :     #Copying essential reaction lists if the packages in this biomasses reaction exactly match those in another biomass reaction
4321 :     my $matches = $bioMgr->get_objects({cofactorPackage=>$bioObj->cofactorPackage(),lipidPackage=>$bioObj->lipidPackage(),cellWallPackage=>$bioObj->cellWallPackage()});
4322 :     my $matchFound = 0;
4323 :