[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.28, Wed Aug 25 20:39:16 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);
1230          $stats->{"Build date"}->[0] = time();          $self->update_model_stats();
         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;  
1231  }  }
1232    
1233  =head3 update_model_stats  =head3 update_model_stats
# Line 974  Line 1241 
1241          #Getting reaction table          #Getting reaction table
1242          my $rxntbl = $self->reaction_table();          my $rxntbl = $self->reaction_table();
1243          if (!defined($rxntbl)) {          if (!defined($rxntbl)) {
1244                  $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;  
1245          }          }
1246          my $cpdtbl = $self->compound_table();          my $cpdtbl = $self->compound_table();
1247    
1248          #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];  
         }  
   
1249          my %GeneHash;          my %GeneHash;
1250          my %NonpegHash;          my %NonpegHash;
1251          my %CompoundHash;          my %CompoundHash;
1252            my $spontaneousReactions = 0;
1253            my $gapFillReactions = 0;
1254            my $biologReactions = 0;
1255            my $transporters = 0;
1256            my $autoCompleteReactions = 0;
1257            my $associatedSubsystemGenes = 0;
1258          for (my $i=0; $i < $rxntbl->size(); $i++) {          for (my $i=0; $i < $rxntbl->size(); $i++) {
1259                  my $Row = $rxntbl->get_row($i);                  my $Row = $rxntbl->get_row($i);
1260                  if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {                  if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) {
# Line 1037  Line 1262 
1262                          if (defined($ReactionRow->{"EQUATION"}->[0])) {                          if (defined($ReactionRow->{"EQUATION"}->[0])) {
1263                                  #Checking for extracellular metabolites which indicate that this is a transporter                                  #Checking for extracellular metabolites which indicate that this is a transporter
1264                                  if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {                                  if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) {
1265                                          $CurrentStats->{"Transport reaction"}->[0]++;                                          $transporters++;
1266                                  }                                  }
1267                          }                          }
1268                          #Identifying spontaneous/biolog/gapfilling/gene associated reactions                          #Identifying spontaneous/biolog/gapfilling/gene associated reactions
1269                          if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {                          if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) {
1270                                  $CurrentStats->{"Biolog gap filling reactions"}->[0]++;                                  $biologReactions++;
1271                            } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GROW/i) {
1272                                    $gapFillReactions++;
1273                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {                          } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) {
1274                                  $CurrentStats->{"Spontaneous"}->[0]++;                                  $spontaneousReactions++;
1275                          } 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) {
1276                                  $CurrentStats->{"Gap filling reactions"}->[0]++;                                  $autoCompleteReactions++;
1277                          } else {                          } else {
1278                                  foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {                                  foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) {
1279                                          my @GeneList = split(/\+/,$GeneSet);                                          $_ = $GeneSet;
1280                                            my @GeneList = /(peg\.\d+)/g;
1281                                          foreach my $Gene (@GeneList) {                                          foreach my $Gene (@GeneList) {
1282                                                  if ($Gene =~ m/(peg\.\d+)/) {                                                  if ($Gene =~ m/(peg\.\d+)/) {
1283                                                          $GeneHash{$1} = 1;                                                          $GeneHash{$1} = 1;
# Line 1063  Line 1291 
1291          }          }
1292          my @genes = keys(%GeneHash);          my @genes = keys(%GeneHash);
1293          my @othergenes = keys(%NonpegHash);          my @othergenes = keys(%NonpegHash);
         $CurrentStats->{"Genes with reactions"}->[0] = @genes + @othergenes;  
1294    
1295          #Updating the stats stored in the table          #Setting the reaction count
1296          $self->figmodel()->database()->update_row("MODEL STATS",$CurrentStats,"Model ID");          $self->{_data}->reactions($rxntbl->size());
1297          $self->{_stats} = $CurrentStats;          #Setting the metabolite count
1298          return $CurrentStats;          $self->{_data}->compounds($cpdtbl->size());
1299            #Setting the gene count
1300            my $geneCount = @genes + @othergenes;
1301            $self->{_data}->associatedGenes($geneCount);
1302            #Setting remaining stats
1303            $self->{_data}->spontaneousReactions($spontaneousReactions);
1304            $self->{_data}->gapFillReactions($gapFillReactions);
1305            $self->{_data}->biologReactions($biologReactions);
1306            $self->{_data}->transporters($transporters);
1307            $self->{_data}->autoCompleteReactions($autoCompleteReactions);
1308            $self->{_data}->associatedSubsystemGenes($associatedSubsystemGenes);
1309            #Setting the model class
1310            my $class = "";
1311            for (my $i=0; $i < @{$self->figmodel()->config("class list")}; $i++) {
1312                    if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i]))) {
1313                            if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->id()})) {
1314                                    $class = $self->figmodel()->config("class list")->[$i];
1315                                    last;
1316                            }
1317                            if ($class eq "" && defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->genome()})) {
1318                                    $class = $self->figmodel()->config("class list")->[$i];
1319                            }
1320                    }
1321            }
1322            if ($class eq "") {
1323                    $class = $self->figmodel()->get_genome_stats($self->genome())->{CLASS}->[0];
1324            }
1325            if ($class eq "") {
1326                    $class = "unknown";
1327            }
1328            $self->{_data}->cellwalltype($class);
1329  }  }
1330    
1331  =head3 GapFillModel  =head3 GapFillModel
# Line 1088  Line 1345 
1345          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
1346          my $StartTime = time();          my $StartTime = time();
1347    
1348          #Archiving the existing model          #Reading original reaction table
1349          $self->ArchiveModel();          my $OriginalRxn = $self->reaction_table();
1350            #Clearing the table
1351            $self->reaction_table(1);
1352    
1353          #Removing any gapfilling reactions that may be currently present in the model          #Removing any gapfilling reactions that may be currently present in the model
1354          if (!defined($donotclear) || $donotclear != 1) {          if (!defined($donotclear) || $donotclear != 1) {
1355                  my $ModelTable = $self->reaction_table();                  my $ModelTable = $self->reaction_table();
1356                  for (my $i=0; $i < $ModelTable->size(); $i++) {                  for (my $i=0; $i < $ModelTable->size(); $i++) {
1357                          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/) {
1358                                  $ModelTable->delete_row($i);                                  $ModelTable->delete_row($i);
1359                                  $i--;                                  $i--;
1360                          }                          }
# Line 1104  Line 1363 
1363          }          }
1364    
1365          #Calling the MFAToolkit to run the gap filling optimization          #Calling the MFAToolkit to run the gap filling optimization
1366          my $MinimalMediaTable = $self->figmodel()->database()->GetDBTable("MINIMAL MEDIA TABLE");          my $Media = $self->autocompleteMedia();
1367          my $Media = "Complete";          if ($Media eq "Complete") {
1368          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));
1369                  $Media = $MinimalMediaTable->get_row_by_key($self->genome(),"Organism")->{"Minimal media"}->[0];          } else {
1370                  #Loading media, changing bounds, saving media as a test media                  #Loading media, changing bounds, saving media as a test media
1371                  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"]);
1372                  for (my $i=0; $i < $MediaTable->size(); $i++) {                  for (my $i=0; $i < $MediaTable->size(); $i++) {
# Line 1119  Line 1378 
1378                          }                          }
1379                  }                  }
1380                  $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");                  $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
1381                  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";
1382                    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));
1383                  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));  
1384          }          }
1385    
1386          #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
1387          my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);          my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename);
1388    
1389          #Clearing the mfatoolkit output and log file          #Clearing the mfatoolkit output and log file
1390          $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");          $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log");
1391          if (!defined($SolutionData) || $SolutionData->size() == 0) {          if (!defined($SolutionData) || $SolutionData->size() == 0) {
# Line 1134  Line 1393 
1393                  $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");                  $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");
1394                  return $self->fail();                  return $self->fail();
1395          }          }
         $SolutionData->save("/home/chenry/Solution.txt");  
1396    
1397          #Looking for the last printed recursive MILP solution          #Looking for the last printed recursive MILP solution
1398          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
1399                  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/) {
1400                          my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);                          my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15);
                         print "Solution:".$AllSolutions."\n";  
