[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.18, Mon May 31 18:05:08 2010 UTC revision 1.28, Wed Aug 25 20:39:16 2010 UTC
# 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 26  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            if (!defined($metagenome) || $metagenome != 1) {
30          my $modelHandle = $self->figmodel()->database()->get_object_manager("model");          my $modelHandle = $self->figmodel()->database()->get_object_manager("model");
31          if (!defined($modelHandle)) {                  if (defined($modelHandle)) {
                 $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!");  
                 return undef;  
         }  
   
32          #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
33          if ($id =~ m/^\d+$/) {          if ($id =~ m/^\d+$/) {
34                  my $objects = $modelHandle->get_objects();                  my $objects = $modelHandle->get_objects();
35                                    if (defined($objects->[$id])) {
36                  $self->{_data} = $objects->[$id];                  $self->{_data} = $objects->[$id];
37                  $self->figmodel()->{_models}->{$id} = $self;                  $self->figmodel()->{_models}->{$id} = $self;
38                                    }
39          } else {          } else {
40                  my $objects = $modelHandle->get_objects({id => $id});                  my $objects = $modelHandle->get_objects({id => $id});
41                  if (!defined($objects->[0]) && $id =~ m/(.+)(V[^V]+)/) {                                  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});                          $objects = $modelHandle->get_objects({id => $1});
45                          if (!defined($objects->[0])) {                                          if (defined($objects->[0])) {
                                 $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!");  
                                 return undef;  
                         }  
46                          $self->{_selectedversion} = $2;                          $self->{_selectedversion} = $2;
47                  } elsif (!defined($objects->[0])) {                                                  $self->{_data} = $objects->[0];
48                          $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!");                                          }
                         return undef;  
49                  }                  }
50                            }
51                    }
52                    if (defined($self->{_data})) {
53                            $self->{_modeltype} = "genome";
54                    }
55            }
56            if (!defined($self->{_data})) {
57                    my $modelHandle = $self->figmodel()->database()->get_object_manager("mgmodel");
58                    if (defined($modelHandle)) {
59                            #If the id is a number, we get the model row by index
60                            if ($id =~ m/^\d+$/) {
61                                    my $objects = $modelHandle->get_objects();
62                                    if (defined($objects->[$id])) {
63                                            $self->{_data} = $objects->[$id];
64                                            $self->figmodel()->{_models}->{"MG".$id} = $self;
65                                    }
66                            } else {
67                                    my $objects = $modelHandle->get_objects({id => $id});
68                                    if (defined($objects->[0])) {
69                  $self->{_data} = $objects->[0];                  $self->{_data} = $objects->[0];
70                                    } elsif (!defined($objects->[0]) && $id =~ m/(.+)(V[^V]+)/) {
71                                            $objects = $modelHandle->get_objects({id => $1});
72                                            if (defined($objects->[0])) {
73                                                    $self->{_selectedversion} = $2;
74                                                    $self->{_data} = $objects->[0];
75                                            }
76                                    }
77                            }
78                    }
79                    if (defined($self->{_data})) {
80                            $self->{_modeltype} = "metagenome";
81                    }
82          }          }
83          if (!defined($self->{_data})) {          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          }          }
87          $self->figmodel()->{_models}->{$self->id()} = $self;          $self->figmodel()->{_models}->{$self->id()} = $self;
   
88          return $self;          return $self;
89  }  }
90    
# Line 75  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 166  Line 201 
201  =cut  =cut
202  sub mgdata {  sub mgdata {
203          my ($self) = @_;          my ($self) = @_;
204          if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) {          if (!defined($self->{_mgdata})) {
205                  require MGRAST;                  my $mgrastdbhandle = $self->figmodel()->database()->get_object_manager("mgjob");
206                  $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } )                  my $objects = $mgrastdbhandle->get_objects( { 'genome_id' => $self->genome() } );
207                    $self->{_mgdata} = $objects->[0];
208          }          }
209          return $self->{_mgdata};          return $self->{_mgdata};
210  }  }
# Line 184  Line 220 
220          return $self->{_data}->id();          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
235  Definition:  Definition:
236          string = FIGMODELmodel->owner();          string = FIGMODELmodel->owner();
# Line 195  Line 242 
242          return $self->{_data}->owner();          return $self->{_data}->owner();
243  }  }
244    
245    =head3 users
246    Definition:
247            string = FIGMODELmodel->users();
248    Description:
249    =cut
250    sub users {
251            my ($self) = @_;
252            return $self->{_data}->users();
253    }
254    
255  =head3 name  =head3 name
256  Definition:  Definition:
257          string = FIGMODELmodel->name();          string = FIGMODELmodel->name();
# Line 248  Line 305 
305                                  $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $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",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $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>";                                  $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",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $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>";                                  $brief_flux ? $class.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $class.="<br>[Flux: ".$min." to ".$max."]<br>";
# Line 295  Line 352 
352                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
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                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $NewClass.="<br>[Flux: ".$max." to ".$min."]<br>";
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                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
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                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$max)." to ".sprintf("%.3g",$min)."]<br>" : $NewClass.="<br>[Flux: ".$max." to ".$min."]<br>";
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                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";                                  $brief_flux ? $NewClass.="<br>[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]<br>" : $NewClass.="<br>[Flux: ".$min." to ".$max."]<br>";
# Line 352  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 387  Line 443 
443  =cut  =cut
444  sub load_model_table {  sub load_model_table {
445          my ($self,$name,$refresh) = @_;          my ($self,$name,$refresh) = @_;
446          if (!defined($refresh)) {          if (defined($refresh) && $refresh == 1) {
                 $refresh = 1;  
         }  
         if ($refresh == 1) {  
447                  delete $self->{"_".$name};                  delete $self->{"_".$name};
448          }          }
449          if (!defined($self->{"_".$name})) {          if (!defined($self->{"_".$name})) {
450                  my $tbldef = $self->figmodel()->config("$name");                  my $tbldef = $self->figmodel()->config($name);
451                  if (!defined($tbldef)) {                  if (!defined($tbldef)) {
452                          return undef;                          return undef;
453                  }                  }
454                  my $filename = $self->directory().$name."-".$self->id().$self->selected_version().".tbl";                  my $itemDelim = "|";
455                  $self->{"_".$name} = $self->figmodel()->database()->load_table($filename,"\t","|",$tbldef->{headingline}->[0],$tbldef->{hashcolumns});                  if (defined($tbldef->{itemdelimiter}->[0])) {
456                  if (!defined($self->{"_".$name})) {                          $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})) {                          if (defined($tbldef->{prefix})) {
484                                  $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},"\t","|",join(@{$tbldef->{prefix}},"\n"));                                  $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim,join(@{$tbldef->{prefix}},"\n"));
485                          } else {                          } else {
486                                  $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},"\t","|");                                  $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim);
487                          }                          }
488                  }                  }
489          }          }
# Line 421  Line 500 
500  =cut  =cut
501  sub create_table_prototype {  sub create_table_prototype {
502          my ($self,$TableName) = @_;          my ($self,$TableName) = @_;
   
503          #Checking if the table definition exists in the FIGMODELconfig file          #Checking if the table definition exists in the FIGMODELconfig file
504          if (!defined($self->figmodel()->config($TableName))) {          my $tbldef = $self->figmodel()->config($TableName);
505            if (!defined($tbldef)) {
506                  $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:Definition not found for ".$TableName);                  $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:Definition not found for ".$TableName);
507                  return undef;                  return undef;
508          }          }
509          #Checking that this is a database table          #Checking that this is a database table
510          if (!defined($self->config($TableName)->{tabletype}) || $self->config($TableName)->{tabletype}->[0] ne "ModelTable") {          if (!defined($tbldef->{tabletype}) || $tbldef->{tabletype}->[0] ne "ModelTable") {
511                  $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:".$TableName." is not a model table!");                  $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:".$TableName." is not a model table!");
512                  return undef;                  return undef;
513          }          }
514          if (!defined($self->config($TableName)->{delimiter})) {          #Setting default values for table parameters
515                  $self->config($TableName)->{delimiter}->[0] = ";";          my $prefix;
516            if (defined($tbldef->{prefix})) {
517                    $prefix = join("\n",@{$self->config($TableName)->{prefix}})."\n";
518          }          }
519          if (!defined($self->config($TableName)->{itemdelimiter})) {          my $itemDelim = "|";
520                  $self->config($TableName)->{itemdelimiter}->[0] = "|";          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          }          }
         my $prefix;  
         if (defined($self->config($TableName)->{prefix})) {  
                 $prefix = join("\n",@{$self->config($TableName)->{prefix}});  
544          }          }
545          my $tbl = FIGMODELTable->new($self->config($TableName)->{columns},$self->directory().$self->config($TableName)->{filename_prefix}->[0]."-".$self->id().$self->selected_version().".txt",$self->config($TableName)->{hashcolumns},$self->config($TableName)->{delimiter}->[0],$self->config($TableName)->{itemdelimiter}->[0],$prefix);          #Creating the table prototype
546            my $tbl = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},$columnDelim,$itemDelim,$prefix);
547          return $tbl;          return $tbl;
548  }  }
549    
# Line 467  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          if (!defined($self->{_reaction_data})) {                  delete $self->{_reaction_table};
574                  $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id());          }
575            if (!defined($self->{_reaction_table})) {
576                    $self->{_reaction_table} = $self->load_model_table("ModelReactions",$clear);
577                  my $classTbl = $self->reaction_class_table();                  my $classTbl = $self->reaction_class_table();
578                    if (defined($classTbl)) {
579                  for (my $i=0; $i < $classTbl->size(); $i++) {                  for (my $i=0; $i < $classTbl->size(); $i++) {
580                          my $row = $classTbl->get_row($i);                          my $row = $classTbl->get_row($i);
581                          my $rxnRow = $self->{_reaction_data}->get_row_by_key($row->{"REACTION"}->[0],"LOAD");                                  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++) {                          for (my $j=0; $j < @{$row->{MEDIA}};$j++) {
585                                  my $class = "Active <=>";                                  my $class = "Active <=>";
586                                  if ($row->{CLASS}->[$j] eq "Positive") {                                  if ($row->{CLASS}->[$j] eq "Positive") {
# Line 496  Line 602 
602                          }                          }
603                  }                  }
604          }          }
605                            }
606                    }
607            }
608            return $self->{_reaction_table};
609    }
610    
611    =head3 essentials_table
612    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 feature_table  =head3 feature_table
# Line 551  Line 682 
682                  }                  }
683              }              }
684              #Loading predictions              #Loading predictions
685              my @files = glob($self->directory()."EssentialGenes-".$self->id()."-*");                  my $esstbl = $self->essentials_table();
             foreach my $file (@files) {  
                 if ($file =~ m/\-([^\-]+)\.tbl/) {  
                         my $media = $1;  
                         my $list = $self->figmodel()->database()->load_single_column_file($file,"");  
                         my $hash;  
                         foreach my $gene (@{$list}) {  
                                 $hash->{$gene} = 1;  
                         }  
686                          for (my $i=0; $i < $self->{_feature_data}->size(); $i++) {                          for (my $i=0; $i < $self->{_feature_data}->size(); $i++) {
687                                  my $Row = $self->{_feature_data}->get_row($i);                                  my $Row = $self->{_feature_data}->get_row($i);
688                                  if ($Row->{ID}->[0] =~ m/(peg\.\d+)/) {                                  if ($Row->{ID}->[0] =~ m/(peg\.\d+)/) {
689                                          my $gene = $1;                                          my $gene = $1;
690                                          if (defined($hash->{$gene})) {                                  my @rows = $esstbl->get_rows_by_key($gene,"ESSENTIAL GENES");
691                                                  push(@{$Row->{$self->id()."PREDICTIONS"}},$media.":essential");                                  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                                                  push(@{$Row->{$self->id()."PREDICTIONS"}},$media.":nonessential");                                                  push(@{$Row->{$self->id()."PREDICTIONS"}},$mediaList[$j].":nonessential");
                                         }  
704                                  }                                  }
705                          }                          }
706                  }                  }
# Line 584  Line 716 
716          Returns FIGMODELTable with the reaction class data, and creates the table file  if it does not exist          Returns FIGMODELTable with the reaction class data, and creates the table file  if it does not exist
717  =cut  =cut
718  sub reaction_class_table {  sub reaction_class_table {
719          my ($self) = @_;          my ($self,$clear) = @_;
720            return $self->load_model_table("ModelReactionClasses",$clear);
         if (!defined($self->{_reaction_class_table})) {  
                 if (-e $self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl") {  
                         $self->{_reaction_class_table} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",";","|",0,["REACTION","CLASS","MEDIA"]);  
                 } else {  
                         $self->{_reaction_class_table} = FIGMODELTable->new(["REACTION","MEDIA","CLASS","MIN","MAX"],$self->directory()."ReactionClassification-".$self->id().$self->selected_version().".tbl",["REACTION","CLASS","MEDIA"],";","|",undef);  
                 }  
         }  
   
         return $self->{_reaction_class_table};  
721  }  }
722    
723  =head3 compound_class_table  =head3 compound_class_table
# Line 604  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);  
                 }  
         }  
   
         return $self->{_compound_class_table};  
