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