[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.2, Fri Dec 11 21:37:02 2009 UTC revision 1.27, Fri Aug 20 14:35:38 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,$figmodel,$id) = @_;          my ($class,$figmodel,$id,$metagenome) = @_;
17          #Error checking first          #Error checking first
18          if (!defined($figmodel)) {          if (!defined($figmodel)) {
19                  print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n";                  print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n";
# Line 25  Line 25 
25                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,undef):id must be defined to create a model object");                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,undef):id must be defined to create a model object");
26                  return undef;                  return undef;
27          }          }
   
28          #Checking that the id exists          #Checking that the id exists
29          my $tbl = $self->figmodel()->database()->GetDBTable("MODELS");          if (!defined($metagenome) || $metagenome != 1) {
30          if (!defined($tbl)) {                  my $modelHandle = $self->figmodel()->database()->get_object_manager("model");
31                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!");                  if (defined($modelHandle)) {
32                  return undef;                          #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          #If the id is a number, we get the model row by index
         my $index = $id;  
60          if ($id =~ m/^\d+$/) {          if ($id =~ m/^\d+$/) {
61                  $self->{_data} = $tbl->get_row($id);                                  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 {          } else {
67                  $self->{_data} = $tbl->get_row_by_key($id,"id");                                  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})) {                  if (defined($self->{_data})) {
80                          $index = $tbl->row_index($self->{_data});                          $self->{_modeltype} = "metagenome";
81                  }                  }
82          }          }
83          if (!defined($self->{_data})) {          if (!defined($self->{_data})) {
84                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find specified id in database!");                  $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!");
85                  return undef;                  return undef;
86          }          }
         $self->{_index} = $index;  
87          $self->figmodel()->{_models}->{$self->id()} = $self;          $self->figmodel()->{_models}->{$self->id()} = $self;
         $self->figmodel()->{_models}->{$self->index()} = $self;  
   
88          return $self;          return $self;
89  }  }
90    
# Line 66  Line 100 
100          return $self->figmodel()->config($key);          return $self->figmodel()->config($key);
101  }  }
102    
103    =head3 error_message
104    Definition:
105            string:message text = FIGMODELmodel->error_message(string::message);
106    Description:
107    =cut
108    sub error_message {
109            my ($self,$message) = @_;
110            return $self->figmodel()->error_message("FIGMODELmodel:".$self->id().":".$message);
111    }
112    
113  =head3 fail  =head3 fail
114  Definition:  Definition:
115          -1 = FIGMODELmodel->fail();          -1 = FIGMODELmodel->fail();
# Line 107  Line 151 
151  =cut  =cut
152  sub fig {  sub fig {
153          my ($self) = @_;          my ($self) = @_;
   
154          if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {          if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) {
155                  if ($self->source() =~ /^RAST/) {                  if ($self->source() =~ /^RAST/) {
156                          $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());                          $self->{"_fig"} = $self->figmodel()->fig();#$self->genome());
# Line 115  Line 158 
158                          $self->{"_fig"} = $self->figmodel()->fig();                          $self->{"_fig"} = $self->figmodel()->fig();
159                  }                  }
160          }          }
   
161          return $self->{"_fig"};          return $self->{"_fig"};
162  }  }
163    
164    =head3 aquireModelLock
165    
166    Definition:
167    
168            FIGMODELmodel->aquireModelLock();
169    
170    Description:
171    
172            Locks the database for alterations relating to the current model object
173    
174    =cut
175    sub aquireModelLock {
176            my ($self) = @_;
177            $self->figmodel()->database()->genericLock($self->id());
178    }
179    
180    =head3 releaseModelLock
181    
182    Definition:
183    
184            FIGMODELmodel->releaseModelLock();
185    
186    Description:
187    
188            Unlocks the database for alterations relating to the current model object
189    
190    =cut
191    sub releaseModelLock {
192            my ($self) = @_;
193            $self->figmodel()->database()->genericUnlock($self->id());
194    }
195    
196  =head3 mgdata  =head3 mgdata
197  Definition:  Definition:
198          FIGMODEL = FIGMODELmodel->mgdata();          FIGMODEL = FIGMODELmodel->mgdata();
# Line 127  Line 201 
201  =cut  =cut
202  sub mgdata {  sub mgdata {
203          my ($self) = @_;          my ($self) = @_;
204            if (!defined($self->{_mgdata})) {
205          if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) {                  my $mgrastdbhandle = $self->figmodel()->database()->get_object_manager("mgjob");
206                  require MGRAST;                  my $objects = $mgrastdbhandle->get_objects( { 'genome_id' => $self->genome() } );
207                  $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } )                  $self->{_mgdata} = $objects->[0];
208          }          }
   
209          return $self->{_mgdata};          return $self->{_mgdata};
210  }  }
211    
# Line 144  Line 217 
217  =cut  =cut
218  sub id {  sub id {
219          my ($self) = @_;          my ($self) = @_;
220          return $self->{_data}->{id}->[0];          return $self->{_data}->id();
221    }
222    
223    =head3 get_model_type
224    Definition:
225            string = FIGMODELmodel->get_model_type();
226    Description:
227            Returns the type of the model
228    =cut
229    sub get_model_type {
230            my ($self) = @_;
231            return $self->{_modeltype};
232  }  }
233    
234  =head3 owner  =head3 owner
# Line 155  Line 239 
239  =cut  =cut
240  sub owner {  sub owner {
241          my ($self) = @_;          my ($self) = @_;
242          return $self->{_data}->{owner}->[0];          return $self->{_data}->owner();
243  }  }
244    
245  =head3 index  =head3 users
246  Definition:  Definition:
247          string = FIGMODELmodel->index();          string = FIGMODELmodel->users();
248  Description:  Description:
         Returns model index  
249  =cut  =cut
250  sub index {  sub users {
251          my ($self) = @_;          my ($self) = @_;
252          return $self->{_index};          return $self->{_data}->users();
253  }  }
254    
255  =head3 name  =head3 name
# Line 185  Line 268 
268                          if (defined($self->mgdata())) {                          if (defined($self->mgdata())) {
269                                  $self->{_name} = $self->mgdata()->genome_name;                                  $self->{_name} = $self->mgdata()->genome_name;
270                          }                          }
                 } elsif (defined($self->stats())) {  
                         $self->{_name} = $self->stats()->{'Organism name'}->[0];  
271                  } elsif ($source !~ /^RAST/) {                  } elsif ($source !~ /^RAST/) {
272                          $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());                          $self->{_name} = $self->fig()->orgname_of_orgid($self->genome());
273                    } else {
274                            $self->{_name} = $self->figmodel()->get_genome_stats($self->genome())->{NAME}->[0];
275                  }                  }
276          }          }
277    
278          return $self->{_name};          return $self->{_name};
279  }  }
280    
 =head3 stats  
 Definition:  
         string = FIGMODELmodel->stats();  
 Description:  
         Returns model stats  
 =cut  
 sub stats {  
         my ($self) = @_;  
   
         if (!defined($self->{_stats})) {  
                 $self->{_stats} = $self->figmodel()->database()->GetDBTable("MODEL STATS")->get_row_by_key($self->id(),"Model ID");  
         }  
         return $self->{_stats};  
 }  
   
281  =head3 get_reaction_class  =head3 get_reaction_class
282  Definition:  Definition:
283          string = FIGMODELmodel->get_reaction_class(string::reaction ID);          string = FIGMODELmodel->get_reaction_class(string::reaction ID);
# Line 217  Line 285 
285          Returns reaction class          Returns reaction class
286  =cut  =cut
287  sub get_reaction_class {  sub get_reaction_class {
288          my ($self,$reaction,$nohtml) = @_;          my ($self,$reaction,$nohtml,$brief_flux) = @_;
289    
290          if (!-e $self->directory()."ReactionClassification-".$self->id().".tbl") {          if (!-e $self->directory()."ReactionClassification-".$self->id().".tbl") {
291                  if (!defined($self->{_reaction_classes})) {                  if (!defined($self->{_reaction_classes})) {
# Line 230  Line 298 
298                  my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");                  my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION");
299                  if (defined($ClassRow) && defined($ClassRow->{CLASS})) {                  if (defined($ClassRow) && defined($ClassRow->{CLASS})) {
300                          my $class;                          my $class;
301                            my $min = $ClassRow->{MIN}->[0];
302                            my $max = $ClassRow->{MAX}->[0];
303                          if ($ClassRow->{CLASS}->[0] eq "Positive") {                          if ($ClassRow->{CLASS}->[0] eq "Positive") {
304                                  $class = "Essential =>";                                  $class = "Essential =>";
305                                    $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
306                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative") {
307                                  $class = "Essential <=";                                  $class = "Essential <=";
308                                    $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $class.="<br>[Flux: ".$max." to ".$min."]<br>";
309                          } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") {
310                                  $class = "Active =>";                                  $class = "Active =>";
311                                    $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
312                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") {
313                                  $class = "Active <=";                                  $class = "Active <=";
314                                    $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $class.="<br>[Flux: ".$max." to ".$min."]<br>";
315                          } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Variable") {
316                                  $class = "Active <=>";                                  $class = "Active <=>";
317                                    $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
318                          } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") {
319                                  $class = "Inactive";                                  $class = "Inactive";
320                          } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {                          } elsif ($ClassRow->{CLASS}->[0] eq "Dead") {
321                                  $class = "Disconnected";                                  $class = "Disconnected";
322                          }                          }
323    
324                          if (!defined($nohtml) || $nohtml ne "1") {                          if (!defined($nohtml) || $nohtml ne "1") {
325                                  $class = "<span title=\"Flux:".$ClassRow->{MIN}->[0]." to ".$ClassRow->{MAX}->[0]."\">".$class."</span>";                                  $class = "<span title=\"Flux:".$min." to ".$max."\">".$class."</span>";
326                          }                          }
327    
328                          return $class;                          return $class;
329                  }                  }
330                  return undef;                  return undef;
# Line 268  Line 345 
345                                  $classstring .= "<br>";                                  $classstring .= "<br>";
346                          }                          }
347                          my $NewClass;                          my $NewClass;
348                            my $min = $ClassRow->{MIN}->[$i];
349                            my $max = $ClassRow->{MAX}->[$i];
350                          if ($ClassRow->{CLASS}->[$i] eq "Positive") {                          if ($ClassRow->{CLASS}->[$i] eq "Positive") {
351                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>";
352                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
353                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") {
354                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <=";
355                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $NewClass.="<br>[Flux: ".$max." to ".$min."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
356                          } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") {
357                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>";
358                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
359                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") {
360                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=";
361                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $NewClass.="<br>[Flux: ".$max." to ".$min."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
362                          } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") {
363                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>";
364                                  if (!defined($nohtml) || $nohtml ne "1") {                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
365                          } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") {
366                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive";
                                 if (!defined($nohtml) || $nohtml ne "1") {  
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
                                 }  
