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