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

Diff of /FigKernelPackages/FIGMODELmodel.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1, Sun Oct 25 21:42:35 2009 UTC revision 1.2, Fri Dec 11 21:37:02 2009 UTC
# Line 13  Line 13 
13          This is the constructor for the FIGMODELmodel object.          This is the constructor for the FIGMODELmodel object.
14  =cut  =cut
15  sub new {  sub new {
16          my ($class,$id,$figmodel) = @_;          my ($class,$figmodel,$id) = @_;
17          my $self->{"_figmodel"}->[0] = $figmodel;          #Error checking first
18          $self->{"_id"}->[0] = $id;          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    
29            #Checking that the id exists
30            my $tbl = $self->figmodel()->database()->GetDBTable("MODELS");
31            if (!defined($tbl)) {
32                    $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!");
33                    return undef;
34            }
35    
36            #If the id is a number, we get the model row by index
37            my $index = $id;
38            if ($id =~ m/^\d+$/) {
39                    $self->{_data} = $tbl->get_row($id);
40            } else {
41                    $self->{_data} = $tbl->get_row_by_key($id,"id");
42                    if (defined($self->{_data})) {
43                            $index = $tbl->row_index($self->{_data});
44                    }
45            }
46            if (!defined($self->{_data})) {
47                    $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find specified id in database!");
48                    return undef;
49            }
50            $self->{_index} = $index;
51            $self->figmodel()->{_models}->{$self->id()} = $self;
52            $self->figmodel()->{_models}->{$self->index()} = $self;
53    
54            return $self;
55    }
56    
57    =head3 config
58    Definition:
59            ref::key value = FIGMODELmodel->config(string::key);
60    Description:
61            Trying to avoid using calls that assume configuration data is stored in a particular manner.
62            Call this function to get file paths etc.
63    =cut
64    sub config {
65            my ($self,$key) = @_;
66            return $self->figmodel()->config($key);
67    }
68    
69    =head3 fail
70    Definition:
71            -1 = FIGMODELmodel->fail();
72    Description:
73            Standard return for failed functions.
74    =cut
75    sub fail {
76            my ($self) = @_;
77            return $self->figmodel()->fail();
78    }
79    
80          return bless $self;  =head3 success
81    Definition:
82            1 = FIGMODELmodel->success();
83    Description:
84            Standard return for successful functions.
85    =cut
86    sub success {
87            my ($self) = @_;
88            return $self->figmodel()->success();
89  }  }
90    
91  =head3 figmodel  =head3 figmodel
# Line 28  Line 96 
96  =cut  =cut
97  sub figmodel {  sub figmodel {
98          my ($self) = @_;          my ($self) = @_;
99          return $self->{"_figmodel"}->[0];          return $self->{"_figmodel"};
100    }
101    
102    =head3 fig
103    Definition:
104            FIGMODEL = FIGMODELmodel->fig();
105    Description:
106            Returns a FIG object
107    =cut
108    sub fig {
109            my ($self) = @_;
110    
111            if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {
112                    if ($self->source() =~ /^RAST/) {
113                            $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());
114                    } else {
115                            $self->{"_fig"} = $self->figmodel()->fig();
116                    }
117            }
118    
119            return $self->{"_fig"};
120    }
121    
122    =head3 mgdata
123    Definition:
124            FIGMODEL = FIGMODELmodel->mgdata();
125    Description:
126            Returns a mgrast database object
127    =cut
128    sub mgdata {
129            my ($self) = @_;
130    
131            if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) {
132                    require MGRAST;
133                    $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } )
134            }
135    
136            return $self->{_mgdata};
137  }  }
138    
139  =head3 id  =head3 id
# Line 39  Line 144 
144  =cut  =cut
145  sub id {  sub id {
146          my ($self) = @_;          my ($self) = @_;
147          return $self->{"_id"}->[0];          return $self->{_data}->{id}->[0];
148    }
149    
150    =head3 owner
151    Definition:
152            string = FIGMODELmodel->owner();
153    Description:
154            Returns the username for the model owner
155    =cut
156    sub owner {
157            my ($self) = @_;
158            return $self->{_data}->{owner}->[0];
159    }
160    
161    =head3 index
162    Definition:
163            string = FIGMODELmodel->index();
164    Description:
165            Returns model index
166    =cut
167    sub index {
168            my ($self) = @_;
169            return $self->{_index};
170    }
171    
172    =head3 name
173    Definition:
174            string = FIGMODELmodel->name();
175    Description:
176            Returns the name of the organism or metagenome sample being modeled
177    =cut
178    sub name {
179            my ($self) = @_;
180    
181            if (!defined($self->{_name})) {
182                    my $source = $self->source();
183                    if ($source =~ /^MGRAST/) {
184                            $self->{_name} = "NA";
185                            if (defined($self->mgdata())) {
186                                    $self->{_name} = $self->mgdata()->genome_name;
187                            }
188                    } elsif (defined($self->stats())) {
189                            $self->{_name} = $self->stats()->{'Organism name'}->[0];
190                    } elsif ($source !~ /^RAST/) {
191                            $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());
192                    }
193            }
194    
195            return $self->{_name};
196    }
197    
198    =head3 stats
199    Definition:
200            string = FIGMODELmodel->stats();
201    Description:
202            Returns model stats
203    =cut
204    sub stats {
205            my ($self) = @_;
206    
207            if (!defined($self->{_stats})) {
208                    $self->{_stats} = $self->figmodel()->database()->GetDBTable("MODEL STATS")->get_row_by_key($self->id(),"Model ID");
209            }
210            return $self->{_stats};
211    }
212    
213    =head3 get_reaction_class
214    Definition:
215            string = FIGMODELmodel->get_reaction_class(string::reaction ID);
216    Description:
217            Returns reaction class
218    =cut
219    sub get_reaction_class {
220            my ($self,$reaction,$nohtml) = @_;
221    
222            if (!-e $self->directory()."ReactionClassification-".$self->id().".tbl") {
223                    if (!defined($self->{_reaction_classes})) {
224                            $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id()."-Complete.tbl",";","|",0,["REACTION"]);
225                            if (!defined($self->{_reaction_classes})) {
226                                    return undef;
227                            }
228                    }
229    
230                    my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
231                    if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
232                            my $class;
233                            if ($ClassRow->{CLASS}->[0] eq "Positive") {
234                                    $class = "Essential =>";
235                            } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {
236                                    $class = "Essential <=";
237                            } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {
238                                    $class = "Active =>";
239                            } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {
240                                    $class = "Active <=";
241                            } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {
242                                    $class = "Active <=>";
243                            } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {
244                                    $class = "Inactive";
245                            } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {
246                                    $class = "Disconnected";
247                            }
248                            if (!defined($nohtml) || $nohtml ne "1") {
249                                    $class = "<span title=\"Flux:".$ClassRow->{MIN}->[0]." to ".$ClassRow->{MAX}->[0]."\">".$class."</span>";
250                            }
251                            return $class;
252                    }
253                    return undef;
254            }
255    
256            if (!defined($self->{_reaction_classes})) {
257                    $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().".tbl",";","|",0,["REACTION"]);
258                    if (!defined($self->{_reaction_classes})) {
259                            return undef;
260                    }
261            }
262    
263            my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
264            my $classstring = "";
265            if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
266                    for (my $i=0; $i < @{$ClassRow->{CLASS}};$i++) {
267                            if (length($classstring) > 0) {
268                                    $classstring .= "<br>";
269                            }
270                            my $NewClass;
271                            if ($ClassRow->{CLASS}->[$i] eq "Positive") {
272                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";
273                                    if (!defined($nohtml) || $nohtml ne "1") {
274                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
275                                    }
276                            } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {
277                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";
278                                    if (!defined($nohtml) || $nohtml ne "1") {
279                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
280                                    }
281                            } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {
282                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";
283                                    if (!defined($nohtml) || $nohtml ne "1") {
284                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
285                                    }
286                            } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {
287                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";
288                                    if (!defined($nohtml) || $nohtml ne "1") {
289                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
290                                    }
291                            } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {
292                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";
293                                    if (!defined($nohtml) || $nohtml ne "1") {
294                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
295                                    }
296                            } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {
297                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";
298                                    if (!defined($nohtml) || $nohtml ne "1") {
299                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
300                                    }
301                            } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {
302                                    $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";
303                                    if (!defined($nohtml) || $nohtml ne "1") {
304                                            $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";
305                                    }
306                            }
307                            $classstring .= $NewClass;
308                    }
309            }
310            return $classstring;
311    }
312    
313    =head3 get_biomass
314    Definition:
315            string = FIGMODELmodel->get_biomass();
316    Description:
317            Returns data for the biomass reaction
318    =cut
319    sub get_biomass {
320            my ($self) = @_;
321    
322            if (!defined($self->{_biomass})) {
323                    my $rxntbl = $self->reaction_table();
324                    for (my $i=0; $i < $rxntbl->size(); $i++) {
325                            if ($rxntbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {
326                                    $self->{_biomass} = $rxntbl->get_row($i)->{"LOAD"}->[0];
327                                    last;
328                            }
329                    }
330            }
331    
332            return $self->get_reaction_data($self->{_biomass});
333    }
334    
335    =head3 get_reaction_data
336    Definition:
337            string = FIGMODELmodel->get_reaction_data(string::reaction ID);
338    Description:
339            Returns model reaction data
340    =cut
341    sub get_reaction_data {
342            my ($self,$reaction) = @_;
343    
344            if (!defined($self->reaction_table())) {
345                    return undef;
346            }
347            if ($reaction =~ m/^\d+$/) {
348                    $self->reaction_table()->get_row($reaction);
349            }
350            return $self->reaction_table()->get_row_by_key($reaction,"LOAD");
351    }
352    
353    =head3 get_reaction_id
354    Definition:
355            string = FIGMODELmodel->get_reaction_id(string::reaction ID);
356    Description:
357            Returns model reaction id
358    =cut
359    sub get_reaction_id {
360            my ($self,$reaction) = @_;
361    
362            my $Data = $self->get_reaction_data($reaction);
363            if (!defined($Data) || !defined($Data->{LOAD}->[0])) {
364                    return undef;
365            }
366            return $Data->{LOAD}->[0];
367    }
368    
369    =head3 get_reaction_number
370    Definition:
371            int = FIGMODELmodel->get_reaction_number();
372    Description:
373            Returns the number of reactions in the model
374    =cut
375    sub get_reaction_number {
376            my ($self) = @_;
377    
378            if (!defined($self->reaction_table())) {
379                    return 0;
380            }
381    
382            return $self->reaction_table()->size();
383    }
384    
385    =head3 reaction_table
386    Definition:
387            FIGMODELTable = FIGMODELmodel->reaction_table();
388    Description:
389            Returns FIGMODELTable with the reaction list for the model
390    =cut
391    sub reaction_table {
392            my ($self) = @_;
393    
394            if (!defined($self->{_reaction_data})) {
395                    $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id());
396            }
397    
398            return $self->{_reaction_data};
399    }
400    
401    =head3 reaction_class_table
402    Definition:
403            FIGMODELTable = FIGMODELmodel->reaction_class_table();
404    Description:
405            Returns FIGMODELTable with the reaction class data, and creates the table file  if it does not exist
406    =cut
407    sub reaction_class_table {
408            my ($self) = @_;
409    
410            if (!defined($self->{_reaction_class_table})) {
411                    if (-e $self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl") {
412                            $self->{_reaction_class_table} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["REACTION","CLASS","MEDIA"]);
413                    } else {
414                            $self->{_reaction_class_table} = FIGMODELTable->new(["REACTION","MEDIA","CLASS","MIN","MAX"],$self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",["REACTION","CLASS","MEDIA"],";","|",undef);
415                    }
416            }
417    
418            return $self->{_reaction_class_table};
419    }
420    
421    =head3 compound_class_table
422    Definition:
423            FIGMODELTable = FIGMODELmodel->compound_class_table();
424    Description:
425            Returns FIGMODELTable with the compound class data, and creates the table file  if it does not exist
426    =cut
427    sub compound_class_table {
428            my ($self) = @_;
429    
430            if (!defined($self->{_compound_class_table})) {
431                    if (-e $self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl") {
432                            $self->{_compound_class_table} = $self->figmodel()->database()->load_table($self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["COMPOUND","CLASS","MEDIA"]);
433                    } else {
434                            $self->{_compound_class_table} = FIGMODELTable->new(["COMPOUND","MEDIA","CLASS","MIN","MAX"],$self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl",["COMPOUND","CLASS","MEDIA"],";","|",undef);
435                    }
436            }
437    
438            return $self->{_compound_class_table};
439    }
440    
441    =head3 compound_table
442    Definition:
443            FIGMODELTable = FIGMODELmodel->compound_table();
444    Description:
445            Returns FIGMODELTable with the compound list for the model
446    =cut
447    sub compound_table {
448            my ($self) = @_;
449    
450            if (!defined($self->{_compound_table})) {
451                    $self->{_compound_table} = $self->figmodel()->database()->GetDBModelCompounds($self->id());
452            }
453    
454            return $self->{_compound_table};
455    }
456    
457    =head3 get_compound_data
458    Definition:
459            string = FIGMODELmodel->get_compound_data(string::compound ID);
460    Description:
461            Returns model compound data
462    =cut
463    sub get_compound_data {
464            my ($self,$compound) = @_;
465    
466            if (!defined($self->compound_table())) {
467                    return undef;
468            }
469            if ($compound =~ m/^\d+$/) {
470                    $self->compound_table()->get_row($compound);
471            }
472            return $self->compound_table()->get_row_by_key($compound,"DATABASE");
473    }
474    
475    =head3 genome
476    Definition:
477            string = FIGMODELmodel->genome();
478    Description:
479            Returns model genome
480    =cut
481    sub genome {
482            my ($self) = @_;
483            return $self->{_data}->{genome}->[0];
484    }
485    
486    =head3 source
487    Definition:
488            string = FIGMODELmodel->source();
489    Description:
490            Returns model source
491    =cut
492    sub source {
493            my ($self) = @_;
494            return $self->{_data}->{source}->[0];
495    }
496    
497    =head3 rights
498    Definition:
499            1/0 = FIGMODELmodel->rights(string::username);
500    Description:
501            Returns 1 if the input user can view the model, and zero otherwise
502    =cut
503    sub rights {
504            my ($self,$username) = @_;
505    
506            if ($self->public()) {
507                    return 1;
508            }
509            if (!defined($username) || $username eq "NONE") {
510                    return 0;
511            }
512            if (!defined($self->{_userrights}->{$username})) {
513                    if (defined($self->{_data}->{master})) {
514                            for (my $i=0; $i < @{$self->{_data}->{master}};$i++) {
515                                    $self->{_userrights}->{$self->{_data}->{master}->[$i]} = 1;
516                            }
517                    }
518                    if (defined($self->{_data}->{users})) {
519                            for (my $i=0; $i < @{$self->{_data}->{users}};$i++) {
520                                    $self->{_userrights}->{$self->{_data}->{users}->[$i]} = 1;
521                            }
522                    }
523            }
524            return $self->{_userrights}->{$username};
525    }
526    
527    =head3 public
528    Definition:
529            1/0 = FIGMODELmodel->public();
530    Description:
531            Returns 1 if the model is public, and zero otherwise
532    =cut
533    sub public {
534            my ($self) = @_;
535    
536            if (!defined($self->{_public})) {
537                    $self->{_public} = 0;
538                    if (defined($self->{_data}->{users}->[0]) && $self->{_data}->{users}->[0] eq "all") {
539                            $self->{_public} = 1;
540                    }
541            }
542            return $self->{_public};
543    }
544    
545    =head3 directory
546    Definition:
547            string = FIGMODELmodel->directory();
548    Description:
549            Returns model directory
550    =cut
551    sub directory {
552            my ($self) = @_;
553    
554            if (!defined($self->{_directory})) {
555                    my $userdirectory = "";
556                    if ($self->owner() ne "master") {
557                            $userdirectory = $self->owner()."/";
558                    }
559                    my $source = $self->source();
560                    if ($source =~ /^MGRAST/) {
561                            $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
562                    } elsif ($source =~ /^RAST/) {
563                            $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
564                    } elsif ($source =~ /^SEED/) {
565                            $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
566                    } elsif ($source =~ /^PM/) {
567                            if (length($userdirectory) == 0) {
568                                    $self->{_directory} = $self->figmodel()->config("imported model directory")->[0].$self->id()."/";
569                            } else {
570                                    $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->id()."/";
571                            }
572                    }
573            }
574    
575            return $self->{_directory};
576    }
577    
578    =head3 filename
579    Definition:
580            string = FIGMODELmodel->filename();
581    Description:
582            Returns model filename
583    =cut
584    sub filename {
585            my ($self) = @_;
586    
587            return $self->directory().$self->id().$self->selected_version().".txt";
588    }
589    
590    =head3 set_metagenome_stats
591    Definition:
592            string = FIGMODELmodel->set_metagenome_stats();
593    Description:
594            Sets the values of many model stats for a metagenome
595    =cut
596    sub set_metagenome_stats {
597            my ($self) = @_;
598    
599            $self->{_total_compounds} = 0;
600            if (defined($self->compound_table())) {
601                    $self->{_total_compounds} = $self->compound_table()->size();
602            }
603            $self->{_gene_reactions} = 0;
604            $self->{_gapfilling_reactions} = 0;
605            $self->{_model_genes} = 0;
606            $self->{_total_reactions} = 0;
607            if (defined($self->reaction_table())) {
608                    $self->{_total_reactions} = $self->reaction_table()->size();
609                    my $tbl = $self->reaction_table();
610                    my $spontaneous = 0;
611                    for (my $i=0; $i < $tbl->size(); $i++) {
612                            my $row = $tbl->get_row($i);
613                            if (!defined($row->{"ASSOCIATED PEG"}->[0]) || $row->{"ASSOCIATED PEG"}->[0] !~ m/peg/) {
614                                    if ($row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/) {
615                                            $spontaneous++;
616                                    } else {
617                                            $self->{_gapfilling_reactions}++;
618                                    }
619                            } else {
620                                    for (my $j=0; $j < @{$row->{"CONFIDENCE"}}; $j++) {
621                                            my @ecores = split(/;/,$row->{"CONFIDENCE"}->[$j]);
622                                            $self->{_model_genes} += @ecores;
623                                    }
624                            }
625                    }
626                    $self->{_gene_reactions} = $tbl->size() - $spontaneous - $self->{_gapfilling_reactions};
627            }
628    }
629    
630    =head3 version
631    Definition:
632            string = FIGMODELmodel->version();
633    Description:
634            Returns the version of the model
635    =cut
636    sub version {
637            my ($self) = @_;
638    
639            if (!defined($self->{_version})) {
640                    if (!defined($self->{_selectedversion})) {
641                            if (defined($self->stats())) {
642                                    $self->{_version} = "V".$self->stats()->{"Version"}->[0].".".$self->stats()->{"Gap fill version"}->[0];
643                            }
644                    } else {
645                            $self->{_version} = $self->{_selectedversion};
646                    }
647            }
648            return $self->{_version};
649    }
650    
651    =head3 selected_version
652    Definition:
653            string = FIGMODELmodel->selected_version();
654    Description:
655            Returns the selected version of the model
656    =cut
657    sub selected_version {
658            my ($self) = @_;
659    
660            if (!defined($self->{_selectedversion})) {
661                    return "";
662            }
663            return $self->{_selectedversion};
664    }
665    
666    =head3 modification_time
667    Definition:
668            string = FIGMODELmodel->modification_time();
669    Description:
670            Returns the selected version of the model
671    =cut
672    sub modification_time {
673            my ($self) = @_;
674            if (!defined($self->{_modification_time})) {
675                    my $stats = $self->stats();
676                    if (defined($stats)) {
677                            $self->{_modification_time} = 0;
678                            if (defined($stats->{"Build date"}->[0]) && $self->{_modification_time} < $stats->{"Build date"}->[0]) {
679                                    $self->{_modification_time} = $stats->{"Build date"}->[0];
680                            } elsif (defined($stats->{"Gap fill date"}->[0]) && $self->{_modification_time} < $stats->{"Gap fill date"}->[0]) {
681                                    $self->{_modification_time} = $stats->{"Gap fill date"}->[0];
682                            }
683                    } else {
684                            $self->{_modification_time} = 0;
685                    }
686            }
687            return $self->{_modification_time};
688    }
689    
690    =head3 gene_reactions
691    Definition:
692            string = FIGMODELmodel->gene_reactions();
693    Description:
694            Returns the number of reactions added by the gap filling
695    =cut
696    sub gene_reactions {
697            my ($self) = @_;
698    
699            if (!defined($self->{_gene_reactions})) {
700                    if ($self->source() =~ /^MGRAST/) {
701                            $self->set_metagenome_stats();
702                    } elsif (defined($self->stats())) {
703                            $self->{_gene_reactions} = $self->total_reactions() - $self->gapfilling_reactions() - $self->stats()->{'Spontaneous'}->[0] - $self->stats()->{'Growmatch reactions'}->[0] - $self->stats()->{'Biolog gap filling reactions'}->[0];
704                    }
705            }
706            return $self->{_gene_reactions};
707    }
708    
709    =head3 total_compounds
710    Definition:
711            string = FIGMODELmodel->total_compounds();
712    Description:
713            Returns the number of compounds in the model
714    =cut
715    sub total_compounds {
716            my ($self) = @_;
717    
718            if (!defined($self->{_total_compounds})) {
719                    if ($self->source() =~ /^MGRAST/) {
720                            $self->set_metagenome_stats();
721                    } elsif (defined($self->stats())) {
722                            $self->{_total_compounds} = $self->stats()->{'Metabolites'}->[0];
723                    }
724            }
725            return $self->{_total_compounds};
726    }
727    
728    =head3 gapfilling_reactions
729    Definition:
730            string = FIGMODELmodel->gapfilling_reactions();
731    Description:
732            Returns the number of reactions added by the gap filling
733    =cut
734    sub gapfilling_reactions {
735            my ($self) = @_;
736    
737            if (!defined($self->{_gapfilling_reactions})) {
738                    if ($self->source() =~ /^MGRAST/) {
739                            $self->set_metagenome_stats();
740                    } elsif (defined($self->stats())) {
741                            $self->{_gapfilling_reactions} = $self->stats()->{'Gap filling reactions'}->[0];
742                    }
743            }
744            return $self->{_gapfilling_reactions};
745    }
746    
747    =head3 total_reactions
748    Definition:
749            string = FIGMODELmodel->total_reactions();
750    Description:
751            Returns the total number of reactions in the model
752    =cut
753    sub total_reactions {
754            my ($self) = @_;
755    
756            if (!defined($self->{_total_reactions})) {
757                    if ($self->source() =~ /^MGRAST/) {
758                            $self->set_metagenome_stats();
759                    } elsif (defined($self->stats())) {
760                            $self->{_total_reactions} = $self->stats()->{'Number of reactions'}->[0];
761                    }
762            }
763            return $self->{_total_reactions};
764    }
765    
766    =head3 model_genes
767    Definition:
768            string = FIGMODELmodel->model_genes();
769    Description:
770            Returns the number of genes mapped to one or more reactions in the model
771    =cut
772    sub model_genes {
773            my ($self) = @_;
774    
775            if (!defined($self->{_model_genes})) {
776                    if ($self->source() =~ /^MGRAST/) {
777                            $self->set_metagenome_stats();
778                    } elsif (defined($self->stats())) {
779                            $self->{_model_genes} = $self->stats()->{'Genes with reactions'}->[0];
780                    }
781            }
782            return $self->{_model_genes};
783    }
784    
785    =head3 class
786    Definition:
787            string = FIGMODELmodel->class();
788    Description:
789            Returns the class of the model: gram positive, gram negative, other
790    =cut
791    sub class {
792            my ($self) = @_;
793    
794            if (!defined($self->{_class})) {
795                    if ($self->source() =~ /^MGRAST/) {
796                            $self->{_class} = "Other";
797                    } elsif (defined($self->stats())) {
798                            $self->{_class} = $self->stats()->{Class}->[0];
799                    }
800            }
801            return $self->{_class};
802    }
803    
804    =head3 taxonomy
805    Definition:
806            string = FIGMODELmodel->taxonomy();
807    Description:
808            Returns model taxonomy or biome if this is an metagenome model
809    =cut
810    sub taxonomy {
811            my ($self) = @_;
812    
813            if (!defined($self->{_taxonomy})) {
814                    my $source = $self->source();
815                    if ($source =~ /^MGRAST/) {
816                            $self->{_taxonomy} = "NA";
817                            my $mgmodeldata = $self->figmodel()->database()->mg_model_data($self->genome());
818                            if (defined($mgmodeldata)) {
819                                    $self->{_taxonomy} = $mgmodeldata->biome();
820                            }
821                    } else {
822                            $self->{_taxonomy} = $self->fig()->taxonomy_of($self->genome());
823                    }
824            }
825    
826            return $self->{_taxonomy};
827    }
828    
829    =head3 genome_size
830    Definition:
831            string = FIGMODELmodel->genome_size();
832    Description:
833            Returns size of the modeled genome in KB
834    =cut
835    sub genome_size {
836            my ($self) = @_;
837    
838            if (!defined($self->{_genome_size})) {
839                    my $source = $self->source();
840                    if ($source =~ /^MGRAST/) {
841                            $self->{_genome_size} = "NA";
842                            if (defined($self->mgdata())) {
843                                    $self->{_genome_size} = $self->mgdata()->size;
844                            }
845                    } else {
846                            $self->{_genome_size} = $self->fig()->genome_szdna($self->genome());
847                    }
848            }
849    
850            return $self->{_genome_size};
851    }
852    
853    =head3 genome_genes
854    Definition:
855            string = FIGMODELmodel->genome_genes();
856    Description:
857            Returns the number of genes in the modeled genome
858    =cut
859    sub genome_genes {
860            my ($self) = @_;
861    
862            if (!defined($self->{_genome_genes})) {
863                    my $source = $self->source();
864                    if ($source =~ /^MGRAST/) {
865                            $self->{_genome_genes} = "NA";
866                            if (defined($self->mgdata())) {
867                                    $self->{_genome_genes} = $self->mgdata()->genome_contig_count;
868                            }
869                    } elsif (defined($self->stats())) {
870                            $self->{_genome_genes} = $self->stats()->{'Total genes'}->[0];
871                    }
872            }
873    
874            return $self->{_genome_genes};
875    }
876    
877    =head3 run_default_model_predictions
878    Definition:
879            0/1::status = FIGMODELmodel->run_default_model_predictions(string::media ID);
880    Description:
881    =cut
882    sub run_default_model_predictions {
883            my ($self,$Media) = @_;
884    
885            #Assuming complete media if none is provided
886            if (!defined($Media)) {
887                    $Media = "Complete";
888            }
889    
890            #Predicting essentiality
891            my $result = $self->figmodel()->RunFBASimulation($self->id(),"SINGLEKO",undef,undef,[$self->id()],[$Media]);
892            #Checking that the table is defined and the output file exists
893            if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {
894                    $self->figmodel()->database()->print_array_to_file($self->directory()."EssentialGenes-".$self->id()."-".$Media.".tbl",[join("\n",@{$result->get_row(0)->{"ESSENTIALGENES"}})]);
895            }
896    
897            #Classifying reactions and compounds
898            my $tbl = $self->classify_model_reactions($Media);
899            $tbl->save();
900    
901            return $self->figmodel()->success();
902    }
903    
904    =head3 save_obsolete_stats
905    Definition:
906            FIGMODELmodel->save_obsolete_stats();
907    Description:
908    =cut
909    sub save_obsolete_stats {
910            my ($self) = @_;
911    
912            #checking if stats exists
913            my $stats = $self->stats();
914            if (defined($stats)) {
915                    $stats->{"Model ID"}->[0] = $self->id()."V".$stats->{"Version"}->[0].".".$stats->{"Gap fill version"}->[0];
916                    $self->figmodel()->database()->update_row("OBSOLETE MODEL STATS",$stats,"Model ID");
917                    $stats->{"Model ID"}->[0] = $self->id();
918            }
919    }
920    
921    =head3 update_stats_for_gap_filling
922    Definition:
923            {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);
924    Description:
925    =cut
926    sub update_stats_for_gap_filling {
927            my ($self,$gapfilltime) = @_;
928    
929            #preserving the stats for the now obselete model
930            $self->save_obsolete_stats();
931            my $stats = $self->update_model_stats(0);
932            $stats->{"Gap filling time"}->[0] = $gapfilltime;
933            $stats->{"Gap fill date"}->[0] = time();
934            if (!defined($stats->{"Gap fill version"}->[0])) {
935                    $stats->{"Gap fill version"}->[0] = 0;
936            }
937            $stats->{"Gap fill version"}->[0]++;
938    
939            #Updating the stats stored in the table
940            $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");
941            return $stats;
942    }
943    
944    =head3 update_stats_for_build
945    Definition:
946            {string => [string]} = FIGMODELmodel->update_stats_for_build();
947    Description:
948    =cut
949    sub update_stats_for_build {
950            my ($self) = @_;
951    
952            #preserving the stats for the now obselete model
953            $self->save_obsolete_stats();
954            my $stats = $self->update_model_stats(0);
955            $stats->{"Build date"}->[0] = time();
956            if (!defined($stats->{"Version"}->[0])) {
957                    $stats->{"Version"}->[0] = 0;
958            }
959            $stats->{"Version"}->[0]++;
960    
961            #Updating the stats stored in the table
962            $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");
963            return $stats;
964    }
965    
966    =head3 update_model_stats
967    Definition:
968            FIGMODELmodel->update_model_stats();
969    Description:
970    =cut
971    sub update_model_stats {
972            my ($self) = @_;
973    
974            #Getting reaction table
975            my $rxntbl = $self->reaction_table();
976            if (!defined($rxntbl)) {
977                    $self->figmodel()->error_message("FIGMODELmodel:update_model_stats:Could not load reaction list for ".$self->id());
978                    return undef;
979            }
980            my $cpdtbl = $self->compound_table();
981    
982            #Creating empty status row
983            my $CurrentStats = {"Genes with reactions" => [0],
984                                                     "Metabolites" => [$cpdtbl->size()],
985                                                     "Growmatch reactions" => [0],
986                                                     "Spontaneous" => [0],
987                                                     "Biolog gap filling reactions" => [0],
988                                                     "Number of reactions" => [$rxntbl->size()],
989                                                     "Transport reaction"=>[0],
990                                                     "Gap filling reactions" => [0],
991                                                     "Model ID" => [$self->id()],
992                                                     "Subsystem genes with reactions" => [0],
993                                                     "Nonsubsystem genes with reactions" => [0],
994                                                     Version => [0],
995                                                     "Gap fill version" => [0],
996                                                     Source => [$self->source()],
997                                                     "Genes with one reaction" => [0],
998                                                     "Genome ID" => [$self->genome()]};
999    
1000            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1001            if (defined($genomestats)) {
1002                    $CurrentStats->{SOURCE}->[0] = $genomestats->{SOURCE}->[0];
1003                    $CurrentStats->{"Organism name"}->[0] = $genomestats->{NAME}->[0];
1004                    $CurrentStats->{"Total genes"}->[0] = $genomestats->{"TOTAL GENES"}->[0];
1005                    $CurrentStats->{"Gram positive genes"}->[0] = $genomestats->{"GRAM POSITIVE GENES"}->[0];
1006                    $CurrentStats->{"Gram negative genes"}->[0] = $genomestats->{"GRAM NEGATIVE GENES"}->[0];
1007                    $CurrentStats->{"Class"}->[0] = $genomestats->{CLASS}->[0];
1008                    $CurrentStats->{"Genes with functions"}->[0] = $genomestats->{"GENES WITH FUNCTIONS"}->[0];
1009                    $CurrentStats->{"Subsystem genes"}->[0] = $genomestats->{"SUBSYSTEM GENES"}->[0];
1010                    $CurrentStats->{"Nonsubsystem genes"}->[0] = $genomestats->{"NON SUBSYSTEM GENES"}->[0];
1011                    if (defined($genomestats->{TAXONOMY})) {
1012                            $CurrentStats->{"Taxonomy 0"}->[0] = $genomestats->{TAXONOMY}->[0];
1013                            $CurrentStats->{"Taxonomy 1"}->[0] = $genomestats->{TAXONOMY}->[1];
1014                            $CurrentStats->{"Taxonomy 2"}->[0] = $genomestats->{TAXONOMY}->[2];
1015                            $CurrentStats->{"Taxonomy 3"}->[0] = $genomestats->{TAXONOMY}->[3];
1016                            $CurrentStats->{"Taxonomy 4"}->[0] = $genomestats->{TAXONOMY}->[4];
1017                            $CurrentStats->{"Taxonomy 5"}->[0] = $genomestats->{TAXONOMY}->[5];
1018                    }
1019            }
1020    
1021            #Transfering build, version, and gap fill data from existing stats
1022            if (defined($self->stats())) {
1023                    $CurrentStats->{Version}->[0] = $self->stats()->{Version}->[0];
1024                    $CurrentStats->{"Gap fill version"}->[0] = $self->stats()->{"Gap fill version"}->[0];
1025                    $CurrentStats->{"Gap filling time"}->[0] = $self->stats()->{"Gap filling time"}->[0];
1026                    $CurrentStats->{"Gap fill date"}->[0] = $self->stats()->{"Gap fill date"}->[0];
1027                    $CurrentStats->{"Build date"}->[0] = $self->stats()->{"Build date"}->[0];
1028            }
1029    
1030            my %GeneHash;
1031            my %NonpegHash;
1032            my %CompoundHash;
1033            for (my $i=0; $i < $rxntbl->size(); $i++) {
1034                    my $Row = $rxntbl->get_row($i);
1035                    if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
1036                            my $ReactionRow = $self->figmodel()->get_reaction($Row->{"LOAD"}->[0]);
1037                            if (defined($ReactionRow->{"EQUATION"}->[0])) {
1038                                    #Checking for extracellular metabolites which indicate that this is a transporter
1039                                    if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {
1040                                            $CurrentStats->{"Transport reaction"}->[0]++;
1041                                    }
1042                            }
1043                            #Identifying spontaneous/biolog/gapfilling/gene associated reactions
1044                            if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {
1045                                    $CurrentStats->{"Biolog gap filling reactions"}->[0]++;
1046                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {
1047                                    $CurrentStats->{"Spontaneous"}->[0]++;
1048                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {
1049                                    $CurrentStats->{"Gap filling reactions"}->[0]++;
1050                            } else {
1051                                    foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
1052                                            my @GeneList = split(/\+/,$GeneSet);
1053                                            foreach my $Gene (@GeneList) {
1054                                                    if ($Gene =~ m/(peg\.\d+)/) {
1055                                                            $GeneHash{$1} = 1;
1056                                                    } else {
1057                                                            $NonpegHash{$Gene} = 1;
1058                                                    }
1059                                            }
1060                                    }
1061                            }
1062                    }
1063            }
1064            my @genes = keys(%GeneHash);
1065            my @othergenes = keys(%NonpegHash);
1066            $CurrentStats->{"Genes with reactions"}->[0] = @genes + @othergenes;
1067    
1068            #Updating the stats stored in the table
1069            $self->figmodel()->database()->update_row("MODEL STATS",$CurrentStats,"Model ID");
1070            $self->{_stats} = $CurrentStats;
1071            return $CurrentStats;
1072    }
1073    
1074    =head3 GapFillModel
1075    Definition:
1076            (success/fail) = FIGMODELmodel->GapFillModel();
1077    Description:
1078            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.
1079            Before running the gap filling, the existing model is backup in the same directory with the current version numbers appended.
1080            If the model has been gap filled previously, the previous gap filling reactions are removed prior to running the gap filling again.
1081    =cut
1082    sub GapFillModel {
1083            my ($self,$donotclear) = @_;
1084    
1085            #Setting status of model to gap filling
1086            my $OrganismID = $self->genome();
1087            $self->set_status(1,"Auto completion running");
1088            my $UniqueFilename = $self->figmodel()->filename();
1089            my $StartTime = time();
1090    
1091            #Archiving the existing model
1092            $self->ArchiveModel();
1093    
1094            #Removing any gapfilling reactions that may be currently present in the model
1095            if (!defined($donotclear) || $donotclear != 1) {
1096                    my $ModelTable = $self->reaction_table();
1097                    for (my $i=0; $i < $ModelTable->size(); $i++) {
1098                            if (!defined($ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0]) || $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/) {
1099                                    $ModelTable->delete_row($i);
1100                                    $i--;
1101                            }
1102                    }
1103                    $ModelTable->save();
1104            }
1105    
1106            #Calling the MFAToolkit to run the gap filling optimization
1107            my $MinimalMediaTable = $self->figmodel()->database()->GetDBTable("MINIMAL MEDIA TABLE");
1108            my $Media = "Complete";
1109            if (defined($MinimalMediaTable->get_row_by_key($self->genome(),"Organism"))) {
1110                    $Media = $MinimalMediaTable->get_row_by_key($self->genome(),"Organism")->{"Minimal media"}->[0];
1111                    #Loading media, changing bounds, saving media as a test media
1112                    my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);
1113                    for (my $i=0; $i < $MediaTable->size(); $i++) {
1114                            if ($MediaTable->get_row($i)->{"Min"}->[0] < 0) {
1115                                    $MediaTable->get_row($i)->{"Min"}->[0] = -10000;
1116                            }
1117                            if ($MediaTable->get_row($i)->{"Max"}->[0] > 0) {
1118                                    $MediaTable->get_row($i)->{"Max"}->[0] = 10000;
1119                            }
1120                    }
1121                    $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1122                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0},"GapFill".$self->id().".log",undef));
1123                    unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1124            } else {
1125                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],undef,"GapFill".$self->id().".log",undef));
1126            }
1127    
1128            #Parse the solutions file for the model and get the reaction list from it
1129            my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);
1130            #Clearing the mfatoolkit output and log file
1131            $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");
1132            if (!defined($SolutionData) || $SolutionData->size() == 0) {
1133                    $self->set_status(1,"Auto completion failed: auto completion solution not found");
1134                    $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");
1135                    return $self->fail();
1136            }
1137            $SolutionData->save("/home/chenry/Solution.txt");
1138    
1139            #Looking for the last printed recursive MILP solution
1140            for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
1141                    if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {
1142                            my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);
1143                            print "Solution:".$AllSolutions."\n";
1144                            my @TempThree = split(/\|/,$AllSolutions);
1145                            if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {
1146                                    my @TempFour = split(/,/,$1);
1147                                    my $DirectionList;
1148                                    my $ReactionList;
1149                                    for (my $j=0; $j < @TempFour; $j++) {
1150                                            my $ID = "";
1151                                            my $Sign = "=>";
1152                                            if ($TempFour[$j] =~ m/\-(rxn\d\d\d\d\d)/) {
1153                                                    $ID = $1;
1154                                                    $Sign = "<=";
1155                                            } elsif ($TempFour[$j] =~ m/\+(rxn\d\d\d\d\d)/) {
1156                                                    $ID = $1;
1157                                            }
1158                                            if ($ID ne "") {
1159                                                    if ($self->figmodel()->reversibility_of_reaction($ID) ne $Sign) {
1160                                                            $Sign = "<=>";
1161                                                    }
1162                                                    push(@{$DirectionList},$Sign);
1163                                                    push(@{$ReactionList},$ID);
1164                                            }
1165                                    }
1166                                    print "Solution:".$TempThree[0]."\n";
1167                                    $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1);
1168                            }
1169                            last;
1170                    }
1171            }
1172            #Updating model stats with gap filling results
1173            my $ElapsedTime = time() - $StartTime;
1174            $self->figmodel()->ClearDBModel($self->id(),1);
1175            #Determining why each gap filling reaction was added
1176            $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);
1177            if (!defined($donotclear) || $donotclear != 1) {
1178                    $self->figmodel()->AddBiologTransporters($self->id());
1179                    if ($self->id() !~ m/MGRast/) {
1180                            $self->update_stats_for_gap_filling($ElapsedTime);
1181                    }
1182            } else {
1183                    $self->update_model_stats();
1184            }
1185            #Printing the updated SBML file
1186            $self->PrintSBMLFile();
1187            $self->PrintModelLPFile($self->id());
1188            $self->set_status(1,"Auto completion successfully finished");
1189            $self->run_default_model_predictions();
1190    
1191            return $self->success();
1192    }
1193    
1194    =head3 datagapfill
1195    Definition:
1196            success()/fail() = FIGMODELmodel->datagapfill();
1197    Description:
1198            Run gapfilling on the input run specifications
1199    =cut
1200    sub datagapfill {
1201            my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;
1202            my $UniqueFilename = $self->figmodel()->filename();
1203            if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {
1204                    print "Gap filling specs:\n".join("\n",@{$GapFillingRunSpecs})."\n\n";
1205                    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));
1206    
1207                    #Checking that the solution exists
1208                    if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {
1209                            $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");
1210                            $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
1211                            return undef;
1212                    }
1213                    my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);
1214                    if (defined($TansferFileSuffix)) {
1215                            system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
1216                            system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingReport.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
1217                    }
1218    
1219                    #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted
1220                    $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
1221                    return $GapFillingRunSpecs;
1222            }
1223            if (defined($TansferFileSuffix)) {
1224                    $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
1225            }
1226            return undef;
1227    }
1228    
1229    =head3 TestSolutions
1230    Definition:
1231            $model->TestSolutions($ModelID,$NumProcessors,$ProcessorIndex,$GapFill);
1232    Description:
1233    Example:
1234    =cut
1235    
1236    sub TestSolutions {
1237            my ($self,$OriginalErrorFilename,$GapFillResultTable) = @_;
1238            #Getting the filename
1239            my $UniqueFilename = $self->figmodel()->filename();
1240            #Reading in the original error matrix which has the headings for the original model simulation
1241            my $OriginalErrorData;
1242            if (!defined($OriginalErrorFilename) || !-e $self->directory().$OriginalErrorFilename) {
1243                    my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
1244                    $OriginalErrorData = [$HeadingVector,$Errorvector];
1245            } else {
1246                    $OriginalErrorData = $self->figmodel()->database()->load_single_column_file($self->directory().$OriginalErrorFilename,"");
1247            }
1248            my $HeadingHash;
1249            my @HeadingArray = split(/;/,$OriginalErrorData->[0]);
1250            my @OrigErrorArray = split(/;/,$OriginalErrorData->[1]);
1251            for (my $i=0; $i < @HeadingArray; $i++) {
1252                    my @SubArray = split(/:/,$HeadingArray[$i]);
1253                    $HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]} = $i;
1254            }
1255            #Scanning through the gap filling solutions
1256            my $TempVersion = "V".$UniqueFilename;
1257            my $ErrorMatrixLines;
1258            for (my $i=0; $i < $GapFillResultTable->size(); $i++) {
1259                    print "Starting problem solving ".$i."\n";
1260                    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];
1261                    #Integrating solution into test model
1262                    my $ReactionArray;
1263                    my $DirectionArray;
1264                    my @ReactionList = split(/,/,$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0]);
1265                    my %SolutionHash;
1266                    for (my $k=0; $k < @ReactionList; $k++) {
1267                            if ($ReactionList[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) {
1268                                    my $Reaction = $2;
1269                                    my $Sign = $1;
1270                                    if (defined($SolutionHash{$Reaction})) {
1271                                            $SolutionHash{$Reaction} = "<=>";
1272                                    } elsif ($Sign eq "-") {
1273                                            $SolutionHash{$Reaction} = "<=";
1274                                    } elsif ($Sign eq "+") {
1275                                            $SolutionHash{$Reaction} = "=>";
1276                                    } else {
1277                                            $SolutionHash{$Reaction} = $Sign;
1278                                    }
1279                            }
1280                    }
1281                    @ReactionList = keys(%SolutionHash);
1282                    for (my $k=0; $k < @ReactionList; $k++) {
1283                            push(@{$ReactionArray},$ReactionList[$k]);
1284                            push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});
1285                    }
1286                    print "Integrating solution!\n";
1287                    $self->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1);
1288                    my $model = $self->get_model($self->id().$TempVersion);
1289                    $model->PrintModelLPFile();
1290                    #Running the model against all available experimental data
1291                    print "Running test model!\n";
1292                    my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
1293    
1294                    @HeadingArray = split(/;/,$HeadingVector);
1295                    my @ErrorArray = @OrigErrorArray;
1296                    my @TempArray = split(/;/,$Errorvector);
1297                    for (my $j=0; $j < @HeadingArray; $j++) {
1298                            my @SubArray = split(/:/,$HeadingArray[$j]);
1299                            $ErrorArray[$HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]}] = $TempArray[$j];
1300                    }
1301                    $ErrorLine .= ";".$FalsePostives."/".$FalseNegatives.";".join(";",@ErrorArray);
1302                    push(@{$ErrorMatrixLines},$ErrorLine);
1303                    print "Finishing problem solving ".$i."\n";
1304            }
1305            #Clearing out the test model
1306            if (-e $self->directory().$self->id().$TempVersion.".txt") {
1307                    unlink($self->directory().$self->id().$TempVersion.".txt");
1308                    unlink($self->directory()."SimulationOutput".$self->id().$TempVersion.".txt");
1309            }
1310            return $ErrorMatrixLines;
1311    }
1312    
1313    =head3 status
1314    Definition:
1315            int::model status = FIGMODELmodel->status();
1316    Description:
1317            Returns the current status of the SEED model associated with the input genome ID.
1318            model status = 1: model exists
1319            model status = 0: model is being built
1320            model status = -1: model does not exist
1321            model status = -2: model build failed
1322    =cut
1323    sub status {
1324            my ($self) = @_;
1325            return $self->{_data}->{status}->[0];
1326    }
1327    
1328    =head3 message
1329    Definition:
1330            string::model message = FIGMODELmodel->message();
1331    Description:
1332            Returns a message associated with the models current status
1333    =cut
1334    sub message {
1335            my ($self) = @_;
1336            return $self->{_data}->{message}->[0];
1337    }
1338    
1339    =head3 set_status
1340    Definition:
1341            (success/fail) = FIGMODELmodel->set_status(int::new status,string::status message);
1342    Description:
1343            Changes the current status of the SEED model
1344            new status = 1: model exists
1345            new status = 0: model is being built
1346            new status = -1: model does not exist
1347            new status = -2: model build failed
1348    =cut
1349    sub set_status {
1350            my ($self,$NewStatus,$Message) = @_;
1351    
1352            #Getting the model row from the MODELS table
1353            $self->{_data}->{status}->[0] = $NewStatus;
1354            $self->{_data}->{message}->[0] = $Message;
1355            $self->figmodel()->database()->update_row("MODELS",$self->{_data},"id");
1356            return $self->config("SUCCESS")->[0];
1357    }
1358    
1359    =head3 CreateSingleGenomeReactionList
1360    Definition:
1361            FIGMODEL->CreateSingleGenomeReactionList();
1362    Description:
1363            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.
1364    Example:
1365    =cut
1366    
1367    sub CreateSingleGenomeReactionList {
1368            my ($self,$RunGapFilling) = @_;
1369    
1370            #Getting genome stats
1371            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1372            my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
1373            if (!defined($FeatureTable)) {
1374                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome features could not be accessed!");
1375                    return $self->fail();
1376            }
1377            #Checking that the number of genes exceeds the minimum size
1378            if ($FeatureTable->size() < $self->config("minimum genome size for modeling")->[0]) {
1379                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome rejected as too small for modeling!");
1380                    return $self->fail();
1381            }
1382            #Setting up needed variables
1383            my $OriginalModelTable = undef;
1384            #Locking model status table
1385            my $ModelTable = $self->figmodel()->database()->LockDBTable("MODELS");
1386            my $Row = $ModelTable->get_row_by_key($self->id(),"id");
1387            if (!defined($Row) || !defined($Row->{status}->[0])) {
1388                    $self->figmodel()->database()->UnlockDBTable("MODELS");
1389                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." has no row in models table!");
1390                    return $self->fail();
1391            } elsif ($Row->{status}->[0] == 0) {
1392                    $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList:Model is already being built. Canceling current build.");
1393                    $self->figmodel()->database()->UnlockDBTable("MODELS");
1394                    return $self->fail();
1395            }elsif ($Row->{status}->[0] == 1) {
1396                    $OriginalModelTable = $self->reaction_table();
1397                    $self->ArchiveModel();
1398                    $Row->{status}->[0] = 0;
1399                    $Row->{message}->[0] = "Rebuilding preliminary reconstruction";
1400            } else {
1401                    $Row->{status}->[0] = 0;
1402                    $Row->{message}->[0] = "Preliminary reconstruction";
1403            }
1404            #Updating the status table
1405            $self->figmodel()->database()->save_table($ModelTable);
1406            $self->figmodel()->database()->UnlockDBTable("MODELS");
1407            #Sorting GenomeData by gene location on the chromosome
1408            $FeatureTable->sort_rows("MIN LOCATION");
1409            my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);
1410            my %LastGenePosition;
1411            my $GeneRoles;
1412            for (my $j=0; $j < $FeatureTable->size(); $j++) {
1413                    my $CurrentRow = $FeatureTable->get_row($j);
1414                    #"ID","ALIASES","MIN LOCATION","MAX LOCATION","ROLES","SUBSYSTEMS","SUBSYSTEM CLASS"
1415                    if (defined($CurrentRow)) {
1416                            my $GeneID = $CurrentRow->{"ID"}->[0];
1417                            if ($GeneID =~ m/(peg\.\d+)/) {
1418                                    $GeneID = $1;
1419                            }
1420                            foreach my $Role (@{$CurrentRow->{"ROLES"}}) {
1421                                    if ($self->figmodel()->role_is_valid($Role) != 0) {
1422                                            push(@{$GeneRoles->{$GeneID}},$Role);
1423                                            my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);
1424                                            if ($ReactionHashArrayRef != 0) {
1425                                                    foreach my $Mapping (@{$ReactionHashArrayRef}) {
1426                                                            if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
1427                                                                    if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) {
1428                                                                            if ($Mapping->{"MASTER"}->[0] eq 1) {
1429                                                                                    #Creating a complex if consecutive genes have been assigned to the same reaction
1430                                                                                    $ComplexHash->{$Mapping->{"REACTION"}->[0]}->{$Mapping->{"COMPLEX"}->[0]}->{$Role}->{$GeneID} = 1;
1431                                                                                    if (!defined($LastGenePosition{$Mapping->{"REACTION"}->[0]})) {
1432                                                                                            $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j;
1433                                                                                            push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID);
1434                                                                                    } elsif (($j-$LastGenePosition{$Mapping->{"REACTION"}->[0]}) < 3 && $LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) {
1435                                                                                            my $CurrentComplex = pop(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}});
1436                                                                                            push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$CurrentComplex."+".$GeneID);
1437                                                                                    } elsif ($LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) {
1438                                                                                            push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID);
1439                                                                                    }
1440                                                                                    $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j;
1441                                                                                    #Adding a subsystem for the reaction
1442                                                                                    if ($self->figmodel()->subsystem_is_valid($Mapping->{"SUBSYSTEM"}->[0]) == 1) {
1443                                                                                            ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},$Mapping->{"SUBSYSTEM"}->[0]);
1444                                                                                            if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 1) {
1445                                                                                                    if ($Mapping->{"SOURCE"}->[0] =~ m/Hope\sFiles/) {
1446                                                                                                            $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 1;
1447                                                                                                    } elsif ($Mapping->{"SOURCE"}->[0] =~ m/SEED/) {
1448                                                                                                            $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 2;
1449                                                                                                    } elsif (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) {
1450                                                                                                            $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3;
1451                                                                                                    }
1452                                                                                            }
1453                                                                                    }
1454                                                                                    #Handling confidence
1455                                                                                    if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) {
1456                                                                                            if ($Mapping->{"SOURCE"}->[0] =~ m/MATT/) {
1457                                                                                                    $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3;
1458                                                                                            } elsif ($Mapping->{"SOURCE"}->[0] =~ m/CHRIS/) {
1459                                                                                                    $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 4;
1460                                                                                            } else {
1461                                                                                                    $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 5;
1462                                                                                            }
1463                                                                                    }
1464                                                                                    #Parsing sources
1465                                                                                    ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},split(/\|/,$Mapping->{"SOURCE"}->[0]));
1466                                                                            } else {
1467                                                                                    push(@{$SuggestedMappings},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role);
1468                                                                            }
1469                                                                    } else {
1470                                                                            push(@{$UnrecognizedReactions},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role);
1471                                                                    }
1472                                                            }
1473                                                    }
1474                                            }
1475                                    }
1476                            }
1477                    }
1478            }
1479    
1480            #Creating nonadjacent complexes
1481            my @ReactionList = keys(%{$ReactionHash});
1482            foreach my $Reaction (@ReactionList) {
1483                    #If multiple genes are assigned to the reaction, we check if they should should be in a complex
1484                    if (@{$ReactionHash->{$Reaction}->{"GENES"}} > 0 && defined($ComplexHash->{$Reaction})) {
1485                            my $GeneArray;
1486                            foreach my $Complex (keys(%{$ComplexHash->{$Reaction}})) {
1487                                    my %ComplexComponents;
1488                                    foreach my $CurrentGeneSet (@{$ReactionHash->{$Reaction}->{"GENES"}}) {
1489                                            my @GeneList = split(/\+/,$CurrentGeneSet);
1490                                            my %RoleHash;
1491                                            foreach my $Gene (@GeneList) {
1492                                                    foreach my $Role (@{$GeneRoles->{$Gene}}) {
1493                                                            if (defined($ComplexHash->{$Reaction}->{$Complex}->{$Role})) {
1494                                                                    $RoleHash{$Role} = 1;
1495                                                            }
1496                                                    }
1497                                            }
1498                                            if (keys(%RoleHash) > 0) {
1499                                                    if (!defined($ComplexComponents{join("|",sort(keys(%RoleHash)))})) {
1500                                                            my @RoleList = keys(%RoleHash);
1501                                                            my @ComplexList = keys(%ComplexComponents);
1502                                                            foreach my $ComplexSet (@ComplexList) {
1503                                                                    my @RoleList = split(/\|/,$ComplexSet);
1504                                                                    my $Match = 0;
1505                                                                    foreach my $SingleRole (@RoleList) {
1506                                                                            if (defined($RoleHash{$SingleRole})) {
1507                                                                                    $Match = 1;
1508                                                                                    last;
1509                                                                            }
1510                                                                    }
1511                                                                    if ($Match == 1) {
1512                                                                            foreach my $SingleRole (@RoleList) {
1513                                                                                    $RoleHash{$SingleRole} = 1
1514                                                                            }
1515                                                                            push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},@{$ComplexComponents{$ComplexSet}});
1516                                                                            delete $ComplexComponents{$ComplexSet};
1517                                                                    }
1518                                                            }
1519                                                    }
1520                                                    push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},$CurrentGeneSet);
1521                                            }
1522                                    }
1523                                    my @Position;
1524                                    my @Options;
1525                                    my $Count = 0;
1526                                    foreach my $RoleSet (keys(%ComplexComponents)) {
1527                                            push(@Position,0);
1528                                            push(@{$Options[$Count]},@{$ComplexComponents{$RoleSet}});
1529                                            $Count++;
1530                                    }
1531                                    my $Done = 0;
1532                                    $Count = 0;
1533                                    my $NewRelationship;
1534                                    while($Done == 0) {
1535                                            #Creating complex with current indecies
1536                                            $NewRelationship->[$Count] = $Options[0]->[$Position[0]];
1537                                            for (my $i=1; $i < @Position; $i++) {
1538                                                    $NewRelationship->[$Count] .= "+".$Options[$i]->[$Position[$i]];
1539                                            }
1540                                            $NewRelationship->[$Count] = join("+",$self->figmodel()->remove_duplicates(split(/\+/,$NewRelationship->[$Count])));
1541                                            $Count++;
1542                                            #Iterating indecies
1543                                            my $CurrentIndex = 0;
1544                                            while($CurrentIndex >= 0) {
1545                                                    if ($CurrentIndex >= @Position) {
1546                                                            $CurrentIndex = -1000;
1547                                                    } elsif ($Position[$CurrentIndex]+1 == @{$Options[$CurrentIndex]}) {
1548                                                            $Position[$CurrentIndex] = -1;
1549                                                            $CurrentIndex++;
1550                                                    } else {
1551                                                            $Position[$CurrentIndex]++;
1552                                                            $CurrentIndex--;
1553                                                    }
1554                                            }
1555                                            if ($CurrentIndex == -1000) {
1556                                                    $Done = 1;
1557                                            }
1558                                    }
1559                                    push(@{$GeneArray},@{$NewRelationship});
1560                            }
1561                            @{$ReactionHash->{$Reaction}->{"GENES"}} = $self->figmodel()->remove_duplicates(@{$GeneArray});
1562                    }
1563            }
1564    
1565            #Getting the reaction table
1566            my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");
1567    
1568            #Creating the model reaction table
1569            my $NewModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n");
1570            @ReactionList = keys(%{$ReactionHash});
1571            foreach my $ReactionID (@ReactionList) {
1572                    #Getting the thermodynamic reversibility from the database
1573                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1574                    my $Subsystem = "NONE";
1575                    if (defined($ReactionHash->{$ReactionID}->{"SUBSYSTEMS"})) {
1576                            $Subsystem = join("|",@{$ReactionHash->{$ReactionID}->{"SUBSYSTEMS"}});
1577                    }
1578                    my $Source = "NONE";
1579                    if (defined($ReactionHash->{$ReactionID}->{"SOURCE"})) {
1580                            $Source = join("|",@{$ReactionHash->{$ReactionID}->{"SOURCE"}});
1581                    }
1582                    $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"]});
1583            }
1584    
1585            #Adding the spontaneous and universal reactions
1586            foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
1587                    #Getting the thermodynamic reversibility from the database
1588                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1589                    #Checking if the reaction is already in the model
1590                    if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1591                            $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]});
1592                    }
1593            }
1594            foreach my $ReactionID (@{$self->config("universal reactions")}) {
1595                    #Getting the thermodynamic reversibility from the database
1596                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1597                    #Checking if the reaction is already in the model
1598                    if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1599                            $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]});
1600                    }
1601            }
1602    
1603            #Checking if a biomass reaction already exists
1604            my $BiomassReactionRow = $self->figmodel()->database()->get_row_by_key("BIOMASS TABLE",$self->id(),"MODELS");
1605            if (!defined($BiomassReactionRow)) {
1606                    $BiomassReactionRow = $self->BuildSpecificBiomassReaction();
1607                    if (!defined($BiomassReactionRow)) {
1608                            $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
1609                            $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");
1610                            return $self->fail();
1611                    }
1612            }
1613            my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};
1614            push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);
1615    
1616            #Adding biomass reactions to the model table
1617            foreach my $BOFReaction (@{$ReactionList}) {
1618                    #Getting the thermodynamic reversibility from the database
1619                    my $Directionality = $self->figmodel()->reversibility_of_reaction($BOFReaction);
1620                    #Checking if the reaction is already in the model
1621                    if (!defined($NewModelTable->get_row_by_key($BOFReaction,"LOAD"))) {
1622                            if ($BOFReaction =~ m/bio/) {
1623                                    $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["BOF"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [1],"REFERENCE" => ["Biomass objective function"],"NOTES" => ["NONE"]});
1624                            } else {
1625                                    $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["INITIAL GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Initial gap filling"],"NOTES" => ["NONE"]});
1626                            }
1627                    }
1628            }
1629    
1630            #Completing any incomplete reactions sets
1631            my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS");
1632            for (my $i=0; $i < $ReactionSetTable->size(); $i++) {
1633                    if (defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1634                            foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) {
1635                                    if (!defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1636                                            #Getting the thermodynamic reversibility from the database
1637                                            my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction);
1638                                            $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"]});
1639                                    }
1640                            }
1641                    }
1642            }
1643    
1644            #Now we compare the model to the previous model to determine if any differences exist that aren't gap filling reactions
1645            if (defined($OriginalModelTable)) {
1646                    my $PerfectMatch = 1;
1647                    my $ReactionCount = 0;
1648                    for (my $i=0; $i < $OriginalModelTable->size(); $i++) {
1649                            #We only check that nongapfilling reactions exist in the new model
1650                            if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] !~ m/GAP/ || $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "INITIAL GAP FILLING") {
1651                                    $ReactionCount++;
1652                                    my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD");
1653                                    if (defined($Row)) {
1654                                            #We check that the reaction directionality is identical
1655                                            if ($Row->{"DIRECTIONALITY"}->[0] ne $OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]) {
1656                                                    if (defined($OriginalModelTable->get_row($i)->{"NOTES"}->[0]) && $OriginalModelTable->get_row($i)->{"NOTES"}->[0] =~ m/Directionality\sswitched\sfrom\s([^\s])/) {
1657                                                            if ($1 ne $Row->{"DIRECTIONALITY"}->[0]) {
1658                                                                    print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$1." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n";
1659                                                                    $PerfectMatch = 0;
1660                                                                    last;
1661                                                            }
1662                                                    } else {
1663                                                            print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n";
1664                                                            $PerfectMatch = 0;
1665                                                            last;
1666                                                    }
1667                                            }
1668                                            #We check that the genes assigned to the reaction are identical
1669                                            if ($PerfectMatch == 1 && @{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}} != @{$Row->{"ASSOCIATED PEG"}}) {
1670                                                    print "Gene associatation mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}}." vs ".@{$Row->{"ASSOCIATED PEG"}}."\n";
1671                                                    $PerfectMatch = 0;
1672                                                    last;
1673                                            }
1674                                            if ($PerfectMatch == 1) {
1675                                                    my @GeneSetOne = sort(@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}});
1676                                                    my @GeneSetTwo = sort(@{$Row->{"ASSOCIATED PEG"}});
1677                                                    for (my $j=0; $j < @GeneSetOne; $j++) {
1678                                                            if ($GeneSetOne[$j] ne $GeneSetTwo[$j]) {
1679                                                                    print "Gene mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$GeneSetOne[$j]." vs ".$GeneSetTwo[$j]."\n";
1680                                                                    $PerfectMatch = 0;
1681                                                                    $i = $OriginalModelTable->size();
1682                                                                    last;
1683                                                            }
1684                                                    }
1685                                            }
1686                                    } else {
1687                                            print "Original model contains an extra reaction:".$OriginalModelTable->get_row($i)->{"LOAD"}->[0]."\n";
1688                                            $PerfectMatch = 0;
1689                                            last;
1690                                    }
1691                            }
1692                    }
1693                    if ($PerfectMatch == 1 && $ReactionCount == $NewModelTable->size()) {
1694                            #Bailing out of function as the model has not changed
1695                            $self->set_status(1,"rebuild canceled because model has not changed");
1696                            return $self->success();
1697                    }
1698            }
1699    
1700            #Saving the new model to file
1701            $self->set_status(1,"Preliminary reconstruction complete");
1702            $self->figmodel()->database()->save_table($NewModelTable);
1703            $self->{_reaction_data} = $NewModelTable;
1704            #Clearing the previous model from the cache
1705            $self->figmodel()->ClearDBModel($self->id(),1);
1706            #Updating the model stats table
1707            $self->update_stats_for_build();
1708            $self->PrintSBMLFile();
1709    
1710            #Adding model to gapfilling queue
1711            if (defined($RunGapFilling) && $RunGapFilling == 1) {
1712                    $self->set_status(1,"Autocompletion queued");
1713                    $self->figmodel()->add_job_to_queue("gapfillmodel?".$self->id(),"QSUB","cplex",$self->owner(),"BACK");
1714            }
1715            return $self->success();
1716    }
1717    
1718    =head3 CreateMetaGenomeReactionList
1719    Definition:
1720            (success/fail) = FIGMODELmodel->CreateMetaGenomeReactionList();
1721    Description:
1722            This is the code called to create or update the reaction list for a metgenome model
1723    =cut
1724    
1725    sub CreateMetaGenomeReactionList {
1726            my ($self) = @_;
1727    
1728            #Checking if the metagenome file exists
1729            if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {
1730                    $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
1731            }
1732            #Loading metagenome data
1733            my $MGRASTData = $self->figmodel()->database()->load_multiple_column_file($self->config("raw MGRAST directory")->[0].$self->genome().".summary","\t");
1734            if (!defined($MGRASTData)) {
1735                    $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
1736            }
1737    
1738            #Setting up needed variables
1739            my $OriginalModelTable = undef;
1740    
1741            #Checking status
1742            if ($self->status() < 0) {
1743                    $self->set_status(0,"Preliminary reconstruction");
1744            } elsif ($self->status() == 0) {
1745                    $self->error_message("FIGMODEL->CreateModelReactionList:Model is already being built. Canceling current build.");
1746                    return $self->fail();
1747            } else {
1748                    $OriginalModelTable = $self->reaction_table();
1749                    $self->ArchiveModel();
1750                    $self->set_status(0,"Rebuilding preliminary reconstruction");
1751            }
1752    
1753            #Getting the reaction table
1754            my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");
1755            #Creating model table
1756            my $ModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n");
1757            for (my $i=0; $i < @{$MGRASTData};$i++) {
1758                    #MD5,PEG,number of sims,role,sim e-scores
1759                    my $Role = $MGRASTData->[$i]->[3];
1760                    my $MD5 = $MGRASTData->[$i]->[0];
1761                    my $peg = $MGRASTData->[$i]->[1];
1762                    my $sims = $MGRASTData->[$i]->[4];
1763                    $sims =~ s/;/,/g;
1764    
1765                    #Checking for subsystems
1766                    my $GeneSubsystems = $self->figmodel()->subsystems_of_role($Role);
1767                    #Checking if there are reactions associated with this role
1768                    my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);
1769                    if ($ReactionHashArrayRef != 0) {
1770                            foreach my $Mapping (@{$ReactionHashArrayRef}) {
1771                                    if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
1772                                            if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) {
1773                                                    if ($Mapping->{"MASTER"}->[0] eq 1) {
1774                                                            #Checking if the reaction is already in the model
1775                                                            my $ReactionRow = $ModelTable->get_row_by_key($Mapping->{"REACTION"}->[0],"LOAD");
1776                                                            if (!defined($ReactionRow)) {
1777                                                                    $ReactionRow = {"LOAD" => [$Mapping->{"REACTION"}->[0]],"DIRECTIONALITY" => [$self->figmodel()->reversibility_of_reaction($Mapping->{"REACTION"}->[0])],"COMPARTMENT" => ["c"]};
1778                                                                    $ModelTable->add_row($ReactionRow);
1779                                                            }
1780                                                            push(@{$ReactionRow->{"ASSOCIATED PEG"}},substr($peg,4));
1781                                                            push(@{$ReactionRow->{"REFERENCE"}},$MD5.":".$Role);
1782                                                            push(@{$ReactionRow->{"CONFIDENCE"}},$sims);
1783                                                            if (defined($GeneSubsystems)) {
1784                                                                    push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});
1785                                                            }
1786                                                    }
1787                                            }
1788                                    }
1789                            }
1790                    }
1791            }
1792    
1793            #Adding the spontaneous and universal reactions
1794            foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
1795                    #Getting the thermodynamic reversibility from the database
1796                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1797                    #Checking if the reaction is already in the model
1798                    if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1799                            $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]});
1800                    }
1801            }
1802            foreach my $ReactionID (@{$self->config("universal reactions")}) {
1803                    #Getting the thermodynamic reversibility from the database
1804                    my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID);
1805                    #Checking if the reaction is already in the model
1806                    if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) {
1807                            $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]});
1808                    }
1809            }
1810    
1811            #Completing any incomplete reactions sets
1812            my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS");
1813            for (my $i=0; $i < $ReactionSetTable->size(); $i++) {
1814                    if (defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1815                            foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) {
1816                                    if (!defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) {
1817                                            #Getting the thermodynamic reversibility from the database
1818                                            my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction);
1819                                            $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"]});
1820                                    }
1821                            }
1822                    }
1823            }
1824    
1825            #Clearing the previous model from the cache
1826            $self->figmodel()->ClearDBModel($self->id(),1);
1827            $ModelTable->save();
1828    
1829            return $self->success();
1830    }
1831    
1832    =head3 ArchiveModel
1833    Definition:
1834            (success/fail) = FIGMODELmodel->ArchiveModel();
1835    Description:
1836            This function archives the specified model in the model directory with the current version numbers appended.
1837            This function is used to preserve old versions of models prior to overwriting so new versions may be compared with old versions.
1838    =cut
1839    sub ArchiveModel {
1840            my ($self) = @_;
1841    
1842            #Making sure the model exists
1843            if (!defined($self->stats())) {
1844                    $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Could not find specified ".$self->id()." in database!");
1845                    return $self->fail();
1846            }
1847    
1848            #Checking that the model file exists
1849            if (!(-e $self->filename())) {
1850                    $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");
1851                    return $self->fail();
1852            }
1853    
1854            #Copying the model file
1855            system("cp ".$self->filename()." ".$self->directory().$self->id().$self->version().".txt");
1856    }
1857    
1858    =head3 PrintModelDataToFile
1859    Definition:
1860            (success/fail) = FIGMODELmodel->PrintModelDataToFile();
1861    Description:
1862            This function uses the MFAToolkit to print out all of the compound and reaction data for the input model.
1863            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
1864            function. The LoadModel function for example would not give you this data.
1865    =cut
1866    sub PrintModelDataToFile {
1867            my($self) = @_;
1868    
1869            #Running the MFAToolkit on the model file
1870            my $OutputIndex = $self->figmodel()->filename();
1871            my $Command = $self->config("MFAToolkit executable")->[0]." parameterfile ../Parameters/Printing.txt resetparameter output_folder ".$OutputIndex.'/ LoadCentralSystem "'.$self->filename().'"';
1872            system($Command);
1873    
1874            #Copying the model file printed by the toolkit out of the output directory and into the model directory
1875            if (!-e $self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().".txt") {
1876                    $self->figmodel()->error_message("New model file not created due to an error. Check that the input modelfile exists.");
1877                    $self->figmodel()->cleardirectory($OutputIndex);
1878                    return $self->fail();
1879            }
1880    
1881            $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().'.txt" "'.$self->directory().$self->id().$self->selected_version().'Data.txt"';
1882            system($Command);
1883            $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex.'/ErrorLog0.txt" "'.$self->directory().'ModelErrors.txt"';
1884            system($Command);
1885            $self->figmodel()->cleardirectory($OutputIndex);
1886            return $self->success();
1887    }
1888    
1889    =head2 Analysis Functions
1890    
1891    =head3 run_microarray_analysis
1892    Definition:
1893            int::status = FIGMODEL->run_microarray_analysis(string::media,string::job id,string::gene calls);
1894    Description:
1895            Runs microarray analysis attempting to turn off genes that are inactive in the microarray
1896    =cut
1897    sub run_microarray_analysis {
1898            my ($self,$media,$jobid,$index,$genecall) = @_;
1899            $genecall =~ s/_/:/g;
1900            $genecall =~ s/\//;/g;
1901            #print "\n\n".$genecall."\n\n";
1902            my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($jobid,$self->id(),$media,["MFA","MicroarrayAssertions"],{"Microarray assertions" => $self->id().";".$index.";".$genecall,"MFASolver" => "CPLEX","Network output location" => "/scratch/"},"MicroarrayAnalysis-".$jobid.".txt",undef,$self->selected_version());
1903            #print $command."\n";
1904            system($command);
1905            #system("/home/chenry/Software/MFAToolkitRepository/Linux/mfatoolkit resetparameter \"user bounds filename\" \"Carbon-D-Glucose.txt\" resetparameter output_folder \"32749.354149.0/\" resetparameter \"Microarray assertions\" \"Seed83333.1;0;peg.3130:-1;peg.4035:-1\" resetparameter \"MFASolver\" \"CPLEX\" resetparameter \"Network output location\" \"/scratch/\" LoadCentralSystem \"/vol/model-dev/MODEL_DEV_DB/Models/83333.1/Seed83333.1.txt\" > \"/vol/model-dev/MODEL_DEV_DB/ReactionDB/log/MicroarrayAnalysis-32749.354149.0.txt\"");
1906    }
1907    
1908    =head3 find_minimal_pathways
1909    Definition:
1910            int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
1911    Description:
1912            Runs microarray analysis attempting to turn off genes that are inactive in the microarray
1913    =cut
1914    sub find_minimal_pathways {
1915            my ($self,$media,$objective,$solutionnum) = @_;
1916    
1917            #Setting default media
1918            if (!defined($media)) {
1919                    $media = "Complete";
1920            }
1921    
1922            #Setting default solution number
1923            if (!defined($solutionnum)) {
1924                    $solutionnum = "5";
1925            }
1926    
1927            #Translating objective
1928            my $objectivestring;
1929            if ($objective eq "ALL") {
1930                    #Getting the list of universal building blocks
1931                    my $buildingblocks = $self->config("universal building blocks");
1932                    my @objectives = keys(%{$buildingblocks});
1933                    #Getting the nonuniversal building blocks
1934                    my $otherbuildingblocks = $self->config("nonuniversal building blocks");
1935                    my @array = keys(%{$otherbuildingblocks});
1936                    if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
1937                            my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
1938                            if (defined($equation)) {
1939                                    for (my $i=0; $i < @array; $i++) {
1940                                            if (CORE::index($equation,$array[$i]) > 0) {
1941                                                    push(@objectives,$array[$i]);
1942                                            }
1943                                    }
1944                            }
1945                    }
1946                    for (my $i=0; $i < @objectives; $i++) {
1947                            $self->find_minimal_pathways($media,$objectives[$i])
1948                    }
1949            } elsif ($objective eq "ENERGY") {
1950                    $objectivestring = "MAX;FLUX;rxn00062;c;1";
1951            } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
1952                    $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
1953            }
1954    
1955            #Setting additional exchange fluxes
1956            my $additionalexchange = $self->config("default exchange fluxes")->[0];
1957            if ($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") {
1958                    $additionalexchange .= ";cpd15666[c]:0:100";
1959            } elsif ($objective eq "cpd11493") {
1960                    $additionalexchange .= ";cpd12370[c]:0:100";
1961            } elsif ($objective eq "cpd00166") {
1962                    $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
1963            }
1964    
1965            #Running MFAToolkit
1966            my $filename = $self->figmodel()->filename();
1967            my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["MFA"],{"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());
1968            system($command);
1969    
1970            #Loading problem report
1971            my $results = $self->figmodel()->LoadProblemReport($filename);
1972            #Clearing output
1973            $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
1974            if (!defined($results)) {
1975                    return;
1976            }
1977    
1978            #Parsing output
1979            my @Array;
1980            my $row = $results->get_row(1);
1981            if (defined($row->{"Notes"}->[0])) {
1982                    $_ = $row->{"Notes"}->[0];
1983                    @Array = /\d+:([^\|]+)\|/g;
1984            }
1985    
1986            # Storing data in a figmodel table
1987            my $TableObject;
1988            if (-e $self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt") {
1989                    $TableObject->load_table($self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",";","|",0,["OBJECTIVE"]);
1990            } else {
1991                    $TableObject = FIGMODELTable->new(["OBJECTIVE","REACTIONS"],$self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",["OBJECTIVE"],";","|",undef);
1992            }
1993            my $tablerow = $TableObject->get_row_by_key($objective,"OBJECTIVE",1);
1994            push(@{$tablerow->{"REACTIONS"}},@Array);
1995            $TableObject->save();
1996    }
1997    
1998    =head3 calculate_growth
1999    Definition:
2000            string::growth = FIGMODELmodel->calculate_growth(string:media);
2001    Description:
2002            Calculating growth in the input media
2003    =cut
2004    sub calculate_growth {
2005            my ($self,$Media) = @_;
2006            my $UniqueFilename = $self->figmodel()->filename();
2007            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],{"optimize metabolite production if objective is zero" => 1},$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version()));
2008            my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);
2009            my $Result;
2010            if (defined($ProblemReport)) {
2011                    my $Row = $ProblemReport->get_row(0);
2012                    if (defined($Row) && defined($Row->{"Objective"}->[0])) {
2013                            if ($Row->{"Objective"}->[0] < 0.00000001) {
2014                                    $Result = "NOGROWTH:".$Row->{"Individual metabolites with zero production"}->[0];
2015                            } else {
2016                                    $Result = $Row->{"Objective"}->[0];
2017                            }
2018                    }
2019            }
2020            return $Result;
2021    }
2022    
2023    =head3 classify_model_reactions
2024    Definition:
2025            (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);
2026    Description:
2027            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.
2028            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".
2029            Possible values for "Class" include:
2030            1.) Positive: these reactions are essential in the forward direction.
2031            2.) Negative: these reactions are essential in the reverse direction.
2032            3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.
2033            4.) Negative variable: these reactions are nonessential, but they only ever proceed in the reverse direction.
2034            5.) Variable: these reactions are nonessential and proceed in the forward or reverse direction.
2035            6.) Blocked: these reactions never carry any flux at all in the media condition tested.
2036            7.) Dead: these reactions are disconnected from the network.
2037    =cut
2038    sub classify_model_reactions {
2039            my ($self,$Media) = @_;
2040    
2041            #Getting unique file for printing model output
2042            my $UniqueFilename = $self->figmodel()->filename();
2043            #Running the MFAToolkit
2044            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()));
2045            #Reading in the output bounds file
2046            my ($ReactionTB,$CompoundTB,$DeadCompounds,$DeadEndCompounds,$DeadReactions);
2047            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") {
2048                    $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]);
2049            }
2050            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt") {
2051                    $CompoundTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt",";","|",1,["DATABASE ID"]);
2052            }
2053            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt") {
2054                    $DeadReactions = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt",""));
2055            }
2056            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt") {
2057                    $DeadCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt",""));
2058            }
2059            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt") {
2060                    $DeadEndCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt",""));
2061            }
2062            if (!defined($ReactionTB) && !defined($CompoundTB)) {
2063                    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";
2064                    return (undef,undef);
2065            }
2066    
2067            #Clearing output
2068            $self->figmodel()->clearing_output($UniqueFilename,"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
2069            #Creating the table objects that will hold the results of the reaction classification
2070            my $rxnclasstable = $self->reaction_class_table();
2071            my $cpdclasstable = $self->compound_class_table();
2072            #Loading the compound table
2073            if (defined($CompoundTB)) {
2074                    for (my $i=0; $i < $CompoundTB->size(); $i++) {
2075                            my $Row = $CompoundTB->get_row($i);
2076                            if (defined($Row->{"DATABASE ID"})) {
2077                                    #Getting the compound row
2078                                    my $CpdRow = $cpdclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0].$Row->{COMPARTMENT}->[0],"COMPOUND",1);
2079                                    #Setting row values
2080                                    my $Max = 0;
2081                                    my $Min = 0;
2082                                    my $Class = "Unknown";
2083                                    if (defined($DeadCompounds) && defined($DeadCompounds->{$Row->{"DATABASE ID"}->[0]})) {
2084                                            $Class = "Dead";
2085                                    } elsif (defined($DeadEndCompounds) && defined($DeadEndCompounds->{$Row->{"DATABASE ID"}->[0]})) {
2086                                            $Class = "Deadend";
2087                                    } elsif (defined($Row->{"Min DRAIN_FLUX"}) && defined($Row->{"Max DRAIN_FLUX"}) && $Row->{"Min DRAIN_FLUX"}->[0] ne "1e+07") {
2088                                            $Max = $Row->{"Max DRAIN_FLUX"}->[0];
2089                                            $Min = $Row->{"Min DRAIN_FLUX"}->[0];
2090                                            if ($Row->{"Min DRAIN_FLUX"}->[0] > 0.00000001) {
2091                                                    $Class = "Positive";
2092                                            } elsif ($Row->{"Max DRAIN_FLUX"}->[0] < -0.00000001) {
2093                                                    $Class = "Negative";
2094                                            } elsif ($Row->{"Min DRAIN_FLUX"}->[0] < -0.00000001) {
2095                                                    if ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) {
2096                                                            $Class = "Variable";
2097                                                    } else {
2098                                                            $Max = 0;
2099                                                            $Class = "Negative variable";
2100                                                    }
2101                                            } elsif ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) {
2102                                                    $Min = 0;
2103                                                    $Class = "Positive variable";
2104                                            } else {
2105                                                    $Min = 0;
2106                                                    $Max = 0;
2107                                                    $Class = "Blocked";
2108                                            }
2109                                    }
2110                                    my $index = 0;
2111                                    if (defined($CpdRow->{MEDIA})) {
2112                                            for (my $i=0; $i < @{$CpdRow->{MEDIA}};$i++) {
2113                                                    $index++;
2114                                                    if ($CpdRow->{MEDIA}->[$i] eq $Media) {
2115                                                            $index = $i;
2116                                                            last;
2117                                                    }
2118                                            }
2119                                    }
2120                                    $CpdRow->{MIN}->[$index] = $Min;
2121                                    $CpdRow->{MAX}->[$index] = $Max;
2122                                    $CpdRow->{CLASS}->[$index] = $Class;
2123                                    $CpdRow->{MEDIA}->[$index] = $Media;
2124                            }
2125                    }
2126                    $cpdclasstable->save();
2127            }
2128            if (defined($ReactionTB)) {
2129                    for (my $i=0; $i < $ReactionTB->size(); $i++) {
2130                            my $Row = $ReactionTB->get_row($i);
2131                            if (defined($Row->{"DATABASE ID"})) {
2132                                    #Getting the compound row
2133                                    my $Compartment = "c";
2134                                    if (defined($Row->{COMPARTMENT}->[0])) {
2135                                            $Compartment = $Row->{COMPARTMENT}->[0];
2136                                    }
2137                                    my $RxnRow = $rxnclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0],"REACTION",1);
2138                                    my $Max = 0;
2139                                    my $Min = 0;
2140                                    my $Class = "Unknown";
2141                                    if (defined($DeadReactions) && defined($DeadReactions->{$Row->{"DATABASE ID"}->[0]})) {
2142                                            $Class = "Dead";
2143                                    } elsif (defined($Row->{"Min FLUX"}) && defined($Row->{"Max FLUX"})) {
2144                                            $Max = $Row->{"Max FLUX"}->[0];
2145                                            $Min = $Row->{"Min FLUX"}->[0];
2146                                            if ($Row->{"Min FLUX"}->[0] > 0.00000001) {
2147                                                    $Class = "Positive";
2148                                            } elsif ($Row->{"Max FLUX"}->[0] < -0.00000001) {
2149                                                    $Class = "Negative";
2150                                            } elsif ($Row->{"Min FLUX"}->[0] < -0.00000001) {
2151                                                    if ($Row->{"Max FLUX"}->[0] > 0.00000001) {
2152                                                            $Class = "Variable";
2153                                                    } else {
2154                                                            $Max = 0;
2155                                                            $Class = "Negative variable";
2156                                                    }
2157                                            } elsif ($Row->{"Max FLUX"}->[0] > 0.00000001) {
2158                                                    $Min = 0;
2159                                                    $Class = "Positive variable";
2160                                            } else {
2161                                                    $Min = 0;
2162                                                    $Max = 0;
2163                                                    $Class = "Blocked";
2164                                            }
2165                                    }
2166                                    my $index = 0;
2167                                    if (defined($RxnRow->{MEDIA})) {
2168                                            for (my $i=0; $i < @{$RxnRow->{MEDIA}};$i++) {
2169                                                    $index++;
2170                                                    if ($RxnRow->{MEDIA}->[$i] eq $Media) {
2171                                                            $index = $i;
2172                                                            last;
2173                                                    }
2174                                            }
2175                                    }
2176                                    $RxnRow->{MIN}->[$index] = $Min;
2177                                    $RxnRow->{MAX}->[$index] = $Max;
2178                                    $RxnRow->{CLASS}->[$index] = $Class;
2179                                    $RxnRow->{MEDIA}->[$index] = $Media;
2180                            }
2181                    }
2182                    $rxnclasstable->save();
2183            }
2184            return ($rxnclasstable,$cpdclasstable);
2185    }
2186    
2187    =head3 RunAllStudiesWithDataFast
2188    Definition:
2189            (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);
2190    Description:
2191            Simulates every experimental condition currently available for the model.
2192    =cut
2193    
2194    sub RunAllStudiesWithDataFast {
2195            my ($self,$Experiment,$PrintResults) = @_;
2196    
2197            #Printing lp and key file for model
2198            if (!-e $self->directory()."FBA-".$self->id().$self->selected_version().".lp") {
2199                    $self->PrintModelLPFile();
2200            }
2201            my $UniqueFilename = $self->figmodel()->filename();
2202    
2203            #Determing the simulations that need to be run
2204            my $ExperimentalDataTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment);
2205            #Creating the table of jobs to submit
2206            my $JobArray = $self->GetSimulationJobTable($ExperimentalDataTable,$Experiment,$UniqueFilename);
2207            #Printing the job file
2208            if (!-d $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/") {
2209                    system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/");
2210            }
2211            $JobArray->save();
2212    
2213            #Running simulations
2214            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");
2215            print "simulation complete for ".$self->id().$self->selected_version()."\n";
2216            #Parsing the results
2217            my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
2218            if (!defined($Results)) {
2219                    $self->figmodel()->error_message("FIGMODELmodel:RunAllStudiesWithDataFast:Could not find simulation results: ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
2220                    return undef;
2221            }
2222            my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector,$SimulationResults) = $self->EvaluateSimulationResults($Results,$ExperimentalDataTable);
2223            #Printing results to file
2224            $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");
2225            $self->figmodel()->clearing_output($UniqueFilename);
2226    
2227            return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector);
2228    }
2229    
2230    =head3 GetSimulationJobTable
2231    Definition:
2232            my $JobTable = $model->GetSimulationJobTable($Experiment,$PrintResults,$Version);
2233    Description:
2234    =cut
2235    
2236    sub GetSimulationJobTable {
2237            my ($self,$SimulationTable,$Experiment,$Folder) = @_;
2238    
2239            #Determing the simulations that need to be run
2240            if (!defined($SimulationTable)) {
2241                    $SimulationTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment);
2242                    if (!defined($SimulationTable)) {
2243                            return undef;
2244                    }
2245            }
2246    
2247            #Creating the job table
2248            my $JobTable = $self->figmodel()->CreateJobTable($Folder);
2249            for (my $i=0; $i < $SimulationTable->size(); $i++) {
2250                    if ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Gene\sKO/) {
2251                            my $Row = $JobTable->get_row_by_key("Gene KO","LABEL",1);
2252                            $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
2253                    } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Media\sgrowth/) {
2254                            my $Row = $JobTable->get_row_by_key("Growth phenotype","LABEL",1);
2255                            $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
2256                    } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Interval\sKO/) {
2257                            my $Row = $JobTable->get_row_by_key($SimulationTable->get_row($i)->{"Heading"}->[0],"LABEL",1);
2258                            $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1);
2259                            $JobTable->add_data($Row,"GENE KO",$SimulationTable->get_row($i)->{"Experiment type"}->[0],1);
2260                    }
2261            }
2262    
2263            #Filling in model specific elements of the job table
2264            for (my $i=0; $i < $JobTable->size(); $i++) {
2265                    if ($JobTable->get_row($i)->{"LABEL"}->[0] =~ m/Gene\sKO/) {
2266                            $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "SINGLEKO";
2267                            $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 1;
2268                    } else {
2269                            $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "GROWTH";
2270                            $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 0;
2271                    }
2272                    $JobTable->get_row($i)->{"LP FILE"}->[0] = $self->directory()."FBA-".$self->id().$self->selected_version();
2273                    $JobTable->get_row($i)->{"MODEL"}->[0] = $self->directory().$self->id().$self->selected_version().".txt";
2274                    $JobTable->get_row($i)->{"SAVE FLUXES"}->[0] = 0;
2275            }
2276    
2277            return $JobTable;
2278    }
2279    
2280    =head3 EvaluateSimulationResults
2281    Definition:
2282            (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);
2283    Description:
2284            Compares simulation results with experimental data to produce a table indicating where predictions are incorrect.
2285    =cut
2286    
2287    sub EvaluateSimulationResults {
2288            my ($self,$Results,$ExperimentalDataTable) = @_;
2289    
2290            #Comparing experimental results with simulation results
2291            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);
2292            my $FalsePostives = 0;
2293            my $FalseNegatives = 0;
2294            my $CorrectNegatives = 0;
2295            my $CorrectPositives = 0;
2296            my @Errorvector;
2297            my @HeadingVector;
2298            my $ReactionKOWithGeneHash;
2299            for (my $i=0; $i < $Results->size(); $i++) {
2300                    if ($Results->get_row($i)->{"LABEL"}->[0] eq "Gene KO") {
2301                            if (defined($Results->get_row($i)->{"REACTION KO WITH GENES"})) {
2302                                    for (my $j=0; $j < @{$Results->get_row($i)->{"REACTION KO WITH GENES"}}; $j++) {
2303                                            my @Temp = split(/:/,$Results->get_row($i)->{"REACTION KO WITH GENES"}->[$j]);
2304                                            if (defined($Temp[1]) && length($Temp[1]) > 0) {
2305                                                    $ReactionKOWithGeneHash->{$Temp[0]} = $Temp[1];
2306                                            }
2307                                    }
2308                            }
2309                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] == 0) {
2310                                    for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) {
2311                                            my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading");
2312                                            if (defined($Row)) {
2313                                                    my $KOReactions = "none";
2314                                                    if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) {
2315                                                            $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]};
2316                                                    }
2317                                                    push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
2318                                                    my $Status = "Unknown";
2319                                                    if ($Row->{"Growth"}->[0] > 0) {
2320                                                            $Status = "False negative";
2321                                                            $FalseNegatives++;
2322                                                            push(@Errorvector,3);
2323                                                    } else {
2324                                                            $Status = "False positive";
2325                                                            $FalsePostives++;
2326                                                            push(@Errorvector,2);
2327                                                    }
2328                                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
2329                                            }
2330                                    }
2331                            } else {
2332                                    for (my $j=0; $j < @{$Results->get_row($i)->{"ESSENTIALGENES"}}; $j++) {
2333                                            #print $j."\t".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]."\n";
2334                                            my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j],"Heading");
2335                                            if (defined($Row)) {
2336                                                    my $KOReactions = "none";
2337                                                    if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]})) {
2338                                                            $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]};
2339                                                    }
2340                                                    push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
2341                                                    my $Status = "Unknown";
2342                                                    if ($Row->{"Growth"}->[0] > 0) {
2343                                                            $Status = "False negative";
2344                                                            $FalseNegatives++;
2345                                                            push(@Errorvector,3);
2346                                                    } else {
2347                                                            $Status = "Correct negative";
2348                                                            $CorrectNegatives++;
2349                                                            push(@Errorvector,1);
2350                                                    }
2351                                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
2352                                            }
2353                                    }
2354                                    for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) {
2355                                            my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading");
2356                                            if (defined($Row)) {
2357                                                    my $KOReactions = "none";
2358                                                    if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) {
2359                                                            $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]};
2360                                                    }
2361                                                    push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions);
2362                                                    my $Status = "Unknown";
2363                                                    if ($Row->{"Growth"}->[0] > 0) {
2364                                                            $Status = "Correct positive";
2365                                                            $CorrectPositives++;
2366                                                            push(@Errorvector,0);
2367                                                    } else {
2368                                                            $Status = "False positive";
2369                                                            $FalsePostives++;
2370                                                            push(@Errorvector,2);
2371                                                    }
2372                                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]});
2373                                            }
2374                                    }
2375                            }
2376                    } elsif ($Results->get_row($i)->{"LABEL"}->[0] eq "Growth phenotype") {
2377                            my $Row = $ExperimentalDataTable->get_row_by_key("Media growth:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"MEDIA"}->[0],"Heading");
2378                            if (defined($Row)) {
2379                                    push(@HeadingVector,$Row->{"Heading"}->[0].":none");
2380                                    my $Status = "Unknown";
2381                                    if ($Row->{"Growth"}->[0] > 0) {
2382                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2383                                                    $Status = "Correct positive";
2384                                                    $CorrectPositives++;
2385                                                    push(@Errorvector,0);
2386                                            } else {
2387                                                    $Status = "False negative";
2388                                                    $FalseNegatives++;
2389                                                    push(@Errorvector,3);
2390                                            }
2391                                    } else {
2392                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2393                                                    $Status = "False positive";
2394                                                    $FalsePostives++;
2395                                                    push(@Errorvector,2);
2396                                            } else {
2397                                                    $Status = "Correct negative";
2398                                                    $CorrectNegatives++;
2399                                                    push(@Errorvector,1);
2400                                            }
2401                                    }
2402                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Media growth"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Media"}->[0]],"Reactions knocked out" => ["none"]});
2403                            }
2404                    } elsif ($Results->get_row($i)->{"LABEL"}->[0] =~ m/Interval\sKO/ && defined($Results->get_row($i)->{"KOGENES"}->[0])) {
2405                            my $Row = $ExperimentalDataTable->get_row_by_key($Results->get_row($i)->{"LABEL"}->[0],"Heading");
2406                            if (defined($Row)) {
2407                                    my $Status = "Unknown";
2408                                    if ($Row->{"Growth"}->[0] > 0) {
2409                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2410                                                    $Status = "Correct positive";
2411                                                    $CorrectPositives++;
2412                                                    push(@Errorvector,0);
2413                                            } else {
2414                                                    $Status = "False negative";
2415                                                    $FalseNegatives++;
2416                                                    push(@Errorvector,3);
2417                                            }
2418                                    } else {
2419                                            if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) {
2420                                                    $Status = "False positive";
2421                                                    $FalsePostives++;
2422                                                    push(@Errorvector,2);
2423                                            } else {
2424                                                    $Status = "Correct negative";
2425                                                    $CorrectNegatives++;
2426                                                    push(@Errorvector,1);
2427                                            }
2428                                    }
2429                                    $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Interval KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => ["none"]});
2430                            }
2431                    }
2432            }
2433    
2434            return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,join(";",@Errorvector),join(";",@HeadingVector),$SimulationResults);
2435    }
2436    
2437    =head3 InspectSolution
2438    Definition:
2439            $model->InspectSolution(string::gene knocked out,string::media condition,[string]::list of reactions);
2440    Description:
2441    =cut
2442    
2443    sub InspectSolution {
2444            my ($self,$GeneKO,$Media,$ReactionList) = @_;
2445    
2446            #Getting a directory for the results
2447            my $UniqueFilename = $self->figmodel()->filename();
2448            system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/");
2449            my $TempVersion = "V".$UniqueFilename;
2450    
2451            #Setting gene ko to none if no genes are to be knocked out
2452            if ($GeneKO !~ m/^peg\./) {
2453                    $GeneKO = "none";
2454            }
2455    
2456            #Implementing the input solution in the test model
2457            my $ReactionArray;
2458            my $DirectionArray;
2459            my %SolutionHash;
2460            for (my $k=0; $k < @{$ReactionList}; $k++) {
2461                    if ($ReactionList->[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) {
2462                            my $Reaction = $2;
2463                            my $Sign = $1;
2464                            if (defined($SolutionHash{$Reaction})) {
2465                                    $SolutionHash{$Reaction} = "<=>";
2466                            } elsif ($Sign eq "-") {
2467                                    $SolutionHash{$Reaction} = "<=";
2468                            } elsif ($Sign eq "+") {
2469                                    $SolutionHash{$Reaction} = "=>";
2470                            } else {
2471                                    $SolutionHash{$Reaction} = $Sign;
2472                            }
2473                    }
2474            }
2475            my @TempList = keys(%SolutionHash);
2476            for (my $k=0; $k < @TempList; $k++) {
2477                    push(@{$ReactionArray},$TempList[$k]);
2478                    push(@{$DirectionArray},$SolutionHash{$TempList[$k]});
2479            }
2480    
2481            print "Integrating solution!\n";
2482            $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"SolutionInspection",1,1);
2483    
2484            #Printing lp and key file for model
2485            $self->PrintModelLPFile();
2486    
2487            #Running FBA on the test model
2488            my $JobTable = $self->figmodel()->CreateJobTable($UniqueFilename);
2489            $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]});
2490            $JobTable->save();
2491    
2492            #Running simulations
2493            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");
2494    
2495            #Parsing the results
2496            my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
2497            if (!defined($Results)) {
2498                    $self->figmodel()->error_message("FIGMODELmodel:InspectSolution:Could not load problem report ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt");
2499                    return undef;
2500            }
2501    
2502            #Making sure that the model grew with all reactions present
2503            my $Found = 0;
2504            for (my $i=0; $i < $Results->size(); $i++) {
2505                    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) {
2506                            $Found = 1;
2507                    }
2508            }
2509            if ($Found == 0) {
2510                    print "Solution no longer valid\n";
2511                    return undef;
2512            }
2513    
2514            #Making sure all of the reactions added are still necessary
2515            my $FinalReactionList;
2516            for (my $k=0; $k < $Results->size(); $k++) {
2517                    if (defined($Results->get_row($k)->{"KOGENES"}->[0]) && $Results->get_row($k)->{"KOGENES"}->[0] eq $GeneKO) {
2518                            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) {
2519                                    push(@{$FinalReactionList},$Results->get_row($k)->{"KOREACTIONS"}->[0]);
2520                            }
2521                    }
2522            }
2523    
2524            #Deleting extra files created
2525            unlink($self->directory()."FBA-".$self->id().$TempVersion.".lp");
2526            unlink($self->directory()."FBA-".$self->id().$TempVersion.".key");
2527            unlink($self->directory().$self->id().$TempVersion.".txt");
2528    
2529            #Deleting the test model and the MFA folder
2530            $self->figmodel()->clearing_output($UniqueFilename);
2531    
2532            return $FinalReactionList;
2533    }
2534    
2535    =head3 GapFillingAlgorithm
2536    
2537    Definition:
2538            FIGMODELmodel->GapFillingAlgorithm();
2539    
2540    Description:
2541            This is a wrapper for running the gap filling algorithm on any model in the database.
2542            The algorithm performs a gap filling for any false negative prediction of the avialable experimental data.
2543            This function is threaded to improve efficiency: one thread does nothing but using the MFAToolkit to fill gaps for every false negative prediction.
2544            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.
2545            This function prints two important output files in the Model directory:
2546            1.) GapFillingOutput.txt: this is a summary of the results of the gap filling analysis
2547            2.) GapFillingErrorMatrix.txt: this lists the correct and incorrect predictions for each gapfilling solution implemented in a test model.
2548    =cut
2549    
2550    sub GapFillingAlgorithm {
2551            my ($self) = @_;
2552    
2553            #First the input model version and model filename should be simulated and the false negatives identified
2554            my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
2555    
2556            #Getting the filename
2557            my $UniqueFilename = $self->figmodel()->filename();
2558    
2559            #Printing the original performance vector
2560            $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-OPEM".".txt",[$HeadingVector,$Errorvector]);
2561    
2562            my $PreviousGapFilling;
2563            if (-e $self->directory().$self->id().$self->selected_version()."-GFS.txt") {
2564                    #Backing up the old solution file
2565                    system("cp ".$self->directory().$self->id().$self->selected_version()."-GFS.txt ".$self->directory().$self->id().$self->selected_version()."-OldGFS.txt");
2566                    unlink($self->directory().$self->id().$self->selected_version()."-GFS.txt");
2567            }
2568            if (-e $self->directory().$self->id().$self->selected_version()."-OldGFS.txt") {
2569                    #Reading in the solution file from the previous gap filling if it exists
2570                    $PreviousGapFilling = $self->figmodel()->database()->load_table($self->directory().$self->id().$self->selected_version()."-OldGFS.txt",";",",",0,["Experiment"]);
2571            }
2572    
2573            #Now we use the simulation output to make the gap filling run data
2574            my @Errors = split(/;/,$Errorvector);
2575            my @Headings = split(/;/,$HeadingVector);
2576            my $GapFillingRunSpecs = "";
2577            my $Count = 0;
2578            my $RescuedPreviousResults;
2579            my $RunCount = 0;
2580            my $SolutionExistedCount = 0;
2581            my $AcceptedSolutions = 0;
2582            my $RejectedSolutions = 0;
2583            my $NoExistingSolutions = 0;
2584            for (my $i=0; $i < @Errors; $i++) {
2585                    if ($Errors[$i] == 3) {
2586                            my @HeadingDataArray = split(/:/,$Headings[$i]);
2587                            if ($HeadingDataArray[2] !~ m/^peg\./ || $HeadingDataArray[3] ne "none") {
2588                                    my $SolutionFound = 0;
2589                                    if (defined($PreviousGapFilling) && defined($PreviousGapFilling->get_row_by_key($HeadingDataArray[2],"Experiment"))) {
2590                                            my @Rows = $PreviousGapFilling->get_rows_by_key($HeadingDataArray[2],"Experiment");
2591                                            for (my $j=0; $j < @Rows; $j++) {
2592                                                    if ($HeadingDataArray[2] =~ m/^peg\./) {
2593                                                            my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"});
2594                                                            if (defined($ReactionList)) {
2595                                                                    print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n";
2596                                                                    $SolutionFound++;
2597                                                                    push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList}));
2598                                                                    $AcceptedSolutions++;
2599                                                            } else {
2600                                                                    $RejectedSolutions++;
2601                                                            }
2602                                                    } else {
2603                                                            my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"});
2604                                                            if (defined($ReactionList)) {
2605                                                                    print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n";
2606                                                                    $SolutionFound++;
2607                                                                    push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList}));
2608                                                                    $AcceptedSolutions++;
2609                                                            } else {
2610                                                                    $RejectedSolutions++;
2611                                                            }
2612                                                    }
2613                                            }
2614                                    } else {
2615                                            $NoExistingSolutions++;
2616                                    }
2617                                    if ($SolutionFound == 0) {
2618                                            $RunCount++;
2619                                            if (length($GapFillingRunSpecs) > 0) {
2620                                                    $GapFillingRunSpecs .= ";";
2621                                            }
2622                                            $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].".txt:".$HeadingDataArray[3];
2623                                    } else {
2624                                            $SolutionExistedCount++;
2625                                    }
2626                            }
2627                            $Count++;
2628                    }
2629            }
2630    
2631            #Updating the growmatch progress table
2632            my $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1);
2633            $Row->{"INITIAL FP"}->[0] = $FalsePostives;
2634            $Row->{"INITIAL FN"}->[0] = $FalseNegatives;
2635            $Row->{"GF TIMING"}->[0] = time()."-";
2636            $Row->{"FN WITH SOL"}->[0] = $FalseNegatives-$NoExistingSolutions;
2637            $Row->{"FN WITH ACCEPTED SOL"}->[0] = $SolutionExistedCount;
2638            $Row->{"TOTAL ACCEPTED GF SOL"}->[0] = $AcceptedSolutions;
2639            $Row->{"TOTAL REJECTED GF SOL"}->[0] = $RejectedSolutions;
2640            $Row->{"FN WITH NO SOL"}->[0] = $NoExistingSolutions+$RejectedSolutions;
2641            $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM");
2642    
2643            #Running the gap filling once to correct all false negative errors
2644            my $SolutionsFound = 0;
2645            my $GapFillingArray;
2646            push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));
2647            $self->datagapfill($GapFillingArray);
2648    
2649            if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {
2650                    #Printing previous solutions to GFS file
2651                    $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",$RescuedPreviousResults,1);
2652                    $SolutionsFound = 1;
2653            }
2654    
2655            #Recording the finishing of the gapfilling
2656            $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1);
2657            $Row->{"GF TIMING"}->[0] .= time();
2658            $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM");
2659    
2660            if ($SolutionsFound == 1) {
2661                    #Scheduling solution testing
2662                    $self->figmodel()->add_job_to_queue("testsolutions?".$self->id().$self->selected_version()."?-1?GF","QSUB","fast","master","BACK");
2663            } else {
2664                    $self->figmodel()->error_message("No false negative predictions found. Data gap filling not necessary!");
2665            }
2666    
2667            return $self->success();
2668    }
2669    
2670    =head3 BuildSpecificBiomassReaction
2671    Definition:
2672            FIGMODELmodel->BuildSpecificBiomassReaction();
2673    Description:
2674    =cut
2675    sub BuildSpecificBiomassReaction {
2676            my ($self) = @_;
2677    
2678            my $biomassrxn;
2679            my $OrganismID = $self->genome();
2680            #Checking for a biomass override
2681            if (defined($self->config("biomass reaction override")->{$OrganismID})) {
2682                    $biomassrxn = $self->config("biomass reaction override")->{$OrganismID};
2683                    print "Overriding biomass template and selecting ".$biomassrxn." for ".$OrganismID.".\n";
2684            } else {#Creating biomass reaction from the template
2685                    #Getting the genome stats
2686                    my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
2687                    my $Class = $genomestats->{CLASS}->[0];
2688                    my $Name = $genomestats->{NAME}->[0];
2689    
2690                    #Checking for phoenix variants
2691                    my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");
2692                    my $Phoenix = 0;
2693                    my @Rows = $PhoenixVariantTable->get_rows_by_key($OrganismID,"GENOME");
2694                    my $VariantHash;
2695                    for (my $i=0; $i < @Rows; $i++) {
2696                            $Phoenix = 1;
2697                            if (defined($Rows[$i]->{"SUBSYSTEM"}) && defined($Rows[$i]->{"VARIANT"})) {
2698                                    $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];
2699                            }
2700                    }
2701    
2702                    #Collecting genome data
2703                    my $RoleHash;
2704                    my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
2705                    for (my $i=0; $i < $FeatureTable->size(); $i++) {
2706                            for (my $j=0; $j < @{$FeatureTable->get_row($i)->{"ROLES"}}; $j++) {
2707                                    $RoleHash->{$FeatureTable->get_row($i)->{"ROLES"}->[$j]} = 1;
2708                                    my $Subsystems = $self->figmodel()->subsystems_of_role($FeatureTable->get_row($i)->{"ROLES"}->[$j]);
2709                                    if (defined($Subsystems)) {
2710                                            for (my $k=0; $k < @{$Subsystems}; $k++) {
2711                                                    if ($Phoenix == 0) {
2712                                                            $VariantHash->{$Subsystems->[$k]} = 1;
2713                                                    }
2714                                            }
2715                                    }
2716                            }
2717                    }
2718    
2719                    #Scanning through the template item by item and determinine which biomass components should be added
2720                    my $ComponentTypes;
2721                    my $EquationData;
2722                    my $BiomassReactionTemplateTable = $self->figmodel()->database()->GetDBTable("BIOMASS TEMPLATE");
2723                    for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {
2724                            my $Row = $BiomassReactionTemplateTable->get_row($i);
2725                            if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {
2726                                    push(@{$EquationData},$Row);
2727                                    $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;
2728                            } else {
2729                                    my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];
2730                                    my $End = 0;
2731                                    while ($End == 0) {
2732                                            if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {
2733                                                    print $Criteria."\n";
2734                                                    my $Start = "";
2735                                                    my $End = "";
2736                                                    my $Condition = $1;
2737                                                    my $Data = $2;
2738                                                    if ($1 ne "AND" && $1 ne "OR") {
2739                                                            $Start = $1;
2740                                                            $End = $4;
2741                                                            $Condition = $2;
2742                                                            $Data = $3;
2743                                                    }
2744                                                    my $Result = "YES";
2745                                                    if ($Condition eq "OR") {
2746                                                            $Result = "NO";
2747                                                    }
2748                                                    my @Array = split(/\|/,$Data);
2749                                                    for (my $j=0; $j < @Array; $j++) {
2750                                                            if ($Array[$j] eq "YES" && $Condition eq "OR") {
2751                                                                    $Result = "YES";
2752                                                                    last;
2753                                                            } elsif ($Array[$j] eq "NO" && $Condition eq "AND") {
2754                                                                    $Result = "NO";
2755                                                                    last;
2756                                                            } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {
2757                                                                    my $Match = 0;
2758                                                                    for (my $k=0; $k < @{$EquationData}; $k++) {
2759                                                                            if ($EquationData->[$k]->{"ID"}->[0] eq $1) {
2760                                                                                    $Match = 1;
2761                                                                                    last;
2762                                                                            }
2763                                                                    }
2764                                                                    if ($Match == 1 && $Condition eq "OR") {
2765                                                                            $Result = "YES";
2766                                                                            last;
2767                                                                    } elsif ($Match != 1 && $Condition eq "AND") {
2768                                                                            $Result = "NO";
2769                                                                            last;
2770                                                                    }
2771                                                            } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {
2772                                                                    my $Match = 0;
2773                                                                    for (my $k=0; $k < @{$EquationData}; $k++) {
2774                                                                            if ($EquationData->[$k]->{"ID"}->[0] eq $1) {
2775                                                                                    $Match = 1;
2776                                                                                    last;
2777                                                                            }
2778                                                                    }
2779                                                                    if ($Match != 1 && $Condition eq "OR") {
2780                                                                            $Result = "YES";
2781                                                                            last;
2782                                                                    } elsif ($Match == 1 && $Condition eq "AND") {
2783                                                                            $Result = "NO";
2784                                                                            last;
2785                                                                    }
2786                                                            } elsif ($Array[$j] =~ m/^NAME:(.+)/) {
2787                                                                    my $Comparison = $1;
2788                                                                    if ((!defined($Comparison) || !defined($Name) || $Name =~ m/$Comparison/) && $Condition eq "OR") {
2789                                                                            $Result = "YES";
2790                                                                            last;
2791                                                                    } elsif (defined($Comparison) && defined($Name) && $Name !~ m/$Comparison/ && $Condition eq "AND") {
2792                                                                            $Result = "NO";
2793                                                                            last;
2794                                                                    }
2795                                                            } elsif ($Array[$j] =~ m/^!NAME:(.+)/) {
2796                                                                    my $Comparison = $1;
2797                                                                    if ((!defined($Comparison) || !defined($Name) || $Name !~ m/$Comparison/) && $Condition eq "OR") {
2798                                                                            $Result = "YES";
2799                                                                            last;
2800                                                                    } elsif (defined($Comparison) && defined($Name) && $Name =~ m/$Comparison/ && $Condition eq "AND") {
2801                                                                            $Result = "NO";
2802                                                                            last;
2803                                                                    }
2804                                                            } elsif ($Array[$j] =~ m/^SUBSYSTEM:(.+)/) {
2805                                                                    my @SubsystemArray = split(/`/,$1);
2806                                                                    if (@SubsystemArray == 1) {
2807                                                                            if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "OR") {
2808                                                                                    $Result = "YES";
2809                                                                                    last;
2810                                                                            } elsif ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "AND") {
2811                                                                                    $Result = "NO";
2812                                                                                    last;
2813                                                                            }
2814                                                                    } else {
2815                                                                            my $Match = 0;
2816                                                                            for (my $k=1; $k < @SubsystemArray; $k++) {
2817                                                                                    if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) {
2818                                                                                            $Match = 1;
2819                                                                                            last;
2820                                                                                    }
2821                                                                            }
2822                                                                            if ($Match == 1 && $Condition eq "OR") {
2823                                                                                    $Result = "YES";
2824                                                                                    last;
2825                                                                            } elsif ($Match != 1 && $Condition eq "AND") {
2826                                                                                    $Result = "NO";
2827                                                                                    last;
2828                                                                            }
2829                                                                    }
2830                                                            } elsif ($Array[$j] =~ m/^!SUBSYSTEM:(.+)/) {
2831                                                                    my @SubsystemArray = split(/`/,$1);
2832                                                                    if (@SubsystemArray == 1) {
2833                                                                            if ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "OR") {
2834                                                                                    $Result = "YES";
2835                                                                                    last;
2836                                                                            } elsif (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "AND") {
2837                                                                                    $Result = "NO";
2838                                                                                    last;
2839                                                                            }
2840                                                                    } else {
2841                                                                            my $Match = 0;
2842                                                                            for (my $k=1; $k < @SubsystemArray; $k++) {
2843                                                                                    if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) {
2844                                                                                            $Match = 1;
2845                                                                                            last;
2846                                                                                    }
2847                                                                            }
2848                                                                            if ($Match != 1 && $Condition eq "OR") {
2849                                                                                    $Result = "YES";
2850                                                                                    last;
2851                                                                            } elsif ($Match == 1 && $Condition eq "AND") {
2852                                                                                    $Result = "NO";
2853                                                                                    last;
2854                                                                            }
2855                                                                    }
2856                                                            } elsif ($Array[$j] =~ m/^ROLE:(.+)/) {
2857                                                                    if (defined($RoleHash->{$1}) && $Condition eq "OR") {
2858                                                                            $Result = "YES";
2859                                                                            last;
2860                                                                    } elsif (!defined($RoleHash->{$1}) && $Condition eq "AND") {
2861                                                                            $Result = "NO";
2862                                                                            last;
2863                                                                    }
2864                                                            } elsif ($Array[$j] =~ m/^!ROLE:(.+)/) {
2865                                                                    if (!defined($RoleHash->{$1}) && $Condition eq "OR") {
2866                                                                            $Result = "YES";
2867                                                                            last;
2868                                                                    } elsif (defined($RoleHash->{$1}) && $Condition eq "AND") {
2869                                                                            $Result = "NO";
2870                                                                            last;
2871                                                                    }
2872                                                            } elsif ($Array[$j] =~ m/^CLASS:(.+)/) {
2873                                                                    if ($Class eq $1 && $Condition eq "OR") {
2874                                                                            $Result = "YES";
2875                                                                            last;
2876                                                                    } elsif ($Class ne $1 && $Condition eq "AND") {
2877                                                                            $Result = "NO";
2878                                                                            last;
2879                                                                    }
2880                                                            } elsif ($Array[$j] =~ m/^!CLASS:(.+)/) {
2881                                                                    if ($Class ne $1 && $Condition eq "OR") {
2882                                                                            $Result = "YES";
2883                                                                            last;
2884                                                                    } elsif ($Class eq $1 && $Condition eq "AND") {
2885                                                                            $Result = "NO";
2886                                                                            last;
2887                                                                    }
2888                                                            }
2889                                                    }
2890                                                    $Criteria = $Start.$Result.$End;
2891                                                    print $Criteria."\n";
2892                                            } else {
2893                                                    $End = 1;
2894                                                    last;
2895                                            }
2896                                    }
2897                                    if ($Criteria eq "YES") {
2898                                            push(@{$EquationData},$Row);
2899                                            $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;
2900                                    }
2901                            }
2902                    }
2903    
2904                    #Building biomass equation
2905                    my %Reactants;
2906                    my %Products;
2907                    foreach my $EquationRow (@{$EquationData}) {
2908                            #First determine what the coefficient should be
2909                            my $Coefficient;
2910                            if (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^[0-9\.]+$/) {
2911                                    $Coefficient = $EquationRow->{"COEFFICIENT"}->[0];
2912                            } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/cpd\d\d\d\d\d/) {
2913                                    $Coefficient = 0;
2914                                    my @CompoundList = split(/,/,$EquationRow->{"COEFFICIENT"}->[0]);
2915                                    foreach my $Compound (@CompoundList) {
2916                                            if (defined($Reactants{$Compound})) {
2917                                                    $Coefficient += $Reactants{$Compound};
2918                                            }
2919                                    }
2920                            } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^([0-9\.]+)\/(.+)$/) {
2921                                    my @Keys = keys(%{$ComponentTypes->{$2}});
2922                                    my $MW = 1;
2923                                    my $CompoundData = $self->figmodel()->LoadObject($EquationRow->{"ID"}->[0]);
2924                                    if (defined($CompoundData->{"MASS"})) {
2925                                            $MW = $CompoundData->{"MASS"}->[0];
2926                                    }
2927                                    $Coefficient = $1/@Keys/$MW;
2928                            }
2929                            if (defined($EquationRow->{"REACTANT"}) && $EquationRow->{"REACTANT"}->[0] eq "YES") {
2930                                    if (defined($Reactants{$EquationRow->{"ID"}->[0]})) {
2931                                            $Reactants{$EquationRow->{"ID"}->[0]} += $Coefficient;
2932                                    } elsif (defined($Products{$EquationRow->{"ID"}->[0]}) && $Products{$EquationRow->{"ID"}->[0]} > $Coefficient) {
2933                                            $Products{$EquationRow->{"ID"}->[0]} -= $Coefficient;
2934                                    } elsif (defined($Products{$EquationRow->{"ID"}->[0]}) && $Products{$EquationRow->{"ID"}->[0]} < $Coefficient) {
2935                                            $Reactants{$EquationRow->{"ID"}->[0]} = $Coefficient - $Products{$EquationRow->{"ID"}->[0]};
2936                                            delete $Products{$EquationRow->{"ID"}->[0]};
2937                                    } else {
2938                                            $Reactants{$EquationRow->{"ID"}->[0]} = $Coefficient;
2939                                    }
2940                            } else {
2941                                    if (defined($Products{$EquationRow->{"ID"}->[0]})) {
2942                                            $Products{$EquationRow->{"ID"}->[0]} += $Coefficient;
2943                                    } elsif (defined($Reactants{$EquationRow->{"ID"}->[0]}) && $Reactants{$EquationRow->{"ID"}->[0]} > $Coefficient) {
2944                                            $Reactants{$EquationRow->{"ID"}->[0]} -= $Coefficient;
2945                                    } elsif (defined($Reactants{$EquationRow->{"ID"}->[0]}) && $Reactants{$EquationRow->{"ID"}->[0]} < $Coefficient) {
2946                                            $Products{$EquationRow->{"ID"}->[0]} = $Coefficient - $Reactants{$EquationRow->{"ID"}->[0]};
2947                                            delete $Reactants{$EquationRow->{"ID"}->[0]};
2948                                    } else {
2949                                            $Products{$EquationRow->{"ID"}->[0]} = $Coefficient;
2950                                    }
2951                            }
2952                    }
2953                    my $Equation = "";
2954                    my @ReactantList = sort(keys(%Reactants));
2955                    for (my $i=0; $i < @ReactantList; $i++) {
2956                            if (length($Equation) > 0) {
2957                                    $Equation .= " + ";