367                          } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {                          } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") {
368                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";                                  $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected";
                                 if (!defined($nohtml) || $nohtml ne "1") {  
                                         $NewClass = "<span title=\"Flux:".$ClassRow->{MIN}->[$i]." to ".$ClassRow->{MAX}->[$i]."\">".$NewClass."</span>";  
369                                  }                                  }
370    
371                            if (!defined($nohtml) || $nohtml ne "1") {
372                                    $NewClass = "<span title=\"Flux:".$min." to ".$max."\">".$NewClass."</span>";
373                          }                          }
374                          $classstring .= $NewClass;                          $classstring .= $NewClass;
375                  }                  }
# Line 321  Line 388 
388    
389          if (!defined($self->{_biomass})) {          if (!defined($self->{_biomass})) {
390                  my $rxntbl = $self->reaction_table();                  my $rxntbl = $self->reaction_table();
391                    if (defined($rxntbl)) {
392                  for (my $i=0; $i < $rxntbl->size(); $i++) {                  for (my $i=0; $i < $rxntbl->size(); $i++) {
393                          if ($rxntbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {                          if ($rxntbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {
394                                  $self->{_biomass} = $rxntbl->get_row($i)->{"LOAD"}->[0];                                  $self->{_biomass} = $rxntbl->get_row($i)->{"LOAD"}->[0];
# Line 328  Line 396 
396                          }                          }
397                  }                  }
398          }          }
399            }
400    
401          return $self->get_reaction_data($self->{_biomass});          return $self->get_reaction_data($self->{_biomass});
402  }  }
# Line 340  Line 409 
409  =cut  =cut
410  sub get_reaction_data {  sub get_reaction_data {
411          my ($self,$reaction) = @_;          my ($self,$reaction) = @_;
   
412          if (!defined($self->reaction_table())) {          if (!defined($self->reaction_table())) {
413                  return undef;                  return undef;
414          }          }
415          if ($reaction =~ m/^\d+$/) {          if ($reaction =~ m/^\d+$/) {
416                  $self->reaction_table()->get_row($reaction);                  return $self->reaction_table()->get_row($reaction);
417          }          }
418          return $self->reaction_table()->get_row_by_key($reaction,"LOAD");          return $self->reaction_table()->get_row_by_key($reaction,"LOAD");
419  }  }
# Line 366  Line 434 
434          return $Data->{LOAD}->[0];          return $Data->{LOAD}->[0];
435  }  }
436    
437    =head3 load_model_table
438    
439    Definition: FIGMODELTable = FIGMODELmodel->load_model_table(string:table name,0/1:refresh the table));
440    
441    Description: Returns the table specified by the input filename. Table will be stored in a file in the model directory.
442    
443    =cut
444    sub load_model_table {
445            my ($self,$name,$refresh) = @_;
446            if (defined($refresh) && $refresh == 1) {
447                    delete $self->{"_".$name};
448            }
449            if (!defined($self->{"_".$name})) {
450                    my $tbldef = $self->figmodel()->config($name);
451                    if (!defined($tbldef)) {
452                            return undef;
453                    }
454                    my $itemDelim = "|";
455                    if (defined($tbldef->{itemdelimiter}->[0])) {
456                            $itemDelim = $tbldef->{itemdelimiter}->[0];
457                            if ($itemDelim eq "SC") {
458                                    $itemDelim = ";";
459                            }
460                    }
461                    my $columnDelim = "\t";
462                    if (defined($tbldef->{columndelimiter}->[0])) {
463                            $columnDelim = $tbldef->{columndelimiter}->[0];
464                            if ($columnDelim eq "SC") {
465                                    $columnDelim = ";";
466                            }
467                    }
468                    my $suffix = ".tbl";
469                    if (defined($tbldef->{filename_suffix}->[0])) {
470                            $suffix = $tbldef->{filename_suffix}->[0];
471                    }
472                    my $filename = $self->directory().$name."-".$self->id().$self->selected_version().$suffix;
473                    if (defined($tbldef->{filename_prefix}->[0])) {
474                            if ($tbldef->{filename_prefix}->[0] eq "NONE") {
475                                    $filename = $self->directory().$self->id().$self->selected_version().$suffix;
476                            } else {
477                                    $filename = $self->directory().$tbldef->{filename_prefix}->[0]."-".$self->id().$self->selected_version().$suffix;
478                            }
479                    }
480                    if (-e $filename) {
481                            $self->{"_".$name} = $self->figmodel()->database()->load_table($filename,$columnDelim,$itemDelim,$tbldef->{headingline}->[0],$tbldef->{hashcolumns});
482                    } else {
483                            if (defined($tbldef->{prefix})) {
484                                    $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim,join(@{$tbldef->{prefix}},"\n"));
485                            } else {
486                                    $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim);
487                            }
488                    }
489            }
490            return $self->{"_".$name};
491    }
492    
493    =head3 create_table_prototype
494    
495    Definition:
496            FIGMODELTable::table = FIGMODELmodel->create_table_prototype(string::table);
497    Description:
498            Returns a empty FIGMODELTable with all the metadata associated with the input table name
499    
500    =cut
501    sub create_table_prototype {
502            my ($self,$TableName) = @_;
503            #Checking if the table definition exists in the FIGMODELconfig file
504            my $tbldef = $self->figmodel()->config($TableName);
505            if (!defined($tbldef)) {
506                    $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:Definition not found for ".$TableName);
507                    return undef;
508            }
509            #Checking that this is a database table
510            if (!defined($tbldef->{tabletype}) || $tbldef->{tabletype}->[0] ne "ModelTable") {
511                    $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:".$TableName." is not a model table!");
512                    return undef;
513            }
514            #Setting default values for table parameters
515            my $prefix;
516            if (defined($tbldef->{prefix})) {
517                    $prefix = join("\n",@{$self->config($TableName)->{prefix}})."\n";
518            }
519            my $itemDelim = "|";
520            if (defined($tbldef->{itemdelimiter}->[0])) {
521                    $itemDelim = $tbldef->{itemdelimiter}->[0];
522                    if ($itemDelim eq "SC") {
523                            $itemDelim = ";";
524                    }
525            }
526            my $columnDelim = "\t";
527            if (defined($tbldef->{columndelimiter}->[0])) {
528                    $columnDelim = $tbldef->{columndelimiter}->[0];
529                    if ($columnDelim eq "SC") {
530                            $columnDelim = ";";
531                    }
532            }
533            my $suffix = ".tbl";
534            if (defined($tbldef->{filename_suffix}->[0])) {
535                    $suffix = $tbldef->{filename_suffix}->[0];
536            }
537            my $filename = $self->directory().$TableName."-".$self->id().$self->selected_version().$suffix;
538            if (defined($tbldef->{filename_prefix}->[0])) {
539                    if ($tbldef->{filename_prefix}->[0] eq "NONE") {
540                            $filename = $self->directory().$self->id().$self->selected_version().$suffix;
541                    } else {
542                            $filename = $self->directory().$tbldef->{filename_prefix}->[0]."-".$self->id().$self->selected_version().$suffix;
543                    }
544            }
545            #Creating the table prototype
546            my $tbl = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim,$prefix);
547            return $tbl;
548    }
549    
550  =head3 get_reaction_number  =head3 get_reaction_number
551  Definition:  Definition:
552          int = FIGMODELmodel->get_reaction_number();          int = FIGMODELmodel->get_reaction_number();
# Line 374  Line 555 
555  =cut  =cut
556  sub get_reaction_number {  sub get_reaction_number {
557          my ($self) = @_;          my ($self) = @_;
   
558          if (!defined($self->reaction_table())) {          if (!defined($self->reaction_table())) {
559                  return 0;                  return 0;
560          }          }
   
561          return $self->reaction_table()->size();          return $self->reaction_table()->size();
562  }  }
563    
# Line 389  Line 568 
568          Returns FIGMODELTable with the reaction list for the model          Returns FIGMODELTable with the reaction list for the model
569  =cut  =cut
570  sub reaction_table {  sub reaction_table {
571          my ($self) = @_;          my ($self,$clear) = @_;
572            if (defined($clear) && $clear == 1) {
573                    delete $self->{_reaction_table};
574            }
575            if (!defined($self->{_reaction_table})) {
576                    $self->{_reaction_table} = $self->load_model_table("ModelReactions",$clear);
577                    my $classTbl = $self->reaction_class_table();
578                    if (defined($classTbl)) {
579                            for (my $i=0; $i < $classTbl->size(); $i++) {
580                                    my $row = $classTbl->get_row($i);
581                                    if (defined($row->{REACTION})) {
582                                            my $rxnRow = $self->{_reaction_table}->get_row_by_key($row->{"REACTION"}->[0],"LOAD");
583                                            if (defined($row->{MEDIA})) {
584                                                    for (my $j=0; $j < @{$row->{MEDIA}};$j++) {
585                                                            my $class = "Active <=>";
586                                                            if ($row->{CLASS}->[$j] eq "Positive") {
587                                                                    $class = "Essential =>";
588                                                            } elsif ($row->{CLASS}->[$j] eq "Negative") {
589                                                                    $class = "Essential <=";
590                                                            } elsif ($row->{CLASS}->[$j] eq "Blocked") {
591                                                                    $class = "Inactive";
592                                                            } elsif ($row->{CLASS}->[$j] eq "Positive variable") {
593                                                                    $class = "Active =>";
594                                                            } elsif ($row->{CLASS}->[$j] eq "Negative variable") {
595                                                                    $class = "Active <=";
596                                                            } elsif ($row->{CLASS}->[$j] eq "Variable") {
597                                                                    $class = "Active <=>";
598                                                            } elsif ($row->{CLASS}->[$j] eq "Dead") {
599                                                                    $class = "Dead";
600                                                            }
601                                                            push(@{$rxnRow->{PREDICTIONS}},$row->{MEDIA}->[$j].":".$class);
602                                                    }
603                                            }
604                                    }
605                            }
606                    }
607            }
608            return $self->{_reaction_table};
609    }
610    
611          if (!defined($self->{_reaction_data})) {  =head3 essentials_table
612                  $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id());  Definition:
613            FIGMODELTable = FIGMODELmodel->essentials_table();
614    Description:
615            Returns FIGMODELTable with the essential genes for the model
616    =cut
617    sub essentials_table {
618            my ($self,$clear) = @_;
619            my $tbl = $self->load_model_table("ModelEssentialGenes",$clear);
620            return $tbl;
621          }          }
622    
623          return $self->{_reaction_data};  =head3 model_history
624    Definition:
625            FIGMODELTable = FIGMODELmodel->model_history();
626    Description:
627            Returns FIGMODELTable with the history of model changes
628    =cut
629    sub model_history {
630            my ($self,$clear) = @_;
631            return $self->load_model_table("ModelHistory",$clear);
632  }  }
633    
634  =head3 reaction_class_table  =head3 feature_table
635  Definition:  Definition:
636          FIGMODELTable = FIGMODELmodel->reaction_class_table();          FIGMODELTable = FIGMODELmodel->feature_table();
637  Description:  Description:
638          Returns FIGMODELTable with the reaction class data, and creates the table file  if it does not exist          Returns FIGMODELTable with the feature list for the model
639  =cut  =cut
640  sub reaction_class_table {  sub feature_table {
641          my ($self) = @_;          my ($self) = @_;
642    
643          if (!defined($self->{_reaction_class_table})) {          if (!defined($self->{_feature_data})) {
644                  if (-e $self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl") {                  #Getting the genome feature list
645                          $self->{_reaction_class_table} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["REACTION","CLASS","MEDIA"]);                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
646                    if (!defined($FeatureTable)) {
647                            print STDERR "FIGMODELmodel:feature_table:Could not get features for genome ".$self->genome()." in database!";
648                            return undef;
649                    }
650                    #Getting the reaction table for the model
651                    my $rxnTable = $self->reaction_table();
652                    if (!defined($rxnTable)) {
653                            print STDERR "FIGMODELmodel:feature_table:Could not get reaction table for model ".$self->id()." in database!";
654                            return undef;
655                    }
656                    #Cloning the feature table
657                    $self->{_feature_data} = $FeatureTable->clone_table_def();
658                    $self->{_feature_data}->add_headings(($self->id()."REACTIONS",$self->id()."PREDICTIONS"));
659                    for (my $i=0; $i < $rxnTable->size(); $i++) {
660                            my $Row = $rxnTable->get_row($i);
661                            if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
662                                    foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
663                                            my $temp = $GeneSet;
664                                            $temp =~ s/\+/|/g;
665                                            $temp =~ s/\sAND\s/|/gi;
666                                            $temp =~ s/\sOR\s/|/gi;
667                                            $temp =~ s/[\(\)\s]//g;
668                                            my @GeneList = split(/\|/,$temp);
669                                            foreach my $Gene (@GeneList) {
670                                                    my $FeatureRow = $self->{_feature_data}->get_row_by_key("fig|".$self->genome().".".$Gene,"ID");
671                                                    if (!defined($FeatureRow)) {
672                                                            $FeatureRow = $FeatureTable->get_row_by_key("fig|".$self->genome().".".$Gene,"ID");
673                                                            if (defined($FeatureRow)) {
674                                                                    $self->{_feature_data}->add_row($FeatureRow);
675                                                            }
676                                                    }
677                                                    if (defined($FeatureRow)) {
678                                                            $self->{_feature_data}->add_data($FeatureRow,$self->id()."REACTIONS",$Row->{"LOAD"}->[0],1);
679                                                    }
680                                            }
681                                    }
682                            }
683                    }
684                    #Loading predictions
685                    my $esstbl = $self->essentials_table();
686                    for (my $i=0; $i < $self->{_feature_data}->size(); $i++) {
687                            my $Row = $self->{_feature_data}->get_row($i);
688                            if ($Row->{ID}->[0] =~ m/(peg\.\d+)/) {
689                                    my $gene = $1;
690                                    my @rows = $esstbl->get_rows_by_key($gene,"ESSENTIAL GENES");
691                                    my $mediahash;
692                                    for (my $j=0; $j < $esstbl->size(); $j++) {
693                                            $mediahash->{$esstbl->get_row($j)->{MEDIA}->[0]} = 0;
694                                    }
695                                    for (my $j=0; $j < @rows; $j++) {
696                                            $mediahash->{$rows[$j]->{MEDIA}->[0]} = 1;
697                                    }
698                                    my @mediaList = keys(%{$mediahash});
699                                    for (my $j=0; $j < @mediaList; $j++) {
700                                            if ($mediahash->{$mediaList[$j]} == 1) {
701                                                    push(@{$Row->{$self->id()."PREDICTIONS"}},$mediaList[$j].":essential");
702                  } else {                  } else {
703                          $self->{_reaction_class_table} = FIGMODELTable->new(["REACTION","MEDIA","CLASS","MIN","MAX"],$self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",["REACTION","CLASS","MEDIA"],";","|",undef);                                                  push(@{$Row->{$self->id()."PREDICTIONS"}},$mediaList[$j].":nonessential");
704                                            }
705                  }                  }
706          }          }
707                    }
708            }
709            return $self->{_feature_data};
710    }
711    
712          return $self->{_reaction_class_table};  =head3 reaction_class_table
713    Definition:
714            FIGMODELTable = FIGMODELmodel->reaction_class_table();
715    Description:
716            Returns FIGMODELTable with the reaction class data, and creates the table file  if it does not exist
717    =cut
718    sub reaction_class_table {
719            my ($self,$clear) = @_;
720            return $self->load_model_table("ModelReactionClasses",$clear);
721  }  }
722    
723  =head3 compound_class_table  =head3 compound_class_table
# Line 425  Line 727 
727          Returns FIGMODELTable with the compound class data, and creates the table file  if it does not exist          Returns FIGMODELTable with the compound class data, and creates the table file  if it does not exist
728  =cut  =cut
729  sub compound_class_table {  sub compound_class_table {
730          my ($self) = @_;          my ($self,$clear) = @_;
731            return $self->load_model_table("ModelCompoundClasses",$clear);
         if (!defined($self->{_compound_class_table})) {  
                 if (-e $self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl") {  
                         $self->{_compound_class_table} = $self->figmodel()->database()->load_table($self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["COMPOUND","CLASS","MEDIA"]);  
                 } else {  
                         $self->{_compound_class_table} = FIGMODELTable->new(["COMPOUND","MEDIA","CLASS","MIN","MAX"],$self->directory()."CompoundClassification-".$self->id().$self->selected_version().".tbl",["COMPOUND","CLASS","MEDIA"],";","|",undef);  
                 }  
732          }          }
733    
734          return $self->{_compound_class_table};  =head3 get_essential_genes
735    Definition:
736            [string::peg ID] = FIGMODELmodel->get_essential_genes(string::media condition);
737    Description:
738            Returns an reference to an array of the predicted essential genes during growth in the input media condition
739    =cut
740    sub get_essential_genes {
741            my ($self,$media) = @_;
742            my $tbl = $self->essentials_table();
743            my $row = $tbl->get_row_by_key($media,"MEDIA");
744            if (defined($row)) {
745                    return $row->{"ESSENTIAL GENES"};
746            }
747            return undef;
748  }  }
749    
750  =head3 compound_table  =head3 compound_table
# Line 448  Line 757 
757          my ($self) = @_;          my ($self) = @_;
758    
759          if (!defined($self->{_compound_table})) {          if (!defined($self->{_compound_table})) {
760                  $self->{_compound_table} = $self->figmodel()->database()->GetDBModelCompounds($self->id());                  $self->{_compound_table} = $self->create_table_prototype("ModelCompounds");
761                    #Loading the reactions
762                    my $ReactionTable = $self->figmodel()->database()->get_table("REACTIONS");
763                    my $BiomassTable = $self->figmodel()->database()->get_table("BIOMASS");
764                    #Loading the model
765                    my $ModelTable = $self->reaction_table();
766                    #Checking that the tables were loaded
767                    if (!defined($ModelTable) || !defined($ReactionTable)) {
768                            return undef;
769                    }
770                    #Finding the biomass reaction
771                    for (my $i=0; $i < $ModelTable->size(); $i++) {
772                            my $ID = $ModelTable->get_row($i)->{"LOAD"}->[0];
773                            my $Row = $ReactionTable->get_row_by_key($ID,"DATABASE");
774                            my $IsBiomass = 0;
775                            if (!defined($Row)) {
776                                    $Row = $BiomassTable->get_row_by_key($ID,"DATABASE");
777                                    $IsBiomass = 1;
778                            }
779                            if (defined($Row->{"EQUATION"}->[0])) {
780                                    $_ = $Row->{"EQUATION"}->[0];
781                                    my @OriginalArray = /(cpd\d\d\d\d\d[\[\w]*)/g;
782                                    foreach my $Compound (@OriginalArray) {
783                                            my $ID = substr($Compound,0,8);
784                                            my $NewRow = $self->{_compound_table}->get_row_by_key($ID,"DATABASE",1);
785                                            if ($IsBiomass == 1) {
786                                                    $self->{_compound_table}->add_data($NewRow,"BIOMASS",$Row->{"DATABASE"}->[0],1);
787                                            }
788                                            if (length($Compound) > 8) {
789                                                    #print $Compound."\t".$Row->{"EQUATION"}->[0]."\t".$Row->{"DATABASE"}->[0]."\n";
790                                                    my $Compartment = substr($Compound,8,1);
791                                                    $self->{_compound_table}->add_data($NewRow,"COMPARTMENTS",$Compartment,1);
792                                                    $self->{_compound_table}->add_data($NewRow,"TRANSPORTERS",$Row->{"DATABASE"}->[0],1);
793                                            }
794                                    }
795                            }
796                    }
797          }          }
798    
799          return $self->{_compound_table};          return $self->{_compound_table};
# Line 456  Line 801 
801    
802  =head3 get_compound_data  =head3 get_compound_data
803  Definition:  Definition:
804          string = FIGMODELmodel->get_compound_data(string::compound ID);          {string:key=>[string]:values} = FIGMODELmodel->get_compound_data(string::compound ID);
805  Description:  Description:
806          Returns model compound data          Returns model compound data
807  =cut  =cut
# Line 467  Line 812 
812                  return undef;                  return undef;
813          }          }
814          if ($compound =~ m/^\d+$/) {          if ($compound =~ m/^\d+$/) {
815                  $self->compound_table()->get_row($compound);                  return $self->compound_table()->get_row($compound);
816          }          }
817          return $self->compound_table()->get_row_by_key($compound,"DATABASE");          return $self->compound_table()->get_row_by_key($compound,"DATABASE");
818  }  }
819    
820    =head3 get_feature_data
821    Definition:
822            {string:key=>[string]:values} = FIGMODELmodel->get_feature_data(string::feature ID);
823    Description:
824            Returns model feature data
825    =cut
826    sub get_feature_data {
827            my ($self,$feature) = @_;
828            if (!defined($self->feature_table())) {
829                    return undef;
830            }
831            if ($feature =~ m/^\d+$/) {
832                    return $self->feature_table()->get_row($feature);
833            }
834            if ($feature =~ m/(peg\.\d+)/) {
835                    $feature = $1;
836            }
837            return $self->feature_table()->get_row_by_key("fig|".$self->genome().".".$feature,"ID");
838    }
839    
840  =head3 genome  =head3 genome
841  Definition:  Definition:
842          string = FIGMODELmodel->genome();          string = FIGMODELmodel->genome();
# Line 480  Line 845 
845  =cut  =cut
846  sub genome {  sub genome {
847          my ($self) = @_;          my ($self) = @_;
848          return $self->{_data}->{genome}->[0];          return $self->{_data}->genome();
849  }  }
850    
851  =head3 source  =head3 source
# Line 491  Line 856 
856  =cut  =cut
857  sub source {  sub source {
858          my ($self) = @_;          my ($self) = @_;
859          return $self->{_data}->{source}->[0];          return $self->{_data}->source();
860  }  }
861    
862  =head3 rights  =head3 rights
# Line 502  Line 867 
867  =cut  =cut
868  sub rights {  sub rights {
869          my ($self,$username) = @_;          my ($self,$username) = @_;
   
870          if ($self->public()) {          if ($self->public()) {
871                  return 1;                  return 1;
872          }          }
# Line 510  Line 874 
874                  return 0;                  return 0;
875          }          }
876          if (!defined($self->{_userrights}->{$username})) {          if (!defined($self->{_userrights}->{$username})) {
877                  if (defined($self->{_data}->{master})) {                  $self->{_userrights}->{$self->{_data}->owner()} = 1;
878                          for (my $i=0; $i < @{$self->{_data}->{master}};$i++) {                  my @users = split(/\|/,$self->{_data}->users());
879                                  $self->{_userrights}->{$self->{_data}->{master}->[$i]} = 1;                  for (my $i=0; $i < @users;$i++) {
880                          }                          $self->{_userrights}->{$users[$i]} = 1;
                 }  
                 if (defined($self->{_data}->{users})) {  
                         for (my $i=0; $i < @{$self->{_data}->{users}};$i++) {  
                                 $self->{_userrights}->{$self->{_data}->{users}->[$i]} = 1;  
                         }  
881                  }                  }
882          }          }
883          return $self->{_userrights}->{$username};          return $self->{_userrights}->{$username};
# Line 532  Line 891 
891  =cut  =cut
892  sub public {  sub public {
893          my ($self) = @_;          my ($self) = @_;
894            if ($self->{_data}->users() eq "all") {
895          if (!defined($self->{_public})) {                  return 1;
                 $self->{_public} = 0;  
                 if (defined($self->{_data}->{users}->[0]) && $self->{_data}->{users}->[0] eq "all") {  
                         $self->{_public} = 1;  
                 }  
896          }          }
897          return $self->{_public};          return 0;
898  }  }
899    
900  =head3 directory  =head3 directory
# Line 558  Line 913 
913                  }                  }
914                  my $source = $self->source();                  my $source = $self->source();
915                  if ($source =~ /^MGRAST/) {                  if ($source =~ /^MGRAST/) {
916                          $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";                          $self->{_directory} = $self->figmodel()->config("mgrast model directory")->[0].$userdirectory.$self->genome()."/";
917                  } elsif ($source =~ /^RAST/) {                  } elsif ($source =~ /^RAST/) {
918                          $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";                          $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/";
919                  } elsif ($source =~ /^SEED/) {                  } elsif ($source =~ /^SEED/) {
# Line 587  Line 942 
942          return $self->directory().$self->id().$self->selected_version().".txt";          return $self->directory().$self->id().$self->selected_version().".txt";
943  }  }
944    
 =head3 set_metagenome_stats  
 Definition:  
         string = FIGMODELmodel->set_metagenome_stats();  
 Description:  
         Sets the values of many model stats for a metagenome  
 =cut  
 sub set_metagenome_stats {  
         my ($self) = @_;  
   
         $self->{_total_compounds} = 0;  
         if (defined($self->compound_table())) {  
                 $self->{_total_compounds} = $self->compound_table()->size();  
         }  
         $self->{_gene_reactions} = 0;  
         $self->{_gapfilling_reactions} = 0;  
         $self->{_model_genes} = 0;  
         $self->{_total_reactions} = 0;  
         if (defined($self->reaction_table())) {  
                 $self->{_total_reactions} = $self->reaction_table()->size();  
                 my $tbl = $self->reaction_table();  
                 my $spontaneous = 0;  
                 for (my $i=0; $i < $tbl->size(); $i++) {  
                         my $row = $tbl->get_row($i);  
                         if (!defined($row->{"ASSOCIATED PEG"}->[0]) || $row->{"ASSOCIATED PEG"}->[0] !~ m/peg/) {  
                                 if ($row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/) {  
                                         $spontaneous++;  
                                 } else {  
                                         $self->{_gapfilling_reactions}++;  
                                 }  
                         } else {  
                                 for (my $j=0; $j < @{$row->{"CONFIDENCE"}}; $j++) {  
                                         my @ecores = split(/;/,$row->{"CONFIDENCE"}->[$j]);  
                                         $self->{_model_genes} += @ecores;  
                                 }  
                         }  
                 }  
                 $self->{_gene_reactions} = $tbl->size() - $spontaneous - $self->{_gapfilling_reactions};  
         }  
 }  
   
