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