1401                          my @TempThree = split(/\|/,$AllSolutions);                          my @TempThree = split(/\|/,$AllSolutions);
1402                          if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {                          if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) {
1403                                  my @TempFour = split(/,/,$1);                                  my @TempFour = split(/,/,$1);
# Line 1163  Line 1420 
1420                                                  push(@{$ReactionList},$ID);                                                  push(@{$ReactionList},$ID);
1421                                          }                                          }
1422                                  }                                  }
1423                                  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);  
1424                          }                          }
1425                          last;                          last;
1426                  }                  }
1427          }          }
1428    
1429          #Updating model stats with gap filling results          #Updating model stats with gap filling results
1430          my $ElapsedTime = time() - $StartTime;          my $ElapsedTime = time() - $StartTime;
1431          $self->figmodel()->ClearDBModel($self->id(),1);          $self->reaction_table(1);
1432            $self->calculate_model_changes($OriginalRxn,"AUTOCOMPLETION");
1433    
1434          #Determining why each gap filling reaction was added          #Determining why each gap filling reaction was added
1435          $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);          $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media);
1436          if (!defined($donotclear) || $donotclear != 1) {          if (!defined($donotclear) || $donotclear != 1) {
                 $self->figmodel()->AddBiologTransporters($self->id());  
1437                  if ($self->id() !~ m/MGRast/) {                  if ($self->id() !~ m/MGRast/) {
1438                          $self->update_stats_for_gap_filling($ElapsedTime);                          $self->update_stats_for_gap_filling($ElapsedTime);
1439                  }                  }
# Line 1191  Line 1449 
1449          return $self->success();          return $self->success();
1450  }  }
1451    
1452    =head3 calculate_model_changes
1453    Definition:
1454            FIGMODELmodel->calculate_model_changes(FIGMODELTable:original reaction table,string:modification cause);
1455    Description:
1456    
1457    =cut
1458    
1459    sub calculate_model_changes {
1460            my ($self,$originalReactions,$cause,$tbl,$version) = @_;
1461            my $modTime = time();
1462            if (!defined($version)) {
1463                    $version = $self->selected_version();
1464            }
1465            my $user = $self->figmodel()->user();
1466            #Getting the history table
1467            my $histTbl = $self->model_history();
1468            #Checking for differences
1469            if (!defined($tbl)) {
1470                    $tbl = $self->reaction_table();
1471            }
1472            for (my $i=0; $i < $tbl->size(); $i++) {
1473                    my $row = $tbl->get_row($i);
1474                    my $orgRow = $originalReactions->get_row_by_key($row->{LOAD}->[0],"LOAD");
1475                    if (!defined($orgRow)) {
1476                            $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => $row->{DIRECTIONALITY}, GeneChange => $row->{"ASSOCIATED PEG"}, Action => ["ADDED"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1477                    } else {
1478                            my $geneChanges;
1479                            my $directionChange;
1480                            if ($orgRow->{"DIRECTIONALITY"}->[0] ne $row->{"DIRECTIONALITY"}->[0]) {
1481                                    $directionChange = $orgRow->{"DIRECTIONALITY"}->[0]." to ".$row->{"DIRECTIONALITY"}->[0];
1482                            }
1483                            for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
1484                                    my $match = 0;
1485                                    if (defined($orgRow->{"ASSOCIATED PEG"})) {
1486                                            for (my $k=0; $k < @{$orgRow->{"ASSOCIATED PEG"}}; $k++) {
1487                                                    if ($row->{"ASSOCIATED PEG"}->[$j] eq $orgRow->{"ASSOCIATED PEG"}->[$k]) {
1488                                                            $match = 1;
1489                                                    }
1490                                            }
1491                                    }
1492                                    if ($match == 0) {
1493                                            push(@{$geneChanges},"Added ".$row->{"ASSOCIATED PEG"}->[$j]);
1494                                    }
1495                            }
1496                            if (defined($orgRow->{"ASSOCIATED PEG"})) {
1497                                    for (my $k=0; $k < @{$orgRow->{"ASSOCIATED PEG"}}; $k++) {
1498                                            my $match = 0;
1499                                            if (defined($row->{"ASSOCIATED PEG"})) {
1500                                                    for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
1501                                                            if ($row->{"ASSOCIATED PEG"}->[$j] eq $orgRow->{"ASSOCIATED PEG"}->[$k]) {
1502                                                                    $match = 1;
1503                                                            }
1504                                                    }
1505                                            }
1506                                            if ($match == 0) {
1507                                                    push(@{$geneChanges},"Removed ".$orgRow->{"ASSOCIATED PEG"}->[$k]);
1508                                            }
1509                                    }
1510                            }
1511                            if ((defined($directionChange) && length($directionChange) > 0) || defined($geneChanges) && @{$geneChanges} > 0) {
1512                                    $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => [$directionChange], GeneChange => $geneChanges, Action => ["CHANGE"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1513                            }
1514                    }
1515            }
1516            #Looking for removed reactions
1517            for (my $i=0; $i < $originalReactions->size(); $i++) {
1518                    my $row = $originalReactions->get_row($i);
1519                    my $orgRow = $tbl->get_row_by_key($row->{LOAD}->[0],"LOAD");
1520                    if (!defined($orgRow)) {
1521                            $histTbl->add_row({Reaction => [$row->{LOAD}->[0]], DirectionChange => $row->{DIRECTIONALITY}, GeneChange => $row->{"ASSOCIATED PEG"}, Action => ["REMOVED"], ModificationTime => [$modTime], ModifcationCause => [$cause], User => [$user], Version => [$version]});
1522                    }
1523            }
1524            $histTbl->save();
1525    }
1526    
1527    =head3 GapGenModel
1528    Definition:
1529            FIGMODELmodel->GapGenModel();
1530    Description:
1531            Runs the gap generation algorithm to correct a single false positive prediction. Results are loaded into a table.
1532    =cut
1533    
1534    sub GapGenModel {
1535            my ($self,$Media,$KOList,$NoKOList,$Experiment,$SolutionLimit) = @_;
1536    
1537            #Enforcing nonoptional arguments
1538            if (!defined($Media)) {
1539                    return undef;
1540            }
1541            if (!defined($KOList)) {
1542                    $KOList->[0] = "none";
1543            }
1544            if (!defined($NoKOList)) {
1545                    $NoKOList->[0] = "none";
1546            }
1547            if (!defined($Experiment)) {
1548                    $Experiment= "ReactionKO";
1549            }
1550            if (!defined($SolutionLimit)) {
1551                    $SolutionLimit = "10";
1552            }
1553    
1554            #Translating the KO lists into arrays
1555            if (ref($KOList) ne "ARRAY") {
1556                    my $temp = $KOList;
1557                    $KOList = ();
1558                    push(@{$KOList},split(/[,;]/,$temp));
1559            }
1560            my $noKOHash;
1561            if (defined($NoKOList) && ref($NoKOList) ne "ARRAY") {
1562                    my $temp = $NoKOList;
1563                    $NoKOList = ();
1564                    push(@{$NoKOList},split(/[,;]/,$temp));
1565                    foreach my $rxn (@{$NoKOList}) {
1566                            $noKOHash->{$rxn} = 1;
1567                    }
1568            }
1569    
1570            #Checking if solutions exist for the input parameters
1571            $self->aquireModelLock();
1572            my $tbl = $self->load_model_table("GapGenSolutions");
1573            my $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1574            my $solutions;
1575            if (defined($solutionRow)) {
1576                    #Checking if any solutions conform to the no KO list
1577                    foreach my $solution (@{$solutionRow->{Solutions}}) {
1578                            my @reactions = split(/,/,$solution);
1579                            my $include = 1;
1580                            foreach my $rxn (@reactions) {
1581                                    if ($rxn =~ m/(rxn\d\d\d\d\d)/) {
1582                                            if (defined($noKOHash->{$1})) {
1583                                                    $include = 0;
1584                                            }
1585                                    }
1586                            }
1587                            if ($include == 1) {
1588                                    push(@{$solutions},$solution);
1589                            }
1590                    }
1591            } else {
1592                    $solutionRow = {Media => [$Media],Experiment => [$Experiment],KOlist => [join(",",@{$KOList})]};
1593                    $tbl->add_row($solutionRow);
1594                    $self->figmodel()->database()->save_table($tbl);
1595            }
1596            $self->releaseModelLock();
1597    
1598            #Returning solution list of solutions were found
1599            if (defined($solutions) && @{$solutions} > 0) {
1600                    return $solutions;
1601            }
1602    
1603            #Getting unique filename
1604            my $Filename = $self->figmodel()->filename();
1605    
1606            #Running the gap generation
1607            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));
1608            my $ProblemReport = $self->figmodel()->LoadProblemReport($Filename);
1609            if (!defined($ProblemReport)) {
1610                    $self->figmodel()->error_message("FIGMODEL:GapGenerationAlgorithm;No problem report;".$Filename.";".$self->id().$self->selected_version().";".$Media.";".$KOList.";".$NoKOList);
1611                    return undef;
1612            }
1613    
1614            #Clearing the output folder and log file
1615            $self->figmodel()->clearing_output($Filename,"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log");
1616    
1617            #Saving the solution
1618            $self->aquireModelLock();
1619            $tbl = $self->load_model_table("GapGenSolutions");
1620            $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist");
1621            for (my $j=0; $j < $ProblemReport->size(); $j++) {
1622                    if ($ProblemReport->get_row($j)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^)]+)/) {
1623                            my @SolutionList = split(/\|/,$1);
1624                            for (my $k=0; $k < @SolutionList; $k++) {
1625                                    if ($SolutionList[$k] =~ m/(\d+):(.+)/) {
1626                                            push(@{$solutionRow->{Solutions}},$2);
1627                                            push(@{$solutions},$2);
1628                                    }
1629                            }
1630                    }
1631            }
1632            $self->figmodel()->database()->save_table($tbl);
1633            $self->releaseModelLock();
1634    
1635            return $solutions;
1636    }
1637    
1638  =head3 datagapfill  =head3 datagapfill
1639  Definition:  Definition:
1640          success()/fail() = FIGMODELmodel->datagapfill();          success()/fail() = FIGMODELmodel->datagapfill();
# Line 1201  Line 1645 
1645          my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;          my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_;
1646          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
1647          if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {          if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) {
                 print "Gap filling specs:\n".join("\n",@{$GapFillingRunSpecs})."\n\n";  
1648                  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));
   
1649                  #Checking that the solution exists                  #Checking that the solution exists
1650                  if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {                  if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") {
1651                          $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 1655 
1655                  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);
1656                  if (defined($TansferFileSuffix)) {                  if (defined($TansferFileSuffix)) {
1657                          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");  
1658                  }                  }
   
1659                  #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
1660                  $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");
1661                  return $GapFillingRunSpecs;                  return $GapFillResultTable;
1662          }          }
1663          if (defined($TansferFileSuffix)) {          if (defined($TansferFileSuffix)) {
1664                  $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 1724 
1724                          push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});                          push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]});
1725                  }                  }
1726                  print "Integrating solution!\n";                  print "Integrating solution!\n";
1727                  $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);
1728                  my $model = $self->get_model($self->id().$TempVersion);                  $self->PrintModelLPFile();
                 $model->PrintModelLPFile();  
1729                  #Running the model against all available experimental data                  #Running the model against all available experimental data
1730                  print "Running test model!\n";                  print "Running test model!\n";
1731                  my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");                  my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All");
# Line 1322  Line 1761 
1761  =cut  =cut
1762  sub status {  sub status {
1763          my ($self) = @_;          my ($self) = @_;
1764          return $self->{_data}->{status}->[0];          return $self->{_data}->status();
1765  }  }
1766    
1767  =head3 message  =head3 message
# Line 1333  Line 1772 
1772  =cut  =cut
1773  sub message {  sub message {
1774          my ($self) = @_;          my ($self) = @_;
1775          return $self->{_data}->{message}->[0];          return $self->{_data}->message();
1776  }  }
1777    
1778  =head3 set_status  =head3 set_status
# Line 1348  Line 1787 
1787  =cut  =cut
1788  sub set_status {  sub set_status {
1789          my ($self,$NewStatus,$Message) = @_;          my ($self,$NewStatus,$Message) = @_;
1790            $self->{_data}->status($NewStatus);
1791          #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");  
1792          return $self->config("SUCCESS")->[0];          return $self->config("SUCCESS")->[0];
1793  }  }
1794    
# Line 1367  Line 1803 
1803  sub CreateSingleGenomeReactionList {  sub CreateSingleGenomeReactionList {
1804          my ($self,$RunGapFilling) = @_;          my ($self,$RunGapFilling) = @_;
1805    
1806            #Creating directory
1807            if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/") {
1808                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/");
1809            } elsif ($self->owner() eq "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->genome()."/") {
1810                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->genome()."/");
1811            }
1812            if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/") {
1813                    system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/");
1814            }
1815    
1816          #Getting genome stats          #Getting genome stats
1817          my $genomestats = $self->figmodel()->get_genome_stats($self->genome());          my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
1818          my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());          my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
# Line 1381  Line 1827 
1827          }          }
1828          #Setting up needed variables          #Setting up needed variables
1829          my $OriginalModelTable = undef;          my $OriginalModelTable = undef;
1830          #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) {  
1831                  $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");  
1832                  return $self->fail();                  return $self->fail();
1833          }elsif ($Row->{status}->[0] == 1) {          }elsif ($self->status() == 1) {
1834                  $OriginalModelTable = $self->reaction_table();                  $OriginalModelTable = $self->reaction_table();
1835                  $self->ArchiveModel();                  $self->set_status(0,"Rebuilding preliminary reconstruction");
                 $Row->{status}->[0] = 0;  
                 $Row->{message}->[0] = "Rebuilding preliminary reconstruction";  
1836          } else {          } else {
1837                  $Row->{status}->[0] = 0;                  $self->set_status(0,"Preliminary reconstruction");
                 $Row->{message}->[0] = "Preliminary reconstruction";  
1838          }          }
         #Updating the status table  
         $self->figmodel()->database()->save_table($ModelTable);  
         $self->figmodel()->database()->UnlockDBTable("MODELS");  