945  =head3 version  =head3 version
946  Definition:  Definition:
947          string = FIGMODELmodel->version();          string = FIGMODELmodel->version();
# Line 638  Line 953 
953    
954          if (!defined($self->{_version})) {          if (!defined($self->{_version})) {
955                  if (!defined($self->{_selectedversion})) {                  if (!defined($self->{_selectedversion})) {
956                          if (defined($self->stats())) {                          $self->{_version} = "V".$self->{_data}->version().".".$self->{_data}->autocompleteVersion();
                                 $self->{_version} = "V".$self->stats()->{"Version"}->[0].".".$self->stats()->{"Gap fill version"}->[0];  
                         }  
957                  } else {                  } else {
958                          $self->{_version} = $self->{_selectedversion};                          $self->{_version} = $self->{_selectedversion};
959                  }                  }
# Line 671  Line 984 
984  =cut  =cut
985  sub modification_time {  sub modification_time {
986          my ($self) = @_;          my ($self) = @_;
987          if (!defined($self->{_modification_time})) {          return $self->{_data}->modificationDate();
                 my $stats = $self->stats();  
                 if (defined($stats)) {  
                         $self->{_modification_time} = 0;  
                         if (defined($stats->{"Build date"}->[0]) && $self->{_modification_time} < $stats->{"Build date"}->[0]) {  
                                 $self->{_modification_time} = $stats->{"Build date"}->[0];  
                         } elsif (defined($stats->{"Gap fill date"}->[0]) && $self->{_modification_time} < $stats->{"Gap fill date"}->[0]) {  
                                 $self->{_modification_time} = $stats->{"Gap fill date"}->[0];  
                         }  
                 } else {  
                         $self->{_modification_time} = 0;  
                 }  
         }  
         return $self->{_modification_time};  
988  }  }
989    
990  =head3 gene_reactions  =head3 gene_reactions
# Line 695  Line 995 
995  =cut  =cut
996  sub gene_reactions {  sub gene_reactions {
997          my ($self) = @_;          my ($self) = @_;
998            return ($self->{_data}->reactions() - $self->{_data}->autoCompleteReactions() - $self->{_data}->spontaneousReactions() - $self->{_data}->gapFillReactions());
         if (!defined($self->{_gene_reactions})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_gene_reactions} = $self->total_reactions() - $self->gapfilling_reactions() - $self->stats()->{'Spontaneous'}->[0] - $self->stats()->{'Growmatch reactions'}->[0] - $self->stats()->{'Biolog gap filling reactions'}->[0];  
                 }  
         }  
         return $self->{_gene_reactions};  
999  }  }
1000    
1001  =head3 total_compounds  =head3 total_compounds
# Line 714  Line 1006 
1006  =cut  =cut
1007  sub total_compounds {  sub total_compounds {
1008          my ($self) = @_;          my ($self) = @_;
1009            return $self->{_data}->compounds();
         if (!defined($self->{_total_compounds})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_total_compounds} = $self->stats()->{'Metabolites'}->[0];  
                 }  
         }  
         return $self->{_total_compounds};  
1010  }  }
1011    
1012  =head3 gapfilling_reactions  =head3 gapfilling_reactions
# Line 733  Line 1017 
1017  =cut  =cut
1018  sub gapfilling_reactions {  sub gapfilling_reactions {
1019          my ($self) = @_;          my ($self) = @_;
1020            return ($self->{_data}->autoCompleteReactions()+$self->{_data}->gapFillReactions());
         if (!defined($self->{_gapfilling_reactions})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_gapfilling_reactions} = $self->stats()->{'Gap filling reactions'}->[0];  
                 }  
         }  
         return $self->{_gapfilling_reactions};  
1021  }  }
1022    
1023  =head3 total_reactions  =head3 total_reactions
# Line 752  Line 1028 
1028  =cut  =cut
1029  sub total_reactions {  sub total_reactions {
1030          my ($self) = @_;          my ($self) = @_;
1031            return $self->{_data}->reactions();
         if (!defined($self->{_total_reactions})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_total_reactions} = $self->stats()->{'Number of reactions'}->[0];  
                 }  
         }  
         return $self->{_total_reactions};  
1032  }  }
1033    
1034  =head3 model_genes  =head3 model_genes
# Line 771  Line 1039 
1039  =cut  =cut
1040  sub model_genes {  sub model_genes {
1041          my ($self) = @_;          my ($self) = @_;
1042            return $self->{_data}->associatedGenes();
         if (!defined($self->{_model_genes})) {  
                 if ($self->source() =~ /^MGRAST/) {  
                         $self->set_metagenome_stats();  
                 } elsif (defined($self->stats())) {  
                         $self->{_model_genes} = $self->stats()->{'Genes with reactions'}->[0];  
                 }  
         }  
         return $self->{_model_genes};  
1043  }  }
1044    
1045  =head3 class  =head3 class
# Line 790  Line 1050 
1050  =cut  =cut
1051  sub class {  sub class {
1052          my ($self) = @_;          my ($self) = @_;
1053            return $self->{_data}->cellwalltype();
1054    }
1055    
1056    sub autocompleteMedia {
1057            my ($self,$newMedia) = @_;
1058            if (defined($newMedia)) {
1059                    return $self->{_data}->autoCompleteMedia($newMedia);
1060            }
1061            return $self->{_data}->autoCompleteMedia();
1062    }
1063    
1064          if (!defined($self->{_class})) {  sub biomassReaction {
1065                  if ($self->source() =~ /^MGRAST/) {          my ($self,$newBiomass) = @_;
1066                          $self->{_class} = "Other";          if (!defined($newBiomass)) {
1067                  } elsif (defined($self->stats())) {                  return $self->{_data}->biomassReaction();
1068                          $self->{_class} = $self->stats()->{Class}->[0];          } else {
1069                    #Figuring out what the old biomass is
1070                    my $oldBiomass = $self->{_data}->biomassReaction();
1071                    $self->{_data}->biomassReaction($newBiomass);
1072                    #Changing the biomass reaction in the model file
1073                    my $rxnTbl = $self->reaction_table();
1074                    for (my $i=0; $i < $rxnTbl->size(); $i++) {
1075                            my $row = $rxnTbl->get_row($i);
1076                            if ($row->{LOAD}->[0] =~ m/^bio/) {
1077                                    $row->{LOAD}->[0] = $newBiomass;
1078                            }
1079                  }                  }
1080                    $rxnTbl->save();
1081                    if ($newBiomass ne $oldBiomass) {
1082                            #Figuring out if the new biomass exists
1083                            my $handle = $self->figmodel()->database()->get_object_manager("bof");
1084                            my $objects = $handle->get_objects({id=>$newBiomass});
1085                            if (!defined($objects) || !defined($objects->[0])) {
1086                                    print STDERR "Could not find new biomass reaction ".$newBiomass."\n";
1087                                    return $oldBiomass;
1088                            }
1089                    }
1090                    return $newBiomass;
1091          }          }
         return $self->{_class};  
1092  }  }
1093    
1094  =head3 taxonomy  =head3 taxonomy
# Line 866  Line 1156 
1156                          if (defined($self->mgdata())) {                          if (defined($self->mgdata())) {
1157                                  $self->{_genome_genes} = $self->mgdata()->genome_contig_count;                                  $self->{_genome_genes} = $self->mgdata()->genome_contig_count;
1158                          }                          }
1159                  } elsif (defined($self->stats())) {                  } else {
1160                          $self->{_genome_genes} = $self->stats()->{'Total genes'}->[0];                          $self->{_genome_genes} = $self->figmodel()->get_genome_stats($self->genome())->{"TOTAL GENES"}->[0];
1161                  }                  }
1162          }          }
1163    
# Line 884  Line 1174 
1174    
1175          #Assuming complete media if none is provided          #Assuming complete media if none is provided
1176          if (!defined($Media)) {          if (!defined($Media)) {
1177                  $Media = "Complete";                  $Media = $self->autocompleteMedia();
1178          }          }
1179    
1180          #Predicting essentiality          #Predicting essentiality
1181          my $result = $self->figmodel()->RunFBASimulation($self->id(),"SINGLEKO",undef,undef,[$self->id()],[$Media]);          my $result = $self->figmodel()->RunFBASimulation($self->id(),"SINGLEKO",undef,undef,[$self->id()],[$Media]);
1182          #Checking that the table is defined and the output file exists          #Checking that the table is defined and the output file exists
1183          if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {          if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) {
1184                  $self->figmodel()->database()->print_array_to_file($self->directory()."EssentialGenes-".$self->id()."-".$Media.".tbl",[join("\n",@{$result->get_row(0)->{"ESSENTIALGENES"}})]);                  my $tbl = $self->essentials_table();
1185                    my $row = $tbl->get_row_by_key($Media,"MEDIA",1);
1186                    $row->{"ESSENTIAL GENES"} = $result->get_row(0)->{"ESSENTIALGENES"};
1187                    $tbl->save();
1188            } else {
1189                    $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not identify essential reactions for model ".$self->id().$self->selected_version().".");
1190                    return $self->figmodel()->fail();
1191          }          }
1192    
1193          #Classifying reactions and compounds          #Classifying reactions and compounds
1194          my $tbl = $self->classify_model_reactions($Media);          my $tbl = $self->classify_model_reactions($Media);
1195            if (!defined($tbl)) {
1196                    $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not classify reactions for model ".$self->id().$self->selected_version().".");
1197                    return $self->figmodel()->fail();
1198            }
1199          $tbl->save();          $tbl->save();
1200    
1201          return $self->figmodel()->success();          return $self->figmodel()->success();
1202  }  }
1203    
 =head3 save_obsolete_stats  
 Definition:  
         FIGMODELmodel->save_obsolete_stats();  
 Description:  
 =cut  
 sub save_obsolete_stats {  
         my ($self) = @_;  
   
         #checking if stats exists  
         my $stats = $self->stats();  
         if (defined($stats)) {  
                 $stats->{"Model ID"}->[0] = $self->id()."V".$stats->{"Version"}->[0].".".$stats->{"Gap fill version"}->[0];  
                 $self->figmodel()->database()->update_row("OBSOLETE MODEL STATS",$stats,"Model ID");  
                 $stats->{"Model ID"}->[0] = $self->id();  
         }  
 }  
   
1204  =head3 update_stats_for_gap_filling  =head3 update_stats_for_gap_filling
1205  Definition:  Definition:
1206          {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);          {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time);
# Line 925  Line 1208 
1208  =cut  =cut
1209  sub update_stats_for_gap_filling {  sub update_stats_for_gap_filling {
1210          my ($self,$gapfilltime) = @_;          my ($self,$gapfilltime) = @_;
1211            $self->{_data}->autoCompleteTime($gapfilltime);
1212          #preserving the stats for the now obselete model          $self->{_data}->autocompleteDate(time());
1213          $self->save_obsolete_stats();          $self->{_data}->modificationDate(time());
1214          my $stats = $self->update_model_stats(0);          my $version = $self->{_data}->autocompleteVersion();
1215          $stats->{"Gap filling time"}->[0] = $gapfilltime;          $self->{_data}->autocompleteVersion($version+1);
1216          $stats->{"Gap fill date"}->[0] = time();          $self->update_model_stats();
         if (!defined($stats->{"Gap fill version"}->[0])) {  
                 $stats->{"Gap fill version"}->[0] = 0;  
         }  
         $stats->{"Gap fill version"}->[0]++;  
   
         #Updating the stats stored in the table  
         $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");  
         return $stats;  
1217  }  }
1218    
1219  =head3 update_stats_for_build  =head3 update_stats_for_build
# Line 948  Line 1223 
1223  =cut  =cut
1224  sub update_stats_for_build {  sub update_stats_for_build {
1225          my ($self) = @_;          my ($self) = @_;
1226            $self->{_data}->builtDate(time());
1227          #preserving the stats for the now obselete model          $self->{_data}->modificationDate(time());
1228          $self->save_obsolete_stats();          my $version = $self->{_data}->version();
1229          my $stats = $self->update_model_stats(0);          $self->{_data}->version($version+1);
         $stats->{"Build date"}->[0] = time();  
         if (!defined($stats->{"Version"}->[0])) {  
                 $stats->{"Version"}->[0] = 0;  
         }  
         $stats->{"Version"}->[0]++;  
   
         #Updating the stats stored in the table  
         $self->figmodel()->database()->update_row("MODEL STATS",$stats,"Model ID");  
         return $stats;  
1230  }  }
1231    
1232  =head3 update_model_stats  =head3 update_model_stats
# Line 974  Line 1240 
1240          #Getting reaction table          #Getting reaction table
1241          my $rxntbl = $self->reaction_table();          my $rxntbl = $self->reaction_table();
1242          if (!defined($rxntbl)) {          if (!defined($rxntbl)) {
1243                  $self->figmodel()->error_message("FIGMODELmodel:update_model_stats:Could not load reaction list for ".$self->id());                  die $self->error_message("update_model_stats:Could not load reaction list!");
                 return undef;  
1244          }          }
1245          my $cpdtbl = $self->compound_table();          my $cpdtbl = $self->compound_table();
1246    
1247          #Creating empty status row          #Calculating all necessary stats
         my $CurrentStats = {"Genes with reactions" => [0],  
                                                  "Metabolites" => [$cpdtbl->size()],  
                                                  "Growmatch reactions" => [0],  
                                                  "Spontaneous" => [0],  
                                                  "Biolog gap filling reactions" => [0],  
                                                  "Number of reactions" => [$rxntbl->size()],  
                                                  "Transport reaction"=>[0],  
                                                  "Gap filling reactions" => [0],  
                                                  "Model ID" => [$self->id()],  
                                                  "Subsystem genes with reactions" => [0],  
                                                  "Nonsubsystem genes with reactions" => [0],  
                                                  Version => [0],  
                                                  "Gap fill version" => [0],  
                                                  Source => [$self->source()],  
                                                  "Genes with one reaction" => [0],  
                                                  "Genome ID" => [$self->genome()]};  
   
         my $genomestats = $self->figmodel()->get_genome_stats($self->genome());  
         if (defined($genomestats)) {  
                 $CurrentStats->{SOURCE}->[0] = $genomestats->{SOURCE}->[0];  
                 $CurrentStats->{"Organism name"}->[0] = $genomestats->{NAME}->[0];  
                 $CurrentStats->{"Total genes"}->[0] = $genomestats->{"TOTAL GENES"}->[0];  
                 $CurrentStats->{"Gram positive genes"}->[0] = $genomestats->{"GRAM POSITIVE GENES"}->[0];  
                 $CurrentStats->{"Gram negative genes"}->[0] = $genomestats->{"GRAM NEGATIVE GENES"}->[0];  
                 $CurrentStats->{"Class"}->[0] = $genomestats->{CLASS}->[0];  
                 $CurrentStats->{"Genes with functions"}->[0] = $genomestats->{"GENES WITH FUNCTIONS"}->[0];  
                 $CurrentStats->{"Subsystem genes"}->[0] = $genomestats->{"SUBSYSTEM GENES"}->[0];  
                 $CurrentStats->{"Nonsubsystem genes"}->[0] = $genomestats->{"NON SUBSYSTEM GENES"}->[0];  
                 if (defined($genomestats->{TAXONOMY})) {  
                         $CurrentStats->{"Taxonomy 0"}->[0] = $genomestats->{TAXONOMY}->[0];  
                         $CurrentStats->{"Taxonomy 1"}->[0] = $genomestats->{TAXONOMY}->[1];  
                         $CurrentStats->{"Taxonomy 2"}->[0] = $genomestats->{TAXONOMY}->[2];  
                         $CurrentStats->{"Taxonomy 3"}->[0] = $genomestats->{TAXONOMY}->[3];  
                         $CurrentStats->{"Taxonomy 4"}->[0] = $genomestats->{TAXONOMY}->[4];  
                         $CurrentStats->{"Taxonomy 5"}->[0] = $genomestats->{TAXONOMY}->[5];  
                 }  
         }  
   
         #Transfering build, version, and gap fill data from existing stats  
         if (defined($self->stats())) {  
                 $CurrentStats->{Version}->[0] = $self->stats()->{Version}->[0];  
                 $CurrentStats->{"Gap fill version"}->[0] = $self->stats()->{"Gap fill version"}->[0];  
                 $CurrentStats->{"Gap filling time"}->[0] = $self->stats()->{"Gap filling time"}->[0];  
                 $CurrentStats->{"Gap fill date"}->[0] = $self->stats()->{"Gap fill date"}->[0];  
                 $CurrentStats->{"Build date"}->[0] = $self->stats()->{"Build date"}->[0];  
         }  
   