732  }  }
733    
734  =head3 get_essential_genes  =head3 get_essential_genes
# Line 625  Line 739 
739  =cut  =cut
740  sub get_essential_genes {  sub get_essential_genes {
741          my ($self,$media) = @_;          my ($self,$media) = @_;
742            my $tbl = $self->essentials_table();
743          if (!defined($media)) {          my $row = $tbl->get_row_by_key($media,"MEDIA");
744                  $media = "Complete";          if (defined($row)) {
745                    return $row->{"ESSENTIAL GENES"};
746          }          }
747          if (!defined($self->{_essential_genes}->{$media})) {          return undef;
                 $self->{_essential_genes}->{$media} = undef;  
                 if (-e $self->directory()."EssentialGenes-".$self->id().$self->selected_version()."-".$media.".tbl") {  
                         $self->{_essential_genes}->{$media} = $self->figmodel()->database()->load_single_column_file($self->directory()."EssentialGenes-".$self->id().$self->selected_version()."-".$media.".tbl","");  
                 }  
         }  
   
         return $self->{_essential_genes}->{$media};  
748  }  }
749    
750  =head3 compound_table  =head3 compound_table
# Line 649  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 769  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 909  Line 1053 
1053          return $self->{_data}->cellwalltype();          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    sub biomassReaction {
1065            my ($self,$newBiomass) = @_;
1066            if (!defined($newBiomass)) {
1067                    return $self->{_data}->biomassReaction();
1068            } 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            }
1092    }
1093    
1094  =head3 taxonomy  =head3 taxonomy
1095  Definition:  Definition:
1096          string = FIGMODELmodel->taxonomy();          string = FIGMODELmodel->taxonomy();
# Line 992  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 {          } else {
1189                  $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not identify essential reactions for model ".$self->id().$self->selected_version().".");                  $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();                  return $self->figmodel()->fail();
# Line 1028  Line 1213 
1213          $self->{_data}->modificationDate(time());          $self->{_data}->modificationDate(time());
1214          my $version = $self->{_data}->autocompleteVersion();          my $version = $self->{_data}->autocompleteVersion();
1215          $self->{_data}->autocompleteVersion($version+1);          $self->{_data}->autocompleteVersion($version+1);
1216            $self->update_model_stats();
1217  }  }
1218    
1219  =head3 update_stats_for_build  =head3 update_stats_for_build
# Line 1041  Line 1227 
1227          $self->{_data}->modificationDate(time());          $self->{_data}->modificationDate(time());
1228          my $version = $self->{_data}->version();          my $version = $self->{_data}->version();
1229          $self->{_data}->version($version+1);          $self->{_data}->version($version+1);
1230            $self->update_model_stats();
1231  }  }
1232    
1233  =head3 update_model_stats  =head3 update_model_stats
# Line 1054  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    
# Line 1109  Line 1295 
1295          #Setting the reaction count          #Setting the reaction count
1296          $self->{_data}->reactions($rxntbl->size());          $self->{_data}->reactions($rxntbl->size());
1297          #Setting the metabolite count          #Setting the metabolite count
1298          $self->{_data}->compounds($rxntbl->size());          $self->{_data}->compounds($cpdtbl->size());
1299          #Setting the gene count          #Setting the gene count
1300          my $geneCount = @genes + @othergenes;          my $geneCount = @genes + @othergenes;
1301          $self->{_data}->associatedGenes($geneCount);          $self->{_data}->associatedGenes($geneCount);
# Line 1159  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 1175  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 1190  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                    #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));                  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"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"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 1205  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--) {
# Line 1233  Line 1420 
1420                                                  push(@{$ReactionList},$ID);                                                  push(@{$ReactionList},$ID);
1421                                          }                                          }
1422                                  }                                  }
1423                                  $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1);                                  $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"AUTOCOMPLETION",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()->database()->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 1260  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  =head3 GapGenModel
1528  Definition:  Definition:
1529          FIGMODELmodel->GapGenModel();          FIGMODELmodel->GapGenModel();
# Line 1497  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 1508  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 1523  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 1566  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 1661  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 1767  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 1785  Line 2059 
2059                  }                  }
2060          }          }
2061    
2062          #Checking if a biomass reaction already exists          #Creating biomass reaction for model
2063          my $BiomassReactionRow = $self->BuildSpecificBiomassReaction();          my $biomassID = $self->BuildSpecificBiomassReaction();
2064          if (!defined($BiomassReactionRow)) {          if ($biomassID !~ m/bio\d\d\d\d\d/) {
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 1823  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 1882  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()->database()->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 1906  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 1931  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 1959  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 2387  Line 2723 
2723          Calculating growth in the input media          Calculating growth in the input media
2724  =cut  =cut
2725  sub calculate_growth {  sub calculate_growth {
2726          my ($self,$Media) = @_;          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();          my $UniqueFilename = $self->figmodel()->filename();
2742          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()));          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);          my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename);
2755          my $Result;          my $Result;
2756          if (defined($ProblemReport)) {          if (defined($ProblemReport)) {
2757                  my $Row = $ProblemReport->get_row(0);                  my $Row = $ProblemReport->get_row(0);
2758                  if (defined($Row) && defined($Row->{"Objective"}->[0])) {                  if (defined($Row) && defined($Row->{"Objective"}->[0])) {
2759                          if ($Row->{"Objective"}->[0] < 0.00000001) {                          if ($Row->{"Objective"}->[0] < 0.00000001 || $Row->{"Objective"}->[0] == 1e7) {
2760                                  $Result = "NOGROWTH:".$Row->{"Individual metabolites with zero production"}->[0];                                  $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 {                          } 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];                                  $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;          return $Result;
2780  }  }
2781    
# Line 3413  Line 3787 
3787          return 1;          return 1;
3788  }  }
3789    
3790  =head3 BuildSpecificBiomassReaction  =head3 DetermineCofactorLipidCellWallComponents
3791  Definition:  Definition:
3792          FIGMODELmodel->BuildSpecificBiomassReaction();          {cofactor=>{string:compound id=>float:coefficient},lipid=>...cellWall=>} = FIGMODELmodel->DetermineCofactorLipidCellWallComponents();
3793  Description:  Description:
3794  =cut  =cut
3795  sub BuildSpecificBiomassReaction {  sub DetermineCofactorLipidCellWallComponents {
3796          my ($self) = @_;          my ($self) = @_;
3797            my $templateResults;
         my $biomassrxn;  
         my $OrganismID = $self->genome();  
         #Checking for a biomass override  
         if (defined($self->config("biomass reaction override")->{$OrganismID})) {  
                 my $biomassID = $self->config("biomass reaction override")->{$OrganismID};  
                 my $tbl = $self->database()->get_table("BIOMASS",1);  
                 $biomassrxn = $tbl->get_row_by_key($biomassID,"DATABASE");  
                 print "Overriding biomass template and selecting ".$biomassID." for ".$OrganismID.".\n";  
         } else {#Creating biomass reaction from the template  
                 #Getting the genome stats  
3798                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
3799                  my $Class = $genomestats->{CLASS}->[0];          my $Class = $self->{_data}->cellwalltype();
3800                  my $Name = $genomestats->{NAME}->[0];          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 3446  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 3463  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 3502  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 3643  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") {
4190                                    $fraction = $fraction/$netMass;
4191                            }
4192                            if ($compounds->{$types->[$i]}->{$compounds[$j]} eq 1e-5) {
4193                                    $fraction = 1;
4194                            }
4195                            $equationCompounds->{$compounds[$j]}->{"coef"} += $fraction*$compounds->{$types->[$i]}->{$compounds[$j]};
4196                    }
4197                    chop($coefString);
4198                    if (length($coefString) == 0) {
4199                            $coefString = "NONE";
4200                    }
4201                    my $function = $types->[$i]."Coef";
4202                    if ($types->[$i] ne "energy") {
4203                            $bioObj->$function($coefString);
4204                    }
4205            }
4206            #Adding biomass to compound list
4207            $equationCompounds->{cpd11416}->{coef} = 1;
4208            $equationCompounds->{cpd11416}->{type} = "macromolecule";
4209            #Building equation from hash and populating compound biomass table
4210            my @compoundList = keys(%{$equationCompounds});
4211            my ($reactants,$products);
4212            #Deleting existing Biomass Compound info
4213            my $cpdbofMgr = $self->figmodel()->database()->get_object_manager("cpdbof");
4214            my $matchingObjs = $cpdbofMgr->get_objects({BIOMASS=>$biomassID});
4215            for (my $i=0; $i < @{$matchingObjs}; $i++) {
4216                    $matchingObjs->[$i]->delete();
4217            }
4218            my $typeCategories = {"macromolecule"=>"M","RNA"=>"R","DNA"=>"D","protein"=>"P","lipid"=>"L","cellWall"=>"W","cofactor"=>"C","energy"=>"E"};
4219            my $productmass = 0;
4220            my $reactantmass = 0;
4221            my $totalmass = 0;
4222            foreach my $compound (@compoundList) {
4223                    if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4224                            $totalmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4225                    }
4226                    if ($equationCompounds->{$compound}->{coef} < 0) {
4227                            if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4228                                    $reactantmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4229                            }
4230                            $reactants->{$compound} = $self->figmodel()->format_coefficient(-1*$equationCompounds->{$compound}->{coef});
4231                    } else {
4232                            if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4233                                    $productmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4234                            }
4235                            $products->{$compound} = $self->figmodel()->format_coefficient($equationCompounds->{$compound}->{coef});
4236                    }
4237                    #Adding biomass reaction compounds to the biomass compound table
4238                    $cpdbofMgr = $self->figmodel()->database()->get_object_manager("cpdbof");
4239                    $cpdbofMgr->create({COMPOUND=>$compound,BIOMASS=>$biomassID,coefficient=>$equationCompounds->{$compound}->{coef},compartment=>"c",category=>$typeCategories->{$equationCompounds->{$compound}->{type}}});
4240                  }                  }
4241            print "Total mass = ".$totalmass.", Reactant mass = ".$reactantmass.", Product mass = ".$productmass."\n";
4242                  my $Equation = "";                  my $Equation = "";
4243                  my @ReactantList = sort(keys(%Reactants));          my @ReactantList = sort(keys(%{$reactants}));
4244                  for (my $i=0; $i < @ReactantList; $i++) {                  for (my $i=0; $i < @ReactantList; $i++) {
4245                          if (length($Equation) > 0) {                          if (length($Equation) > 0) {
4246                                  $Equation .= " + ";                                  $Equation .= " + ";
4247                          }                          }
4248                          $Equation .= $self->figmodel()->format_coefficient($Reactants{$ReactantList[$i]})." ".$ReactantList[$i];                  $Equation .= "(".$reactants->{$ReactantList[$i]}.") ".$ReactantList[$i];
4249                  }                  }
4250                  $Equation .= " => ";                  $Equation .= " => ";
4251                  my $First = 1;                  my $First = 1;
4252                  @ReactantList = sort(keys(%Products));          @ReactantList = sort(keys(%{$products}));
4253                  for (my $i=0; $i < @ReactantList; $i++) {                  for (my $i=0; $i < @ReactantList; $i++) {
4254                          if ($First == 0) {                          if ($First == 0) {
4255                                  $Equation .= " + ";                                  $Equation .= " + ";
4256                          }                          }
4257                          $First = 0;                          $First = 0;
4258                          $Equation .= $self->figmodel()->format_coefficient($Products{$ReactantList[$i]})." ".$ReactantList[$i];                  $Equation .= "(".$products->{$ReactantList[$i]}.") ".$ReactantList[$i];
4259            }
4260            $bioObj->equation($Equation);
4261            #Setting the biomass reaction of this model
4262            $self->biomassReaction($biomassID);
4263            $self->figmodel()->print_biomass_reaction_file($biomassID);
4264            #Checking if the biomass reaction remained unchanged
4265            if ($originalPackages ne "" && $originalPackages eq $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage()) {
4266                    print "UNCHANGED!\n";
4267                    $bioObj->essentialRxn($originalEssReactions);
4268            } else {
4269                    #Copying essential reaction lists if the packages in this biomasses reaction exactly match those in another biomass reaction
4270                    my $matches = $bioMgr->get_objects({cofactorPackage=>$bioObj->cofactorPackage(),lipidPackage=>$bioObj->lipidPackage(),cellWallPackage=>$bioObj->cellWallPackage()});
4271                    my $matchFound = 0;
4272                    for (my $i=0; $i < @{$matches}; $i++) {
4273                            if ($matches->[$i]->id() ne $biomassID && defined($matches->[$i]->essentialRxn()) && length($matches->[$i]->essentialRxn())) {
4274                                    $bioObj->essentialRxn($matches->[$i]->essentialRxn());
4275                                    print "MATCH!\n";
4276                                    $matchFound = 1;
4277                                    last;
4278                            }
4279                    }
4280                    #Otherwise, we calculate essential reactions
4281                    if ($matchFound == 0) {
4282                            print "NOMATCH!\n";
4283                            $self->figmodel()->determine_biomass_essential_reactions($biomassID);
4284                  }                  }
   
                 #Adding the biomass equation to the biomass table  
                 my $NewRow = $self->figmodel()->add_biomass_reaction($Equation,$self->id(),"Template:".$self->genome());  
                 $biomassrxn = $NewRow;  
4285          }          }
4286          return $biomassrxn;          return $biomassID;
4287  }  }
4288    
4289  =head3 PrintSBMLFile  =head3 PrintSBMLFile
# Line 3732  Line 4294 
4294  =cut  =cut
4295  sub PrintSBMLFile {  sub PrintSBMLFile {
4296          my($self) = @_;          my($self) = @_;
   
4297          #Opening the SBML file for printing          #Opening the SBML file for printing
4298          my $Filename = $self->config("SBML files")->[0].$self->id().".xml";          my $Filename = $self->directory().$self->id().".xml";
         if ($self->owner() ne "master") {  
                 if (!-d $self->config("SBML files")->[0].$self->owner()."/") {  
                         system("mkdir ".$self->config("SBML files")->[0].$self->owner()."/");  
                 }  
                 $Filename = $self->config("SBML files")->[0].$self->owner()."/".$self->id().".xml";  
         }  
4299          if (!open (SBMLOUTPUT, ">$Filename")) {          if (!open (SBMLOUTPUT, ">$Filename")) {
4300                  return;                  return;
4301          }          }
   
4302          #Loading and parsing the model data          #Loading and parsing the model data
4303          my $ModelTable = $self->reaction_table();          my $mdlTbl = $self->reaction_table();
4304          if (!defined($ModelTable) || !defined($ModelTable->{"array"})) {          if (!defined($mdlTbl) || !defined($mdlTbl->{"array"})) {
4305                  print "Failed to load ".$self->id()."\n";                  return $self->fail();
                 return;  
4306          }          }
4307            my $rxnMgr = $self->figmodel()->database()->get_object_manager("reaction");
4308            my $cmpTbl = $self->figmodel()->database()->get_table("COMPARTMENTS");
4309            my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4310            my $bioMgr = $self->figmodel()->database()->get_object_manager("bof");
4311          #Adding intracellular metabolites that also need exchange fluxes to the exchange hash          #Adding intracellular metabolites that also need exchange fluxes to the exchange hash
4312          my $ExchangeHash = {"cpd11416" => "c"};          my $ExchangeHash = {"cpd11416" => "c"};
   
4313          my %CompartmentsPresent;          my %CompartmentsPresent;
4314          $CompartmentsPresent{"c"} = 1;          $CompartmentsPresent{"c"} = 1;
4315          my %CompoundList;          my %CompoundList;
4316          my @ReactionList;          my @ReactionList;
4317          my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS");          for (my $i=0; $i < $mdlTbl->size(); $i++) {
4318          for (my $i=0; $i < $ModelTable->size(); $i++) {                  my $rxnObj;
4319                  my $Reaction = $ModelTable->get_row($i)->{"LOAD"}->[0];                  if ($mdlTbl->get_row($i)->{"LOAD"}->[0] =~ m/rxn\d\d\d\d\d/) {
4320                  if (defined($ReactionTable->get_row_by_key($Reaction,"DATABASE")) && defined($ReactionTable->get_row_by_key($Reaction,"DATABASE")->{"EQUATION"}->[0])) {                          $rxnObj = $rxnMgr->get_objects({id=>$mdlTbl->get_row($i)->{"LOAD"}->[0]})->[0];
4321                          push(@ReactionList,$Reaction);                  } elsif ($mdlTbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {
4322                          $_ = $ReactionTable->get_row_by_key($Reaction,"DATABASE")->{"EQUATION"}->[0];                          $rxnObj = $bioMgr->get_objects({id=>$mdlTbl->get_row($i)->{"LOAD"}->[0]})->[0];
4323                    }
4324                    if (!defined($rxnObj)) {
4325                            next;
4326                    }
4327                    push(@ReactionList,$rxnObj);
4328                    $_ = $rxnObj->equation();
4329                          my @MatchArray = /(cpd\d\d\d\d\d)/g;                          my @MatchArray = /(cpd\d\d\d\d\d)/g;
4330                          for (my $j=0; $j < @MatchArray; $j++) {                          for (my $j=0; $j < @MatchArray; $j++) {
4331                                  $CompoundList{$MatchArray[$j]}->{"c"} = 1;                                  $CompoundList{$MatchArray[$j]}->{"c"} = 1;
4332                          }                          }
4333                          $_ = $ReactionTable->get_row_by_key($Reaction,"DATABASE")->{"EQUATION"}->[0];                  $_ = $rxnObj->equation();
4334                          @MatchArray = /(cpd\d\d\d\d\d\[\D\])/g;                          @MatchArray = /(cpd\d\d\d\d\d\[\D\])/g;
4335                          for (my $j=0; $j < @MatchArray; $j++) {                          for (my $j=0; $j < @MatchArray; $j++) {
4336                                  if ($MatchArray[$j] =~ m/(cpd\d\d\d\d\d)\[(\D)\]/) {                                  if ($MatchArray[$j] =~ m/(cpd\d\d\d\d\d)\[(\D)\]/) {
# Line 3778  Line 4339 
4339                                  }                                  }
4340                          }                          }
4341                  }                  }
         }  
4342    
4343          #Printing header to SBML file          #Printing header to SBML file
4344          my $ModelName = $self->id();          my $ModelName = $self->id().$self->selected_version();
4345          $ModelName =~ s/\./_/;          $ModelName =~ s/\./_/;
4346          print SBMLOUTPUT '<?xml version="1.0" encoding="UTF-8"?>'."\n";          print SBMLOUTPUT '<?xml version="1.0" encoding="UTF-8"?>'."\n";
4347      print SBMLOUTPUT '<sbml xmlns="http://www.sbml.org/sbml/level2" level="2" version="1" xmlns:html="http://www.w3.org/1999/xhtml">' . "\n";      print SBMLOUTPUT '<sbml xmlns="http://www.sbml.org/sbml/level2" level="2" version="1" xmlns:html="http://www.w3.org/1999/xhtml">' . "\n";
4348      if (defined($self->figmodel()->database()->GetDBTable("MODEL STATS")->{$self->id()}->[0]->{"Organism name"}->[0])) {          if (defined($self->name())) {
4349          print SBMLOUTPUT '<model id="'.$ModelName.'" name="'.$self->figmodel()->database()->GetDBTable("MODEL STATS")->{$self->id()}->[0]->{"Organism name"}->[0].' SEED model">'."\n";                  print SBMLOUTPUT '<model id="'.$ModelName.'" name="'.$self->name().' SEED model">'."\n";
4350      } else {      } else {
4351          print SBMLOUTPUT '<model id="'.$ModelName.'" name="'.$self->id().' SEED model">'."\n";                  print SBMLOUTPUT '<model id="'.$ModelName.'" name="'.$self->id().$self->selected_version().' SEED model">'."\n";
4352      }      }
4353    
4354          #Printing the unit data          #Printing the unit data
# Line 3805  Line 4365 
4365      #Printing compartments for SBML file      #Printing compartments for SBML file
4366      print SBMLOUTPUT '<listOfCompartments>'."\n";      print SBMLOUTPUT '<listOfCompartments>'."\n";
4367      foreach my $Compartment (keys(%CompartmentsPresent)) {      foreach my $Compartment (keys(%CompartmentsPresent)) {
4368          if (defined($self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Name"}->[0])) {                  my $row = $cmpTbl->get_row_by_key($Compartment,"Abbreviation");
4369                  my @OutsideList = split(/\//,$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Outside"}->[0]);                  if (!defined($row) && !defined($row->{"Name"}->[0])) {
4370                            next;
4371                    }
4372                    my @OutsideList = split(/\//,$row->{"Outside"}->[0]);
4373                  my $Printed = 0;                  my $Printed = 0;
4374                  foreach my $Outside (@OutsideList) {                  foreach my $Outside (@OutsideList) {
4375                                  if (defined($CompartmentsPresent{$Outside}) && defined($self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Outside}->[0]->{"Name"}->[0])) {                          if (defined($CompartmentsPresent{$Outside}) && defined($row->{"Name"}->[0])) {
4376                                  print SBMLOUTPUT '<compartment id="'.$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Name"}->[0].'" outside="'.$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Outside}->[0]->{"Name"}->[0].'"/>'."\n";                                  print SBMLOUTPUT '<compartment id="'.$row->{"Name"}->[0].'" outside="'.$row->{"Name"}->[0].'"/>'."\n";
4377                                  $Printed = 1;                                  $Printed = 1;
4378                                  last;                                  last;
4379                          }                          }
4380                  }                  }
4381                  if ($Printed eq 0) {                  if ($Printed eq 0) {
4382                          print SBMLOUTPUT '<compartment id="'.$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Name"}->[0].'"/>'."\n";                          print SBMLOUTPUT '<compartment id="'.$row->{"Name"}->[0].'"/>'."\n";
                         }  
4383          }          }
4384      }      }
4385      print SBMLOUTPUT '</listOfCompartments>'."\n";      print SBMLOUTPUT '</listOfCompartments>'."\n";
# Line 3825  Line 4387 
4387      #Printing the list of metabolites involved in the model      #Printing the list of metabolites involved in the model
4388          print SBMLOUTPUT '<listOfSpecies>'."\n";          print SBMLOUTPUT '<listOfSpecies>'."\n";
4389      foreach my $Compound (keys(%CompoundList)) {      foreach my $Compound (keys(%CompoundList)) {
4390          my $Name = $Compound;                  my $cpdObj = $cpdMgr->get_objects({id=>$Compound})->[0];
4391                    if (!defined($cpdObj)) {
4392                            next;
4393                    }
4394                  my $Formula = "";                  my $Formula = "";
4395                  if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0])) {                  if (defined($cpdObj->formula())) {
4396                          $Formula = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0];                          $Formula = $cpdObj->formula();
4397                  }                  }
4398                  if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0])) {                  my $Name = $cpdObj->name();
                         $Name = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0];  