1839          #Sorting GenomeData by gene location on the chromosome          #Sorting GenomeData by gene location on the chromosome
1840            my $ftrTbl = $self->figmodel()->database()->get_table("ROLERXNMAPPING");
1841          $FeatureTable->sort_rows("MIN LOCATION");          $FeatureTable->sort_rows("MIN LOCATION");
1842          my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);          my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash);
1843          my %LastGenePosition;          my %LastGenePosition;
# Line 1476  Line 1909 
1909                          }                          }
1910                  }                  }
1911          }          }
   
1912          #Creating nonadjacent complexes          #Creating nonadjacent complexes
1913          my @ReactionList = keys(%{$ReactionHash});          my @ReactionList = keys(%{$ReactionHash});
1914          foreach my $Reaction (@ReactionList) {          foreach my $Reaction (@ReactionList) {
# Line 1582  Line 2014 
2014                  $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"]});
2015          }          }
2016    
2017            #Getting feature rows for features that are lumped
2018            my @rows = $ftrTbl->get_rows_by_key("2","MASTER");
2019            for (my $i=0; $i < @rows; $i++) {
2020                    my $rxn = $rows[$i]->{REACTION}->[0];
2021                    my $role = $rows[$i]->{ROLE}->[0];
2022                    my @orgrows = $FeatureTable->get_rows_by_key($role,"ROLES");
2023                    my $genes;
2024                    for (my $j=0; $j < @orgrows; $j++) {
2025                            if ($orgrows[$j]->{ID}->[0] =~ m/(peg\.\d+)/) {
2026                                    push(@{$genes},$1);
2027                            }
2028                    }
2029                    if (defined($genes) && @{$genes} > 0) {
2030                            my $row = $NewModelTable->get_row_by_key($rxn,"LOAD");
2031                            my $newGeneAssociations;
2032                            for (my $k=0; $k < @{$genes}; $k++) {
2033                                    for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) {
2034                                            my $peg = $genes->[$k];
2035                                            if ($row->{"ASSOCIATED PEG"}->[$j] !~ m/$peg/) {
2036                                                    push(@{$newGeneAssociations},$row->{"ASSOCIATED PEG"}->[$j]."+".$peg);
2037                                            }
2038                                    }
2039                            }
2040                            $row->{"ASSOCIATED PEG"} = $newGeneAssociations;
2041                    }
2042            }
2043    
2044          #Adding the spontaneous and universal reactions          #Adding the spontaneous and universal reactions
2045          foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {          foreach my $ReactionID (@{$self->config("spontaneous reactions")}) {
2046                  #Getting the thermodynamic reversibility from the database                  #Getting the thermodynamic reversibility from the database
# Line 1600  Line 2059 
2059                  }                  }
2060          }          }
2061    
2062          #Checking if a biomass reaction already exists          #Creating biomass reaction for model
2063          my $BiomassReactionRow = $self->figmodel()->database()->get_row_by_key("BIOMASS TABLE",$self->id(),"MODELS");          my $biomassID = $self->BuildSpecificBiomassReaction();
2064          if (!defined($BiomassReactionRow)) {          if ($biomassID !~ m/bio\d\d\d\d\d/) {
                 $BiomassReactionRow = $self->BuildSpecificBiomassReaction();  
                 if (!defined($BiomassReactionRow)) {  
2065                          $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");                          $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
2066                          $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");                  die $self->error_message("CreateSingleGenomeReactionList: Could not generate biomass reaction!");
2067                          return $self->fail();          }
2068            #Getting the biomass reaction PPO object
2069            my $bioRxn = $self->figmodel()->database()->get_object("bof",{id=>$biomassID});
2070            if (!defined($bioRxn)) {
2071                    die $self->error_message("CreateSingleGenomeReactionList: Could not find biomass reaction ".$biomassID."!");
2072            }
2073            #Getting the list of essential reactions for biomass reaction
2074            my $ReactionList;
2075            my $essentialReactions = $bioRxn->essentialRxn();
2076            if (defined($essentialReactions) && $essentialReactions =~ m/rxn\d\d\d\d\d/) {
2077                    push(@{$ReactionList},split(/\|/,$essentialReactions));
2078                    if ($essentialReactions !~ m/$biomassID/) {
2079                            push(@{$ReactionList},$biomassID);
2080                  }                  }
2081          }          }
         my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};  
         push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);  