1248          my %GeneHash;          my %GeneHash;
1249          my %NonpegHash;          my %NonpegHash;
1250          my %CompoundHash;          my %CompoundHash;
1251            my $spontaneousReactions = 0;
1252            my $gapFillReactions = 0;
1253            my $biologReactions = 0;
1254            my $transporters = 0;
1255            my $autoCompleteReactions = 0;
1256            my $associatedSubsystemGenes = 0;
1257          for (my $i=0; $i < $rxntbl->size(); $i++) {          for (my $i=0; $i < $rxntbl->size(); $i++) {
1258                  my $Row = $rxntbl->get_row($i);                  my $Row = $rxntbl->get_row($i);
1259                  if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {                  if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
# Line 1037  Line 1261 
1261                          if (defined($ReactionRow->{"EQUATION"}->[0])) {                          if (defined($ReactionRow->{"EQUATION"}->[0])) {
1262                                  #Checking for extracellular metabolites which indicate that this is a transporter                                  #Checking for extracellular metabolites which indicate that this is a transporter
1263                                  if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {                                  if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {
1264                                          $CurrentStats->{"Transport reaction"}->[0]++;                                          $transporters++;
1265                                  }                                  }
1266                          }                          }
1267                          #Identifying spontaneous/biolog/gapfilling/gene associated reactions                          #Identifying spontaneous/biolog/gapfilling/gene associated reactions
1268                          if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {                          if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {
1269                                  $CurrentStats->{"Biolog gap filling reactions"}->[0]++;                                  $biologReactions++;
1270                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GROW/i) {
1271                                    $gapFillReactions++;
1272                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {
1273                                  $CurrentStats->{"Spontaneous"}->[0]++;                                  $spontaneousReactions++;
1274                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) {
1275                                  $CurrentStats->{"Gap filling reactions"}->[0]++;                                  $autoCompleteReactions++;
1276                          } else {                          } else {
1277                                  foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {                                  foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
1278                                          my @GeneList = split(/\+/,$GeneSet);                                          $_ = $GeneSet;
1279                                            my @GeneList = /(peg\.\d+)/g;
1280                                          foreach my $Gene (@GeneList) {                                          foreach my $Gene (@GeneList) {
1281                                                  if ($Gene =~ m/(peg\.\d+)/) {                                                  if ($Gene =~ m/(peg\.\d+)/) {
1282                                                          $GeneHash{$1} = 1;                                                          $GeneHash{$1} = 1;
# Line 1063  Line 1290 
1290          }          }
1291          my @genes = keys(%GeneHash);          my @genes = keys(%GeneHash);
1292          my @othergenes = keys(%NonpegHash);          my @othergenes = keys(%NonpegHash);
         $CurrentStats->{"Genes with reactions"}->[0] = @genes + @othergenes;  
1293    
1294          #Updating the stats stored in the table          #Setting the reaction count
1295          $self->figmodel()->database()->update_row("MODEL STATS",$CurrentStats,"Model ID");          $self->{_data}->reactions($rxntbl->size());
1296          $self->{_stats} = $CurrentStats;          #Setting the metabolite count
1297          return $CurrentStats;          $self->{_data}->compounds($cpdtbl->size());
1298            #Setting the gene count
1299            my $geneCount = @genes + @othergenes;
1300            $self->{_data}->associatedGenes($geneCount);
1301            #Setting remaining stats
1302            $self->{_data}->spontaneousReactions($spontaneousReactions);
1303            $self->{_data}->gapFillReactions($gapFillReactions);
1304            $self->{_data}->biologReactions($biologReactions);
1305            $self->{_data}->transporters($transporters);
1306            $self->{_data}->autoCompleteReactions($autoCompleteReactions);
1307            $self->{_data}->associatedSubsystemGenes($associatedSubsystemGenes);
1308            #Setting the model class
1309            my $class = "";
1310            for (my $i=0; $i < @{$self->figmodel()->config("class list")}; $i++) {
1311                    if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i]))) {
1312                            if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->id()})) {
1313                                    $class = $self->figmodel()->config("class list")->[$i];
1314                                    last;
1315                            }
1316                            if ($class eq "" && defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->genome()})) {
1317                                    $class = $self->figmodel()->config("class list")->[$i];
1318                            }
1319                    }
1320            }
1321            if ($class eq "") {
1322                    $class = $self->figmodel()->get_genome_stats($self->genome())->{CLASS}->[0];
1323            }
1324            if ($class eq "") {
1325                    $class = "unknown";
1326            }
1327            $self->{_data}->cellwalltype($class);
1328  }  }
1329    
1330  =head3 GapFillModel  =head3 GapFillModel
# Line 1088  Line 1344 
1344          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
1345          my $StartTime = time();          my $StartTime = time();
1346    
1347          #Archiving the existing model          #Reading original reaction table
1348          $self->ArchiveModel();          my $OriginalRxn = $self->reaction_table();
1349            #Clearing the table
1350            $self->reaction_table(1);
1351    
1352          #Removing any gapfilling reactions that may be currently present in the model          #Removing any gapfilling reactions that may be currently present in the model
1353          if (!defined($donotclear) || $donotclear != 1) {          if (!defined($donotclear) || $donotclear != 1) {
1354                  my $ModelTable = $self->reaction_table();                  my $ModelTable = $self->reaction_table();
1355                  for (my $i=0; $i < $ModelTable->size(); $i++) {                  for (my $i=0; $i < $ModelTable->size(); $i++) {
1356                          if (!defined($ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0]) || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "GAP FILLING" || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/ || $ModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/GROWMATCH/) {                          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/) {
1357                                  $ModelTable->delete_row($i);                                  $ModelTable->delete_row($i);
1358                                  $i--;                                  $i--;
1359                          }                          }
# Line 1104  Line 1362 
1362          }          }
1363    
1364          #Calling the MFAToolkit to run the gap filling optimization          #Calling the MFAToolkit to run the gap filling optimization
1365          my $MinimalMediaTable = $self->figmodel()->database()->GetDBTable("MINIMAL MEDIA TABLE");          my $Media = $self->autocompleteMedia();
1366          my $Media = "Complete";          if ($Media eq "Complete") {
1367          if (defined($MinimalMediaTable->get_row_by_key($self->genome(),"Organism"))) {                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1368                  $Media = $MinimalMediaTable->get_row_by_key($self->genome(),"Organism")->{"Minimal media"}->[0];          } else {
1369                  #Loading media, changing bounds, saving media as a test media                  #Loading media, changing bounds, saving media as a test media
1370                  my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);                  my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);
1371                  for (my $i=0; $i < $MediaTable->size(); $i++) {                  for (my $i=0; $i < $MediaTable->size(); $i++) {
# Line 1119  Line 1377 
1377                          }                          }
1378                  }                  }
1379                  $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");                  $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1380                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0},"GapFill".$self->id().".log",undef));                  #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";
1381                    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));
1382                  unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");                  unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
         } else {  
                 system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],undef,"GapFill".$self->id().".log",undef));  
1383          }          }
1384    
1385          #Parse the solutions file for the model and get the reaction list from it          #Parse the solutions file for the model and get the reaction list from it
1386          my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);          my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);
1387    
1388          #Clearing the mfatoolkit output and log file          #Clearing the mfatoolkit output and log file
1389          $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");          $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");
1390          if (!defined($SolutionData) || $SolutionData->size() == 0) {          if (!defined($SolutionData) || $SolutionData->size() == 0) {
# Line 1134  Line 1392 
1392                  $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");                  $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");
1393                  return $self->fail();                  return $self->fail();
1394          }          }
         $SolutionData->save("/home/chenry/Solution.txt");  
1395    
1396          #Looking for the last printed recursive MILP solution          #Looking for the last printed recursive MILP solution
1397          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
1398                  if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {                  if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) {
1399                          my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);                          my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);
                         print "Solution:".$AllSolutions."\n";  
1400                          my @TempThree = split(/\|/,$AllSolutions);                          my @TempThree = split(/\|/,$AllSolutions);
1401                          if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {                          if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {
1402                                  my @TempFour = split(/,/,$1);                                  my @TempFour = split(/,/,$1);
# Line 1163  Line 1419 
1419                                                  push(@{$ReactionList},$ID);                                                  push(@{$ReactionList},$ID);
1420                                          }                                          }
1421                                  }                                  }
1422                                  print "Solution:".$TempThree[0]."\n";                                  $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"AUTOCOMPLETION",0,1);
                                 $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1);  
1423                          }                          }
1424                          last;                          last;
1425                  }                  }
1426          }          }
1427    
1428          #Updating model stats with gap filling results          #Updating model stats with gap filling results
1429          my $ElapsedTime = time() - $StartTime;          my $ElapsedTime = time() - $StartTime;
1430          $self->figmodel()->ClearDBModel($self->id(),1);          $self->reaction_table(1);
1431            $self->calculate_model_changes($OriginalRxn,"AUTOCOMPLETION");
1432    
1433          #Determining why each gap filling reaction was added          #Determining why each gap filling reaction was added
1434          $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);          $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);
1435          if (!defined($donotclear) || $donotclear != 1) {          if (!defined($donotclear) || $donotclear != 1) {
                 $self->figmodel()->AddBiologTransporters($self->id());  
1436                  if ($self->id() !~ m/MGRast/) {                  if ($self->id() !~ m/MGRast/) {
1437                          $self->update_stats_for_gap_filling($ElapsedTime);                          $self->update_stats_for_gap_filling($ElapsedTime);
1438                  }                  }
# Line 1191  Line 1448 
1448          return $self->success();          return $self->success();
1449  }  }
1450    
1451    =head3 calculate_model_changes
1452    Definition:
1453            FIGMODELmodel->calculate_model_changes(FIGMODELTable:original reaction table,string:modification cause);
1454    Description:
1455    
1456    =cut
1457    
1458    sub calculate_model_changes {
1459            my ($self,$originalReactions,$cause,$tbl,$version) = @_;
1460            my $modTime = time();
1461            if (!defined($version)) {
1462                    $version = $self->selected_version();
1463            }
1464            my $user = $self->figmodel()->user();
1465            #Getting the history table
1466            my $histTbl = $self->model_history();
1467            #Checking for differences
1468            if (!defined($tbl)) {
1469                    $tbl = $self->reaction_table();
1470            }
1471            for (my $i=0; $i < $tbl->size(); $i++) {
1472                    my $row = $tbl->get_row($i);
1473                    my $orgRow = $originalReactions->get_row_by_key($row->{LOAD}->[0],"LOAD");
1474                    if (!defined($orgRow)) {
1475                            $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => $row->{DIRECTIONALITY}, GeneChange => $row->{"ASSOCIATED PEG"}, Action => ["ADDED"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1476                    } else {
1477                            my $geneChanges;
1478                            my $directionChange;
1479                            if ($orgRow->{"DIRECTIONALITY"}->[0] ne $row->{"DIRECTIONALITY"}->[0]) {
1480                                    $directionChange = $orgRow->{"DIRECTIONALITY"}->[0]." to ".$row->{"DIRECTIONALITY"}->[0];
1481                            }
1482                            for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
1483                                    my $match = 0;
1484                                    if (defined($orgRow->{"ASSOCIATED PEG"})) {
1485                                            for (my $k=0; $k < @{$orgRow->{"ASSOCIATED PEG"}}; $k++) {
1486                                                    if ($row->{"ASSOCIATED PEG"}->[$j] eq $orgRow->{"ASSOCIATED PEG"}->[$k]) {
1487                                                            $match = 1;
1488                                                    }
1489                                            }
1490                                    }
1491                                    if ($match == 0) {
1492                                            push(@{$geneChanges},"Added ".$row->{"ASSOCIATED PEG"}->[$j]);
1493                                    }
1494                            }
1495                            if (defined($orgRow->{"ASSOCIATED PEG"})) {
1496                                    for (my $k=0; $k < @{$orgRow->{"ASSOCIATED PEG"}}; $k++) {
1497                                            my $match = 0;
1498                                            if (defined($row->{"ASSOCIATED PEG"})) {
1499                                                    for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
1500                                                            if ($row->{"ASSOCIATED PEG"}->[$j] eq $orgRow->{"ASSOCIATED PEG"}->[$k]) {
1501                                                                    $match = 1;
1502                                                            }
1503                                                    }
1504                                            }
1505                                            if ($match == 0) {
1506                                                    push(@{$geneChanges},"Removed ".$orgRow->{"ASSOCIATED PEG"}->[$k]);
1507                                            }
1508                                    }
1509                            }
1510                            if ((defined($directionChange) && length($directionChange) > 0) || defined($geneChanges) && @{$geneChanges} > 0) {
1511                                    $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => [$directionChange], GeneChange => $geneChanges, Action => ["CHANGE"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1512                            }
1513                    }
1514            }
1515            #Looking for removed reactions
1516            for (my $i=0; $i < $originalReactions->size(); $i++) {
1517                    my $row = $originalReactions->get_row($i);
1518                    my $orgRow = $tbl->get_row_by_key($row->{LOAD}->[0],"LOAD");
1519                    if (!defined($orgRow)) {
1520                            $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => $row->{DIRECTIONALITY}, GeneChange => $row->{"ASSOCIATED PEG"}, Action => ["REMOVED"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1521                    }
1522            }
1523            $histTbl->save();
1524    }
1525    
1526    =head3 GapGenModel
1527    Definition:
1528            FIGMODELmodel->GapGenModel();
1529    Description:
1530            Runs the gap generation algorithm to correct a single false positive prediction. Results are loaded into a table.
1531    =cut
1532    
1533    sub GapGenModel {
1534            my ($self,$Media,$KOList,$NoKOList,$Experiment,$SolutionLimit) = @_;
1535    
1536            #Enforcing nonoptional arguments
1537            if (!defined($Media)) {
1538                    return undef;
1539            }
1540            if (!defined($KOList)) {
1541                    $KOList->[0] = "none";
1542            }
1543            if (!defined($NoKOList)) {
1544                    $NoKOList->[0] = "none";
1545            }
1546            if (!defined($Experiment)) {
1547                    $Experiment= "ReactionKO";
1548            }
1549            if (!defined($SolutionLimit)) {
1550                    $SolutionLimit = "10";
1551            }
1552    
1553            #Translating the KO lists into arrays
1554            if (ref($KOList) ne "ARRAY") {
1555                    my $temp = $KOList;
1556                    $KOList = ();
1557                    push(@{$KOList},split(/[,;]/,$temp));
1558            }
1559            my $noKOHash;
1560            if (defined($NoKOList) && ref($NoKOList) ne "ARRAY") {
1561                    my $temp = $NoKOList;
1562                    $NoKOList = ();
1563                    push(@{$NoKOList},split(/[,;]/,$temp));
1564                    foreach my $rxn (@{$NoKOList}) {
1565                            $noKOHash->{$rxn} = 1;
1566                    }
1567            }
1568    
1569            #Checking if solutions exist for the input parameters
1570            $self->aquireModelLock();
1571            my $tbl = $self->load_model_table("GapGenSolutions");
1572            my $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1573            my $solutions;
1574            if (defined($solutionRow)) {
1575                    #Checking if any solutions conform to the no KO list
1576                    foreach my $solution (@{$solutionRow->{Solutions}}) {
1577                            my @reactions = split(/,/,$solution);
1578                            my $include = 1;
1579                            foreach my $rxn (@reactions) {
1580                                    if ($rxn =~ m/(rxn\d\d\d\d\d)/) {
1581                                            if (defined($noKOHash->{$1})) {
1582                                                    $include = 0;
1583                                            }
1584                                    }
1585                            }
1586                            if ($include == 1) {
1587                                    push(@{$solutions},$solution);
1588                            }
1589                    }
1590            } else {
1591                    $solutionRow = {Media => [$Media],Experiment => [$Experiment],KOlist => [join(",",@{$KOList})]};
1592                    $tbl->add_row($solutionRow);
1593                    $self->figmodel()->database()->save_table($tbl);
1594            }
1595            $self->releaseModelLock();
1596    
1597            #Returning solution list of solutions were found
1598            if (defined($solutions) && @{$solutions} > 0) {
1599                    return $solutions;
1600            }
1601    
1602            #Getting unique filename
1603            my $Filename = $self->figmodel()->filename();
1604    
1605            #Running the gap generation
1606            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));
1607            my $ProblemReport = $self->figmodel()->LoadProblemReport($Filename);
1608            if (!defined($ProblemReport)) {
1609                    $self->figmodel()->error_message("FIGMODEL:GapGenerationAlgorithm;No problem report;".$Filename.";".$self->id().$self->selected_version().";".$Media.";".$KOList.";".$NoKOList);
1610                    return undef;
1611            }
1612    
1613            #Clearing the output folder and log file
1614            $self->figmodel()->clearing_output($Filename,"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log");
1615    
1616            #Saving the solution
1617            $self->aquireModelLock();
1618            $tbl = $self->load_model_table("GapGenSolutions");
1619            $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1620            for (my $j=0; $j < $ProblemReport->size(); $j++) {
1621                    if ($ProblemReport->get_row($j)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^)]+)/) {
1622                            my @SolutionList = split(/\|/,$1);
1623                            for (my $k=0; $k < @SolutionList; $k++) {
1624                                    if ($SolutionList[$k] =~ m/(\d+):(.+)/) {
1625                                            push(@{$solutionRow->{Solutions}},$2);
1626                                            push(@{$solutions},$2);
1627                                    }
1628                            }
1629                    }
1630            }
1631            $self->figmodel()->database()->save_table($tbl);
1632            $self->releaseModelLock();
1633    
1634            return $solutions;
1635    }
1636    
1637  =head3 datagapfill  =head3 datagapfill
1638  Definition:  Definition:
1639          success()/fail() = FIGMODELmodel->datagapfill();          success()/fail() = FIGMODELmodel->datagapfill();
# Line 1201  Line 1644 
1644          my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;          my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;
1645          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
1646          if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {          if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {
                 print "Gap filling specs:\n".join("\n",@{$GapFillingRunSpecs})."\n\n";  
1647                  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));                  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));
   
1648                  #Checking that the solution exists                  #Checking that the solution exists
1649                  if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {                  if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {
1650                          $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");                          $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!");
# Line 1213  Line 1654 
1654                  my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);                  my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef);
1655                  if (defined($TansferFileSuffix)) {                  if (defined($TansferFileSuffix)) {
1656                          system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");                          system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");
                         system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingReport.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt");  
1657                  }                  }
   
1658                  #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted                  #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted
1659                  $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");                  $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
1660                  return $GapFillingRunSpecs;                  return $GapFillResultTable;
1661          }          }
1662          if (defined($TansferFileSuffix)) {          if (defined($TansferFileSuffix)) {
1663                  $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);                  $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]);
# Line 1284  Line 1723 
1723                          push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});                          push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});
1724                  }                  }
1725                  print "Integrating solution!\n";                  print "Integrating solution!\n";
1726                  $self->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1);                  $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);
1727                  my $model = $self->get_model($self->id().$TempVersion);                  $self->PrintModelLPFile();
                 $model->PrintModelLPFile();  
1728                  #Running the model against all available experimental data                  #Running the model against all available experimental data
1729                  print "Running test model!\n";                  print "Running test model!\n";
1730                  my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");                  my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
# Line 1322  Line 1760 
1760  =cut  =cut
1761  sub status {  sub status {
1762          my ($self) = @_;          my ($self) = @_;
1763          return $self->{_data}->{status}->[0];          return $self->{_data}->status();
1764  }  }
1765    
1766  =head3 message  =head3 message
# Line 1333  Line 1771 
1771  =cut  =cut
1772  sub message {  sub message {
1773          my ($self) = @_;          my ($self) = @_;
1774          return $self->{_data}->{message}->[0];          return $self->{_data}->message();
1775  }  }
1776    
1777  =head3 set_status  =head3 set_status
# Line 1348  Line 1786 
1786  =cut  =cut
1787  sub set_status {  sub set_status {
1788          my ($self,$NewStatus,$Message) = @_;          my ($self,$NewStatus,$Message) = @_;
1789            $self->{_data}->status($NewStatus);
1790          #Getting the model row from the MODELS table          $self->{_data}->message($Message);
         $self->{_data}->{status}->[0] = $NewStatus;  
         $self->{_data}->{message}->[0] = $Message;  
         $self->figmodel()->database()->update_row("MODELS",$self->{_data},"id");  
1791          return $self->config("SUCCESS")->[0];          return $self->config("SUCCESS")->[0];
1792  }  }
1793    
# Line 1367  Line 1802 
1802  sub CreateSingleGenomeReactionList {  sub CreateSingleGenomeReactionList {
1803          my ($self,$RunGapFilling) = @_;          my ($self,$RunGapFilling) = @_;
1804    
1805            #Creating directory
1806            if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/") {
1807                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/");
1808            } elsif ($self->owner() eq "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->genome()."/") {
1809                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->genome()."/");
1810            }
1811            if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/") {
1812                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/");
1813            }
1814    
1815          #Getting genome stats          #Getting genome stats
1816          my $genomestats = $self->figmodel()->get_genome_stats($self->genome());          my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1817          my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());          my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
# Line 1381  Line 1826 
1826          }          }
1827          #Setting up needed variables          #Setting up needed variables
1828          my $OriginalModelTable = undef;          my $OriginalModelTable = undef;
1829          #Locking model status table          if ($self->status() == 0) {
         my $ModelTable = $self->figmodel()->database()->LockDBTable("MODELS");  
         my $Row = $ModelTable->get_row_by_key($self->id(),"id");  
         if (!defined($Row) || !defined($Row->{status}->[0])) {  
                 $self->figmodel()->database()->UnlockDBTable("MODELS");  
                 $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." has no row in models table!");  
                 return $self->fail();  
         } elsif ($Row->{status}->[0] == 0) {  
1830                  $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList:Model is already being built. Canceling current build.");                  $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList:Model is already being built. Canceling current build.");
                 $self->figmodel()->database()->UnlockDBTable("MODELS");  