4399                          $Name =~ s/\s/_/;                          $Name =~ s/\s/_/;
4400                          $Name .= "_".$Formula;                          $Name .= "_".$Formula;
                 }  
4401                  $Name =~ s/[<>;&\*]//;                  $Name =~ s/[<>;&\*]//;
4402                  my $Charge = 0;                  my $Charge = 0;
4403                  if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0])) {                  if (defined($cpdObj->charge())) {
4404                          $Charge = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0];                          $Charge = $cpdObj->charge();
4405                  }                  }
4406                  foreach my $Compartment (keys(%{$CompoundList{$Compound}})) {                  foreach my $Compartment (keys(%{$CompoundList{$Compound}})) {
4407                  if ($Compartment eq "e") {                  if ($Compartment eq "e") {
4408                          $ExchangeHash->{$Compound} = "e";                          $ExchangeHash->{$Compound} = "e";
4409                  }                  }
4410                  print SBMLOUTPUT '<species id="'.$Compound.'_'.$Compartment.'" name="'.$Name.'" compartment="'.$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Name"}->[0].'" charge="'.$Charge.'" boundaryCondition="false"/>'."\n";                          my $cmprow = $cmpTbl->get_row_by_key($Compartment,"Abbreviation");
4411                            print SBMLOUTPUT '<species id="'.$Compound.'_'.$Compartment.'" name="'.$Name.'" compartment="'.$cmprow->{"Name"}->[0].'" charge="'.$Charge.'" boundaryCondition="false"/>'."\n";
4412          }          }
4413      }      }
4414    
4415          #Printing the boundary species          #Printing the boundary species
4416          foreach my $Compound (keys(%{$ExchangeHash})) {          foreach my $Compound (keys(%{$ExchangeHash})) {
4417                  my $Name = $Compound;                  my $cpdObj = $cpdMgr->get_objects({id=>$Compound})->[0];
4418                    if (!defined($cpdObj)) {
4419                            next;
4420                    }
4421                  my $Formula = "";                  my $Formula = "";
4422                  if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0])) {                  if (defined($cpdObj->formula())) {
4423                          $Formula = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0];                          $Formula = $cpdObj->formula();
4424                  }                  }
4425                  if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0])) {                  my $Name = $cpdObj->name();
                         $Name = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0];  