2082    
2083          #Adding biomass reactions to the model table          #Adding biomass reactions to the model table
2084          foreach my $BOFReaction (@{$ReactionList}) {          foreach my $BOFReaction (@{$ReactionList}) {
# Line 1641  Line 2108 
2108                  }                  }
2109          }          }
2110    
2111            #If an original model exists, we copy the gap filling solution from that model
2112            if (defined($OriginalModelTable)) {
2113                    for (my $i=0; $i < $OriginalModelTable->size(); $i++) {
2114                            if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ && $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] ne "INITIAL GAP FILLING") {
2115                                    my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD");
2116                                    if (!defined($Row)) {
2117                                            $NewModelTable->add_row($OriginalModelTable->get_row($i));
2118                                    }
2119                            }
2120                    }
2121            }
2122    
2123          #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
2124          if (defined($OriginalModelTable)) {          if (defined($OriginalModelTable)) {
2125                  my $PerfectMatch = 1;                  my $PerfectMatch = 1;
# Line 1700  Line 2179 
2179          #Saving the new model to file          #Saving the new model to file
2180          $self->set_status(1,"Preliminary reconstruction complete");          $self->set_status(1,"Preliminary reconstruction complete");
2181          $self->figmodel()->database()->save_table($NewModelTable);          $self->figmodel()->database()->save_table($NewModelTable);
2182          $self->{_reaction_data} = $NewModelTable;          $self->reaction_table(1);
         #Clearing the previous model from the cache  
         $self->figmodel()->ClearDBModel($self->id(),1);  
2183          #Updating the model stats table          #Updating the model stats table
2184          $self->update_stats_for_build();          $self->update_stats_for_build();
2185          $self->PrintSBMLFile();          $self->PrintSBMLFile();
2186            $self->PrintModelLPFile();
2187            $self->PrintModelSimpleReactionTable();
2188            if (defined($OriginalModelTable)) {
2189                    $self->calculate_model_changes($OriginalModelTable,"REBUILD");
2190            }
2191    
2192          #Adding model to gapfilling queue          #Adding model to gapfilling queue
2193          if (defined($RunGapFilling) && $RunGapFilling == 1) {          if (defined($RunGapFilling) && $RunGapFilling == 1) {
# Line 1724  Line 2206 
2206    
2207  sub CreateMetaGenomeReactionList {  sub CreateMetaGenomeReactionList {
2208          my ($self) = @_;          my ($self) = @_;
   
2209          #Checking if the metagenome file exists          #Checking if the metagenome file exists
2210          if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {          if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") {
2211                  $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());
2212                    return $self->fail();
2213          }          }
2214          #Loading metagenome data          #Loading metagenome data
2215          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");
2216          if (!defined($MGRASTData)) {          if (!defined($MGRASTData)) {
2217                  $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());
2218                    return $self->fail();
2219          }          }
   
2220          #Setting up needed variables          #Setting up needed variables
2221          my $OriginalModelTable = undef;          my $OriginalModelTable = undef;
   
2222          #Checking status          #Checking status
2223          if ($self->status() < 0) {          if ($self->status() < 0) {
2224                  $self->set_status(0,"Preliminary reconstruction");                  $self->set_status(0,"Preliminary reconstruction");
# Line 1749  Line 2230 
2230                  $self->ArchiveModel();                  $self->ArchiveModel();
2231                  $self->set_status(0,"Rebuilding preliminary reconstruction");                  $self->set_status(0,"Rebuilding preliminary reconstruction");
2232          }          }
2233            #Creating a hash of escores and pegs associated with each role
2234            my $rolePegHash;
2235            my $roleEscores;
2236            for (my $i=0; $i < @{$MGRASTData};$i++) {
2237                    #MD5,PEG,number of sims,role,sim e-scores,max escore,min escore,ave escore,stdev escore,ave exponent,stddev exponent
2238                    $rolePegHash->{$MGRASTData->[$i]->[3]}->{substr($MGRASTData->[$i]->[1],4)} = 1;
2239                    push(@{$roleEscores->{$MGRASTData->[$i]->[3]}},split(/;/,$MGRASTData->[$i]->[4]));
2240            }
2241          #Getting the reaction table          #Getting the reaction table
2242          my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");          my $ReactionTable = $self->figmodel()->database()->get_table("REACTIONS");
2243          #Creating model table          #Creating model table
2244          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");
2245          for (my $i=0; $i < @{$MGRASTData};$i++) {          print $ModelTable->filename();
2246                  #MD5,PEG,number of sims,role,sim e-scores          my @roles = keys(%{$rolePegHash});
2247                  my $Role = $MGRASTData->[$i]->[3];          for (my $i=0; $i < @roles; $i++) {
2248                  my $MD5 = $MGRASTData->[$i]->[0];                  my $min = -1;
2249                  my $peg = $MGRASTData->[$i]->[1];                  my $max = -1;
2250                  my $sims = $MGRASTData->[$i]->[4];                  my $count = @{$roleEscores->{$roles[$i]}};
2251                  $sims =~ s/;/,/g;                  my $ave = 0;
2252                    my $stdev = 0;
2253                    my $aveexp = 0;
2254                    my $stdevexp = 0;
2255                    for (my $j=0; $j < @{$roleEscores->{$roles[$i]}}; $j++) {
2256                            if ($roleEscores->{$roles[$i]} < $min || $min == -1) {
2257                                    $min = $roleEscores->{$roles[$i]};
2258                            }
2259                            if ($roleEscores->{$roles[$i]} > $max || $max == -1) {
2260                                    $max = $roleEscores->{$roles[$i]};
2261                            }
2262                            $ave += $roleEscores->{$roles[$i]}->[$j];
2263                            if ($roleEscores->{$roles[$i]}->[$j] =~ m/e(-\d+$)/) {
2264                                    $aveexp += $1;
2265                            }
2266                    }
2267                    $ave = $ave/$count;
2268                    $aveexp = $aveexp/$count;
2269                    for (my $j=0; $j < @{$roleEscores->{$roles[$i]}}; $j++) {
2270                            $stdev += ($roleEscores->{$roles[$i]}->[$j]-$ave)*($roleEscores->{$roles[$i]}->[$j]-$ave);
2271                            if ($roleEscores->{$roles[$i]}->[$j] =~ m/e(-\d+$)/) {
2272                                    $stdevexp += ($1-$aveexp)*($1-$aveexp);
2273                            }
2274                    }
2275                    $stdev = sqrt($stdev/$count);
2276                    $stdevexp = sqrt($stdevexp/$count);
2277                  #Checking for subsystems                  #Checking for subsystems
2278                  my $GeneSubsystems = $self->figmodel()->subsystems_of_role($Role);                  my $GeneSubsystems = $self->figmodel()->subsystems_of_role($roles[$i]);
2279                  #Checking if there are reactions associated with this role                  #Checking if there are reactions associated with this role
2280                  my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role);                  my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($roles[$i]);
2281                  if ($ReactionHashArrayRef != 0) {                  if ($ReactionHashArrayRef != 0) {
2282                          foreach my $Mapping (@{$ReactionHashArrayRef}) {                          foreach my $Mapping (@{$ReactionHashArrayRef}) {
2283                                  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 2289 
2289                                                                  $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"]};
2290                                                                  $ModelTable->add_row($ReactionRow);                                                                  $ModelTable->add_row($ReactionRow);
2291                                                          }                                                          }
2292                                                          push(@{$ReactionRow->{"ASSOCIATED PEG"}},substr($peg,4));                                                          my %pegHash = %{$rolePegHash->{$roles[$i]}};
2293                                                          push(@{$ReactionRow->{"REFERENCE"}},$MD5.":".$Role);                                                          if (defined($ReactionRow->{"ASSOCIATED PEG"})) {
2294                                                          push(@{$ReactionRow->{"CONFIDENCE"}},$sims);                                                                  for (my $j=0; $j < @{$ReactionRow->{"ASSOCIATED PEG"}}; $j++) {
2295                                                                            $pegHash{$ReactionRow->{"ASSOCIATED PEG"}->[$j]} = 1;
2296                                                                    }
2297                                                            }
2298                                                            delete $ReactionRow->{"ASSOCIATED PEG"};
2299                                                            push(@{$ReactionRow->{"ASSOCIATED PEG"}},keys(%pegHash));
2300                                                            push(@{$ReactionRow->{"REFERENCE"}},$count.":".$ave.":".$stdev.":".$aveexp.":".$stdevexp.":".$min.":".$max);
2301                                                          if (defined($GeneSubsystems)) {                                                          if (defined($GeneSubsystems)) {
2302                                                                  push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});                                                                  push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems});
2303                                                          }                                                          }
# Line 1823  Line 2341 
2341          }          }
2342    
2343          #Clearing the previous model from the cache          #Clearing the previous model from the cache
2344          $self->figmodel()->ClearDBModel($self->id(),1);          $self->figmodel()->database()->ClearDBModel($self->id(),1);
2345          $ModelTable->save();          $ModelTable->save();
2346    
2347          return $self->success();          return $self->success();
# Line 1839  Line 2357 
2357  sub ArchiveModel {  sub ArchiveModel {
2358          my ($self) = @_;          my ($self) = @_;
2359    
         #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();  
         }  
   
2360          #Checking that the model file exists          #Checking that the model file exists
2361          if (!(-e $self->filename())) {          if (!(-e $self->filename())) {
2362                  $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 2407 
2407          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
2408  =cut  =cut
2409  sub run_microarray_analysis {  sub run_microarray_analysis {
2410          my ($self,$media,$jobid,$index,$genecall) = @_;          my ($self,$media,$label,$index,$genecall) = @_;
2411          $genecall =~ s/_/:/g;          $genecall =~ s/_/:/g;
2412          $genecall =~ s/\//;/g;          $genecall =~ s/\//;/g;
2413          #print "\n\n".$genecall."\n\n";          my $uniqueFilename = $self->figmodel()->filename();
2414          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";  
2415          system($command);          system($command);
2416          #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";
2417            if (-e $filename) {
2418                    my $output = $self->figmodel()->database()->load_single_column_file($filename);
2419                    if (defined($output->[0])) {
2420                            my @array = split(/;/,$output->[0]);
2421                            $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2422                            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]);
2423                    }
2424                    print STDERR $filename." is empty!";
2425            }
2426            print STDERR $filename." not found!";
2427            $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt");
2428    
2429            return undef;
2430  }  }
2431    
2432  =head3 find_minimal_pathways  =head3 find_minimal_pathways
# Line 1912  Line 2436 
2436          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
2437  =cut  =cut
2438  sub find_minimal_pathways {  sub find_minimal_pathways {
2439          my ($self,$media,$objective,$solutionnum) = @_;          my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2440    
2441          #Setting default media          #Setting default media
2442          if (!defined($media)) {          if (!defined($media)) {
# Line 1924  Line 2448 
2448                  $solutionnum = "5";                  $solutionnum = "5";
2449          }          }
2450    
2451            #Setting additional exchange fluxes
2452            if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2453                    if ($self->id() eq "iAF1260") {
2454                            $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2455                    } else {
2456                            $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2457                    }
2458            }
2459    
2460          #Translating objective          #Translating objective
2461          my $objectivestring;          my $objectivestring;
2462          if ($objective eq "ALL") {          if ($objective eq "ALL") {
# Line 1944  Line 2477 
2477                          }                          }
2478                  }                  }
2479                  for (my $i=0; $i < @objectives; $i++) {                  for (my $i=0; $i < @objectives; $i++) {
2480                          $self->find_minimal_pathways($media,$objectives[$i])                          $self->find_minimal_pathways($media,$objectives[$i]);
2481                  }                  }
2482                    return;
2483          } elsif ($objective eq "ENERGY") {          } elsif ($objective eq "ENERGY") {
2484                  $objectivestring = "MAX;FLUX;rxn00062;c;1";                  $objectivestring = "MAX;FLUX;rxn00062;c;1";
2485          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {          } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2486                    if ($objective =~ m/\[(\w)\]/) {
2487                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2488                            $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2489                    } else {
2490                  $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";                  $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2491                            $additionalexchange .= ";".$objective."[c]:-100:0";
2492                    }
2493            } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2494                    my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2495                    for (my $i=0; $i < @{$Products};$i++) {
2496                            my $temp = $Products->[$i]->{"DATABASE"}->[0];
2497                            if ($additionalexchange !~ m/$temp/) {
2498                                    #$additionalexchange .= ";".$temp."[c]:-100:0";
2499                            }
2500                    }
2501                    for (my $i=0; $i < @{$Reactants};$i++) {
2502                            print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2503                            $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2504                            print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2505                    }
2506                    return;
2507          }          }
2508    
2509          #Setting additional exchange fluxes          #Adding additional drains
2510          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") {  
2511                  $additionalexchange .= ";cpd15666[c]:0:100";                  $additionalexchange .= ";cpd15666[c]:0:100";
2512          } elsif ($objective eq "cpd11493") {          } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2513                  $additionalexchange .= ";cpd12370[c]:0:100";                  $additionalexchange .= ";cpd12370[c]:0:100";
2514          } elsif ($objective eq "cpd00166") {          } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2515                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2516          }          }
2517    
2518          #Running MFAToolkit          #Running MFAToolkit
2519          my $filename = $self->figmodel()->filename();          my $filename = $self->figmodel()->filename();
2520          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;
2521            if (defined($AllReversible) && $AllReversible == 1) {
2522                    $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());
2523            } else {
2524                    $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());
2525            }
2526          system($command);          system($command);
2527    
2528          #Loading problem report          #Loading problem report
# Line 1972  Line 2530 
2530          #Clearing output          #Clearing output
2531          $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");          $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2532          if (!defined($results)) {          if (!defined($results)) {
2533                    print STDERR $objective." pathway results not found!\n";
2534                  return;                  return;
2535          }          }
2536    
# Line 1983  Line 2542 
2542                  @Array = /\d+:([^\|]+)\|/g;                  @Array = /\d+:([^\|]+)\|/g;
2543          }          }
2544    
2545          # Storing data in a figmodel table          #Writing output to file
2546          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();  
2547  }  }
2548    
2549  =head3 calculate_growth  =head3 find_minimal_pathways
2550  Definition:  Definition:
2551          string::growth = FIGMODELmodel->calculate_growth(string:media);          int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective);
2552  Description:  Description:
2553          Calculating growth in the input media          Runs microarray analysis attempting to turn off genes that are inactive in the microarray
2554  =cut  =cut
2555  sub calculate_growth {  sub find_minimal_pathways_two {
2556          my ($self,$Media) = @_;          my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_;
2557          my $UniqueFilename = $self->figmodel()->filename();  
2558          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
2559          my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);          if (!defined($media)) {
2560          my $Result;                  $media = "Complete";
2561          if (defined($ProblemReport)) {          }
2562                  my $Row = $ProblemReport->get_row(0);  
2563                  if (defined($Row) && defined($Row->{"Objective"}->[0])) {          #Setting default solution number
2564                          if ($Row->{"Objective"}->[0] < 0.00000001) {          if (!defined($solutionnum)) {
2565                                  $Result = "NOGROWTH:".$Row->{"Individual metabolites with zero production"}->[0];                  $solutionnum = "5";
2566            }
2567    
2568            #Setting additional exchange fluxes
2569            if (!defined($additionalexchange) || length($additionalexchange) == 0) {
2570                    if ($self->id() eq "iAF1260") {
2571                            $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0";
2572                          } else {                          } else {
2573                                  $Result = $Row->{"Objective"}->[0];                          $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0];
2574                    }
2575            }
2576    
2577            #Translating objective
2578            my $objectivestring;
2579            if ($objective eq "ALL") {
2580                    #Getting the list of universal building blocks
2581                    my $buildingblocks = $self->config("universal building blocks");
2582                    my @objectives = keys(%{$buildingblocks});
2583                    #Getting the nonuniversal building blocks
2584                    my $otherbuildingblocks = $self->config("nonuniversal building blocks");
2585                    my @array = keys(%{$otherbuildingblocks});
2586                    if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) {
2587                            my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0];
2588                            if (defined($equation)) {
2589                                    for (my $i=0; $i < @array; $i++) {
2590                                            if (CORE::index($equation,$array[$i]) > 0) {
2591                                                    push(@objectives,$array[$i]);
2592                          }                          }
2593                  }                  }
2594          }          }
2595          return $Result;                  }
2596                    for (my $i=0; $i < @objectives; $i++) {
2597                            $self->find_minimal_pathways($media,$objectives[$i]);
2598                    }
2599                    return;
2600            } elsif ($objective eq "ENERGY") {
2601                    $objectivestring = "MAX;FLUX;rxn00062;c;1";
2602            } elsif ($objective =~ m/cpd\d\d\d\d\d/) {
2603                    if ($objective =~ m/\[(\w)\]/) {
2604                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1";
2605                            $additionalexchange .= ";".$objective."[".$1."]:-100:0";
2606                    } else {
2607                            $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1";
2608                            $additionalexchange .= ";".$objective."[c]:-100:0";
2609                    }
2610            } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) {
2611                    my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective);
2612                    for (my $i=0; $i < @{$Products};$i++) {
2613                            my $temp = $Products->[$i]->{"DATABASE"}->[0];
2614                            if ($additionalexchange !~ m/$temp/) {
2615                                    #$additionalexchange .= ";".$temp."[c]:-100:0";
2616                            }
2617                    }
2618                    for (my $i=0; $i < @{$Reactants};$i++) {
2619                            print $Reactants->[$i]->{"DATABASE"}->[0]." started\n";
2620                            $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange);
2621                            print $Reactants->[$i]->{"DATABASE"}->[0]." done\n";
2622                    }
2623                    return;
2624  }  }
2625    
2626  =head3 classify_model_reactions          #Adding additional drains
2627  Definition:          if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) {
2628          (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);                  $additionalexchange .= ";cpd15666[c]:0:100";
2629  Description:          } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) {
2630          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.                  $additionalexchange .= ";cpd12370[c]:0:100";
2631          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".          } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) {
2632          Possible values for "Class" include:                  $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100";
2633            }
2634    
2635            #Running MFAToolkit
2636            my $filename = $self->figmodel()->filename();
2637            my $command;
2638            if (defined($AllReversible) && $AllReversible == 1) {
2639                    $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());
2640            } else {
2641                    $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());
2642            }
2643            print $command."\n";
2644            system($command);
2645    
2646            #Loading problem report
2647            my $results = $self->figmodel()->LoadProblemReport($filename);
2648            #Clearing output
2649            $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt");
2650            if (!defined($results)) {
2651                    print STDERR $objective." pathway results not found!\n";
2652                    return;
2653            }
2654    
2655            #Parsing output
2656            my @Array;
2657            my $row = $results->get_row(1);
2658            if (defined($row->{"Notes"}->[0])) {
2659                    $_ = $row->{"Notes"}->[0];
2660                    @Array = /\d+:([^\|]+)\|/g;
2661            }
2662    
2663            #Writing output to file
2664            $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]);
2665    }
2666    
2667    sub combine_minimal_pathways {
2668            my ($self) = @_;
2669    
2670            my $tbl;
2671            if (-e $self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl") {
2672                    $tbl = FIGMODELTable::load_table($self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",";","|",0,["Objective","Media","Reversible"]);
2673            } else {
2674                    $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"],";","|");
2675            }
2676            my @files = glob($self->directory()."MinimalPathways-*");
2677            for (my $i=0; $i < @files;$i++) {
2678                    if ($files[$i] =~ m/MinimalPathways\-(\S+)\-(cpd\d\d\d\d\d)\-(\w+)\-(\d)\-/ || $files[$i] =~ m/MinimalPathways\-(\S+)\-(ENERGY)\-(\w+)\-(\d)\-/) {
2679                            my $reactions = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2680                            if (defined($reactions) && @{$reactions} > 0 && length($reactions->[0]) > 0) {
2681                                    my $newrow = {"Objective"=>[$2],"Media"=>[$1],"Reversible"=>[$4]};
2682                                    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");
2683                                    if (!defined($row)) {
2684                                            $row = $tbl->add_row($newrow);
2685                                    }
2686                                    $row->{Reactions} = $self->figmodel()->database()->load_single_column_file($files[$i],"");
2687                                    delete($row->{"Shortest path"});
2688                                    delete($row->{"Number of essentials"});
2689                                    delete($row->{"Essentials"});
2690                                    delete($row->{"Length"});
2691                                    for (my $j=0; $j < @{$row->{Reactions}}; $j++) {
2692                                            my @array = split(/,/,$row->{Reactions}->[$j]);
2693                                            $row->{"Length"}->[$j] = @array;
2694                                            if (!defined($row->{"Shortest path"}->[0]) || $row->{"Length"}->[$j] < $row->{"Shortest path"}->[0]) {
2695                                                    $row->{"Shortest path"}->[0] = $row->{"Length"}->[$j];
2696                                            }
2697                                            $row->{"Number of essentials"}->[0] = 0;
2698                                            for (my $k=0; $k < @array;$k++) {
2699                                                    if ($array[$k] =~ m/(rxn\d\d\d\d\d)/) {
2700                                                            my $class = $self->get_reaction_class($1,1);
2701                                                            my $temp = $row->{Media}->[0].":Essential";
2702                                                            if ($class =~ m/$temp/) {
2703                                                                    $row->{"Number of essentials"}->[$j]++;
2704                                                                    if (!defined($row->{"Essentials"}->[$j]) && length($row->{"Essentials"}->[$j]) > 0) {
2705                                                                            $row->{"Essentials"}->[$j] = $array[$k];
2706                                                                    } else {
2707                                                                            $row->{"Essentials"}->[$j] .= ",".$array[$k];
2708                                                                    }
2709                                                            }
2710                                                    }
2711                                            }
2712                                    }
2713                            }
2714                    }
2715            }
2716            $tbl->save();
2717    }
2718    
2719    =head3 calculate_growth
2720    Definition:
2721            string::growth = FIGMODELmodel->calculate_growth(string:media);
2722    Description:
2723            Calculating growth in the input media
2724    =cut
2725    sub calculate_growth {
2726            my ($self,$Media,$outputDirectory,$InParameters,$saveLPFile) = @_;
2727            #Setting the Media
2728            if (!defined($Media) || length($Media) == 0) {
2729                    $Media = $self->autocompleteMedia();
2730            }
2731            #Setting parameters for the run
2732            my $DefaultParameters = $self->figmodel()->defaultParameters();
2733            if (defined($InParameters)) {
2734                    my @parameters = keys(%{$InParameters});
2735                    for (my $i=0; $i < @parameters; $i++) {
2736                            $DefaultParameters->{$parameters[$i]} = $InParameters->{$parameters[$i]};
2737                    }
2738            }
2739            $DefaultParameters->{"optimize metabolite production if objective is zero"} = 1;
2740            #Setting filenames
2741            my $UniqueFilename = $self->figmodel()->filename();
2742            if (!defined($outputDirectory)) {
2743                    $outputDirectory = $self->config("database message file directory")->[0];
2744            }
2745            my $fluxFilename = $outputDirectory."Fluxes-".$self->id()."-".$Media.".txt";
2746            my $cpdFluxFilename = $outputDirectory."CompoundFluxes-".$self->id()."-".$Media.".txt";
2747            #Running FBA
2748            #print $self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],$DefaultParameters,$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version())."\n";
2749            system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],$DefaultParameters,$self->id()."-".$Media."-GrowthTest.txt",undef,$self->selected_version()));
2750            #Saving LP file if requested
2751            if (defined($saveLPFile) && $saveLPFile == 1 && -e $self->figmodel()->{"MFAToolkit output directory"}->[0].$UniqueFilename."/CurrentProblem.lp") {
2752                    system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/CurrentProblem.lp ".$self->directory().$self->id().".lp");
2753            }
2754            my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);
2755            my $Result;
2756            if (defined($ProblemReport)) {
2757                    my $Row = $ProblemReport->get_row(0);
2758                    if (defined($Row) && defined($Row->{"Objective"}->[0])) {
2759                            if ($Row->{"Objective"}->[0] < 0.00000001 || $Row->{"Objective"}->[0] == 1e7) {
2760                                    $Result = "NOGROWTH";
2761                                    if (defined($Row->{"Individual metabolites with zero production"}->[0]) && $Row->{"Individual metabolites with zero production"}->[0] =~ m/cpd\d\d\d\d\d/) {
2762                                            $Result .= ":".$Row->{"Individual metabolites with zero production"}->[0];
2763                                    }
2764                            } else {
2765                                    if (-e $self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionReactionData0.txt") {
2766                                            system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionReactionData0.txt ".$fluxFilename);
2767                                            system("cp ".$self->figmodel()->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/SolutionCompoundData0.txt ".$cpdFluxFilename);
2768                                    }
2769                                    $Result = $Row->{"Objective"}->[0];
2770                            }
2771                    }
2772            }
2773            #Deleting files if necessary
2774            if ($self->figmodel()->config("preserve all log files")->[0] ne "yes") {
2775                    $self->figmodel()->cleardirectory($UniqueFilename);
2776                    unlink($self->figmodel()->config("database message file directory")->[0].$self->id()."-".$Media."-GrowthTest.txt");
2777            }
2778            #Returning result
2779            return $Result;
2780    }
2781    
2782    =head3 classify_model_reactions
2783    Definition:
2784            (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media);
2785    Description:
2786            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.
2787            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".
2788            Possible values for "Class" include:
2789          1.) Positive: these reactions are essential in the forward direction.          1.) Positive: these reactions are essential in the forward direction.
2790          2.) Negative: these reactions are essential in the reverse direction.          2.) Negative: these reactions are essential in the reverse direction.
2791          3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.          3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction.
# Line 2036  Line 2795 
2795          7.) Dead: these reactions are disconnected from the network.          7.) Dead: these reactions are disconnected from the network.
2796  =cut  =cut
2797  sub classify_model_reactions {  sub classify_model_reactions {
2798          my ($self,$Media) = @_;          my ($self,$Media,$SaveChanges) = @_;
2799    
2800          #Getting unique file for printing model output          #Getting unique file for printing model output
2801          my $UniqueFilename = $self->figmodel()->filename();          my $UniqueFilename = $self->figmodel()->filename();
# Line 2123  Line 2882 
2882                                  $CpdRow->{MEDIA}->[$index] = $Media;                                  $CpdRow->{MEDIA}->[$index] = $Media;
2883                          }                          }
2884                  }                  }
2885                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2886                  $cpdclasstable->save();                  $cpdclasstable->save();
2887          }          }
2888            }
2889          if (defined($ReactionTB)) {          if (defined($ReactionTB)) {
2890                  for (my $i=0; $i < $ReactionTB->size(); $i++) {                  for (my $i=0; $i < $ReactionTB->size(); $i++) {
2891                          my $Row = $ReactionTB->get_row($i);                          my $Row = $ReactionTB->get_row($i);
# Line 2179  Line 2940 
2940                                  $RxnRow->{MEDIA}->[$index] = $Media;                                  $RxnRow->{MEDIA}->[$index] = $Media;
2941                          }                          }
2942                  }                  }
2943                    if (!defined($SaveChanges) || $SaveChanges == 1) {
2944                  $rxnclasstable->save();                  $rxnclasstable->save();
2945          }          }
2946            }
2947          return ($rxnclasstable,$cpdclasstable);          return ($rxnclasstable,$cpdclasstable);
2948  }  }
2949    
# Line 2212  Line 2975 
2975    
2976          #Running simulations          #Running simulations
2977          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";  
2978          #Parsing the results          #Parsing the results
2979          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);
2980          if (!defined($Results)) {          if (!defined($Results)) {
# Line 2619  Line 3381 
3381                                          if (length($GapFillingRunSpecs) > 0) {                                          if (length($GapFillingRunSpecs) > 0) {
3382                                                  $GapFillingRunSpecs .= ";";                                                  $GapFillingRunSpecs .= ";";
3383                                          }                                          }
3384                                          $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].".txt:".$HeadingDataArray[3];                                          $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].":".$HeadingDataArray[3];
3385                                  } else {                                  } else {
3386                                          $SolutionExistedCount++;                                          $SolutionExistedCount++;
3387                                  }                                  }
# Line 2644  Line 3406 
3406          my $SolutionsFound = 0;          my $SolutionsFound = 0;
3407          my $GapFillingArray;          my $GapFillingArray;
3408          push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));          push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs));
3409          $self->datagapfill($GapFillingArray);          my $GapFillingResults = $self->datagapfill($GapFillingArray,"GFS");
3410            if (defined($GapFillingResults)) {
3411                    $SolutionsFound = 1;
3412            }
3413    
3414          if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {          if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) {
3415                  #Printing previous solutions to GFS file                  #Printing previous solutions to GFS file
# Line 2667  Line 3432 
3432          return $self->success();          return $self->success();
3433  }  }
3434    
3435  =head3 BuildSpecificBiomassReaction  =head3 SolutionReconciliation
3436  Definition:  Definition:
3437          FIGMODELmodel->BuildSpecificBiomassReaction();          FIGMODELmodel->SolutionReconciliation();
3438  Description:  Description:
3439            This is a wrapper for running the solution reconciliation algorithm on any model in the database.
3440            The algorithm performs a reconciliation of any gap filling solutions to identify the combination of solutions that results in the optimal model.
3441            This function prints out one output file in the Model directory: ReconciliationOutput.txt: this is a summary of the results of the reconciliation analysis
3442  =cut  =cut
 sub BuildSpecificBiomassReaction {  
         my ($self) = @_;  
3443    
3444          my $biomassrxn;  sub SolutionReconciliation {
3445          my $OrganismID = $self->genome();          my ($self,$GapFill,$Stage) = @_;
3446          #Checking for a biomass override  
3447          if (defined($self->config("biomass reaction override")->{$OrganismID})) {          #Setting the output filenames
3448                  $biomassrxn = $self->config("biomass reaction override")->{$OrganismID};          my $OutputFilename;
3449                  print "Overriding biomass template and selecting ".$biomassrxn." for ".$OrganismID.".\n";          my $OutputFilenameTwo;
3450          } else {#Creating biomass reaction from the template          if ($GapFill == 1) {
3451                  #Getting the genome stats                  $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GFReconciliation.txt";
3452                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());                  $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GFSRS.txt";
3453                  my $Class = $genomestats->{CLASS}->[0];          } else {
3454                  my $Name = $genomestats->{NAME}->[0];                  $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GGReconciliation.txt";
3455                    $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GGSRS.txt";
3456            }
3457    
3458            #In stage one, we run the reconciliation and create a test file to check combined solution performance
3459            if (!defined($Stage) || $Stage == 1) {
3460                    my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3461                    my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3462                    $Row->{"GF RECONCILATION TIMING"}->[0] = time()."-";
3463                    $GrowMatchTable->save();
3464                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3465    
3466                    #Getting a unique filename
3467                    my $UniqueFilename = $self->figmodel()->filename();
3468    
3469                    #Copying over the necessary files
3470                    if ($GapFill == 1) {
3471                            if (!-e $self->directory().$self->id().$self->selected_version()."-GFEM.txt") {
3472                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GFEM.txt file not found. Could not reconcile!";
3473                                    return 0;
3474                            }
3475                            if (!-e $self->directory().$self->id().$self->selected_version()."-OPEM.txt") {
3476                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-OPEM.txt file not found. Could not reconcile!";
3477                                    return 0;
3478                            }
3479                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GFEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GFEM.txt");
3480                            system("cp ".$self->directory().$self->id().$self->selected_version()."-OPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3481                            #Backing up and deleting the existing reconciliation file
3482                            if (-e $OutputFilename) {
3483                                    system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGFReconciliation.txt");
3484                                    unlink($OutputFilename);
3485                            }
3486                    } else {
3487                            if (!-e $self->directory().$self->id().$self->selected_version()."-GGEM.txt") {
3488                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGEM.txt file not found. Could not reconcile!";
3489                                    return 0;
3490                            }
3491                            if (!-e $self->directory().$self->id().$self->selected_version()."-GGOPEM.txt") {
3492                                    print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt file not found. Could not reconcile!";
3493                                    return 0;
3494                            }
3495                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GGEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GGEM.txt");
3496                            system("cp ".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt");
3497                            #Backing up and deleting the existing reconciliation file
3498                            if (-e $OutputFilename) {
3499                                    system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGGReconciliation.txt");
3500                                    unlink($OutputFilename);
3501                            }
3502                    }
3503    
3504                    #Running the reconciliation
3505                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),"NONE",["SolutionReconciliation"],{"Solution data for model optimization" => $UniqueFilename},"Reconciliation".$UniqueFilename.".log",undef,$self->selected_version()));
3506                    $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3507                    $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3508                    $Row->{"GF RECONCILATION TIMING"}->[0] .= time();
3509                    $GrowMatchTable->save();
3510                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3511    
3512                    #Loading the problem report from the reconciliation run
3513                    my $ReconciliatonOutput = $self->figmodel()->LoadProblemReport($UniqueFilename);
3514                    print $UniqueFilename."\n";
3515                    #Clearing output files
3516                    $self->figmodel()->clearing_output($UniqueFilename,"Reconciliation".$UniqueFilename.".log");
3517                    $ReconciliatonOutput->save("/home/chenry/Test.txt");
3518    
3519                    #Checking the a problem report was found and was loaded
3520                    if (!defined($ReconciliatonOutput) || $ReconciliatonOutput->size() < 1 || !defined($ReconciliatonOutput->get_row(0)->{"Notes"}->[0])) {
3521                            print STDERR "FIGMODEL:SolutionReconciliation: MFAToolkit output from SolutionReconciliation of ".$self->id()." not found!\n\n";
3522                            return 0;
3523                    }
3524    
3525                    #Processing the solutions
3526                    my $SolutionCount = 0;
3527                    my $ReactionSetHash;
3528                    my $SingleReactionHash;
3529                    my $ReactionDataHash;
3530                    for (my $n=0; $n < $ReconciliatonOutput->size(); $n++) {
3531                            if (defined($ReconciliatonOutput->get_row($n)->{"Notes"}->[0]) && $ReconciliatonOutput->get_row($n)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^;]+)/) {
3532                                    #Breaking up the solution into reaction sets
3533                                    my @ReactionSets = split(/\|/,$1);
3534                                    #Creating reaction lists for each set
3535                                    my $SolutionHash;
3536                                    for (my $i=0; $i < @ReactionSets; $i++) {
3537                                            if (length($ReactionSets[$i]) > 0) {
3538                                                    my @Alternatives = split(/:/,$ReactionSets[$i]);
3539                                                    for (my $j=1; $j < @Alternatives; $j++) {
3540                                                            if (length($Alternatives[$j]) > 0) {
3541                                                                    push(@{$SolutionHash->{$Alternatives[$j]}},$Alternatives[0]);
3542                                                            }
3543                                                    }
3544                                                    if (@Alternatives == 1) {
3545                                                            $SingleReactionHash->{$Alternatives[0]}->{$SolutionCount} = 1;
3546                                                            if (!defined($SingleReactionHash->{$Alternatives[0]}->{"COUNT"})) {
3547                                                                    $SingleReactionHash->{$Alternatives[0]}->{"COUNT"} = 0;
3548                                                            }
3549                                                            $SingleReactionHash->{$Alternatives[0]}->{"COUNT"}++;
3550                                                    }
3551                                            }
3552                                    }
3553                                    #Identifying reactions sets and storing the sets in the reactions set hash
3554                                    foreach my $Solution (keys(%{$SolutionHash})) {
3555                                            my $SetKey = join(",",sort(@{$SolutionHash->{$Solution}}));
3556                                            if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount})) {
3557                                                    $ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount} = 1;
3558                                                    if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"})) {
3559                                                            $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"} = 0;
3560                                                    }
3561                                                    $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"}++;
3562                                            }
3563                                            $ReactionSetHash->{$SetKey}->{$Solution}->{$SolutionCount} = 1;
3564                                            if (!defined($ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"})) {
3565                                                    $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"} = 0;
3566                                            }
3567                                            $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"}++;
3568                                    }
3569                                    $SolutionCount++;
3570                            }
3571                    }
3572    
3573                    #Handling the scenario where no solutions were found
3574                    if ($SolutionCount == 0) {
3575                            print STDERR "FIGMODEL:SolutionReconciliation: Reconciliation unsuccessful. No solution found.\n\n";
3576                            return 0;
3577                    }
3578    
3579                    #Printing results without solution performance figures. Also printing solution test file
3580                    open (RECONCILIATION, ">$OutputFilename");
3581                    #Printing the file heading
3582                    print RECONCILIATION "DATABASE;DEFINITION;REVERSIBLITY;DELTAG;DIRECTION;NUMBER OF SOLUTIONS";
3583                    for (my $i=0; $i < $SolutionCount; $i++) {
3584                            print RECONCILIATION ";Solution ".$i;
3585                    }
3586                    print RECONCILIATION "\n";
3587                    #Printing the singlet reactions first
3588                    my $Solutions;
3589                    print RECONCILIATION "SINGLET REACTIONS\n";
3590                    my @SingletReactions = keys(%{$SingleReactionHash});
3591                    for (my $j=0; $j < $SolutionCount; $j++) {
3592                            $Solutions->[$j]->{"BASE"} = $j;
3593                    }
3594                    for (my $i=0; $i < @SingletReactions; $i++) {
3595                            my $ReactionData;
3596                            if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3597                                    $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3598                            } else {
3599                                    my $Direction = substr($SingletReactions[$i],0,1);
3600                                    if ($Direction eq "+") {
3601                                            $Direction = "=>";
3602                                    } else {
3603                                            $Direction = "<=";
3604                                    }
3605                                    my $Reaction = substr($SingletReactions[$i],1);
3606                                    $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3607                                    $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3608                                    $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3609                                    if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3610                                            $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3611                                    }
3612                                    if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3613                                            $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3614                                    }
3615                                    if (!defined($ReactionData->{"DELTAG"}->[0])) {
3616                                            $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3617                                    }
3618                                    $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3619                            }
3620                            print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$SingleReactionHash->{$SingletReactions[$i]}->{"COUNT"};
3621                            for (my $j=0; $j < $SolutionCount; $j++) {
3622                                    print RECONCILIATION ";";
3623                                    if (defined($SingleReactionHash->{$SingletReactions[$i]}->{$j})) {
3624                                            $Solutions->[$j]->{$SingletReactions[$i]} = 1;
3625                                            $Solutions->[$j]->{"BASE"} = $j;
3626                                            print RECONCILIATION "|".$j."|";
3627                                    }
3628                            }
3629                            print RECONCILIATION "\n";
3630                    }
3631                    #Printing the reaction sets with alternatives
3632                    print RECONCILIATION "Reaction sets with alternatives\n";
3633                    my @ReactionSets = keys(%{$ReactionSetHash});
3634                    foreach my $ReactionSet (@ReactionSets) {
3635                            my $NewSolutions;
3636                            my $BaseReactions;
3637                            my $AltList = [$ReactionSet];
3638                            push(@{$AltList},keys(%{$ReactionSetHash->{$ReactionSet}}));
3639                            for (my $j=0; $j < @{$AltList}; $j++) {
3640                                    my $CurrentNewSolutions;
3641                                    my $Index;
3642                                    if ($j == 0) {
3643                                            print RECONCILIATION "NEW SET\n";
3644                                    } elsif ($AltList->[$j] ne $ReactionSet) {
3645                                            print RECONCILIATION "ALTERNATIVE SET\n";
3646                                            #For each base solution in which this set is represented, we copy the base solution to the new solution
3647                                            my $NewSolutionCount = 0;
3648                                            for (my $k=0; $k < $SolutionCount; $k++) {
3649                                                    if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3650                                                            if (defined($Solutions)) {
3651                                                                    $Index->{$k} = @{$Solutions} + $NewSolutionCount;
3652                                                            } else {
3653                                                                    $Index->{$k} = $NewSolutionCount;
3654                                                            }
3655                                                            if (defined($NewSolutions) && @{$NewSolutions} > 0) {
3656                                                                    $Index->{$k} += @{$NewSolutions};
3657                                                            }
3658                                                            $CurrentNewSolutions->[$NewSolutionCount] = {};
3659                                                            foreach my $Reaction (keys(%{$Solutions->[$k]})) {
3660                                                                    $CurrentNewSolutions->[$NewSolutionCount]->{$Reaction} = $Solutions->[$k]->{$Reaction};
3661                                                            }
3662                                                            $NewSolutionCount++;
3663                                                    }
3664                                            }
3665                                    }
3666                                    if ($j == 0 || $AltList->[$j] ne $ReactionSet) {
3667                                            my @SingletReactions = split(/,/,$AltList->[$j]);
3668                                            for (my $i=0; $i < @SingletReactions; $i++) {
3669                                                    #Adding base reactions to base solutions and set reactions the new solutions
3670                                                    if ($j == 0) {
3671                                                            push(@{$BaseReactions},$SingletReactions[$i]);
3672                                                    } else {
3673                                                            for (my $k=0; $k < @{$CurrentNewSolutions}; $k++) {
3674                                                                    $CurrentNewSolutions->[$k]->{$SingletReactions[$i]} = 1;
3675                                                            }
3676                                                    }
3677                                                    #Getting reaction data and printing reaction in output file
3678                                                    my $ReactionData;
3679                                                    if (defined($ReactionDataHash->{$SingletReactions[$i]})) {
3680                                                            $ReactionData = $ReactionDataHash->{$SingletReactions[$i]};
3681                                                    } else {
3682                                                            my $Direction = substr($SingletReactions[$i],0,1);
3683                                                            if ($Direction eq "+") {
3684                                                                    $Direction = "=>";
3685                                                            } else {
3686                                                                    $Direction = "<=";
3687                                                            }
3688                                                            my $Reaction = substr($SingletReactions[$i],1);
3689                                                            $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t");
3690                                                            $ReactionData->{"DIRECTIONS"}->[0] = $Direction;
3691                                                            $ReactionData->{"REACTIONS"}->[0] = $Reaction;
3692                                                            if (!defined($ReactionData->{"DEFINITION"}->[0])) {
3693                                                                    $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN";
3694                                                            }
3695                                                            if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) {
3696                                                                    $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN";
3697                                                            }
3698                                                            if (!defined($ReactionData->{"DELTAG"}->[0])) {
3699                                                                    $ReactionData->{"DELTAG"}->[0] = "UNKNOWN";
3700                                                            }
3701                                                            $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData;
3702                                                    }
3703                                                    print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{"COUNT"};
3704                                                    for (my $k=0; $k < $SolutionCount; $k++) {
3705                                                            print RECONCILIATION ";";
3706                                                            if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) {
3707                                                                    if ($j == 0) {
3708                                                                            print RECONCILIATION "|".$k."|";
3709                                                                    } else {
3710                                                                            print RECONCILIATION "|".$Index->{$k}."|";
3711                                                                    }
3712                                                            }
3713                                                    }
3714                                                    print RECONCILIATION "\n";
3715                                            }
3716                                            #Adding the current new solutions to the new solutions array
3717                                            if (defined($CurrentNewSolutions) && @{$CurrentNewSolutions} > 0) {
3718                                                    push(@{$NewSolutions},@{$CurrentNewSolutions});
3719                                            }
3720                                    }
3721                            }
3722                            #Adding the base reactions to all existing solutions
3723                            for (my $j=0; $j < @{$Solutions}; $j++) {
3724                                    if (defined($ReactionSetHash->{$ReactionSet}->{$ReactionSet}->{$Solutions->[$j]->{"BASE"}})) {
3725                                            foreach my $SingleReaction (@{$BaseReactions}) {
3726                                                    $Solutions->[$j]->{$SingleReaction} = 1;
3727                                            }
3728                                    }
3729                            }
3730                            #Adding the new solutions to the set of existing solutions
3731                            push(@{$Solutions},@{$NewSolutions});
3732                    }
3733                    close(RECONCILIATION);
3734                    #Now printing a file that defines all of the solutions in a format the testsolutions function understands
3735                    open (RECONCILIATION, ">$OutputFilenameTwo");
3736                    print RECONCILIATION "Experiment;Solution index;Solution cost;Solution reactions\n";
3737                    for (my $i=0; $i < @{$Solutions}; $i++) {
3738                            delete($Solutions->[$i]->{"BASE"});
3739                            print RECONCILIATION "SR".$i.";".$i.";10;".join(",",keys(%{$Solutions->[$i]}))."\n";
3740                    }
3741                    close(RECONCILIATION);
3742    
3743                    $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3744                    $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3745                    $Row->{"GF RECON TESTING TIMING"}->[0] = time()."-";
3746                    $Row->{"GF RECON SOLUTIONS"}->[0] = @{$Solutions};
3747                    $GrowMatchTable->save();
3748                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3749    
3750                    #Scheduling the solution testing
3751                    if ($GapFill == 1) {
3752                            system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GFSR:BACK:fast:QSUB\"");
3753                    } else {
3754                            system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GGSR:BACK:fast:QSUB\"");
3755                    }
3756            } else {
3757                    #Reading in the solution testing results
3758                    my $Data;
3759                    if ($GapFill == 1) {
3760                            $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GFSREM.txt","");
3761                    } else {
3762                            $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GGSREM.txt","");
3763                    }
3764    
3765                    #Reading in the preliminate reconciliation report
3766                    my $OutputData = $self->figmodel()->database()->load_single_column_file($OutputFilename,"");
3767                    #Replacing the file tags with actual performance data
3768                    my $Count = 0;
3769                    for (my $i=0; $i < @{$Data}; $i++) {
3770                            if ($Data->[$i] =~ m/^SR(\d+);.+;(\d+\/\d+);/) {
3771                                    my $Index = $1;
3772                                    my $Performance = $Index."/".$2;
3773                                    for (my $j=0; $j < @{$OutputData}; $j++) {
3774                                            $OutputData->[$j] =~ s/\|$Index\|/$Performance/g;
3775                                    }
3776                            }
3777                    }
3778                    $self->figmodel()->database()->print_array_to_file($OutputFilename,$OutputData);
3779    
3780                    my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE");
3781                    my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1);
3782                    $Row->{"GF RECON TESTING TIMING"}->[0] .= time();
3783                    $GrowMatchTable->save();
3784                    $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE");
3785            }
3786    
3787            return 1;
3788    }
3789    
3790    =head3 DetermineCofactorLipidCellWallComponents
3791    Definition:
3792            {cofactor=>{string:compound id=>float:coefficient},lipid=>...cellWall=>} = FIGMODELmodel->DetermineCofactorLipidCellWallComponents();
3793    Description:
3794    =cut
3795    sub DetermineCofactorLipidCellWallComponents {
3796            my ($self) = @_;
3797            my $templateResults;
3798            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
3799            my $Class = $self->{_data}->cellwalltype();
3800            my $Name = $self->name();
3801            my $translation = {COFACTOR=>"cofactor",LIPIDS=>"lipid","CELL WALL"=>"cellWall"};
3802                  #Checking for phoenix variants                  #Checking for phoenix variants
3803                  my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");                  my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");
3804                  my $Phoenix = 0;                  my $Phoenix = 0;
3805                  my @Rows = $PhoenixVariantTable->get_rows_by_key($OrganismID,"GENOME");          my @Rows = $PhoenixVariantTable->get_rows_by_key($self->genome(),"GENOME");
3806                  my $VariantHash;                  my $VariantHash;
3807                  for (my $i=0; $i < @Rows; $i++) {                  for (my $i=0; $i < @Rows; $i++) {
3808                          $Phoenix = 1;                          $Phoenix = 1;
# Line 2698  Line 3810 
3810                                  $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];                                  $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];
3811                          }                          }
3812                  }                  }
   