1831                  return $self->fail();                  return $self->fail();
1832          }elsif ($Row->{status}->[0] == 1) {          }elsif ($self->status() == 1) {
1833                  $OriginalModelTable = $self->reaction_table();                  $OriginalModelTable = $self->reaction_table();
1834                  $self->ArchiveModel();                  $self->set_status(0,"Rebuilding preliminary reconstruction");
                 $Row->{status}->[0] = 0;  
                 $Row->{message}->[0] = "Rebuilding preliminary reconstruction";  
1835          } else {          } else {
1836                  $Row->{status}->[0] = 0;                  $self->set_status(0,"Preliminary reconstruction");
                 $Row->{message}->[0] = "Preliminary reconstruction";  
1837          }          }
         #Updating the status table  
         $self->figmodel()->database()->save_table($ModelTable);  
         $self->figmodel()->database()->UnlockDBTable("MODELS");  
1838          #Sorting GenomeData by gene location on the chromosome          #Sorting GenomeData by gene location on the chromosome
1839            my $ftrTbl = $self->figmodel()->database()->get_table("ROLERXNMAPPING");
1840          $FeatureTable->sort_rows("MIN LOCATION");          $FeatureTable->sort_rows("MIN LOCATION");
1841          my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);          my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);
1842          my %LastGenePosition;          my %LastGenePosition;
# Line 1476  Line 1908 
1908                          }                          }
1909                  }                  }
1910          }          }
   
1911          #Creating nonadjacent complexes          #Creating nonadjacent complexes
1912          my @ReactionList = keys(%{$ReactionHash});          my @ReactionList = keys(%{$ReactionHash});
1913          foreach my $Reaction (@ReactionList) {          foreach my $Reaction (@ReactionList) {
# Line 1582  Line 2013 
2013                  $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"]});                  $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"]});
2014          }          }
2015    
2016            #Getting feature rows for features that are lumped
2017            my @rows = $ftrTbl->get_rows_by_key("2","MASTER");
2018            for (my $i=0; $i < @rows; $i++) {
2019                    my $rxn = $rows[$i]->{REACTION}->[0];
2020                    my $role = $rows[$i]->{ROLE}->[0];
2021                    my @orgrows = $FeatureTable->get_rows_by_key($role,"ROLES");
2022                    my $genes;
2023                    for (my $j=0; $j < @orgrows; $j++) {
2024                            if ($orgrows[$j]->{ID}->[0] =~ m/(peg\.\d+)/) {
2025                                    push(@{$genes},$1);
2026                            }
2027                    }
2028                    if (defined($genes) && @{$genes} > 0) {
2029                            my $row = $NewModelTable->get_row_by_key($rxn,"LOAD");
2030                            my $newGeneAssociations;
2031                            for (my $k=0; $k < @{$genes}; $k++) {
2032                                    for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
2033                                            my $peg = $genes->[$k];
2034                                            if ($row->{"ASSOCIATED PEG"}->[$j] !~ m/$peg/) {
2035                                                    push(@{$newGeneAssociations},$row->{"ASSOCIATED PEG"}->[$j]."+".$peg);
2036                                            }
2037                                    }
2038                            }
2039                            $row->{"ASSOCIATED PEG"} = $newGeneAssociations;
2040                    }
2041            }
2042    
2043          #Adding the spontaneous and universal reactions          #Adding the spontaneous and universal reactions
2044          foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {          foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
2045                  #Getting the thermodynamic reversibility from the database                  #Getting the thermodynamic reversibility from the database
# Line 1600  Line 2058 
2058                  }                  }
2059          }          }
2060    
2061          #Checking if a biomass reaction already exists          #Creating biomass reaction for model
2062          my $BiomassReactionRow = $self->figmodel()->database()->get_row_by_key("BIOMASS TABLE",$self->id(),"MODELS");          my $biomassID = $self->BuildSpecificBiomassReaction();
2063          if (!defined($BiomassReactionRow)) {          if ($biomassID !~ m/bio\d\d\d\d\d/) {
                 $BiomassReactionRow = $self->BuildSpecificBiomassReaction();  
                 if (!defined($BiomassReactionRow)) {  
2064                          $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");                          $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
2065                          $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");                  die $self->error_message("CreateSingleGenomeReactionList: Could not generate biomass reaction!");
2066                          return $self->fail();          }
2067            #Getting the biomass reaction PPO object
2068            my $bioRxn = $self->figmodel()->database()->get_object("bof",{id=>$biomassID});
2069            if (!defined($bioRxn)) {
2070                    die $self->error_message("CreateSingleGenomeReactionList: Could not find biomass reaction ".$biomassID."!");
2071            }
2072            #Getting the list of essential reactions for biomass reaction
2073            my $ReactionList;
2074            my $essentialReactions = $bioRxn->essentialRxn();
2075            if (defined($essentialReactions) && $essentialReactions =~ m/rxn\d\d\d\d\d/) {
2076                    push(@{$ReactionList},split(/\|/,$essentialReactions));
2077                    if ($essentialReactions !~ m/$biomassID/) {
2078                            push(@{$ReactionList},$biomassID);
2079                  }                  }
2080          }          }
         my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};  
         push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);  
2081    
2082          #Adding biomass reactions to the model table          #Adding biomass reactions to the model table
2083          foreach my $BOFReaction (@{$ReactionList}) {          foreach my $BOFReaction (@{$ReactionList}) {
# Line 1641  Line 2107 
2107                  }                  }
2108          }          }
2109    
2110            #If an original model exists, we copy the gap filling solution from that model
2111            if (defined($OriginalModelTable)) {
2112                    for (my $i=0; $i < $OriginalModelTable->size(); $i++) {
2113                            if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ && $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] ne "INITIAL GAP FILLING") {
2114                                    my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD");
2115                                    if (!defined($Row)) {
2116                                            $NewModelTable->add_row($OriginalModelTable->get_row($i));
2117                                    }
2118                            }
2119                    }
2120            }
2121    
2122          #Now we compare the model to the previous model to determine if any differences exist that aren't gap filling reactions          #Now we compare the model to the previous model to determine if any differences exist that aren't gap filling reactions
2123          if (defined($OriginalModelTable)) {          if (defined($OriginalModelTable)) {
2124                  my $PerfectMatch = 1;                  my $PerfectMatch = 1;
# Line 1700  Line 2178 
2178          #Saving the new model to file          #Saving the new model to file
2179          $self->set_status(1,"Preliminary reconstruction complete");          $self->set_status(1,"Preliminary reconstruction complete");
2180          $self->figmodel()->database()->save_table($NewModelTable);          $self->figmodel()->database()->save_table($NewModelTable);
2181          $self->{_reaction_data} = $NewModelTable;          $self->reaction_table(1);
         #Clearing the previous model from the cache  
         $self->figmodel()->ClearDBModel($self->id(),1);  
2182          #Updating the model stats table          #Updating the model stats table
2183          $self->update_stats_for_build();          $self->update_stats_for_build();
2184          $self->PrintSBMLFile();          $self->PrintSBMLFile();
2185            if (defined($OriginalModelTable)) {
2186                    $self->calculate_model_changes($OriginalModelTable,"REBUILD");
2187            }
2188    
2189          #Adding model to gapfilling queue          #Adding model to gapfilling queue
2190          if (defined($RunGapFilling) && $RunGapFilling == 1) {          if (defined($RunGapFilling) && $RunGapFilling == 1) {
# Line 1724  Line 2203 
2203    
2204  sub CreateMetaGenomeReactionList {  sub CreateMetaGenomeReactionList {
2205          my ($self) = @_;          my ($self) = @_;
   
2206          #Checking if the metagenome file exists          #Checking if the metagenome file exists
2207          if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {          if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {
2208                  $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());                  $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
2209                    return $self->fail();
2210          }          }
2211          #Loading metagenome data          #Loading metagenome data
2212          my $MGRASTData = $self->figmodel()->database()->load_multiple_column_file($self->config("raw MGRAST directory")->[0].$self->genome().".summary","\t");          my $MGRASTData = $self->figmodel()->database()->load_multiple_column_file($self->config("raw MGRAST directory")->[0].$self->genome().".summary","\t");
2213          if (!defined($MGRASTData)) {          if (!defined($MGRASTData)) {
2214                  $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());                  $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome());
2215                    return $self->fail();
2216          }          }
   
2217          #Setting up needed variables          #Setting up needed variables
2218          my $OriginalModelTable = undef;          my $OriginalModelTable = undef;
   
2219          #Checking status          #Checking status
2220          if ($self->status() < 0) {          if ($self->status() < 0) {
2221                  $self->set_status(0,"Preliminary reconstruction");                  $self->set_status(0,"Preliminary reconstruction");
# Line 1749  Line 2227 
2227                  $self->ArchiveModel();                  $self->ArchiveModel();
2228                  $self->set_status(0,"Rebuilding preliminary reconstruction");                  $self->set_status(0,"Rebuilding preliminary reconstruction");
2229          }          }
2230            #Creating a hash of escores and pegs associated with each role
2231            my $rolePegHash;
2232            my $roleEscores;
2233            for (my $i=0; $i < @{$MGRASTData};$i++) {
2234                    #MD5,PEG,number of sims,role,sim e-scores,max escore,min escore,ave escore,stdev escore,ave exponent,stddev exponent
2235                    $rolePegHash->{$MGRASTData->[$i]->[3]}->{substr($MGRASTData->[$i]->[1],4)} = 1;
2236                    push(@{$roleEscores->{$MGRASTData->[$i]->[3]}},split(/;/,$MGRASTData->[$i]->[4]));
2237            }
2238          #Getting the reaction table          #Getting the reaction table
2239          my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");          my $ReactionTable = $self->figmodel()->database()->get_table("REACTIONS");
2240          #Creating model table          #Creating model table
2241          my $ModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n");          my $ModelTable = $self->create_table_prototype("ModelReactions");
2242          for (my $i=0; $i < @{$MGRASTData};$i++) {          print $ModelTable->filename();
2243                  #MD5,PEG,number of sims,role,sim e-scores          my @roles = keys(%{$rolePegHash});
2244                  my $Role = $MGRASTData->[$i]->[3];          for (my $i=0; $i < @roles; $i++) {
2245                  my $MD5 = $MGRASTData->[$i]->[0];                  my $min = -1;
2246                  my $peg = $MGRASTData->[$i]->[1];                  my $max = -1;
2247                  my $sims = $MGRASTData->[$i]->[4];                  my $count = @{$roleEscores->{$roles[$i]}};
2248                  $sims =~ s/;/,/g;                  my $ave = 0;
2249                    my $stdev = 0;
2250                    my $aveexp = 0;
2251                    my $stdevexp = 0;
2252                    for (my $j=0; $j < @{$roleEscores->{$roles[$i]}}; $j++) {
2253                            if ($roleEscores->{$roles[$i]} < $min || $min == -1) {
2254                                    $min = $roleEscores->{$roles[$i]};
2255                            }
2256                            if ($roleEscores->{$roles[$i]} > $max || $max == -1) {
2257                                    $max = $roleEscores->{$roles[$i]};
2258                            }
2259                            $ave += $roleEscores->{$roles[$i]}->[$j];
2260                            if ($roleEscores->{$roles[$i]}->[$j] =~ m/e(-\d+$)/) {
2261                                    $aveexp += $1;
2262                            }
2263                    }
2264                    $ave = $ave/$count;
2265                    $aveexp = $aveexp/$count;
2266                    for (my $j=0; $j < @{$roleEscores->{$roles[$i]}}; $j++) {
2267                            $stdev += ($roleEscores->{$roles[$i]}->[$j]-$ave)*($roleEscores->{$roles[$i]}->[$j]-$ave);
2268                            if ($roleEscores->{$roles[$i]}->[$j] =~ m/e(-\d+$)/) {
2269                                    $stdevexp += ($1-$aveexp)*($1-$aveexp);
2270                            }
2271                    }
2272                    $stdev = sqrt($stdev/$count);
2273                    $stdevexp = sqrt($stdevexp/$count);
2274                  #Checking for subsystems                  #Checking for subsystems
2275                  my $GeneSubsystems = $self->figmodel()->subsystems_of_role($Role);                  my $GeneSubsystems = $self->figmodel()->subsystems_of_role($roles[$i]);
2276                  #Checking if there are reactions associated with this role                  #Checking if there are reactions associated with this role
2277                  my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);                  my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($roles[$i]);
2278                  if ($ReactionHashArrayRef != 0) {                  if ($ReactionHashArrayRef != 0) {
2279                          foreach my $Mapping (@{$ReactionHashArrayRef}) {                          foreach my $Mapping (@{$ReactionHashArrayRef}) {
2280                                  if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {                                  if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) {
# Line 1777  Line 2286 
2286                                                                  $ReactionRow = {"LOAD" => [$Mapping->{"REACTION"}->[0]],"DIRECTIONALITY" => [$self->figmodel()->reversibility_of_reaction($Mapping->{"REACTION"}->[0])],"COMPARTMENT" => ["c"]};                                                                  $ReactionRow = {"LOAD" => [$Mapping->{"REACTION"}->[0]],"DIRECTIONALITY" => [$self->figmodel()->reversibility_of_reaction($Mapping->{"REACTION"}->[0])],"COMPARTMENT" => ["c"]};
2287                                                                  $ModelTable->add_row($ReactionRow);                                                                  $ModelTable->add_row($ReactionRow);
2288                                                          }                                                          }
2289                                                          push(@{$ReactionRow->{"ASSOCIATED PEG"}},substr($peg,4));                                                          my %pegHash = %{$rolePegHash->{$roles[$i]}};
2290                                                          push(@{$ReactionRow->{"REFERENCE"}},$MD5.":".$Role);                                                          if (defined($ReactionRow->{"ASSOCIATED PEG"})) {
2291                                                          push(@{$ReactionRow->{"CONFIDENCE"}},$sims);                                                                  for (my $j=0; $j < @{$ReactionRow->{"ASSOCIATED PEG"}}; $j++) {
2292                                                                            $pegHash{$ReactionRow->{"ASSOCIATED PEG"}->[$j]} = 1;
2293                                                                    }
2294                                                            }
2295                                                            delete $ReactionRow->{"ASSOCIATED PEG"};
2296                                                            push(@{$ReactionRow->{"ASSOCIATED PEG"}},keys(%pegHash));
2297                                                            push(@{$ReactionRow->{"REFERENCE"}},$count.":".$ave.":".$stdev.":".$aveexp.":".$stdevexp.":".$min.":".$max);
2298                                                          if (defined($GeneSubsystems)) {                                                          if (defined($GeneSubsystems)) {
2299                                                                  push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});                                                                  push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});
2300                                                          }                                                          }
# Line 1823  Line 2338 
2338          }          }
2339    
2340          #Clearing the previous model from the cache          #Clearing the previous model from the cache
2341          $self->figmodel()->ClearDBModel($self->id(),1);          $self->figmodel()->database()->ClearDBModel($self->id(),1);
2342          $ModelTable->save();          $ModelTable->save();
2343    
2344          return $self->success();          return $self->success();
# Line 1839  Line 2354 
2354  sub ArchiveModel {  sub ArchiveModel {
2355          my ($self) = @_;          my ($self) = @_;
2356    
         #Making sure the model exists  
         if (!defined($self->stats())) {  
                 $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Could not find specified ".$self->id()." in database!");  
                 return $self->fail();  
         }  
   
2357          #Checking that the model file exists          #Checking that the model file exists
2358          if (!(-e $self->filename())) {          if (!(-e $self->filename())) {
2359                  $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");                  $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!");
# Line 1895  Line 2404 
2404          Runs microarray analysis attempting to turn off genes that are inactive in the microarray          Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2405  =cut  =cut
2406  sub run_microarray_analysis {  sub run_microarray_analysis {
2407          my ($self,$media,$jobid,$index,$genecall) = @_;          my ($self,$media,$label,$index,$genecall) = @_;
2408          $genecall =~ s/_/:/g;          $genecall =~ s/_/:/g;
2409          $genecall =~ s/\//;/g;          $genecall =~ s/\//;/g;
2410          #print "\n\n".$genecall."\n\n";          my $uniqueFilename = $self->figmodel()->filename();
2411          my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($jobid,$self->id(),$media,["MFA","MicroarrayAssertions"],{"Microarray assertions" => $self->id().";".$index.";".$genecall,"MFASolver" => "CPLEX","Network output location" => "/scratch/"},"MicroarrayAnalysis-".$jobid.".txt",undef,$self->selected_version());          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());
         #print $command."\n";  
2412          system($command);          system($command);
2413          #system("/home/chenry/Software/MFAToolkitRepository/Linux/mfatoolkit resetparameter \"user bounds filename\" \"Carbon-D-Glucose.txt\" resetparameter output_folder \"32749.354149.0/\" resetparameter \"Microarray assertions\" \"Seed83333.1;0;peg.3130:-1;peg.4035:-1\" resetparameter \"MFASolver\" \"CPLEX\" resetparameter \"Network output location\" \"/scratch/\" LoadCentralSystem \"/vol/model-dev/MODEL_DEV_DB/Models/83333.1/Seed83333.1.txt\" > \"/vol/model-dev/MODEL_DEV_DB/ReactionDB/log/MicroarrayAnalysis-32749.354149.0.txt\"");          my $filename = $self->figmodel()->config("MFAToolkit output directory")->[0].$uniqueFilename."/MicroarrayOutput-".$index.".txt";
2414            if (-e $filename) {
2415                    my $output = $self->figmodel()->database()->load_single_column_file($filename);
2416                    if (defined($output->[0])) {
2417                            my @array = split(/;/,$output->[0]);
2418                            $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2419                            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]);
2420                    }
2421                    print STDERR $filename." is empty!";
2422            }
2423            print STDERR $filename." not found!";
2424            $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2425    
2426            return undef;
2427  }  }
2428    
2429  =head3 find_minimal_pathways  =head3 find_minimal_pathways
# Line 1912  Line 2433 
2433          Runs microarray analysis attempting to turn off genes that are inactive in the microarray          Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2434  =cut  =cut
2435  sub find_minimal_pathways {  sub find_minimal_pathways {
2436          my ($self,$media,$objective,$solutionnum) = @_;          my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2437    
2438          #Setting default media          #Setting default media
2439          if (!defined($media)) {          if (!defined($media)) {
# Line 1924  Line 2445 
2445                  $solutionnum = "5";                  $solutionnum = "5";
2446          }          }
2447    
2448            #Setting additional exchange fluxes
2449            if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2450                    if ($self->id() eq "iAF1260") {
2451                            $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2452                    } else {
2453                            $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2454                    }
2455            }
2456    
2457          #Translating objective          #Translating objective
2458          my $objectivestring;          my $objectivestring;
2459          if ($objective eq "ALL") {          if ($objective eq "ALL") {
# Line 1944  Line 2474 
2474                          }                          }
2475                  }                  }
2476                  for (my $i=0; $i < @objectives; $i++) {                  for (my $i=0; $i < @objectives; $i++) {
2477                          $self->find_minimal_pathways($media,$objectives[$i])                          $self->find_minimal_pathways($media,$objectives[$i]);
2478                  }                  }
2479                    return;
2480          } elsif ($objective eq "ENERGY") {          } elsif ($objective eq "ENERGY") {
2481                  $objectivestring = "MAX;FLUX;rxn00062;c;1";                  $objectivestring = "MAX;FLUX;rxn00062;c;1";
2482          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2483                    if ($objective =~ m/\[(\w)\]/) {
2484                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2485                            $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2486                    } else {
2487                  $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";                  $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2488                            $additionalexchange .= ";".$objective."[c]:-100:0";
2489                    }
2490            } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2491                    my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2492                    for (my $i=0; $i < @{$Products};$i++) {
2493                            my $temp = $Products->[$i]->{"DATABASE"}->[0];
2494                            if ($additionalexchange !~ m/$temp/) {
2495                                    #$additionalexchange .= ";".$temp."[c]:-100:0";
2496                            }
2497                    }
2498                    for (my $i=0; $i < @{$Reactants};$i++) {
2499                            print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2500                            $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2501                            print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2502                    }
2503                    return;
2504          }          }
2505    
2506          #Setting additional exchange fluxes          #Adding additional drains
2507          my $additionalexchange = $self->config("default exchange fluxes")->[0];          if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
         if ($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") {  
2508                  $additionalexchange .= ";cpd15666[c]:0:100";                  $additionalexchange .= ";cpd15666[c]:0:100";
2509          } elsif ($objective eq "cpd11493") {          } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2510                  $additionalexchange .= ";cpd12370[c]:0:100";                  $additionalexchange .= ";cpd12370[c]:0:100";
2511          } elsif ($objective eq "cpd00166") {          } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2512                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2513          }          }
2514    
2515          #Running MFAToolkit          #Running MFAToolkit
2516          my $filename = $self->figmodel()->filename();          my $filename = $self->figmodel()->filename();
2517          my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["MFA"],{"Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "CPLEX","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version());          my $command;
2518            if (defined($AllReversible) && $AllReversible == 1) {
2519                    $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());
2520            } else {
2521                    $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());
2522            }
2523          system($command);          system($command);
2524    
2525          #Loading problem report          #Loading problem report
# Line 1972  Line 2527 
2527          #Clearing output          #Clearing output
2528          $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");          $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2529          if (!defined($results)) {          if (!defined($results)) {
2530                    print STDERR $objective." pathway results not found!\n";
2531                  return;                  return;
2532          }          }
2533    
# Line 1983  Line 2539 
2539                  @Array = /\d+:([^\|]+)\|/g;                  @Array = /\d+:([^\|]+)\|/g;
2540          }          }
2541    
2542          # Storing data in a figmodel table          #Writing output to file
2543          my $TableObject;          $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
         if (-e $self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt") {  
                 $TableObject->load_table($self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",";","|",0,["OBJECTIVE"]);  
         } else {  
                 $TableObject = FIGMODELTable->new(["OBJECTIVE","REACTIONS"],$self->directory()."MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",["OBJECTIVE"],";","|",undef);  
         }  
         my $tablerow = $TableObject->get_row_by_key($objective,"OBJECTIVE",1);  
         push(@{$tablerow->{"REACTIONS"}},@Array);  
         $TableObject->save();  
2544  }  }
2545    
2546  =head3 calculate_growth  =head3 find_minimal_pathways
2547  Definition:  Definition:
2548          string::growth = FIGMODELmodel->calculate_growth(string:media);          int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
2549  Description:  Description:
2550          Calculating growth in the input media          Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2551  =cut  =cut
2552  sub calculate_growth {  sub find_minimal_pathways_two {
2553          my ($self,$Media) = @_;          my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2554          my $UniqueFilename = $self->figmodel()->filename();  
2555          system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],{"optimize metabolite production if objective is zero" => 1},$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version()));          #Setting default media
2556          my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);          if (!defined($media)) {
2557          my $Result;                  $media = "Complete";
         if (defined($ProblemReport)) {  
                 my $Row = $ProblemReport->get_row(0);  
                 if (defined($Row) && defined($Row->{"Objective"}->[0])) {  
                         if ($Row->{"Objective"}->[0] < 0.00000001) {  
                                 $Result = "NOGROWTH:".$Row->{"Individual metabolites with zero production"}->[0];  
                         } else {  
                                 $Result = $Row->{"Objective"}->[0];  
2558                          }                          }
2559    
2560            #Setting default solution number
2561            if (!defined($solutionnum)) {
2562                    $solutionnum = "5";
2563                  }                  }
2564    
2565            #Setting additional exchange fluxes
2566            if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2567                    if ($self->id() eq "iAF1260") {
2568                            $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2569                    } else {
2570                            $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2571          }          }
         return $Result;  