4426                          $Name =~ s/\s/_/;                          $Name =~ s/\s/_/;
4427                          $Name .= "_".$Formula;                          $Name .= "_".$Formula;
                 }  
4428                  $Name =~ s/[<>;&\*]//;                  $Name =~ s/[<>;&\*]//;
4429                  my $Charge = 0;                  my $Charge = 0;
4430                  if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0])) {                  if (defined($cpdObj->charge())) {
4431                          $Charge = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0];                          $Charge = $cpdObj->charge();
4432                  }                  }
4433                  print SBMLOUTPUT '<species id="'.$Compound.'_b" name="'.$Name.'" compartment="Extracellular" charge="'.$Charge.'" boundaryCondition="true"/>'."\n";                  print SBMLOUTPUT '<species id="'.$Compound.'_b" name="'.$Name.'" compartment="Extracellular" charge="'.$Charge.'" boundaryCondition="true"/>'."\n";
4434          }          }
# Line 3871  Line 4437 
4437      #Printing the list of reactions involved in the model      #Printing the list of reactions involved in the model
4438          my $ObjectiveCoef;          my $ObjectiveCoef;
4439      print SBMLOUTPUT '<listOfReactions>'."\n";      print SBMLOUTPUT '<listOfReactions>'."\n";
4440          foreach my $Reaction (@ReactionList) {  
4441            foreach my $rxnObj (@ReactionList) {
4442          $ObjectiveCoef = "0.0";          $ObjectiveCoef = "0.0";
4443                  if (defined($self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"EQUATION"}->[0])) {                  my $mdlrow = $mdlTbl->get_row_by_key($rxnObj->id(),"LOAD");
4444                  if ($Reaction =~ m/^bio/) {                  if ($rxnObj->id() =~ m/^bio/) {
4445                                  $ObjectiveCoef = "1.0";                                  $ObjectiveCoef = "1.0";
4446                          }                          }
4447                          my $LowerBound = -10000;                          my $LowerBound = -10000;
4448                  my $UpperBound = 10000;                  my $UpperBound = 10000;
4449                          my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($Reaction);                  my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateDataFromEquation($rxnObj->equation());
4450                  my $Name = $Reaction;                  my $Name = $rxnObj->name();
                 if (defined($self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"NAME"}->[0])) {  
                         $Name = $self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"NAME"}->[0];  
4451                                  $Name =~ s/[<>;&]//g;                                  $Name =~ s/[<>;&]//g;
                 }  