3813                  #Collecting genome data                  #Collecting genome data
3814                  my $RoleHash;                  my $RoleHash;
3815                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
# Line 2715  Line 3826 
3826                                  }                                  }
3827                          }                          }
3828                  }                  }
   
3829                  #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
3830                  my $ComponentTypes;          my $includedHash;
3831                  my $EquationData;          my $BiomassReactionTemplateTable = $self->figmodel()->database()->get_table("BIOMASSTEMPLATE");
                 my $BiomassReactionTemplateTable = $self->figmodel()->database()->GetDBTable("BIOMASS TEMPLATE");  
3832                  for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {                  for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {
3833                          my $Row = $BiomassReactionTemplateTable->get_row($i);                          my $Row = $BiomassReactionTemplateTable->get_row($i);
3834                    if (defined($translation->{$Row->{CLASS}->[0]})) {
3835                            my $coef = -1;
3836                            if ($Row->{"REACTANT"}->[0] eq "NO") {
3837                                    $coef = 1;
3838                                    if ($Row->{"COEFFICIENT"}->[0] =~ m/cpd/) {
3839                                            $coef = $Row->{"COEFFICIENT"}->[0];
3840                                    }
3841                            }
3842                          if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {                          if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {
3843                                  push(@{$EquationData},$Row);                                  $includedHash->{$Row->{"ID"}->[0]} = 1;
3844                                  $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;                                  $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
3845                          } else {                          } elsif (defined($Row->{"INCLUSION CRITERIA"}->[0])) {
3846                                  my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];                                  my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];
3847                                  my $End = 0;                                  my $End = 0;
3848                                  while ($End == 0) {                                  while ($End == 0) {
3849                                          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)\{([^{^}]+)\}$/) {
3850                                                  print $Criteria."\n";                                                  print $Criteria." : ";
3851                                                  my $Start = "";                                                  my $Start = "";
3852                                                  my $End = "";                                                  my $End = "";
3853                                                  my $Condition = $1;                                                  my $Condition = $1;
# Line 2754  Line 3871 
3871                                                                  $Result = "NO";                                                                  $Result = "NO";
3872                                                                  last;                                                                  last;
3873                                                          } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {                                                          } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {
3874                                                                  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") {  
3875                                                                          $Result = "YES";                                                                          $Result = "YES";
3876                                                                          last;                                                                          last;
3877                                                                  } elsif ($Match != 1 && $Condition eq "AND") {                                                                  } elsif (!defined($includedHash->{$1}) && $Condition eq "AND") {
3878                                                                          $Result = "NO";                                                                          $Result = "NO";
3879                                                                          last;                                                                          last;
3880                                                                  }                                                                  }
3881                                                          } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {                                                          } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {
3882                                                                  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") {  
3883                                                                          $Result = "YES";                                                                          $Result = "YES";
3884                                                                          last;                                                                          last;
3885                                                                  } elsif ($Match == 1 && $Condition eq "AND") {                                                                  } elsif (defined($includedHash->{$1}) && $Condition eq "AND") {
3886                                                                          $Result = "NO";                                                                          $Result = "NO";
3887                                                                          last;                                                                          last;
3888                                                                  }                                                                  }
# Line 2895  Line 3998 
3998                                          }                                          }
3999                                  }                                  }
4000                                  if ($Criteria eq "YES") {                                  if ($Criteria eq "YES") {
4001                                          push(@{$EquationData},$Row);                                          $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
4002                                          $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;                                          $includedHash->{$Row->{"ID"}->[0]} = 1;
4003                                  }                                  }
4004                          }                          }
4005                  }                  }
4006            }
4007                  #Building biomass equation          my $types = ["cofactor","lipid","cellWall"];
4008                  my %Reactants;          my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4009                  my %Products;          for (my $i=0; $i < @{$types}; $i++) {
4010                  foreach my $EquationRow (@{$EquationData}) {                  my @list =      keys(%{$templateResults->{$types->[$i]}});
4011                          #First determine what the coefficient should be                  my $entries = 0;
4012                          my $Coefficient;                  for (my $j=0; $j < @list; $j++) {
4013                          if (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^[0-9\.]+$/) {                          if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4014                                  $Coefficient = $EquationRow->{"COEFFICIENT"}->[0];                                  my $objs = $cpdMgr->get_objects({id=>$list[$j]});
4015                          } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/cpd\d\d\d\d\d/) {                                  if (!defined($objs->[0]) || $objs->[0]->mass() == 0) {
4016                                  $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]};  
4017                                  } else {                                  } else {
4018                                          $Products{$EquationRow->{"ID"}->[0]} = $Coefficient;                                          $entries++;
4019                                    }
4020                            }
4021                                  }                                  }
4022                    for (my $j=0; $j < @list; $j++) {
4023                            if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4024                                    $templateResults->{$types->[$i]}->{$list[$j]} = -1/$entries;
4025                            } elsif ($templateResults->{$types->[$i]}->{$list[$j]} =~ m/cpd/) {
4026                                    my $netCoef = 0;
4027                                    my @allcpd = split(/,/,$templateResults->{$types->[$i]}->{$list[$j]});
4028                                    for (my $k=0; $k < @allcpd; $k++) {
4029                                            if (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} ne "-1e-5") {
4030                                                    $netCoef += (1/$entries);
4031                                            } elsif (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} eq "-1e-5") {
4032                                                    $netCoef += 1e-5;
4033                          }                          }
4034                  }                  }
4035                                    $templateResults->{$types->[$i]}->{$list[$j]} = $netCoef;
4036                            }
4037                    }
4038            }
4039            return $templateResults;
4040    }
4041    
4042    =head3 BuildSpecificBiomassReaction
4043    Definition:
4044            FIGMODELmodel->BuildSpecificBiomassReaction();
4045    Description:
4046    =cut
4047    sub BuildSpecificBiomassReaction {
4048            my ($self) = @_;
4049            #Getting the database handle for biomass reactions
4050            my $bioMgr = $self->figmodel()->database()->get_object_manager("bof");
4051            #Checking if the current biomass reaction appears in more than on model, if not, this biomass reaction is conserved for this model
4052            my $biomassID = $self->biomassReaction();
4053            if ($biomassID =~ m/bio\d\d\d\d\d/) {
4054                    my $mdlMgr = $self->figmodel()->database()->get_object_manager("model");
4055                    my $mdlObs = $mdlMgr->get_objects({biomassReaction=>$biomassID});
4056                    if (defined($mdlObs->[1])) {
4057                            $biomassID = "NONE";
4058                    }
4059            }
4060            #If the biomass ID is "NONE", then we create a new biomass reaction for the model
4061            my $bioObj;
4062            my $originalPackages = "";
4063            my $originalEssReactions = "";
4064            if ($biomassID !~ m/bio\d\d\d\d\d/) {
4065                    #Getting the current largest ID
4066                    $biomassID = $self->figmodel()->database()->check_out_new_id("bof");
4067                    $bioObj = $bioMgr->create({
4068                            id=>$biomassID,owner=>$self->owner(),users=>$self->users(),name=>"Biomass",equation=>"NONE",protein=>"0",
4069                            energy=>"0",DNA=>"0",RNA=>"0",lipid=>"0",cellWall=>"0",cofactor=>"0",
4070                            modificationDate=>time(),creationDate=>time(),
4071                            cofactorPackage=>"NONE",lipidPackage=>"NONE",cellWallPackage=>"NONE",
4072                            DNACoef=>"NONE",RNACoef=>"NONE",proteinCoef=>"NONE",lipidCoef=>"NONE",
4073                            cellWallCoef=>"NONE",cofactorCoef=>"NONE",essentialRxn=>"NONE"});
4074                    if (!defined($bioObj)) {
4075                            die $self->error_message("BuildSpecificBiomassReaction():Could not create new biomass reaction ".$biomassID."!");
4076                    }
4077            } else {
4078                    #Getting the biomass DB handler from the database
4079                    my $objs = $bioMgr->get_objects({id=>$biomassID});
4080                    if (!defined($objs->[0])) {
4081                            die $self->error_message("BuildSpecificBiomassReaction():Could not find biomass reaction ".$biomassID." in database!");
4082                    }
4083                    $bioObj = $objs->[0];
4084                    $bioObj->owner($self->owner());
4085                    $bioObj->users($self->users());
4086                    if (defined($bioObj->essentialRxn())) {
4087                            $originalEssReactions = $bioObj->essentialRxn();
4088                            $originalPackages = $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage();
4089                    }
4090            }
4091            #Getting genome stats
4092            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
4093            my $Class = $self->{_data}->cellwalltype();
4094            #Setting global coefficients based on cell wall type
4095            my $biomassCompounds;
4096            my $compounds;
4097            if ($Class eq "Gram positive") {
4098                    $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.323,cpd00052=>-0.199,cpd00062=>-0.215};
4099                    $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};
4100                    $bioObj->protein("0.5284");
4101                    $bioObj->DNA("0.026");
4102                    $bioObj->RNA("0.0655");
4103                    $bioObj->lipid("0.075");
4104                    $bioObj->cellWall("0.25");
4105                    $bioObj->cofactor("0.10");
4106            } else {
4107                    $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.322,cpd00052=>-0.2,cpd00062=>-0.216};
4108                    $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};
4109                    $bioObj->protein("0.563");
4110                    $bioObj->DNA("0.031");
4111                    $bioObj->RNA("0.21");
4112                    $bioObj->lipid("0.093");
4113                    $bioObj->cellWall("0.177");
4114                    $bioObj->cofactor("0.039");
4115            }
4116            #Setting energy coefficient for all reactions
4117            $bioObj->energy("40");
4118            $compounds->{energy} = {cpd00002=>-1,cpd00001=>-1,cpd00008=>1,cpd00009=>1,cpd00067=>1};
4119            #Setting DNA coefficients based on GC content
4120            my $gc = $self->figmodel()->get_genome_gc_content($self->genome());
4121            $compounds->{DNA} = {cpd00012=>1,cpd00115=>0.5*(1-$gc),cpd00241=>0.5*$gc,cpd00356=>0.5*$gc,cpd00357=>0.5*(1-$gc)};
4122            #Setting Lipid,cell wall,and cofactor coefficients based on biomass template
4123            my $templateResults = $self->DetermineCofactorLipidCellWallComponents();
4124            $compounds->{cofactor} = $templateResults->{cofactor};
4125            $compounds->{lipid} = $templateResults->{lipid};
4126            $compounds->{cellWall} = $templateResults->{cellWall};
4127            #Getting package number for cofactor, lipid, and cell wall
4128            my $packages;
4129            my $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4130            my $packageTypes = ["Cofactor","Lipid","CellWall"];
4131            my $translation = {"Cofactor"=>"cofactor","Lipid"=>"lipid","CellWall"=>"cellWall"};
4132            for (my $i=0; $i < @{$packageTypes}; $i++) {
4133                    my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4134                    my $function = $translation->{$packageTypes->[$i]}."Package";
4135                    if (@cpdList == 0) {
4136                            $bioObj->$function("NONE");
4137                    } else {
4138                            my $cpdgrpObs = $cpdgrpMgr->get_objects({type=>$packageTypes->[$i]."Package"});
4139                            for (my $j=0; $j < @{$cpdgrpObs}; $j++) {
4140                                    $packages->{$packageTypes->[$i]}->{$cpdgrpObs->[$j]->grouping()}->{$cpdgrpObs->[$j]->COMPOUND()} = 1;
4141                            }
4142                            my @packageList = keys(%{$packages->{$packageTypes->[$i]}});
4143                            my $packageHash;
4144                            for (my $j=0; $j < @packageList; $j++) {
4145                                    $packageHash->{join("|",sort(keys(%{$packages->{$packageTypes->[$i]}->{$packageList[$j]}})))} = $packageList[$j];
4146                            }
4147                            if (defined($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))})) {
4148                                    $bioObj->$function($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))});
4149                            } else {
4150                                    my $newPackageID = $self->figmodel()->database()->check_out_new_id($packageTypes->[$i]."Pkg");
4151                                    $bioObj->$function($newPackageID);
4152                                    my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4153                                    for (my $j=0; $j < @cpdList; $j++) {
4154                                            $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4155                                            $cpdgrpMgr->create({COMPOUND=>$cpdList[$j],grouping=>$newPackageID,type=>$packageTypes->[$i]."Package"});
4156                                    }
4157                            }
4158                    }
4159            }
4160            #Filling in coefficient terms in database and calculating global reaction coefficients based on classification abundancies
4161            my $equationCompounds;
4162            my $types = ["RNA","DNA","protein","lipid","cellWall","cofactor","energy"];
4163            my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4164            for (my $i=0; $i < @{$types}; $i++) {
4165                    my $coefString = "";
4166                    my @compounds = sort(keys(%{$compounds->{$types->[$i]}}));
4167                    #Building coefficient strings and determining net mass for component types
4168                    my $netMass = 0;
4169                    for (my $j=0; $j < @compounds; $j++) {
4170                            my $objs = $cpdMgr->get_objects({id=>$compounds[$j]});
4171                            my $mass = 0;
4172                            if (defined($objs->[0]) && $objs->[0]->mass() != 0) {
4173                                    $mass = $objs->[0]->mass();
4174                                    $netMass += -$compounds->{$types->[$i]}->{$compounds[$j]}*$objs->[0]->mass();
4175                            }
4176                            if (!defined($equationCompounds->{$compounds[$j]})) {
4177                                    $equationCompounds->{$compounds[$j]}->{"coef"} = 0;
4178                                    $equationCompounds->{$compounds[$j]}->{"type"} = $types->[$i];
4179                                    $equationCompounds->{$compounds[$j]}->{"mass"} = $mass;
4180                            }
4181                            $coefString .= $compounds->{$types->[$i]}->{$compounds[$j]}."|";
4182                    }
4183                    $netMass = 0.001*$netMass;
4184                    #Calculating coefficients for all component compounds
4185                    for (my $j=0; $j < @compounds; $j++) {
4186                            #Normalizing certain type coefficients by mass
4187                            my $function = $types->[$i];
4188                            my $fraction = $bioObj->$function();
4189                            if ($types->[$i] ne "energy") {