2572  }  }
2573    
2574  =head3 classify_model_reactions          #Translating objective
2575  Definition:          my $objectivestring;
2576          (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);          if ($objective eq "ALL") {
2577  Description:                  #Getting the list of universal building blocks
2578          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.                  my $buildingblocks = $self->config("universal building blocks");
2579          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".                  my @objectives = keys(%{$buildingblocks});
2580          Possible values for "Class" include:                  #Getting the nonuniversal building blocks
2581          1.) Positive: these reactions are essential in the forward direction.                  my $otherbuildingblocks = $self->config("nonuniversal building blocks");
2582          2.) Negative: these reactions are essential in the reverse direction.                  my @array = keys(%{$otherbuildingblocks});
2583          3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.                  if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
2584          4.) Negative variable: these reactions are nonessential, but they only ever proceed in the reverse direction.                          my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
2585          5.) Variable: these reactions are nonessential and proceed in the forward or reverse direction.                          if (defined($equation)) {
2586          6.) Blocked: these reactions never carry any flux at all in the media condition tested.                                  for (my $i=0; $i < @array; $i++) {
2587          7.) Dead: these reactions are disconnected from the network.                                          if (CORE::index($equation,$array[$i]) > 0) {
2588  =cut                                                  push(@objectives,$array[$i]);
2589  sub classify_model_reactions {                                          }
2590          my ($self,$Media) = @_;                                  }
2591                            }
2592          #Getting unique file for printing model output                  }
2593          my $UniqueFilename = $self->figmodel()->filename();                  for (my $i=0; $i < @objectives; $i++) {
2594          #Running the MFAToolkit                          $self->find_minimal_pathways($media,$objectives[$i]);
2595          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()));                  }
2596          #Reading in the output bounds file                  return;
2597          my ($ReactionTB,$CompoundTB,$DeadCompounds,$DeadEndCompounds,$DeadReactions);          } elsif ($objective eq "ENERGY") {
2598          if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") {                  $objectivestring = "MAX;FLUX;rxn00062;c;1";
2599                  $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]);          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2600                    if ($objective =~ m/\[(\w)\]/) {
2601                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2602                            $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2603                    } else {
2604                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2605                            $additionalexchange .= ";".$objective."[c]:-100:0";
2606                    }
2607            } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2608                    my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2609                    for (my $i=0; $i < @{$Products};$i++) {
2610                            my $temp = $Products->[$i]->{"DATABASE"}->[0];
2611                            if ($additionalexchange !~ m/$temp/) {
2612                                    #$additionalexchange .= ";".$temp."[c]:-100:0";
2613                            }
2614                    }
2615                    for (my $i=0; $i < @{$Reactants};$i++) {
2616                            print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2617                            $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2618                            print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2619                    }
2620                    return;
2621            }
2622    
2623            #Adding additional drains
2624            if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
2625                    $additionalexchange .= ";cpd15666[c]:0:100";
2626            } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2627                    $additionalexchange .= ";cpd12370[c]:0:100";
2628            } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2629                    $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2630            }
2631    
2632            #Running MFAToolkit
2633            my $filename = $self->figmodel()->filename();
2634            my $command;
2635            if (defined($AllReversible) && $AllReversible == 1) {
2636                    $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());
2637            } else {
2638                    $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());
2639            }
2640            print $command."\n";
2641            system($command);
2642    
2643            #Loading problem report
2644            my $results = $self->figmodel()->LoadProblemReport($filename);
2645            #Clearing output
2646            $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2647            if (!defined($results)) {
2648                    print STDERR $objective." pathway results not found!\n";
2649                    return;
2650            }
2651    
2652            #Parsing output
2653            my @Array;
2654            my $row = $results->get_row(1);
2655            if (defined($row->{"Notes"}->[0])) {
2656                    $_ = $row->{"Notes"}->[0];
2657                    @Array = /\d+:([^\|]+)\|/g;
2658            }
2659    
2660            #Writing output to file
2661            $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
2662    }
2663    
2664    sub combine_minimal_pathways {
2665            my ($self) = @_;
2666    
2667            my $tbl;
2668            if (-e $self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl") {
2669                    $tbl = FIGMODELTable::load_table($self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",";","|",0,["Objective","Media","Reversible"]);
2670            } else {
2671                    $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"],";","|");
2672            }
2673            my @files = glob($self->directory()."MinimalPathways-*");
2674            for (my $i=0; $i < @files;$i++) {
2675                    if ($files[$i] =~ m/MinimalPathways\-(\S+)\-(cpd\d\d\d\d\d)\-(\w+)\-(\d)\-/ || $files[$i] =~ m/MinimalPathways\-(\S+)\-(ENERGY)\-(\w+)\-(\d)\-/) {
2676                            my $reactions = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2677                            if (defined($reactions) && @{$reactions} > 0 && length($reactions->[0]) > 0) {
2678                                    my $newrow = {"Objective"=>[$2],"Media"=>[$1],"Reversible"=>[$4]};
2679                                    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");
2680                                    if (!defined($row)) {
2681                                            $row = $tbl->add_row($newrow);
2682                                    }
2683                                    $row->{Reactions} = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2684                                    delete($row->{"Shortest path"});
2685                                    delete($row->{"Number of essentials"});
2686                                    delete($row->{"Essentials"});
2687                                    delete($row->{"Length"});
2688                                    for (my $j=0; $j < @{$row->{Reactions}}; $j++) {
2689                                            my @array = split(/,/,$row->{Reactions}->[$j]);
2690                                            $row->{"Length"}->[$j] = @array;
2691                                            if (!defined($row->{"Shortest path"}->[0]) || $row->{"Length"}->[$j] < $row->{"Shortest path"}->[0]) {
2692                                                    $row->{"Shortest path"}->[0] = $row->{"Length"}->[$j];
2693                                            }
2694                                            $row->{"Number of essentials"}->[0] = 0;
2695                                            for (my $k=0; $k < @array;$k++) {
2696                                                    if ($array[$k] =~ m/(rxn\d\d\d\d\d)/) {
2697                                                            my $class = $self->get_reaction_class($1,1);
2698                                                            my $temp = $row->{Media}->[0].":Essential";
2699                                                            if ($class =~ m/$temp/) {
2700                                                                    $row->{"Number of essentials"}->[$j]++;
2701                                                                    if (!defined($row->{"Essentials"}->[$j]) && length($row->{"Essentials"}->[$j]) > 0) {
2702                                                                            $row->{"Essentials"}->[$j] = $array[$k];
2703                                                                    } else {
2704                                                                            $row->{"Essentials"}->[$j] .= ",".$array[$k];
2705                                                                    }
2706                                                            }
2707                                                    }
2708                                            }
2709                                    }
2710                            }
2711                    }
2712            }
2713            $tbl->save();
2714    }
2715    
2716    =head3 calculate_growth
2717    Definition:
2718            string::growth = FIGMODELmodel->calculate_growth(string:media);
2719    Description:
2720            Calculating growth in the input media
2721    =cut
2722    sub calculate_growth {
2723            my ($self,$Media,$outputDirectory,$InParameters,$saveLPFile) = @_;
2724            #Setting the Media
2725            if (!defined($Media) || length($Media) == 0) {
2726                    $Media = $self->autocompleteMedia();
2727            }
2728            #Setting parameters for the run
2729            my $DefaultParameters = $self->figmodel()->defaultParameters();
2730            if (defined($InParameters)) {
2731                    my @parameters = keys(%{$InParameters});
2732                    for (my $i=0; $i < @parameters; $i++) {
2733                            $DefaultParameters->{$parameters[$i]} = $InParameters->{$parameters[$i]};
2734                    }
2735            }
2736            $DefaultParameters->{"optimize metabolite production if objective is zero"} = 1;
2737            #Setting filenames
2738            my $UniqueFilename = $self->figmodel()->filename();
2739            if (!defined($outputDirectory)) {
2740                    $outputDirectory = $self->config("database message file directory")->[0];
2741            }
2742            my $fluxFilename = $outputDirectory."Fluxes-".$self->id()."-".$Media.".txt";
2743            my $cpdFluxFilename = $outputDirectory."CompoundFluxes-".$self->id()."-".$Media.".txt";
2744            #Running FBA
2745            #print $self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],$DefaultParameters,$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version())."\n";
2746            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],$DefaultParameters,$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version()));
2747            #Saving LP file if requested
2748            if (defined($saveLPFile) && $saveLPFile == 1 && -e $self->figmodel()->{"MFAToolkit output directory"}->[0].$UniqueFilename."/CurrentProblem.lp") {
2749                    system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/CurrentProblem.lp ".$self->directory().$self->id().".lp");
2750            }
2751            my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);
2752            my $Result;
2753            if (defined($ProblemReport)) {
2754                    my $Row = $ProblemReport->get_row(0);
2755                    if (defined($Row) && defined($Row->{"Objective"}->[0])) {
2756                            if ($Row->{"Objective"}->[0] < 0.00000001 || $Row->{"Objective"}->[0] == 1e7) {
2757                                    $Result = "NOGROWTH";
2758                                    if (defined($Row->{"Individual metabolites with zero production"}->[0]) && $Row->{"Individual metabolites with zero production"}->[0] =~ m/cpd\d\d\d\d\d/) {
2759                                            $Result .= ":".$Row->{"Individual metabolites with zero production"}->[0];
2760                                    }
2761                            } else {
2762                                    if (-e $self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionReactionData0.txt") {
2763                                            system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionReactionData0.txt ".$fluxFilename);
2764                                            system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionCompoundData0.txt ".$cpdFluxFilename);
2765                                    }
2766                                    $Result = $Row->{"Objective"}->[0];
2767                            }
2768                    }
2769            }
2770            #Deleting files if necessary
2771            if ($self->figmodel()->config("preserve all log files")->[0] ne "yes") {
2772                    $self->figmodel()->cleardirectory($UniqueFilename);
2773                    unlink($self->figmodel()->config("database message file directory")->[0].$self->id()."-".$Media."-GrowthTest.txt");
2774            }
2775            #Returning result
2776            return $Result;
2777    }
2778    
2779    =head3 classify_model_reactions
2780    Definition:
2781            (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);
2782    Description:
2783            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.
2784            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".
2785            Possible values for "Class" include:
2786            1.) Positive: these reactions are essential in the forward direction.
2787            2.) Negative: these reactions are essential in the reverse direction.
2788            3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.
2789            4.) Negative variable: these reactions are nonessential, but they only ever proceed in the reverse direction.
2790            5.) Variable: these reactions are nonessential and proceed in the forward or reverse direction.
2791            6.) Blocked: these reactions never carry any flux at all in the media condition tested.
2792            7.) Dead: these reactions are disconnected from the network.
2793    =cut
2794    sub classify_model_reactions {
2795            my ($self,$Media,$SaveChanges) = @_;
2796    
2797            #Getting unique file for printing model output
2798            my $UniqueFilename = $self->figmodel()->filename();
2799            #Running the MFAToolkit
2800            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()));
2801            #Reading in the output bounds file
2802            my ($ReactionTB,$CompoundTB,$DeadCompounds,$DeadEndCompounds,$DeadReactions);
2803            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") {
2804                    $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]);
2805          }          }
2806          if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt") {          if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt") {
2807                  $CompoundTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt",";","|",1,["DATABASE ID"]);                  $CompoundTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt",";","|",1,["DATABASE ID"]);
# Line 2123  Line 2879 
2879                                  $CpdRow->{MEDIA}->[$index] = $Media;                                  $CpdRow->{MEDIA}->[$index] = $Media;
2880                          }                          }
2881                  }                  }
2882                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2883                  $cpdclasstable->save();                  $cpdclasstable->save();
2884          }          }
2885            }
2886          if (defined($ReactionTB)) {          if (defined($ReactionTB)) {
2887                  for (my $i=0; $i < $ReactionTB->size(); $i++) {                  for (my $i=0; $i < $ReactionTB->size(); $i++) {
2888                          my $Row = $ReactionTB->get_row($i);                          my $Row = $ReactionTB->get_row($i);
# Line 2179  Line 2937 
2937                                  $RxnRow->{MEDIA}->[$index] = $Media;                                  $RxnRow->{MEDIA}->[$index] = $Media;
2938                          }                          }
2939                  }                  }
2940                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2941                  $rxnclasstable->save();                  $rxnclasstable->save();
2942          }          }
2943            }
2944          return ($rxnclasstable,$cpdclasstable);          return ($rxnclasstable,$cpdclasstable);
2945  }  }
2946    
# Line 2212  Line 2972 
2972    
2973          #Running simulations          #Running simulations
2974          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");          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");
         print "simulation complete for ".$self->id().$self->selected_version()."\n";  