4452                  my $Reversibility = "true";                  my $Reversibility = "true";
4453                  if (defined($ModelTable->{$Reaction}->[0]->{"DIRECTIONALITY"}->[0])) {                  if (defined($mdlrow->{"DIRECTIONALITY"}->[0])) {
4454                          if ($ModelTable->{$Reaction}->[0]->{"DIRECTIONALITY"}->[0] ne "<=>") {                          if ($mdlrow->{"DIRECTIONALITY"}->[0] ne "<=>") {
4455                                  $LowerBound = 0;                                  $LowerBound = 0;
4456                                          $Reversibility = "false";                                          $Reversibility = "false";
4457                          }                          }
4458                          if ($ModelTable->{$Reaction}->[0]->{"DIRECTIONALITY"}->[0] eq "<=") {                          if ($mdlrow->{"DIRECTIONALITY"}->[0] eq "<=") {
4459                                  my $Temp = $Products;                                  my $Temp = $Products;
4460                                  $Products = $Reactants;                                  $Products = $Reactants;
4461                                  $Reactants = $Temp;                                  $Reactants = $Temp;
4462                          }                          }
4463                  }                  }
4464                          print SBMLOUTPUT '<reaction id="'.$Reaction.'" name="'.$Name.'" reversible="'.$Reversibility.'">'."\n";                  print SBMLOUTPUT '<reaction id="'.$rxnObj->id().'" name="'.$Name.'" reversible="'.$Reversibility.'">'."\n";
4465                          print SBMLOUTPUT "<notes>\n";                          print SBMLOUTPUT "<notes>\n";
4466                          my $ECData = "";                          my $ECData = "";
4467                          if (defined($self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"ENZYME"}->[0])) {                  if ($rxnObj->id() !~ m/^bio/) {
4468                                  $ECData = $self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"ENZYME"}->[0];                          if (defined($rxnObj->enzyme())) {
4469                                    my @ecList = split(/\|/,$rxnObj->enzyme());
4470                                    if (defined($ecList[1])) {
4471                                            $ECData = $ecList[1];
4472                                    }
4473                            }
4474                          }                          }
4475                          my $SubsystemData = "";                          my $SubsystemData = "";
4476                          if (defined($ModelTable->{$Reaction}->[0]->{"SUBSYSTEM"}->[0])) {                  if (defined($mdlrow->{"SUBSYSTEM"}->[0])) {
4477                                  $SubsystemData = $ModelTable->{$Reaction}->[0]->{"SUBSYSTEM"}->[0];                          $SubsystemData = $mdlrow->{"SUBSYSTEM"}->[0];
4478                          }                          }
4479                          my $GeneAssociation = "";                          my $GeneAssociation = "";
4480                          my $ProteinAssociation = "";                          my $ProteinAssociation = "";
4481                          if (defined($ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0])) {                  if (defined($mdlrow->{"ASSOCIATED PEG"}->[0])) {
4482                                  if (@{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}} == 1 && $ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0] !~ m/\+/) {                          if (@{$mdlrow->{"ASSOCIATED PEG"}} == 1 && $mdlrow->{"ASSOCIATED PEG"}->[0] !~ m/\+/) {
4483                                          $GeneAssociation = $ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0];                                  $GeneAssociation = $mdlrow->{"ASSOCIATED PEG"}->[0];
4484                                  } else {                                  } else {
4485                                          if (@{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}} > 1) {                                  if (@{$mdlrow->{"ASSOCIATED PEG"}} > 1) {
4486                                                  $GeneAssociation = "( ";                                                  $GeneAssociation = "( ";
4487                                          }                                          }
4488                                          for (my $i=0; $i < @{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}}; $i++) {                                  for (my $i=0; $i < @{$mdlrow->{"ASSOCIATED PEG"}}; $i++) {
4489                                                  if ($i > 0) {                                                  if ($i > 0) {
4490                                                          $GeneAssociation .= " )  or  ( ";                                                          $GeneAssociation .= " )  or  ( ";
4491                                                  }                                                  }
4492                                                  $GeneAssociation .= $ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[$i];                                          $GeneAssociation .= $mdlrow->{"ASSOCIATED PEG"}->[$i];
4493                                          }                                          }
4494                                          if (@{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}} > 1) {                                  if (@{$mdlrow->{"ASSOCIATED PEG"}} > 1) {
4495                                                  $GeneAssociation .= " )";                                                  $GeneAssociation .= " )";
4496                                          }                                          }
4497                                  }                                  }
# Line 3931  Line 4500 
4500                                          $GeneAssociation = "( ".$GeneAssociation." )";                                          $GeneAssociation = "( ".$GeneAssociation." )";
4501                                  }                                  }
4502                                  $ProteinAssociation = $GeneAssociation;                                  $ProteinAssociation = $GeneAssociation;
4503                                  if (defined($self->figmodel()->database()->GetDBTable("MODEL STATS")->{$self->id()}->[0]->{"Genome ID"}->[0])) {                          if (defined($self->genome())) {
4504                                          $ProteinAssociation = $self->figmodel()->translate_gene_to_protein($ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0],$self->figmodel()->database()->GetDBTable("MODEL STATS")->{$self->id()}->[0]->{"Genome ID"}->[0]);                                  $ProteinAssociation = $self->figmodel()->translate_gene_to_protein($mdlrow->{"ASSOCIATED PEG"}->[0],$self->genome());
4505                                  }                                  }
4506                          }                          }
4507                          print SBMLOUTPUT "<html:p>GENE_ASSOCIATION:".$GeneAssociation."</html:p>\n";                          print SBMLOUTPUT "<html:p>GENE_ASSOCIATION:".$GeneAssociation."</html:p>\n";
# Line 3963  Line 4532 
4532              print SBMLOUTPUT "</kineticLaw>\n";              print SBMLOUTPUT "</kineticLaw>\n";
4533                          print SBMLOUTPUT '</reaction>'."\n";                          print SBMLOUTPUT '</reaction>'."\n";
4534                  }                  }
         }  
