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