2975          #Parsing the results          #Parsing the results
2976          my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);          my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef);
2977          if (!defined($Results)) {          if (!defined($Results)) {
# Line 2619  Line 3378 
3378                                          if (length($GapFillingRunSpecs) > 0) {                                          if (length($GapFillingRunSpecs) > 0) {
3379                                                  $GapFillingRunSpecs .= ";";                                                  $GapFillingRunSpecs .= ";";
3380                                          }                                          }
3381                                          $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].".txt:".$HeadingDataArray[3];                                          $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].":".$HeadingDataArray[3];
3382                                  } else {                                  } else {
3383                                          $SolutionExistedCount++;                                          $SolutionExistedCount++;
3384                                  }                                  }
# Line 2644  Line 3403 
3403          my $SolutionsFound = 0;          my $SolutionsFound = 0;
3404          my $GapFillingArray;          my $GapFillingArray;
3405          push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));          push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));
3406          $self->datagapfill($GapFillingArray);          my $GapFillingResults = $self->datagapfill($GapFillingArray,"GFS");
3407            if (defined($GapFillingResults)) {
3408                    $SolutionsFound = 1;
3409            }
3410    
3411          if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {          if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {
3412                  #Printing previous solutions to GFS file                  #Printing previous solutions to GFS file
# Line 2667  Line 3429 
3429          return $self->success();          return $self->success();
3430  }  }
3431    
3432  =head3 BuildSpecificBiomassReaction  =head3 SolutionReconciliation
3433  Definition:  Definition:
3434          FIGMODELmodel->BuildSpecificBiomassReaction();          FIGMODELmodel->SolutionReconciliation();
3435  Description:  Description:
3436            This is a wrapper for running the solution reconciliation algorithm on any model in the database.
3437            The algorithm performs a reconciliation of any gap filling solutions to identify the combination of solutions that results in the optimal model.
3438            This function prints out one output file in the Model directory: ReconciliationOutput.txt: this is a summary of the results of the reconciliation analysis
3439  =cut  =cut
 sub BuildSpecificBiomassReaction {  
         my ($self) = @_;  
3440    
3441          my $biomassrxn;  sub SolutionReconciliation {
3442          my $OrganismID = $self->genome();          my ($self,$GapFill,$Stage) = @_;
3443          #Checking for a biomass override  
3444          if (defined($self->config("biomass reaction override")->{$OrganismID})) {          #Setting the output filenames
3445                  $biomassrxn = $self->config("biomass reaction override")->{$OrganismID};          my $OutputFilename;
3446                  print "Overriding biomass template and selecting ".$biomassrxn." for ".$OrganismID.".\n";          my $OutputFilenameTwo;
3447          } else {#Creating biomass reaction from the template          if ($GapFill == 1) {
3448                  #Getting the genome stats                  $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GFReconciliation.txt";
3449                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());                  $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GFSRS.txt";
3450                  my $Class = $genomestats->{CLASS}->[0];          } else {
3451                  my $Name = $genomestats->{NAME}->[0];                  $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GGReconciliation.txt";
3452                    $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GGSRS.txt";
3453            }
3454    
3455            #In stage one, we run the reconciliation and create a test file to check combined solution performance
3456            if (!defined($Stage) || $Stage == 1) {
3457                    my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3458                    my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3459                    $Row->{"GF RECONCILATION TIMING"}->[0] = time()."-";
3460                    $GrowMatchTable->save();
3461                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3462    
3463                    #Getting a unique filename
3464                    my $UniqueFilename = $self->figmodel()->filename();
3465    
3466                    #Copying over the necessary files
3467                    if ($GapFill == 1) {
3468                            if (!-e $self->directory().$self->id().$self->selected_version()."-GFEM.txt") {
3469                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GFEM.txt file not found. Could not reconcile!";
3470                                    return 0;
3471                            }
3472                            if (!-e $self->directory().$self->id().$self->selected_version()."-OPEM.txt") {
3473                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-OPEM.txt file not found. Could not reconcile!";
3474                                    return 0;
3475                            }
3476                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GFEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GFEM.txt");
3477                            system("cp ".$self->directory().$self->id().$self->selected_version()."-OPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3478                            #Backing up and deleting the existing reconciliation file
3479                            if (-e $OutputFilename) {
3480                                    system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGFReconciliation.txt");
3481                                    unlink($OutputFilename);
3482                            }
3483                    } else {
3484                            if (!-e $self->directory().$self->id().$self->selected_version()."-GGEM.txt") {
3485                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGEM.txt file not found. Could not reconcile!";
3486                                    return 0;
3487                            }
3488                            if (!-e $self->directory().$self->id().$self->selected_version()."-GGOPEM.txt") {
3489                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt file not found. Could not reconcile!";
3490                                    return 0;
3491                            }
3492                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GGEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GGEM.txt");
3493                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3494                            #Backing up and deleting the existing reconciliation file
3495                            if (-e $OutputFilename) {
3496                                    system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGGReconciliation.txt");
3497                                    unlink($OutputFilename);
3498                            }
3499                    }
3500    
3501                    #Running the reconciliation
3502                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),"NONE",["SolutionReconciliation"],{"Solution data for model optimization" => $UniqueFilename},"Reconciliation".$UniqueFilename.".log",undef,$self->selected_version()));
3503                    $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3504                    $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3505                    $Row->{"GF RECONCILATION TIMING"}->[0] .= time();
3506                    $GrowMatchTable->save();
3507                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3508    
3509                    #Loading the problem report from the reconciliation run
3510                    my $ReconciliatonOutput = $self->figmodel()->LoadProblemReport($UniqueFilename);
3511                    print $UniqueFilename."\n";
3512                    #Clearing output files
3513                    $self->figmodel()->clearing_output($UniqueFilename,"Reconciliation".$UniqueFilename.".log");
3514                    $ReconciliatonOutput->save("/home/chenry/Test.txt");
3515    
3516                    #Checking the a problem report was found and was loaded
3517                    if (!defined($ReconciliatonOutput) || $ReconciliatonOutput->size() < 1 || !defined($ReconciliatonOutput->get_row(0)->{"Notes"}->[0])) {
3518                            print STDERR "FIGMODEL:SolutionReconciliation: MFAToolkit output from SolutionReconciliation of ".$self->id()." not found!\n\n";
3519                            return 0;
3520                    }
3521    
3522                    #Processing the solutions
3523                    my $SolutionCount = 0;
3524                    my $ReactionSetHash;
3525                    my $SingleReactionHash;
3526                    my $ReactionDataHash;
3527                    for (my $n=0; $n < $ReconciliatonOutput->size(); $n++) {
3528                            if (defined($ReconciliatonOutput->get_row($n)->{"Notes"}->[0]) && $ReconciliatonOutput->get_row($n)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^;]+)/) {
3529                                    #Breaking up the solution into reaction sets
3530                                    my @ReactionSets = split(/\|/,$1);
3531                                    #Creating reaction lists for each set
3532                                    my $SolutionHash;
3533                                    for (my $i=0; $i < @ReactionSets; $i++) {
3534                                            if (length($ReactionSets[$i]) > 0) {
3535                                                    my @Alternatives = split(/:/,$ReactionSets[$i]);
3536                                                    for (my $j=1; $j < @Alternatives; $j++) {
3537                                                            if (length($Alternatives[$j]) > 0) {
3538                                                                    push(@{$SolutionHash->{$Alternatives[$j]}},$Alternatives[0]);
3539                                                            }
3540                                                    }
3541                                                    if (@Alternatives == 1) {
3542                                                            $SingleReactionHash->{$Alternatives[0]}->{$SolutionCount} = 1;
3543                                                            if (!defined($SingleReactionHash->{$Alternatives[0]}->{"COUNT"})) {
3544                                                                    $SingleReactionHash->{$Alternatives[0]}->{"COUNT"} = 0;
3545                                                            }
3546                                                            $SingleReactionHash->{$Alternatives[0]}->{"COUNT"}++;
3547                                                    }
3548                                            }
3549                                    }
3550                                    #Identifying reactions sets and storing the sets in the reactions set hash
3551                                    foreach my $Solution (keys(%{$SolutionHash})) {
3552                                            my $SetKey = join(",",sort(@{$SolutionHash->{$Solution}}));
3553                                            if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount})) {
3554                                                    $ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount} = 1;
3555                                                    if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"})) {
3556                                                            $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"} = 0;
3557                                                    }
3558                                                    $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"}++;
3559                                            }
3560                                            $ReactionSetHash->{$SetKey}->{$Solution}->{$SolutionCount} = 1;
3561                                            if (!defined($ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"})) {
3562                                                    $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"} = 0;
3563                                            }
3564                                            $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"}++;
3565                                    }
3566                                    $SolutionCount++;
3567                            }
3568                    }
3569    
3570                    #Handling the scenario where no solutions were found
3571                    if ($SolutionCount == 0) {
3572                            print STDERR "FIGMODEL:SolutionReconciliation: Reconciliation unsuccessful. No solution found.\n\n";
3573                            return 0;
3574                    }
3575    
3576                    #Printing results without solution performance figures. Also printing solution test file
3577                    open (RECONCILIATION, ">$OutputFilename");
3578                    #Printing the file heading
3579                    print RECONCILIATION "DATABASE;DEFINITION;REVERSIBLITY;DELTAG;DIRECTION;NUMBER OF SOLUTIONS";
3580                    for (my $i=0; $i < $SolutionCount; $i++) {
3581                            print RECONCILIATION ";Solution ".$i;
3582                    }
3583                    print RECONCILIATION "\n";
3584                    #Printing the singlet reactions first
3585                    my $Solutions;
3586                    print RECONCILIATION "SINGLET REACTIONS\n";
3587                    my @SingletReactions = keys(%{$SingleReactionHash});
3588                    for (my $j=0; $j < $SolutionCount; $j++) {
3589                            $Solutions->[$j]->{"BASE"} = $j;
3590                    }
3591                    for (my $i=0; $i < @SingletReactions; $i++) {
3592                            my $ReactionData;
3593                            if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3594                                    $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3595                            } else {
3596                                    my $Direction = substr($SingletReactions[$i],0,1);
3597                                    if ($Direction eq "+") {
3598                                            $Direction = "=>";
3599                                    } else {
3600                                            $Direction = "<=";
3601                                    }
3602                                    my $Reaction = substr($SingletReactions[$i],1);
3603                                    $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3604                                    $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3605                                    $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3606                                    if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3607                                            $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3608                                    }
3609                                    if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3610                                            $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3611                                    }
3612                                    if (!defined($ReactionData->{"DELTAG"}->[0])) {
3613                                            $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3614                                    }
3615                                    $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3616                            }
3617                            print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$SingleReactionHash->{$SingletReactions[$i]}->{"COUNT"};
3618                            for (my $j=0; $j < $SolutionCount; $j++) {
3619                                    print RECONCILIATION ";";
3620                                    if (defined($SingleReactionHash->{$SingletReactions[$i]}->{$j})) {
3621                                            $Solutions->[$j]->{$SingletReactions[$i]} = 1;
3622                                            $Solutions->[$j]->{"BASE"} = $j;
3623                                            print RECONCILIATION "|".$j."|";
3624                                    }
3625                            }
3626                            print RECONCILIATION "\n";
3627                    }
3628                    #Printing the reaction sets with alternatives
3629                    print RECONCILIATION "Reaction sets with alternatives\n";
3630                    my @ReactionSets = keys(%{$ReactionSetHash});
3631                    foreach my $ReactionSet (@ReactionSets) {
3632                            my $NewSolutions;
3633                            my $BaseReactions;
3634                            my $AltList = [$ReactionSet];
3635                            push(@{$AltList},keys(%{$ReactionSetHash->{$ReactionSet}}));
3636                            for (my $j=0; $j < @{$AltList}; $j++) {
3637                                    my $CurrentNewSolutions;
3638                                    my $Index;
3639                                    if ($j == 0) {
3640                                            print RECONCILIATION "NEW SET\n";
3641                                    } elsif ($AltList->[$j] ne $ReactionSet) {
3642                                            print RECONCILIATION "ALTERNATIVE SET\n";
3643                                            #For each base solution in which this set is represented, we copy the base solution to the new solution
3644                                            my $NewSolutionCount = 0;
3645                                            for (my $k=0; $k < $SolutionCount; $k++) {
3646                                                    if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3647                                                            if (defined($Solutions)) {
3648                                                                    $Index->{$k} = @{$Solutions} + $NewSolutionCount;
3649                                                            } else {
3650                                                                    $Index->{$k} = $NewSolutionCount;
3651                                                            }
3652                                                            if (defined($NewSolutions) && @{$NewSolutions} > 0) {
3653                                                                    $Index->{$k} += @{$NewSolutions};
3654                                                            }
3655                                                            $CurrentNewSolutions->[$NewSolutionCount] = {};
3656                                                            foreach my $Reaction (keys(%{$Solutions->[$k]})) {
3657                                                                    $CurrentNewSolutions->[$NewSolutionCount]->{$Reaction} = $Solutions->[$k]->{$Reaction};
3658                                                            }
3659                                                            $NewSolutionCount++;
3660                                                    }
3661                                            }
3662                                    }
3663                                    if ($j == 0 || $AltList->[$j] ne $ReactionSet) {
3664                                            my @SingletReactions = split(/,/,$AltList->[$j]);
3665                                            for (my $i=0; $i < @SingletReactions; $i++) {
3666                                                    #Adding base reactions to base solutions and set reactions the new solutions
3667                                                    if ($j == 0) {
3668                                                            push(@{$BaseReactions},$SingletReactions[$i]);
3669                                                    } else {
3670                                                            for (my $k=0; $k < @{$CurrentNewSolutions}; $k++) {
3671                                                                    $CurrentNewSolutions->[$k]->{$SingletReactions[$i]} = 1;
3672                                                            }
3673                                                    }
3674                                                    #Getting reaction data and printing reaction in output file
3675                                                    my $ReactionData;
3676                                                    if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3677                                                            $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3678                                                    } else {
3679                                                            my $Direction = substr($SingletReactions[$i],0,1);
3680                                                            if ($Direction eq "+") {
3681                                                                    $Direction = "=>";
3682                                                            } else {
3683                                                                    $Direction = "<=";
3684                                                            }
3685                                                            my $Reaction = substr($SingletReactions[$i],1);
3686                                                            $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3687                                                            $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3688                                                            $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3689                                                            if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3690                                                                    $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3691                                                            }
3692                                                            if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3693                                                                    $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3694                                                            }
3695                                                            if (!defined($ReactionData->{"DELTAG"}->[0])) {
3696                                                                    $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3697                                                            }
3698                                                            $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3699                                                    }
3700                                                    print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{"COUNT"};
3701                                                    for (my $k=0; $k < $SolutionCount; $k++) {
3702                                                            print RECONCILIATION ";";
3703                                                            if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3704                                                                    if ($j == 0) {
3705                                                                            print RECONCILIATION "|".$k."|";
3706                                                                    } else {
3707                                                                            print RECONCILIATION "|".$Index->{$k}."|";
3708                                                                    }
3709                                                            }
3710                                                    }
3711                                                    print RECONCILIATION "\n";
3712                                            }
3713                                            #Adding the current new solutions to the new solutions array
3714                                            if (defined($CurrentNewSolutions) && @{$CurrentNewSolutions} > 0) {
3715                                                    push(@{$NewSolutions},@{$CurrentNewSolutions});
3716                                            }
3717                                    }
3718                            }
3719                            #Adding the base reactions to all existing solutions
3720                            for (my $j=0; $j < @{$Solutions}; $j++) {
3721                                    if (defined($ReactionSetHash->{$ReactionSet}->{$ReactionSet}->{$Solutions->[$j]->{"BASE"}})) {
3722                                            foreach my $SingleReaction (@{$BaseReactions}) {
3723                                                    $Solutions->[$j]->{$SingleReaction} = 1;
3724                                            }
3725                                    }
3726                            }
3727                            #Adding the new solutions to the set of existing solutions
3728                            push(@{$Solutions},@{$NewSolutions});
3729                    }
3730                    close(RECONCILIATION);
3731                    #Now printing a file that defines all of the solutions in a format the testsolutions function understands
3732                    open (RECONCILIATION, ">$OutputFilenameTwo");
3733                    print RECONCILIATION "Experiment;Solution index;Solution cost;Solution reactions\n";
3734                    for (my $i=0; $i < @{$Solutions}; $i++) {
3735                            delete($Solutions->[$i]->{"BASE"});
3736                            print RECONCILIATION "SR".$i.";".$i.";10;".join(",",keys(%{$Solutions->[$i]}))."\n";
3737                    }
3738                    close(RECONCILIATION);
3739    
3740                    $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3741                    $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3742                    $Row->{"GF RECON TESTING TIMING"}->[0] = time()."-";
3743                    $Row->{"GF RECON SOLUTIONS"}->[0] = @{$Solutions};
3744                    $GrowMatchTable->save();
3745                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3746    
3747                    #Scheduling the solution testing
3748                    if ($GapFill == 1) {
3749                            system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GFSR:BACK:fast:QSUB\"");
3750                    } else {
3751                            system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GGSR:BACK:fast:QSUB\"");
3752                    }
3753            } else {
3754                    #Reading in the solution testing results
3755                    my $Data;
3756                    if ($GapFill == 1) {
3757                            $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GFSREM.txt","");
3758                    } else {
3759                            $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GGSREM.txt","");
3760                    }
3761    
3762                    #Reading in the preliminate reconciliation report
3763                    my $OutputData = $self->figmodel()->database()->load_single_column_file($OutputFilename,"");
3764                    #Replacing the file tags with actual performance data
3765                    my $Count = 0;
3766                    for (my $i=0; $i < @{$Data}; $i++) {
3767                            if ($Data->[$i] =~ m/^SR(\d+);.+;(\d+\/\d+);/) {
3768                                    my $Index = $1;
3769                                    my $Performance = $Index."/".$2;
3770                                    for (my $j=0; $j < @{$OutputData}; $j++) {
3771                                            $OutputData->[$j] =~ s/\|$Index\|/$Performance/g;
3772                                    }
3773                            }
3774                    }
3775                    $self->figmodel()->database()->print_array_to_file($OutputFilename,$OutputData);
3776    
3777                    my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3778                    my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3779                    $Row->{"GF RECON TESTING TIMING"}->[0] .= time();
3780                    $GrowMatchTable->save();
3781                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3782            }
3783    
3784            return 1;
3785    }
3786    
3787    =head3 DetermineCofactorLipidCellWallComponents
3788    Definition:
3789            {cofactor=>{string:compound id=>float:coefficient},lipid=>...cellWall=>} = FIGMODELmodel->DetermineCofactorLipidCellWallComponents();
3790    Description:
3791    =cut
3792    sub DetermineCofactorLipidCellWallComponents {
3793            my ($self) = @_;
3794            my $templateResults;
3795            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
3796            my $Class = $self->{_data}->cellwalltype();
3797            my $Name = $self->name();
3798            my $translation = {COFACTOR=>"cofactor",LIPIDS=>"lipid","CELL WALL"=>"cellWall"};
3799                  #Checking for phoenix variants                  #Checking for phoenix variants
3800                  my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");                  my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");
3801                  my $Phoenix = 0;                  my $Phoenix = 0;
3802                  my @Rows = $PhoenixVariantTable->get_rows_by_key($OrganismID,"GENOME");          my @Rows = $PhoenixVariantTable->get_rows_by_key($self->genome(),"GENOME");
3803                  my $VariantHash;                  my $VariantHash;
3804                  for (my $i=0; $i < @Rows; $i++) {                  for (my $i=0; $i < @Rows; $i++) {
3805                          $Phoenix = 1;                          $Phoenix = 1;
# Line 2698  Line 3807 
3807                                  $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];                                  $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];
3808                          }                          }
3809                  }                  }
   