4535    
4536          my @ExchangeList = keys(%{$ExchangeHash});          my @ExchangeList = keys(%{$ExchangeHash});
4537          foreach my $ExCompound (@ExchangeList) {          foreach my $ExCompound (@ExchangeList) {
4538                  my $ExCompoundName = $ExCompound;                  my $cpdObj = $cpdMgr->get_objects({id=>$ExCompound})->[0];
4539                  my $Row = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->get_row_by_key($ExCompound,"DATABASE");                  if (!defined($cpdObj)) {
4540                  if (defined($Row) && defined($Row->{"NAME"}->[0])) {                          next;
                         $ExCompoundName = $Row->{"NAME"}->[0];  
                         $ExCompoundName =~ s/[<>;&]//g;  
4541                  }                  }
4542                    my $ExCompoundName = $cpdObj->name();
4543                    $ExCompoundName =~ s/[<>;&]//g;
4544                  $ObjectiveCoef = "0.0";                  $ObjectiveCoef = "0.0";
4545                  print SBMLOUTPUT '<reaction id="EX_'.$ExCompound.'_'.$ExchangeHash->{$ExCompound}.'" name="EX_'.$ExCompoundName.'_'.$ExchangeHash->{$ExCompound}.'" reversible="true">'."\n";                  print SBMLOUTPUT '<reaction id="EX_'.$ExCompound.'_'.$ExchangeHash->{$ExCompound}.'" name="EX_'.$ExCompoundName.'_'.$ExchangeHash->{$ExCompound}.'" reversible="true">'."\n";
4546                  print SBMLOUTPUT "\t".'<notes>'."\n";                  print SBMLOUTPUT "\t".'<notes>'."\n";
# Line 4025  Line 4593 
4593                  $tbl->add_row($row);                  $tbl->add_row($row);
4594          }          }
4595          $tbl->save();          $tbl->save();
         system("cp ".$tbl->filename()." ".$self->figmodel()->config("Model table download directory")->[0].$self->figmodel()->config("ModelSimpleReactionTable")->{filename_prefix}->[0]."-".$self->id().".txt");  
