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