3810                  #Collecting genome data                  #Collecting genome data
3811                  my $RoleHash;                  my $RoleHash;
3812                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
# Line 2715  Line 3823 
3823                                  }                                  }
3824                          }                          }
3825                  }                  }
   
3826                  #Scanning through the template item by item and determinine which biomass components should be added                  #Scanning through the template item by item and determinine which biomass components should be added
3827                  my $ComponentTypes;          my $includedHash;
3828                  my $EquationData;          my $BiomassReactionTemplateTable = $self->figmodel()->database()->get_table("BIOMASSTEMPLATE");
                 my $BiomassReactionTemplateTable = $self->figmodel()->database()->GetDBTable("BIOMASS TEMPLATE");  
3829                  for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {                  for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {
3830                          my $Row = $BiomassReactionTemplateTable->get_row($i);                          my $Row = $BiomassReactionTemplateTable->get_row($i);
3831                    if (defined($translation->{$Row->{CLASS}->[0]})) {
3832                            my $coef = -1;
3833                            if ($Row->{"REACTANT"}->[0] eq "NO") {
3834                                    $coef = 1;
3835                                    if ($Row->{"COEFFICIENT"}->[0] =~ m/cpd/) {
3836                                            $coef = $Row->{"COEFFICIENT"}->[0];
3837                                    }
3838                            }
3839                          if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {                          if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {
3840                                  push(@{$EquationData},$Row);                                  $includedHash->{$Row->{"ID"}->[0]} = 1;
3841                                  $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;                                  $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
3842                          } else {                          } elsif (defined($Row->{"INCLUSION CRITERIA"}->[0])) {
3843                                  my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];                                  my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];
3844                                  my $End = 0;                                  my $End = 0;
3845                                  while ($End == 0) {                                  while ($End == 0) {
3846                                          if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {                                          if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {
3847                                                  print $Criteria."\n";                                                  print $Criteria." : ";
3848                                                  my $Start = "";                                                  my $Start = "";
3849                                                  my $End = "";                                                  my $End = "";
3850                                                  my $Condition = $1;                                                  my $Condition = $1;
# Line 2754  Line 3868 
3868                                                                  $Result = "NO";                                                                  $Result = "NO";
3869                                                                  last;                                                                  last;
3870                                                          } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {                                                          } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {
3871                                                                  my $Match = 0;                                                                  if (defined($includedHash->{$1}) && $Condition eq "OR") {
                                                                 for (my $k=0; $k < @{$EquationData}; $k++) {  
                                                                         if ($EquationData->[$k]->{"ID"}->[0] eq $1) {  
                                                                                 $Match = 1;  
                                                                                 last;  
                                                                         }  
                                                                 }  
                                                                 if ($Match == 1 && $Condition eq "OR") {  
3872                                                                          $Result = "YES";                                                                          $Result = "YES";
3873                                                                          last;                                                                          last;
3874                                                                  } elsif ($Match != 1 && $Condition eq "AND") {                                                                  } elsif (!defined($includedHash->{$1}) && $Condition eq "AND") {
3875                                                                          $Result = "NO";                                                                          $Result = "NO";
3876                                                                          last;                                                                          last;
3877                                                                  }                                                                  }
3878                                                          } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {                                                          } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {
3879                                                                  my $Match = 0;                                                                  if (!defined($includedHash->{$1}) && $Condition eq "OR") {
                                                                 for (my $k=0; $k < @{$EquationData}; $k++) {  
                                                                         if ($EquationData->[$k]->{"ID"}->[0] eq $1) {  
                                                                                 $Match = 1;  
                                                                                 last;  
                                                                         }  
                                                                 }  
                                                                 if ($Match != 1 && $Condition eq "OR") {  
3880                                                                          $Result = "YES";                                                                          $Result = "YES";
3881                                                                          last;                                                                          last;
3882                                                                  } elsif ($Match == 1 && $Condition eq "AND") {                                                                  } elsif (defined($includedHash->{$1}) && $Condition eq "AND") {
3883                                                                          $Result = "NO";                                                                          $Result = "NO";
3884                                                                          last;                                                                          last;
3885                                                                  }                                                                  }
# Line 2895  Line 3995 
3995                                          }                                          }
3996                                  }                                  }
3997                                  if ($Criteria eq "YES") {                                  if ($Criteria eq "YES") {
3998                                          push(@{$EquationData},$Row);                                          $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
3999                                          $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;                                          $includedHash->{$Row->{"ID"}->[0]} = 1;
4000                                  }                                  }
4001                          }                          }
4002                  }                  }
4003            }
4004                  #Building biomass equation          my $types = ["cofactor","lipid","cellWall"];
4005                  my %Reactants;          my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4006                  my %Products;          for (my $i=0; $i < @{$types}; $i++) {
4007                  foreach my $EquationRow (@{$EquationData}) {                  my @list =      keys(%{$templateResults->{$types->[$i]}});
4008                          #First determine what the coefficient should be                  my $entries = 0;
4009                          my $Coefficient;                  for (my $j=0; $j < @list; $j++) {
4010                          if (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^[0-9\.]+$/) {                          if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4011                                  $Coefficient = $EquationRow->{"COEFFICIENT"}->[0];                                  my $objs = $cpdMgr->get_objects({id=>$list[$j]});
4012                          } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/cpd\d\d\d\d\d/) {                                  if (!defined($objs->[0]) || $objs->[0]->mass() == 0) {
4013                                  $Coefficient = 0;                                          $templateResults->{$types->[$i]}->{$list[$j]} = -1e-5;
                                 my @CompoundList = split(/,/,$EquationRow->{"COEFFICIENT"}->[0]);  
                                 foreach my $Compound (@CompoundList) {  
                                         if (defined($Reactants{$Compound})) {  
                                                 $Coefficient += $Reactants{$Compound};  
                                         }  
                                 }  
                         } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^([0-9\.]+)\/(.+)$/) {  
                                 my @Keys = keys(%{$ComponentTypes->{$2}});  
                                 my $MW = 1;  
                                 my $CompoundData = $self->figmodel()->LoadObject($EquationRow->{"ID"}->[0]);  
                                 if (defined($CompoundData->{"MASS"})) {  
                                         $MW = $CompoundData->{"MASS"}->[0];  
                                 }  
                                 $Coefficient = $1/@Keys/$MW;  
                         }  
                         if (defined($EquationRow->{"REACTANT"}) && $EquationRow->{"REACTANT"}->[0] eq "YES") {  
                                 if (defined($Reactants{$EquationRow->{"ID"}->[0]})) {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} += $Coefficient;  
                                 } elsif (defined($Products{$EquationRow->{"ID"}->[0]}) && $Products{$EquationRow->{"ID"}->[0]} > $Coefficient) {  
                                         $Products{$EquationRow->{"ID"}->[0]} -= $Coefficient;  
                                 } elsif (defined($Products{$EquationRow->{"ID"}->[0]}) && $Products{$EquationRow->{"ID"}->[0]} < $Coefficient) {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} = $Coefficient - $Products{$EquationRow->{"ID"}->[0]};  
                                         delete $Products{$EquationRow->{"ID"}->[0]};  
                                 } else {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} = $Coefficient;  
                                 }  
                         } else {  
                                 if (defined($Products{$EquationRow->{"ID"}->[0]})) {  
                                         $Products{$EquationRow->{"ID"}->[0]} += $Coefficient;  
                                 } elsif (defined($Reactants{$EquationRow->{"ID"}->[0]}) && $Reactants{$EquationRow->{"ID"}->[0]} > $Coefficient) {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} -= $Coefficient;  
                                 } elsif (defined($Reactants{$EquationRow->{"ID"}->[0]}) && $Reactants{$EquationRow->{"ID"}->[0]} < $Coefficient) {  
                                         $Products{$EquationRow->{"ID"}->[0]} = $Coefficient - $Reactants{$EquationRow->{"ID"}->[0]};  
                                         delete $Reactants{$EquationRow->{"ID"}->[0]};  
4014                                  } else {                                  } else {
4015                                          $Products{$EquationRow->{"ID"}->[0]} = $Coefficient;                                          $entries++;
4016                                    }
4017                            }
4018                    }
4019                    for (my $j=0; $j < @list; $j++) {
4020                            if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4021                                    $templateResults->{$types->[$i]}->{$list[$j]} = -1/$entries;
4022                            } elsif ($templateResults->{$types->[$i]}->{$list[$j]} =~ m/cpd/) {
4023                                    my $netCoef = 0;
4024                                    my @allcpd = split(/,/,$templateResults->{$types->[$i]}->{$list[$j]});
4025                                    for (my $k=0; $k < @allcpd; $k++) {
4026                                            if (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} ne "-1e-5") {
4027                                                    $netCoef += (1/$entries);
4028                                            } elsif (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} eq "-1e-5") {
4029                                                    $netCoef += 1e-5;
4030                                            }
4031                                    }
4032                                    $templateResults->{$types->[$i]}->{$list[$j]} = $netCoef;
4033                                  }                                  }
4034                          }                          }
4035                  }                  }
4036            return $templateResults;
4037    }
4038    
4039    =head3 BuildSpecificBiomassReaction
4040    Definition:
4041            FIGMODELmodel->BuildSpecificBiomassReaction();
4042    Description:
4043    =cut
4044    sub BuildSpecificBiomassReaction {
4045            my ($self) = @_;
4046            #Getting the database handle for biomass reactions
4047            my $bioMgr = $self->figmodel()->database()->get_object_manager("bof");
4048            #Checking if the current biomass reaction appears in more than on model, if not, this biomass reaction is conserved for this model
4049            my $biomassID = $self->biomassReaction();
4050            if ($biomassID =~ m/bio\d\d\d\d\d/) {
4051                    my $mdlMgr = $self->figmodel()->database()->get_object_manager("model");
4052                    my $mdlObs = $mdlMgr->get_objects({biomassReaction=>$biomassID});
4053                    if (defined($mdlObs->[1])) {
4054                            $biomassID = "NONE";
4055                    }
4056            }
4057            #If the biomass ID is "NONE", then we create a new biomass reaction for the model
4058            my $bioObj;
4059            my $originalPackages = "";
4060            my $originalEssReactions = "";
4061            if ($biomassID eq "NONE") {
4062                    #Getting the current largest ID
4063                    $biomassID = $self->figmodel()->database()->check_out_new_id("bof");
4064                    $bioObj = $bioMgr->create({
4065                            id=>$biomassID,owner=>$self->owner(),users=>$self->users(),name=>"Biomass",equation=>"NONE",protein=>"0",
4066                            energy=>"0",DNA=>"0",RNA=>"0",lipid=>"0",cellWall=>"0",cofactor=>"0",
4067                            modificationDate=>time(),creationDate=>time(),
4068                            cofactorPackage=>"NONE",lipidPackage=>"NONE",cellWallPackage=>"NONE",
4069                            DNACoef=>"NONE",RNACoef=>"NONE",proteinCoef=>"NONE",lipidCoef=>"NONE",
4070                            cellWallCoef=>"NONE",cofactorCoef=>"NONE",essentialRxn=>"NONE"});
4071                    if (!defined($bioObj)) {
4072                            die $self->error_message("BuildSpecificBiomassReaction():Could not create new biomass reaction ".$biomassID."!");
4073                    }
4074            } else {
4075                    #Getting the biomass DB handler from the database
4076                    my $objs = $bioMgr->get_objects({id=>$biomassID});
4077                    if (!defined($objs->[0])) {
4078                            die $self->error_message("BuildSpecificBiomassReaction():Could not find biomass reaction ".$biomassID." in database!");
4079                    }
4080                    $bioObj = $objs->[0];
4081                    $bioObj->owner($self->owner());
4082                    $bioObj->users($self->users());
4083                    if (defined($bioObj->essentialRxn())) {
4084                            $originalEssReactions = $bioObj->essentialRxn();
4085                            $originalPackages = $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage();
4086                    }
4087            }
4088            #Getting genome stats
4089            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
4090            my $Class = $self->{_data}->cellwalltype();
4091            #Setting global coefficients based on cell wall type
4092            my $biomassCompounds;
4093            my $compounds;
4094            if ($Class eq "Gram positive") {
4095                    $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.323,cpd00052=>-0.199,cpd00062=>-0.215};
4096                    $compounds->{protein} = {cpd00001=>1,cpd00023=>-0.0637,cpd00033=>-0.0999,cpd00035=>-0.0653,cpd00039=>-0.0790,cpd00041=>-0.0362,cpd00051=>-0.0472,cpd00053=>-0.0637,cpd00054=>-0.0529,cpd00060=>-0.0277,cpd00065=>-0.0133,cpd00066=>-0.0430,cpd00069=>-0.0271,cpd00084=>-0.0139,cpd00107=>-0.0848,cpd00119=>-0.0200,cpd00129=>-0.0393,cpd00132=>-0.0362,cpd00156=>-0.0751,cpd00161=>-0.0456,cpd00322=>-0.0660};
4097                    $bioObj->protein("0.5284");
4098                    $bioObj->DNA("0.026");
4099                    $bioObj->RNA("0.0655");
4100                    $bioObj->lipid("0.075");
4101                    $bioObj->cellWall("0.25");
4102                    $bioObj->cofactor("0.10");
4103            } else {
4104                    $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.322,cpd00052=>-0.2,cpd00062=>-0.216};
4105                    $compounds->{protein} = {cpd00001=>1,cpd00023=>-0.0492,cpd00033=>-0.1145,cpd00035=>-0.0961,cpd00039=>-0.0641,cpd00041=>-0.0451,cpd00051=>-0.0554,cpd00053=>-0.0492,cpd00054=>-0.0403,cpd00060=>-0.0287,cpd00065=>-0.0106,cpd00066=>-0.0347,cpd00069=>-0.0258,cpd00084=>-0.0171,cpd00107=>-0.0843,cpd00119=>-0.0178,cpd00129=>-0.0414,cpd00132=>-0.0451,cpd00156=>-0.0791,cpd00161=>-0.0474,cpd00322=>-0.0543};
4106                    $bioObj->protein("0.563");
4107                    $bioObj->DNA("0.031");
4108                    $bioObj->RNA("0.21");
4109                    $bioObj->lipid("0.093");
4110                    $bioObj->cellWall("0.177");
4111                    $bioObj->cofactor("0.039");
4112            }
4113            #Setting energy coefficient for all reactions
4114            $bioObj->energy("40");
4115            $compounds->{energy} = {cpd00002=>-1,cpd00001=>-1,cpd00008=>1,cpd00009=>1,cpd00067=>1};
4116            #Setting DNA coefficients based on GC content
4117            my $gc = $self->figmodel()->get_genome_gc_content($self->genome());
4118            $compounds->{DNA} = {cpd00012=>1,cpd00115=>0.5*(1-$gc),cpd00241=>0.5*$gc,cpd00356=>0.5*$gc,cpd00357=>0.5*(1-$gc)};
4119            #Setting Lipid,cell wall,and cofactor coefficients based on biomass template
4120            my $templateResults = $self->DetermineCofactorLipidCellWallComponents();
4121            $compounds->{cofactor} = $templateResults->{cofactor};
4122            $compounds->{lipid} = $templateResults->{lipid};
4123            $compounds->{cellWall} = $templateResults->{cellWall};
4124            #Getting package number for cofactor, lipid, and cell wall
4125            my $packages;
4126            my $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4127            my $packageTypes = ["Cofactor","Lipid","CellWall"];
4128            my $translation = {"Cofactor"=>"cofactor","Lipid"=>"lipid","CellWall"=>"cellWall"};
4129            for (my $i=0; $i < @{$packageTypes}; $i++) {
4130                    my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4131                    my $function = $translation->{$packageTypes->[$i]}."Package";
4132                    if (@cpdList == 0) {
4133                            $bioObj->$function("NONE");
4134                    } else {
4135                            my $cpdgrpObs = $cpdgrpMgr->get_objects({type=>$packageTypes->[$i]."Package"});
4136                            for (my $j=0; $j < @{$cpdgrpObs}; $j++) {
4137                                    $packages->{$packageTypes->[$i]}->{$cpdgrpObs->[$j]->grouping()}->{$cpdgrpObs->[$j]->COMPOUND()} = 1;
4138                            }
4139                            my @packageList = keys(%{$packages->{$packageTypes->[$i]}});
4140                            my $packageHash;
4141                            for (my $j=0; $j < @packageList; $j++) {
4142                                    $packageHash->{join("|",sort(keys(%{$packages->{$packageTypes->[$i]}->{$packageList[$j]}})))} = $packageList[$j];
4143                            }
4144                            if (defined($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))})) {
4145                                    $bioObj->$function($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))});
4146                            } else {
4147                                    my $newPackageID = $self->figmodel()->database()->check_out_new_id($packageTypes->[$i]."Pkg");
4148                                    $bioObj->$function($newPackageID);
4149                                    my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4150                                    for (my $j=0; $j < @cpdList; $j++) {
4151                                            $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4152                                            $cpdgrpMgr->create({COMPOUND=>$cpdList[$j],grouping=>$newPackageID,type=>$packageTypes->[$i]."Package"});
4153                                    }
4154                            }
4155                    }
4156            }
4157            #Filling in coefficient terms in database and calculating global reaction coefficients based on classification abundancies
4158            my $equationCompounds;
4159            my $types = ["RNA","DNA","protein","lipid","cellWall","cofactor","energy"];
4160            my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4161            for (my $i=0; $i < @{$types}; $i++) {
4162                    my $coefString = "";
4163                    my @compounds = sort(keys(%{$compounds->{$types->[$i]}}));
4164                    #Building coefficient strings and determining net mass for component types