4596  }  }
4597    
4598  =head3 PrintModelLPFile  =head3 PrintModelLPFile
# Line 4054  Line 4621 
4621    
4622  =head3 patch_model  =head3 patch_model
4623  Definition:  Definition:
4624          FIGMODELTable:patch results = FIGMODELmodel->patch_model(FIGMODELTable:patch table);          FIGMODELmodel->patch_model([] -or- {} of patch arguments);
4625  Description:  Description:
4626  =cut  =cut
4627  sub patch_model {  sub patch_model {
4628          my ($self,$tbl) = @_;          my ($self,$arguments) = @_;
4629            if (($self->owner() eq "master" || $self->owner() eq "chenry") && $self->id() =~ m/Seed/) {
4630          #Instantiating table                  print "Model qualifies.";
4631          my $results = FIGMODELTable->new(["Reactions","New genes","Old genes","Genes","Roles","Status"],$self->directory()."PatchResults-".$self->id().$self->selected_version().".tbl",["Reaction"],"\t",";",undef);                  $self->BuildSpecificBiomassReaction();
4632          #Getting genome annotations                  my $reactionTable = $self->reaction_table();
4633          my $features = $self->figmodel()->database()->get_genome_feature_table($self->genome());                  my $row = $reactionTable->get_row_by_key("rxn05294","LOAD");
4634          #Gettubg reaction table                  if (defined($row)) {
4635          my $reactions = $self->reaction_table();                          $reactionTable->delete_row($reactionTable->row_index($row));
4636          #Checking for patched roles                  }
4637          for (my $i=0; $i < $tbl->size(); $i++) {                  $row = $reactionTable->get_row_by_key("rxn05295","LOAD");
4638                  my $row = $tbl->get_row($i);                  if (defined($row)) {
4639                  my @genes = $features->get_rows_by_key($row->{ROLE}->[0],"ROLES");                          $reactionTable->delete_row($reactionTable->row_index($row));
4640                  if (@genes > 0) {                  }
4641                          for (my $j=0; $j < @{$row->{REACTIONS}};$j++) {                  $row = $reactionTable->get_row_by_key("rxn05296","LOAD");
4642                                  my $resultrxn = $results->get_row_by_key($row->{REACTIONS}->[$j],"Reactions");                  if (defined($row)) {
4643                                  if (!defined($resultrxn)) {                          $reactionTable->delete_row($reactionTable->row_index($row));
4644                                          $resultrxn = $results->add_row({"Reactions"=>[$row->{REACTIONS}->[$j]],"Roles"=>[$row->{ROLE}->[0]]});                  }
4645                                  }                  $reactionTable->save();
4646                                  my $rxnrow = $reactions->get_row_by_key($row->{REACTIONS}->[$j],"LOAD");                  my $growth = $self->calculate_growth();
4647                                  if (defined($rxnrow) && !defined($resultrxn->{"Old genes"})) {                  if ($growth =~ m/^NOGROWTH/ ||  $growth < 1e-5) {
4648                                          $resultrxn->{"Old genes"} = $rxnrow->{"ASSOCIATED PEG"};                          $self->error_message("patch_model:model did not grow, gapfilling!");
4649                                          if ($resultrxn->{"Old genes"}->[0] !~ m/GAP|BOF|UNIVERSAL|SPONTANEOUS/) {                          $self->GapFillModel(1);
4650                                                  push(@{$resultrxn->{"Genes"}},@{$resultrxn->{"Old genes"}});                          $growth = $self->calculate_growth();
4651                                          }                          if ($growth =~ m/^NOGROWTH/ ||  $growth < 1e-5) {
4652                                  }                                  $self->error_message("patch_model:model still did not grow!");
4653                                  delete $resultrxn->{"Current gene set"};                                  print "Patch failed.\n";
4654                                  if (defined($resultrxn->{"Genes"})) {                                  return;
                                         push(@{$resultrxn->{"Current gene set"}},@{$resultrxn->{"Genes"}});  
                                 }  
                                 for (my $k=0; $k < @genes; $k++) {  
                                         if ($genes[$k]->{ID}->[0] =~ m/(peg\.\d+)/) {  
                                                 my $gene = $1;  
                                                 my $addgene = 1;  
                                                 if (defined($resultrxn->{"Old genes"})) {  
                                                         for (my $m=0; $m < @{$resultrxn->{"Old genes"}}; $m++) {  
                                                                 if ($resultrxn->{"Old genes"}->[$m] =~ m/$gene/) {  
                                                                         $addgene = 0;  
                                                                 }  
                                                         }  
                                                 }  
                                                 if ($addgene == 1) {  
                                                         push(@{$resultrxn->{"New genes"}},$gene);  
                                                         if ($row->{COMPLEX}->[0] ne "0" && defined($resultrxn->{"Current gene set"})) {  
                                                                 my $added = 0;  
                                                                 for (my $m=0; $m < @{$resultrxn->{"Current gene set"}}; $m++) {  
                                                                         if ($row->{COMPLEX}->[0] eq "1") {  
                                                                                 $resultrxn->{"Current gene set"}->[$m] = $resultrxn->{"Current gene set"}->[$m]."+".$gene;  
                                                                                 $added = 1;  
                                                                         } else {  
                                                                                 my @geneset = split(/\+/,$resultrxn->{"Current gene set"}->[$m]);  
                                                                                 for (my $n=0; $n < @geneset;$n++) {  
                                                                                         if ($self->figmodel()->colocalized_genes($geneset[$n],$gene,$self->genome()) == 1) {  
                                                                                                 $resultrxn->{"Current gene set"}->[$m] = $resultrxn->{"Current gene set"}->[$m]."+".$gene;  
                                                                                                 $added = 1;  
                                                                                                 last;  
                                                                                         }  
                                                                                 }  
                                                                         }  
                                                                 }  
                                                                 if ($added == 0) {  
                                                                         push(@{$resultrxn->{"Current gene set"}},$gene);  
                                                                 }  
                                                         } else {  
                                                                 push(@{$resultrxn->{"Current gene set"}},$gene);  
                                                         }  
                                                 }  
                                         }  
                                 }  
                                 delete $resultrxn->{"Genes"};  
                                 push(@{$resultrxn->{"Genes"}},@{$resultrxn->{"Current gene set"}});  
                         }  
                 }  
         }  
   
         #Ensuring that the old model is preserved  
         $self->ArchiveModel();  
         #Modifing the reaction list  
         for (my $i=0; $i < $results->size();$i++) {  
                 my $row = $results->get_row($i);  
                 my $rxnrow = $reactions->get_row_by_key($row->{"Reactions"}->[0],"LOAD");  
                 if (defined($rxnrow)) {  
                         $rxnrow->{"ASSOCIATED PEG"} = $row->{"Genes"};  
                 } else {  
                         $reactions->add_row({LOAD=>[$row->{"Reactions"}->[0]],DIRECTIONALITY=>[$self->figmodel()->reversibility_of_reaction($row->{"Reactions"}->[0])],COMPARTMENT=>["c"],"ASSOCIATED PEG"=>$row->{"Genes"},SUBSYSTEM=>["NONE"],CONFIDENCE=>[2],REFERENCE=>["NONE"],NOTES=>["PATCH"]});  
4655                  }                  }
4656          }          }
4657          $reactions->save();                  $self->PrintSBMLFile();
         $results->save();  
         $self->update_model_stats();  
4658          $self->PrintModelLPFile();          $self->PrintModelLPFile();
4659                    $self->PrintModelSimpleReactionTable();
4660          $self->run_default_model_predictions();          $self->run_default_model_predictions();
4661          #Returning results                  $self->update_model_stats();
4662          return $results;                  print "Patch complete.\n";
4663            }
4664  }  }
4665    
4666  =head3 translate_genes  =head3 translate_genes
# Line 4354  Line 4864 
4864          }          }
4865  }  }
4866    
4867    =pod
4868    
4869    =item * [string]:I<list of essential genes> = B<run_geneKO_slow> (string:I<media>,0/1:I<max growth>,0/1:I<save results>);
4870    
4871    =cut
4872    
4873    sub run_geneKO_slow {
4874            my ($self,$media,$maxGrowth,$save) = @_;
4875            my $output;
4876            my $UniqueFilename = $self->figmodel()->filename();
4877            if (defined($maxGrowth) && $maxGrowth == 1) {
4878                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$media,["ProductionMFA"],{"perform single KO experiments" => 1,"MFASolver" => "GLPK","Constrain objective to this fraction of the optimal value" => 0.999},"SlowGeneKO-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
4879            } else {
4880                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$media,["ProductionMFA"],{"perform single KO experiments" => 1,"MFASolver" => "GLPK","Constrain objective to this fraction of the optimal value" => 0.1},"SlowGeneKO-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
4881            }
4882            if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."DeletionStudyResults.txt") {
4883                    print "Deletion study file not found!.\n";
4884                    return undef;
4885            }
4886            my $deltbl = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."DeletionStudyResults.txt",";","|",1,["Experiment"]);
4887            for (my $i=0; $i < $deltbl->size(); $i++) {
4888                    my $row = $deltbl->get_row($i);
4889                    if ($row->{"Insilico growth"}->[0] < 0.0000001) {
4890                            push(@{$output},$row->{Experiment}->[0]);
4891                    }
4892            }
4893            if (defined($output)) {
4894                    if (defined($save) && $save == 1) {
4895                            my $tbl = $self->essentials_table();
4896                            my $row = $tbl->get_row_by_key($media,"MEDIA",1);
4897                            $row->{"ESSENTIAL GENES"} = $output;
4898                            $tbl->save();
4899                    }
4900            }
4901            return $output;
4902    }
4903    
4904    =pod
4905    
4906    =item * [string]:I<list of minimal genes> = B<run_gene_minimization> (string:I<media>,0/1:I<max growth>,0/1:I<save results>);
4907    
4908    =cut
4909    
4910    sub run_gene_minimization {
4911            my ($self,$media,$maxGrowth,$save) = @_;
4912            my $output;
4913    
4914            #Running the MFAToolkit
4915            my $UniqueFilename = $self->figmodel()->filename();
4916            if (defined($maxGrowth) && $maxGrowth == 1) {
4917                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$media,["ProductionMFA"],{"optimize organism genes" => 1,"MFASolver" => "CPLEX","Constrain objective to this fraction of the optimal value" => 0.999},"MinimizeGenes-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
4918            } else {
4919                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$media,["ProductionMFA"],{"optimize organism genes" => 1,"MFASolver" => "CPLEX","Constrain objective to this fraction of the optimal value" => 0.1},"MinimizeGenes-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
4920            }
4921            my $tbl = $self->figmodel()->LoadProblemReport($UniqueFilename);
4922            if (!defined($tbl)) {
4923                    return undef;
4924            }
4925            for (my $i=0; $i < $tbl->size(); $i++) {
4926                    my $row = $tbl->get_row($i);
4927                    if ($row->{Notes}->[0] =~ m/Recursive\sMILP\sGENE_USE\soptimization/) {
4928                            my @array = split(/\|/,$row->{Notes}->[0]);
4929                            my $solution = $array[0];
4930                            $_ = $solution;
4931                            my @OriginalArray = /(peg\.\d+)/g;
4932                            push(@{$output},@OriginalArray);
4933                            last;
4934                    }
4935            }
4936    
4937            if (defined($output)) {
4938                    if (defined($save) && $save == 1) {
4939                            my $tbl = $self->load_model_table("MinimalGenes");
4940                            my $row = $tbl->get_table_by_key("MEDIA",$media)->get_row_by_key("MAXGROWTH",$maxGrowth);
4941                            if (defined($row)) {
4942                                    $row->{GENES} = $output;
4943                            } else {
4944                                    $tbl->add_row({GENES => $output,MEDIA => [$media],MAXGROWTH => [$maxGrowth]});
4945                            }
4946                            $tbl->save();
4947                    }
4948            }
4949            return $output;
4950    }
4951    
4952    =pod
4953    
4954    =item * [string]:I<list of inactive genes> = B<identify_inactive_genes> (string:I<media>,0/1:I<max growth>,0/1:I<save results>);
4955    
4956    =cut
4957    
4958    sub identify_inactive_genes {
4959            my ($self,$media,$maxGrowth,$save) = @_;
4960            my $output;
4961            #Running the MFAToolkit
4962            my $UniqueFilename = $self->figmodel()->filename();
4963            if (defined($maxGrowth) && $maxGrowth == 1) {
4964                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$media,["ProductionMFA"],{"find tight bounds" => 1,"MFASolver" => "GLPK","Constrain objective to this fraction of the optimal value" => 0.999},"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
4965            } else {
4966                    system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$media,["ProductionMFA"],{"find tight bounds" => 1,"MFASolver" => "GLPK","Constrain objective to this fraction of the optimal value" => 0.1},"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version()));
4967            }
4968            #Reading in the output bounds file
4969            my $ReactionTB;
4970            if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") {
4971                    $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]);
4972            }
4973            if (!defined($ReactionTB)) {
4974                    print STDERR "FIGMODEL:ClassifyModelReactions: Classification file not found when classifying reactions in ".$self->id().$self->selected_version()." with ".$media." media. Most likely the model did not grow.\n";
4975                    return undef;
4976            }
4977            #Clearing output
4978            $self->figmodel()->clearing_output($UniqueFilename,"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log");
4979            my $geneHash;
4980            my $activeGeneHash;
4981            for (my $i=0; $i < $ReactionTB->size(); $i++) {
4982                    my $Row = $ReactionTB->get_row($i);
4983                    if (defined($Row->{"Min FLUX"}) && defined($Row->{"Max FLUX"}) && defined($Row->{"DATABASE ID"}) && $Row->{"DATABASE ID"}->[0] =~ m/rxn\d\d\d\d\d/) {
4984                            my $data = $self->get_reaction_data($Row->{"DATABASE ID"}->[0]);
4985                            if (defined($data->{"ASSOCIATED PEG"})) {
4986                                    my $active = 0;
4987                                    if ($Row->{"Min FLUX"}->[0] > 0.00000001 || $Row->{"Max FLUX"}->[0] < -0.00000001 || ($Row->{"Max FLUX"}->[0]-$Row->{"Min FLUX"}->[0]) > 0.00000001) {
4988                                            $active = 1;
4989                                    }
4990                                    for (my $j=0; $j < @{$data->{"ASSOCIATED PEG"}}; $j++) {
4991                                            $_ = $data->{"ASSOCIATED PEG"}->[$j];
4992                                            my @OriginalArray = /(peg\.\d+)/g;
4993                                            for (my $k=0; $k < @OriginalArray; $k++) {
4994                                                    if ($active == 1) {
4995                                                            $activeGeneHash->{$OriginalArray[$k]} = 1;
4996                                                    }
4997                                                    $geneHash->{$OriginalArray[$k]} = 1;
4998                                            }
4999                                    }
5000                            }
5001                    }
5002            }
5003            my @allGenes = keys(%{$geneHash});
5004            for (my $i=0; $i < @allGenes; $i++) {
5005                    if (!defined($activeGeneHash->{$allGenes[$i]})) {
5006                            push(@{$output},$allGenes[$i]);
5007                    }
5008            }
5009            if (defined($output)) {
5010                    if (defined($save) && $save == 1) {
5011                            my $tbl = $self->load_model_table("InactiveGenes");
5012                            my $row = $tbl->get_table_by_key("MEDIA",$media)->get_row_by_key("MAXGROWTH",$maxGrowth);
5013                            if (defined($row)) {
5014                                    $row->{GENES} = $output;
5015                            } else {
5016                                    $tbl->add_row({GENES => $output,MEDIA => [$media],MAXGROWTH => [$maxGrowth]});
5017                            }
5018                            $tbl->save();
5019                    }
5020            }
5021            return $output;
5022    }
5023    
5024    sub ConvertVersionsToHistoryFile {
5025            my ($self) = @_;
5026            my $vone = 0;
5027            my $vtwo = 0;
5028            my $continue = 1;
5029            my $lastTable;
5030            my $currentTable;
5031            my $cause;
5032            my $lastChanged = 0;
5033            my $noHitCount = 0;
5034            while ($continue == 1) {
5035                    $cause = "NONE";
5036                    $currentTable = undef;
5037                    if (-e $self->directory().$self->id()."V".($vone+1).".".$vtwo.".txt") {
5038                            $noHitCount = 0;
5039                            $lastChanged = 0;
5040                            $vone = $vone+1;
5041                            $currentTable = $self->figmodel()->database()->load_table($self->directory().$self->id()."V".$vone.".".$vtwo.".txt",";","|",1,["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG"]);
5042                            $cause = "RECONSTRUCTION";
5043                    } elsif (-e $self->directory().$self->id()."V".$vone.".".($vtwo+1).".txt") {
5044                            $noHitCount = 0;
5045                            $lastChanged = 0;
5046                            $vtwo = $vtwo+1;
5047                            $currentTable = $self->figmodel()->database()->load_table($self->directory().$self->id()."V".$vone.".".$vtwo.".txt",";","|",1,["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG"]);
5048                            $cause = "AUTOCOMPLETION";
5049                    } elsif ($lastChanged == 0) {
5050                            $lastChanged = 1;
5051                            $vone = $vone+1;
5052                            $cause = "RECONSTRUCTION";
5053                    } elsif ($lastChanged == 1) {
5054                            $lastChanged = 2;
5055                            $vone = $vone-1;
5056                            $vtwo = $vtwo+1;
5057                            $cause = "AUTOCOMPLETION";
5058                    } elsif ($lastChanged == 2) {
5059                            $lastChanged = 0;
5060                            $vone = $vone+1;
5061                            $cause = "RECONSTRUCTION";
5062                    }
5063                    if (defined($currentTable)) {
5064                            if (defined($lastTable)) {
5065                                    print $cause."\t".$self->directory().$self->id()."V".$vone.".".$vtwo.".txt\n";
5066                                    $self->calculate_model_changes($lastTable,$cause,$currentTable,"V".$vone.".".$vtwo);
5067                            }
5068                            $lastTable = $currentTable;
5069                    } else {
5070                            $noHitCount++;
5071                            if ($noHitCount >= 40) {
5072                                    last;
5073                            }
5074                    }
5075            }
5076    }
5077    
5078  1;  1;

Legend:
Removed from v.1.18  
changed lines
  Added in v.1.28

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3