[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.21, Sat Jul 10 22:43:52 2010 UTC revision 1.27, Fri Aug 20 14:35:38 2010 UTC
# Line 100  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 191  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 231  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 284  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 331  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 548  Line 569 
569  =cut  =cut
570  sub reaction_table {  sub reaction_table {
571          my ($self,$clear) = @_;          my ($self,$clear) = @_;
572          if (defined($self->{_reaction_table})) {          if (defined($clear) && $clear == 1) {
573                  return $self->{_reaction_table};                  delete $self->{_reaction_table};
574          }          }
575            if (!defined($self->{_reaction_table})) {
576          $self->{_reaction_table} = $self->load_model_table("ModelReactions",$clear);          $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)) {          if (defined($classTbl)) {
# Line 582  Line 604 
604                          }                          }
605                  }                  }
606          }          }
607            }
608          return $self->{_reaction_table};          return $self->{_reaction_table};
609  }  }
610    
# Line 1031  Line 1054 
1054  }  }
1055    
1056  sub autocompleteMedia {  sub autocompleteMedia {
1057          my ($self) = @_;          my ($self,$newMedia) = @_;
1058            if (defined($newMedia)) {
1059                    return $self->{_data}->autoCompleteMedia($newMedia);
1060            }
1061          return $self->{_data}->autoCompleteMedia();          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 1157  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 1183  Line 1240 
1240          #Getting reaction table          #Getting reaction table
1241          my $rxntbl = $self->reaction_table();          my $rxntbl = $self->reaction_table();
1242          if (!defined($rxntbl)) {          if (!defined($rxntbl)) {
1243                  $self->figmodel()->error_message("FIGMODELmodel:update_model_stats:Could not load reaction list for ".$self->id());                  die $self->error_message("update_model_stats:Could not load reaction list!");
                 return undef;  
1244          }          }
1245          my $cpdtbl = $self->compound_table();          my $cpdtbl = $self->compound_table();
1246    
# Line 1306  Line 1362 
1362          }          }
1363    
1364          #Calling the MFAToolkit to run the gap filling optimization          #Calling the MFAToolkit to run the gap filling optimization
1365          my $MinimalMediaTable = $self->figmodel()->database()->GetDBTable("MINIMAL MEDIA TABLE");          my $Media = $self->autocompleteMedia();
1366          my $Media = "Complete";          if ($Media eq "Complete") {
1367          if (defined($MinimalMediaTable->get_row_by_key($self->genome(),"Organism"))) {                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1368                  $Media = $MinimalMediaTable->get_row_by_key($self->genome(),"Organism")->{"Minimal media"}->[0];          } else {
1369                  #Loading media, changing bounds, saving media as a test media                  #Loading media, changing bounds, saving media as a test media
1370                  my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);                  my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]);
1371                  for (my $i=0; $i < $MediaTable->size(); $i++) {                  for (my $i=0; $i < $MediaTable->size(); $i++) {
# Line 1324  Line 1380 
1380                  #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";                  #print $self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0,"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef)."\n";
1381                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0,"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));                  system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$UniqueFilename."TestMedia",["GapFilling"],{"Default max drain flux" => 0,"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));
1382                  unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");                  unlink($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt");
         } else {  
                 #print $self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef)."\n";  
                 system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),undef,["GapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0]},"GapFill".$self->id().".log",undef));  
1383          }          }
1384    
1385          #Parse the solutions file for the model and get the reaction list from it          #Parse the solutions file for the model and get the reaction list from it
# Line 1339  Line 1392 
1392                  $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");                  $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!");
1393                  return $self->fail();                  return $self->fail();
1394          }          }
         $SolutionData->save("/home/chenry/Solution.txt");  
1395    
1396          #Looking for the last printed recursive MILP solution          #Looking for the last printed recursive MILP solution
1397          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {          for (my $i=($SolutionData->size()-1); $i >=0; $i--) {
# Line 2006  Line 2058 
2058                  }                  }
2059          }          }
2060    
2061          #Checking if a biomass reaction already exists          #Creating biomass reaction for model
2062          my $BiomassReactionRow = $self->BuildSpecificBiomassReaction();          my $biomassID = $self->BuildSpecificBiomassReaction();
2063          if (!defined($BiomassReactionRow)) {          if ($biomassID !~ m/bio\d\d\d\d\d/) {
2064                  $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");                  $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction");
2065                  $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!");                  die $self->error_message("CreateSingleGenomeReactionList: Could not generate biomass reaction!");
2066                  return $self->fail();          }
2067            #Getting the biomass reaction PPO object
2068            my $bioRxn = $self->figmodel()->database()->get_object("bof",{id=>$biomassID});
2069            if (!defined($bioRxn)) {
2070                    die $self->error_message("CreateSingleGenomeReactionList: Could not find biomass reaction ".$biomassID."!");
2071            }
2072            #Getting the list of essential reactions for biomass reaction
2073            my $ReactionList;
2074            my $essentialReactions = $bioRxn->essentialRxn();
2075            if (defined($essentialReactions) && $essentialReactions =~ m/rxn\d\d\d\d\d/) {
2076                    push(@{$ReactionList},split(/\|/,$essentialReactions));
2077                    if ($essentialReactions !~ m/$biomassID/) {
2078                            push(@{$ReactionList},$biomassID);
2079                    }
2080          }          }
         my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"};  
         push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]);  
2081    
2082          #Adding biomass reactions to the model table          #Adding biomass reactions to the model table
2083          foreach my $BOFReaction (@{$ReactionList}) {          foreach my $BOFReaction (@{$ReactionList}) {
# Line 3721  Line 3784 
3784          return 1;          return 1;
3785  }  }
3786    
3787  =head3 BuildSpecificBiomassReaction  =head3 DetermineCofactorLipidCellWallComponents
3788  Definition:  Definition:
3789          FIGMODELmodel->BuildSpecificBiomassReaction();          {cofactor=>{string:compound id=>float:coefficient},lipid=>...cellWall=>} = FIGMODELmodel->DetermineCofactorLipidCellWallComponents();
3790  Description:  Description:
3791  =cut  =cut
3792  sub BuildSpecificBiomassReaction {  sub DetermineCofactorLipidCellWallComponents {
3793          my ($self) = @_;          my ($self) = @_;
3794            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  
3795                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());                  my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
3796                  my $Class = $genomestats->{CLASS}->[0];          my $Class = $self->{_data}->cellwalltype();
3797                  my $Name = $genomestats->{NAME}->[0];          my $Name = $self->name();
3798            my $translation = {COFACTOR=>"cofactor",LIPIDS=>"lipid","CELL WALL"=>"cellWall"};
3799                  #Checking for phoenix variants                  #Checking for phoenix variants
3800                  my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");                  my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table");
3801                  my $Phoenix = 0;                  my $Phoenix = 0;
3802                  my @Rows = $PhoenixVariantTable->get_rows_by_key($OrganismID,"GENOME");          my @Rows = $PhoenixVariantTable->get_rows_by_key($self->genome(),"GENOME");
3803                  my $VariantHash;                  my $VariantHash;
3804                  for (my $i=0; $i < @Rows; $i++) {                  for (my $i=0; $i < @Rows; $i++) {
3805                          $Phoenix = 1;                          $Phoenix = 1;
# Line 3754  Line 3807 
3807                                  $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];                                  $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0];
3808                          }                          }
3809                  }                  }
   
3810                  #Collecting genome data                  #Collecting genome data
3811                  my $RoleHash;                  my $RoleHash;
3812                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());                  my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome());
# Line 3771  Line 3823 
3823                                  }                                  }
3824                          }                          }
3825                  }                  }
   
3826                  #Scanning through the template item by item and determinine which biomass components should be added                  #Scanning through the template item by item and determinine which biomass components should be added
3827                  my $ComponentTypes;          my $includedHash;
3828                  my $EquationData;          my $BiomassReactionTemplateTable = $self->figmodel()->database()->get_table("BIOMASSTEMPLATE");
                 my $BiomassReactionTemplateTable = $self->figmodel()->database()->GetDBTable("BIOMASS TEMPLATE");  
3829                  for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {                  for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) {
3830                          my $Row = $BiomassReactionTemplateTable->get_row($i);                          my $Row = $BiomassReactionTemplateTable->get_row($i);
3831                    if (defined($translation->{$Row->{CLASS}->[0]})) {
3832                            my $coef = -1;
3833                            if ($Row->{"REACTANT"}->[0] eq "NO") {
3834                                    $coef = 1;
3835                                    if ($Row->{"COEFFICIENT"}->[0] =~ m/cpd/) {
3836                                            $coef = $Row->{"COEFFICIENT"}->[0];
3837                                    }
3838                            }
3839                          if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {                          if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") {
3840                                  push(@{$EquationData},$Row);                                  $includedHash->{$Row->{"ID"}->[0]} = 1;
3841                                  $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;                                  $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
3842                          } else {                          } elsif (defined($Row->{"INCLUSION CRITERIA"}->[0])) {
3843                                  my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];                                  my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0];
3844                                  my $End = 0;                                  my $End = 0;
3845                                  while ($End == 0) {                                  while ($End == 0) {
3846                                          if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {                                          if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) {
3847                                                  print $Criteria."\n";                                                  print $Criteria." : ";
3848                                                  my $Start = "";                                                  my $Start = "";
3849                                                  my $End = "";                                                  my $End = "";
3850                                                  my $Condition = $1;                                                  my $Condition = $1;
# Line 3810  Line 3868 
3868                                                                  $Result = "NO";                                                                  $Result = "NO";
3869                                                                  last;                                                                  last;
3870                                                          } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {                                                          } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) {
3871                                                                  my $Match = 0;                                                                  if (defined($includedHash->{$1}) && $Condition eq "OR") {
                                                                 for (my $k=0; $k < @{$EquationData}; $k++) {  
                                                                         if ($EquationData->[$k]->{"ID"}->[0] eq $1) {  
                                                                                 $Match = 1;  
                                                                                 last;  
                                                                         }  
                                                                 }  
                                                                 if ($Match == 1 && $Condition eq "OR") {  
3872                                                                          $Result = "YES";                                                                          $Result = "YES";
3873                                                                          last;                                                                          last;
3874                                                                  } elsif ($Match != 1 && $Condition eq "AND") {                                                                  } elsif (!defined($includedHash->{$1}) && $Condition eq "AND") {
3875                                                                          $Result = "NO";                                                                          $Result = "NO";
3876                                                                          last;                                                                          last;
3877                                                                  }                                                                  }
3878                                                          } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {                                                          } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) {
3879                                                                  my $Match = 0;                                                                  if (!defined($includedHash->{$1}) && $Condition eq "OR") {
                                                                 for (my $k=0; $k < @{$EquationData}; $k++) {  
                                                                         if ($EquationData->[$k]->{"ID"}->[0] eq $1) {  
                                                                                 $Match = 1;  
                                                                                 last;  
                                                                         }  
                                                                 }  
                                                                 if ($Match != 1 && $Condition eq "OR") {  
3880                                                                          $Result = "YES";                                                                          $Result = "YES";
3881                                                                          last;                                                                          last;
3882                                                                  } elsif ($Match == 1 && $Condition eq "AND") {                                                                  } elsif (defined($includedHash->{$1}) && $Condition eq "AND") {
3883                                                                          $Result = "NO";                                                                          $Result = "NO";
3884                                                                          last;                                                                          last;
3885                                                                  }                                                                  }
# Line 3951  Line 3995 
3995                                          }                                          }
3996                                  }                                  }
3997                                  if ($Criteria eq "YES") {                                  if ($Criteria eq "YES") {
3998                                          push(@{$EquationData},$Row);                                          $templateResults->{$translation->{$Row->{CLASS}->[0]}}->{$Row->{"ID"}->[0]} = $coef;
3999                                          $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1;                                          $includedHash->{$Row->{"ID"}->[0]} = 1;
4000                                  }                                  }
4001                          }                          }
4002                  }                  }
4003            }
4004                  #Building biomass equation          my $types = ["cofactor","lipid","cellWall"];
4005                  my %Reactants;          my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4006                  my %Products;          for (my $i=0; $i < @{$types}; $i++) {
4007                  foreach my $EquationRow (@{$EquationData}) {                  my @list =      keys(%{$templateResults->{$types->[$i]}});
4008                          #First determine what the coefficient should be                  my $entries = 0;
4009                          my $Coefficient;                  for (my $j=0; $j < @list; $j++) {
4010                          if (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^[0-9\.]+$/) {                          if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4011                                  $Coefficient = $EquationRow->{"COEFFICIENT"}->[0];                                  my $objs = $cpdMgr->get_objects({id=>$list[$j]});
4012                          } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/cpd\d\d\d\d\d/) {                                  if (!defined($objs->[0]) || $objs->[0]->mass() == 0) {
4013                                  $Coefficient = 0;                                          $templateResults->{$types->[$i]}->{$list[$j]} = -1e-5;
                                 my @CompoundList = split(/,/,$EquationRow->{"COEFFICIENT"}->[0]);  
                                 foreach my $Compound (@CompoundList) {  
                                         if (defined($Reactants{$Compound})) {  
                                                 $Coefficient += $Reactants{$Compound};  
                                         }  
                                 }  
                         } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^([0-9\.]+)\/(.+)$/) {  
                                 my @Keys = keys(%{$ComponentTypes->{$2}});  
                                 my $MW = 1;  
                                 my $CompoundData = $self->figmodel()->LoadObject($EquationRow->{"ID"}->[0]);  
                                 if (defined($CompoundData->{"MASS"})) {  
                                         $MW = $CompoundData->{"MASS"}->[0];  
                                 }  
                                 $Coefficient = $1/@Keys/$MW;  
                         }  
                         if (defined($EquationRow->{"REACTANT"}) && $EquationRow->{"REACTANT"}->[0] eq "YES") {  
                                 if (defined($Reactants{$EquationRow->{"ID"}->[0]})) {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} += $Coefficient;  
                                 } elsif (defined($Products{$EquationRow->{"ID"}->[0]}) && $Products{$EquationRow->{"ID"}->[0]} > $Coefficient) {  
                                         $Products{$EquationRow->{"ID"}->[0]} -= $Coefficient;  
                                 } elsif (defined($Products{$EquationRow->{"ID"}->[0]}) && $Products{$EquationRow->{"ID"}->[0]} < $Coefficient) {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} = $Coefficient - $Products{$EquationRow->{"ID"}->[0]};  
                                         delete $Products{$EquationRow->{"ID"}->[0]};  
                                 } else {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} = $Coefficient;  
                                 }  
                         } else {  
                                 if (defined($Products{$EquationRow->{"ID"}->[0]})) {  
                                         $Products{$EquationRow->{"ID"}->[0]} += $Coefficient;  
                                 } elsif (defined($Reactants{$EquationRow->{"ID"}->[0]}) && $Reactants{$EquationRow->{"ID"}->[0]} > $Coefficient) {  
                                         $Reactants{$EquationRow->{"ID"}->[0]} -= $Coefficient;  
                                 } elsif (defined($Reactants{$EquationRow->{"ID"}->[0]}) && $Reactants{$EquationRow->{"ID"}->[0]} < $Coefficient) {  
                                         $Products{$EquationRow->{"ID"}->[0]} = $Coefficient - $Reactants{$EquationRow->{"ID"}->[0]};  
                                         delete $Reactants{$EquationRow->{"ID"}->[0]};  
4014                                  } else {                                  } else {
4015                                          $Products{$EquationRow->{"ID"}->[0]} = $Coefficient;                                          $entries++;
4016                                    }
4017                            }
4018                    }
4019                    for (my $j=0; $j < @list; $j++) {
4020                            if ($templateResults->{$types->[$i]}->{$list[$j]} eq "-1") {
4021                                    $templateResults->{$types->[$i]}->{$list[$j]} = -1/$entries;
4022                            } elsif ($templateResults->{$types->[$i]}->{$list[$j]} =~ m/cpd/) {
4023                                    my $netCoef = 0;
4024                                    my @allcpd = split(/,/,$templateResults->{$types->[$i]}->{$list[$j]});
4025                                    for (my $k=0; $k < @allcpd; $k++) {
4026                                            if (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} ne "-1e-5") {
4027                                                    $netCoef += (1/$entries);
4028                                            } elsif (defined($templateResults->{$types->[$i]}->{$allcpd[$k]}) && $templateResults->{$types->[$i]}->{$allcpd[$k]} eq "-1e-5") {
4029                                                    $netCoef += 1e-5;
4030                                  }                                  }
4031                          }                          }
4032                                    $templateResults->{$types->[$i]}->{$list[$j]} = $netCoef;
4033                  }                  }
4034                    }
4035            }
4036            return $templateResults;
4037    }
4038    
4039    =head3 BuildSpecificBiomassReaction
4040    Definition:
4041            FIGMODELmodel->BuildSpecificBiomassReaction();
4042    Description:
4043    =cut
4044    sub BuildSpecificBiomassReaction {
4045            my ($self) = @_;
4046            #Getting the database handle for biomass reactions
4047            my $bioMgr = $self->figmodel()->database()->get_object_manager("bof");
4048            #Checking if the current biomass reaction appears in more than on model, if not, this biomass reaction is conserved for this model
4049            my $biomassID = $self->biomassReaction();
4050            if ($biomassID =~ m/bio\d\d\d\d\d/) {
4051                    my $mdlMgr = $self->figmodel()->database()->get_object_manager("model");
4052                    my $mdlObs = $mdlMgr->get_objects({biomassReaction=>$biomassID});
4053                    if (defined($mdlObs->[1])) {
4054                            $biomassID = "NONE";
4055                    }
4056            }
4057            #If the biomass ID is "NONE", then we create a new biomass reaction for the model
4058            my $bioObj;
4059            my $originalPackages = "";
4060            my $originalEssReactions = "";
4061            if ($biomassID eq "NONE") {
4062                    #Getting the current largest ID
4063                    $biomassID = $self->figmodel()->database()->check_out_new_id("bof");
4064                    $bioObj = $bioMgr->create({
4065                            id=>$biomassID,owner=>$self->owner(),users=>$self->users(),name=>"Biomass",equation=>"NONE",protein=>"0",
4066                            energy=>"0",DNA=>"0",RNA=>"0",lipid=>"0",cellWall=>"0",cofactor=>"0",
4067                            modificationDate=>time(),creationDate=>time(),
4068                            cofactorPackage=>"NONE",lipidPackage=>"NONE",cellWallPackage=>"NONE",
4069                            DNACoef=>"NONE",RNACoef=>"NONE",proteinCoef=>"NONE",lipidCoef=>"NONE",
4070                            cellWallCoef=>"NONE",cofactorCoef=>"NONE",essentialRxn=>"NONE"});
4071                    if (!defined($bioObj)) {
4072                            die $self->error_message("BuildSpecificBiomassReaction():Could not create new biomass reaction ".$biomassID."!");
4073                    }
4074            } else {
4075                    #Getting the biomass DB handler from the database
4076                    my $objs = $bioMgr->get_objects({id=>$biomassID});
4077                    if (!defined($objs->[0])) {
4078                            die $self->error_message("BuildSpecificBiomassReaction():Could not find biomass reaction ".$biomassID." in database!");
4079                    }
4080                    $bioObj = $objs->[0];
4081                    $bioObj->owner($self->owner());
4082                    $bioObj->users($self->users());
4083                    if (defined($bioObj->essentialRxn())) {
4084                            $originalEssReactions = $bioObj->essentialRxn();
4085                            $originalPackages = $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage();
4086                    }
4087            }
4088            #Getting genome stats
4089            my $genomestats = $self->figmodel()->get_genome_stats($self->genome());
4090            my $Class = $self->{_data}->cellwalltype();
4091            #Setting global coefficients based on cell wall type
4092            my $biomassCompounds;
4093            my $compounds;
4094            if ($Class eq "Gram positive") {
4095                    $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.323,cpd00052=>-0.199,cpd00062=>-0.215};
4096                    $compounds->{protein} = {cpd00001=>1,cpd00023=>-0.0637,cpd00033=>-0.0999,cpd00035=>-0.0653,cpd00039=>-0.0790,cpd00041=>-0.0362,cpd00051=>-0.0472,cpd00053=>-0.0637,cpd00054=>-0.0529,cpd00060=>-0.0277,cpd00065=>-0.0133,cpd00066=>-0.0430,cpd00069=>-0.0271,cpd00084=>-0.0139,cpd00107=>-0.0848,cpd00119=>-0.0200,cpd00129=>-0.0393,cpd00132=>-0.0362,cpd00156=>-0.0751,cpd00161=>-0.0456,cpd00322=>-0.0660};
4097                    $bioObj->protein("0.5284");
4098                    $bioObj->DNA("0.026");
4099                    $bioObj->RNA("0.0655");
4100                    $bioObj->lipid("0.075");
4101                    $bioObj->cellWall("0.25");
4102                    $bioObj->cofactor("0.10");
4103            } else {
4104                    $compounds->{RNA} = {cpd00002=>-0.262,cpd00012=>1,cpd00038=>-0.322,cpd00052=>-0.2,cpd00062=>-0.216};
4105                    $compounds->{protein} = {cpd00001=>1,cpd00023=>-0.0492,cpd00033=>-0.1145,cpd00035=>-0.0961,cpd00039=>-0.0641,cpd00041=>-0.0451,cpd00051=>-0.0554,cpd00053=>-0.0492,cpd00054=>-0.0403,cpd00060=>-0.0287,cpd00065=>-0.0106,cpd00066=>-0.0347,cpd00069=>-0.0258,cpd00084=>-0.0171,cpd00107=>-0.0843,cpd00119=>-0.0178,cpd00129=>-0.0414,cpd00132=>-0.0451,cpd00156=>-0.0791,cpd00161=>-0.0474,cpd00322=>-0.0543};
4106                    $bioObj->protein("0.563");
4107                    $bioObj->DNA("0.031");
4108                    $bioObj->RNA("0.21");
4109                    $bioObj->lipid("0.093");
4110                    $bioObj->cellWall("0.177");
4111                    $bioObj->cofactor("0.039");
4112            }
4113            #Setting energy coefficient for all reactions
4114            $bioObj->energy("40");
4115            $compounds->{energy} = {cpd00002=>-1,cpd00001=>-1,cpd00008=>1,cpd00009=>1,cpd00067=>1};
4116            #Setting DNA coefficients based on GC content
4117            my $gc = $self->figmodel()->get_genome_gc_content($self->genome());
4118            $compounds->{DNA} = {cpd00012=>1,cpd00115=>0.5*(1-$gc),cpd00241=>0.5*$gc,cpd00356=>0.5*$gc,cpd00357=>0.5*(1-$gc)};
4119            #Setting Lipid,cell wall,and cofactor coefficients based on biomass template
4120            my $templateResults = $self->DetermineCofactorLipidCellWallComponents();
4121            $compounds->{cofactor} = $templateResults->{cofactor};
4122            $compounds->{lipid} = $templateResults->{lipid};
4123            $compounds->{cellWall} = $templateResults->{cellWall};
4124            #Getting package number for cofactor, lipid, and cell wall
4125            my $packages;
4126            my $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4127            my $packageTypes = ["Cofactor","Lipid","CellWall"];
4128            my $translation = {"Cofactor"=>"cofactor","Lipid"=>"lipid","CellWall"=>"cellWall"};
4129            for (my $i=0; $i < @{$packageTypes}; $i++) {
4130                    my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4131                    my $function = $translation->{$packageTypes->[$i]}."Package";
4132                    if (@cpdList == 0) {
4133                            $bioObj->$function("NONE");
4134                    } else {
4135                            my $cpdgrpObs = $cpdgrpMgr->get_objects({type=>$packageTypes->[$i]."Package"});
4136                            for (my $j=0; $j < @{$cpdgrpObs}; $j++) {
4137                                    $packages->{$packageTypes->[$i]}->{$cpdgrpObs->[$j]->grouping()}->{$cpdgrpObs->[$j]->COMPOUND()} = 1;
4138                            }
4139                            my @packageList = keys(%{$packages->{$packageTypes->[$i]}});
4140                            my $packageHash;
4141                            for (my $j=0; $j < @packageList; $j++) {
4142                                    $packageHash->{join("|",sort(keys(%{$packages->{$packageTypes->[$i]}->{$packageList[$j]}})))} = $packageList[$j];
4143                            }
4144                            if (defined($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))})) {
4145                                    $bioObj->$function($packageHash->{join("|",sort(keys(%{$compounds->{$translation->{$packageTypes->[$i]}}})))});
4146                            } else {
4147                                    my $newPackageID = $self->figmodel()->database()->check_out_new_id($packageTypes->[$i]."Pkg");
4148                                    $bioObj->$function($newPackageID);
4149                                    my @cpdList = keys(%{$compounds->{$translation->{$packageTypes->[$i]}}});
4150                                    for (my $j=0; $j < @cpdList; $j++) {
4151                                            $cpdgrpMgr = $self->figmodel()->database()->get_object_manager("cpdgrp");
4152                                            $cpdgrpMgr->create({COMPOUND=>$cpdList[$j],grouping=>$newPackageID,type=>$packageTypes->[$i]."Package"});
4153                                    }
4154                            }
4155                    }
4156            }
4157            #Filling in coefficient terms in database and calculating global reaction coefficients based on classification abundancies
4158            my $equationCompounds;
4159            my $types = ["RNA","DNA","protein","lipid","cellWall","cofactor","energy"];
4160            my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4161            for (my $i=0; $i < @{$types}; $i++) {
4162                    my $coefString = "";
4163                    my @compounds = sort(keys(%{$compounds->{$types->[$i]}}));
4164                    #Building coefficient strings and determining net mass for component types
4165                    my $netMass = 0;
4166                    for (my $j=0; $j < @compounds; $j++) {
4167                            my $objs = $cpdMgr->get_objects({id=>$compounds[$j]});
4168                            my $mass = 0;
4169                            if (defined($objs->[0]) && $objs->[0]->mass() != 0) {
4170                                    $mass = $objs->[0]->mass();
4171                                    $netMass += -$compounds->{$types->[$i]}->{$compounds[$j]}*$objs->[0]->mass();
4172                            }
4173                            if (!defined($equationCompounds->{$compounds[$j]})) {
4174                                    $equationCompounds->{$compounds[$j]}->{"coef"} = 0;
4175                                    $equationCompounds->{$compounds[$j]}->{"type"} = $types->[$i];
4176                                    $equationCompounds->{$compounds[$j]}->{"mass"} = $mass;
4177                            }
4178                            $coefString .= $compounds->{$types->[$i]}->{$compounds[$j]}."|";
4179                    }
4180                    $netMass = 0.001*$netMass;
4181                    #Calculating coefficients for all component compounds
4182                    for (my $j=0; $j < @compounds; $j++) {
4183                            #Normalizing certain type coefficients by mass
4184                            my $function = $types->[$i];
4185                            my $fraction = $bioObj->$function();
4186                            if ($types->[$i] ne "energy") {
4187                                    $fraction = $fraction/$netMass;
4188                            }
4189                            if ($compounds->{$types->[$i]}->{$compounds[$j]} eq 1e-5) {
4190                                    $fraction = 1;
4191                            }
4192                            $equationCompounds->{$compounds[$j]}->{"coef"} += $fraction*$compounds->{$types->[$i]}->{$compounds[$j]};
4193                    }
4194                    chop($coefString);
4195                    if (length($coefString) == 0) {
4196                            $coefString = "NONE";
4197                    }
4198                    my $function = $types->[$i]."Coef";
4199                    if ($types->[$i] ne "energy") {
4200                            $bioObj->$function($coefString);
4201                    }
4202            }
4203            #Adding biomass to compound list
4204            $equationCompounds->{cpd11416}->{coef} = 1;
4205            $equationCompounds->{cpd11416}->{type} = "macromolecule";
4206            #Building equation from hash and populating compound biomass table
4207            my @compoundList = keys(%{$equationCompounds});
4208            my ($reactants,$products);
4209            #Deleting existing Biomass Compound info
4210            my $cpdbofMgr = $self->figmodel()->database()->get_object_manager("cpdbof");
4211            my $matchingObjs = $cpdbofMgr->get_objects({BIOMASS=>$biomassID});
4212            for (my $i=0; $i < @{$matchingObjs}; $i++) {
4213                    $matchingObjs->[$i]->delete();
4214            }
4215            my $typeCategories = {"macromolecule"=>"M","RNA"=>"R","DNA"=>"D","protein"=>"P","lipid"=>"L","cellWall"=>"W","cofactor"=>"C","energy"=>"E"};
4216            my $productmass = 0;
4217            my $reactantmass = 0;
4218            my $totalmass = 0;
4219            foreach my $compound (@compoundList) {
4220                    if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4221                            $totalmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4222                    }
4223                    if ($equationCompounds->{$compound}->{coef} < 0) {
4224                            if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4225                                    $reactantmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4226                            }
4227                            $reactants->{$compound} = $self->figmodel()->format_coefficient(-1*$equationCompounds->{$compound}->{coef});
4228                    } else {
4229                            if (defined($equationCompounds->{$compound}->{coef}) && defined($equationCompounds->{$compound}->{mass})) {
4230                                    $productmass += $equationCompounds->{$compound}->{coef}*0.001*$equationCompounds->{$compound}->{mass};
4231                            }
4232                            $products->{$compound} = $self->figmodel()->format_coefficient($equationCompounds->{$compound}->{coef});
4233                    }
4234                    #Adding biomass reaction compounds to the biomass compound table
4235                    $cpdbofMgr = $self->figmodel()->database()->get_object_manager("cpdbof");
4236                    $cpdbofMgr->create({COMPOUND=>$compound,BIOMASS=>$biomassID,coefficient=>$equationCompounds->{$compound}->{coef},compartment=>"c",category=>$typeCategories->{$equationCompounds->{$compound}->{type}}});
4237            }
4238            print "Total mass = ".$totalmass.", Reactant mass = ".$reactantmass.", Product mass = ".$productmass."\n";
4239                  my $Equation = "";                  my $Equation = "";
4240                  my @ReactantList = sort(keys(%Reactants));          my @ReactantList = sort(keys(%{$reactants}));
4241                  for (my $i=0; $i < @ReactantList; $i++) {                  for (my $i=0; $i < @ReactantList; $i++) {
4242                          if (length($Equation) > 0) {                          if (length($Equation) > 0) {
4243                                  $Equation .= " + ";                                  $Equation .= " + ";
4244                          }                          }
4245                          $Equation .= $self->figmodel()->format_coefficient($Reactants{$ReactantList[$i]})." ".$ReactantList[$i];                  $Equation .= "(".$reactants->{$ReactantList[$i]}.") ".$ReactantList[$i];
4246                  }                  }
4247                  $Equation .= " => ";                  $Equation .= " => ";
4248                  my $First = 1;                  my $First = 1;
4249                  @ReactantList = sort(keys(%Products));          @ReactantList = sort(keys(%{$products}));
4250                  for (my $i=0; $i < @ReactantList; $i++) {                  for (my $i=0; $i < @ReactantList; $i++) {
4251                          if ($First == 0) {                          if ($First == 0) {
4252                                  $Equation .= " + ";                                  $Equation .= " + ";
4253                          }                          }
4254                          $First = 0;                          $First = 0;
4255                          $Equation .= $self->figmodel()->format_coefficient($Products{$ReactantList[$i]})." ".$ReactantList[$i];                  $Equation .= "(".$products->{$ReactantList[$i]}.") ".$ReactantList[$i];
4256            }
4257            $bioObj->equation($Equation);
4258            #Setting the biomass reaction of this model
4259            $self->biomassReaction($biomassID);
4260            $self->figmodel()->print_biomass_reaction_file($biomassID);
4261            #Checking if the biomass reaction remained unchanged
4262            if ($originalPackages ne "" && $originalPackages eq $bioObj->cofactorPackage().$bioObj->lipidPackage().$bioObj->cellWallPackage()) {
4263                    print "UNCHANGED!\n";
4264                    $bioObj->essentialRxn($originalEssReactions);
4265            } else {
4266                    #Copying essential reaction lists if the packages in this biomasses reaction exactly match those in another biomass reaction
4267                    my $matches = $bioMgr->get_objects({cofactorPackage=>$bioObj->cofactorPackage(),lipidPackage=>$bioObj->lipidPackage(),cellWallPackage=>$bioObj->cellWallPackage()});
4268                    my $matchFound = 0;
4269                    for (my $i=0; $i < @{$matches}; $i++) {
4270                            if ($matches->[$i]->id() ne $biomassID && defined($matches->[$i]->essentialRxn()) && length($matches->[$i]->essentialRxn())) {
4271                                    $bioObj->essentialRxn($matches->[$i]->essentialRxn());
4272                                    print "MATCH!\n";
4273                                    $matchFound = 1;
4274                                    last;
4275                  }                  }
   
                 #Adding the biomass equation to the biomass table  
                 my $NewRow = $self->figmodel()->add_biomass_reaction($Equation,$self->id(),"Template:".$self->genome());  
                 $biomassrxn = $NewRow;  
4276          }          }
4277          return $biomassrxn;                  #Otherwise, we calculate essential reactions
4278                    if ($matchFound == 0) {
4279                            print "NOMATCH!\n";
4280                            $self->figmodel()->determine_biomass_essential_reactions($biomassID);
4281                    }
4282            }
4283            return $biomassID;
4284  }  }
4285    
4286  =head3 PrintSBMLFile  =head3 PrintSBMLFile
# Line 4040  Line 4291 
4291  =cut  =cut
4292  sub PrintSBMLFile {  sub PrintSBMLFile {
4293          my($self) = @_;          my($self) = @_;
   
4294          #Opening the SBML file for printing          #Opening the SBML file for printing
4295          my $Filename = $self->directory().$self->id().".xml";          my $Filename = $self->directory().$self->id().".xml";
4296          if (!open (SBMLOUTPUT, ">$Filename")) {          if (!open (SBMLOUTPUT, ">$Filename")) {
4297                  return;                  return;
4298          }          }
   
4299          #Loading and parsing the model data          #Loading and parsing the model data
4300          my $mdlTbl = $self->reaction_table();          my $mdlTbl = $self->reaction_table();
4301          if (!defined($mdlTbl) || !defined($mdlTbl->{"array"})) {          if (!defined($mdlTbl) || !defined($mdlTbl->{"array"})) {
4302                  return $self->fail();                  return $self->fail();
4303          }          }
4304          my $rxnTbl = $self->figmodel()->database()->get_table("REACTIONS");          my $rxnMgr = $self->figmodel()->database()->get_object_manager("reaction");
         my $bioTbl = $self->figmodel()->database()->get_table("BIOMASS");  
         my $cpdTbl = $self->figmodel()->database()->get_table("COMPOUNDS");  
4305          my $cmpTbl = $self->figmodel()->database()->get_table("COMPARTMENTS");          my $cmpTbl = $self->figmodel()->database()->get_table("COMPARTMENTS");
4306            my $cpdMgr = $self->figmodel()->database()->get_object_manager("compound");
4307            my $bioMgr = $self->figmodel()->database()->get_object_manager("bof");
4308          #Adding intracellular metabolites that also need exchange fluxes to the exchange hash          #Adding intracellular metabolites that also need exchange fluxes to the exchange hash
4309          my $ExchangeHash = {"cpd11416" => "c"};          my $ExchangeHash = {"cpd11416" => "c"};
4310          my %CompartmentsPresent;          my %CompartmentsPresent;
# Line 4064  Line 4312 
4312          my %CompoundList;          my %CompoundList;
4313          my @ReactionList;          my @ReactionList;
4314          for (my $i=0; $i < $mdlTbl->size(); $i++) {          for (my $i=0; $i < $mdlTbl->size(); $i++) {
4315                  my $Reaction = $mdlTbl->get_row($i)->{"LOAD"}->[0];                  my $rxnObj;
4316                  my $row = $rxnTbl->get_row_by_key($Reaction,"DATABASE");                  if ($mdlTbl->get_row($i)->{"LOAD"}->[0] =~ m/rxn\d\d\d\d\d/) {
4317                  if (!defined($row)) {                          $rxnObj = $rxnMgr->get_objects({id=>$mdlTbl->get_row($i)->{"LOAD"}->[0]})->[0];
4318                          $row = $bioTbl->get_row_by_key($Reaction,"DATABASE");                  } elsif ($mdlTbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) {
4319                          if (!defined($row)) {                          $rxnObj = $bioMgr->get_objects({id=>$mdlTbl->get_row($i)->{"LOAD"}->[0]})->[0];
                                 next;  
                         }  
4320                  }                  }
4321                  if (!defined($row->{"EQUATION"}->[0])) {                  if (!defined($rxnObj)) {
4322                          next;                          next;
4323                  }                  }
4324                  push(@ReactionList,$Reaction);                  push(@ReactionList,$rxnObj);
4325                  $_ = $row->{"EQUATION"}->[0];                  $_ = $rxnObj->equation();
4326                  my @MatchArray = /(cpd\d\d\d\d\d)/g;                  my @MatchArray = /(cpd\d\d\d\d\d)/g;
4327                  for (my $j=0; $j < @MatchArray; $j++) {                  for (my $j=0; $j < @MatchArray; $j++) {
4328                          $CompoundList{$MatchArray[$j]}->{"c"} = 1;                          $CompoundList{$MatchArray[$j]}->{"c"} = 1;
4329                  }                  }
4330                  $_ = $row->{"EQUATION"}->[0];                  $_ = $rxnObj->equation();
4331                  @MatchArray = /(cpd\d\d\d\d\d\[\D\])/g;                  @MatchArray = /(cpd\d\d\d\d\d\[\D\])/g;
4332                  for (my $j=0; $j < @MatchArray; $j++) {                  for (my $j=0; $j < @MatchArray; $j++) {
4333                          if ($MatchArray[$j] =~ m/(cpd\d\d\d\d\d)\[(\D)\]/) {                          if ($MatchArray[$j] =~ m/(cpd\d\d\d\d\d)\[(\D)\]/) {
# Line 4138  Line 4384 
4384          #Printing the list of metabolites involved in the model          #Printing the list of metabolites involved in the model
4385          print SBMLOUTPUT '<listOfSpecies>'."\n";          print SBMLOUTPUT '<listOfSpecies>'."\n";
4386          foreach my $Compound (keys(%CompoundList)) {          foreach my $Compound (keys(%CompoundList)) {
4387                  my $row = $cpdTbl->get_row_by_key($Compound,"DATABASE");                  my $cpdObj = $cpdMgr->get_objects({id=>$Compound})->[0];
4388                  if (!defined($row)) {                  if (!defined($cpdObj)) {
4389                          next;                          next;
4390                  }                  }
                 my $Name = $Compound;  
4391                  my $Formula = "";                  my $Formula = "";
4392                  if (defined($row->{"FORMULA"}->[0])) {                  if (defined($cpdObj->formula())) {
4393                          $Formula = $row->{"FORMULA"}->[0];                          $Formula = $cpdObj->formula();
4394                  }                  }
4395                  if (defined($row->{"NAME"}->[0])) {                  my $Name = $cpdObj->name();
                         $Name = $row->{"NAME"}->[0];  
4396                          $Name =~ s/\s/_/;                          $Name =~ s/\s/_/;
4397                          $Name .= "_".$Formula;                          $Name .= "_".$Formula;
                 }  
4398                  $Name =~ s/[<>;&\*]//;                  $Name =~ s/[<>;&\*]//;
4399                  my $Charge = 0;                  my $Charge = 0;
4400                  if (defined($row->{"CHARGE"}->[0])) {                  if (defined($cpdObj->charge())) {
4401                          $Charge = $row->{"CHARGE"}->[0];                          $Charge = $cpdObj->charge();
4402                  }                  }
4403                  foreach my $Compartment (keys(%{$CompoundList{$Compound}})) {                  foreach my $Compartment (keys(%{$CompoundList{$Compound}})) {
4404                          if ($Compartment eq "e") {                          if ($Compartment eq "e") {
# Line 4168  Line 4411 
4411    
4412          #Printing the boundary species          #Printing the boundary species
4413          foreach my $Compound (keys(%{$ExchangeHash})) {          foreach my $Compound (keys(%{$ExchangeHash})) {
4414                  my $Name = $Compound;                  my $cpdObj = $cpdMgr->get_objects({id=>$Compound})->[0];
4415                  my $Formula = "";                  if (!defined($cpdObj)) {
                 my $row = $cpdTbl->get_row_by_key($Compound,"DATABASE");  
                 if (!defined($row)) {  
4416                          next;                          next;
4417                  }                  }
4418                  if (defined($row->{"FORMULA"}->[0])) {                  my $Formula = "";
4419                          $Formula = $row->{"FORMULA"}->[0];                  if (defined($cpdObj->formula())) {
4420                            $Formula = $cpdObj->formula();
4421                  }                  }
4422                  if (defined($row->{"NAME"}->[0])) {                  my $Name = $cpdObj->name();
                         $Name = $row->{"NAME"}->[0];  
4423                          $Name =~ s/\s/_/;                          $Name =~ s/\s/_/;
4424                          $Name .= "_".$Formula;                          $Name .= "_".$Formula;
                 }  
4425                  $Name =~ s/[<>;&\*]//;                  $Name =~ s/[<>;&\*]//;
4426                  my $Charge = 0;                  my $Charge = 0;
4427                  if (defined($row->{"CHARGE"}->[0])) {                  if (defined($cpdObj->charge())) {
4428                          $Charge = $row->{"CHARGE"}->[0];                          $Charge = $cpdObj->charge();
4429                  }                  }
4430                  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";
4431          }          }
# Line 4195  Line 4435 
4435          my $ObjectiveCoef;          my $ObjectiveCoef;
4436          print SBMLOUTPUT '<listOfReactions>'."\n";          print SBMLOUTPUT '<listOfReactions>'."\n";
4437    
4438          foreach my $Reaction (@ReactionList) {          foreach my $rxnObj (@ReactionList) {
4439                  $ObjectiveCoef = "0.0";                  $ObjectiveCoef = "0.0";
4440                  my $mdlrow = $mdlTbl->get_row_by_key($Reaction,"LOAD");                  my $mdlrow = $mdlTbl->get_row_by_key($rxnObj->id(),"LOAD");
4441                  my $Reaction = $mdlrow->{"LOAD"}->[0];                  if ($rxnObj->id() =~ m/^bio/) {
                 my $row = $rxnTbl->get_row_by_key($Reaction,"DATABASE");  
                 if (!defined($row)) {  
                         $row = $bioTbl->get_row_by_key($Reaction,"DATABASE");  
                         if (!defined($row)) {  
                                 next;  
                         }  
                 }  
                 if (!defined($row->{"EQUATION"}->[0])) {  
                         next;  
                 }  
                 if ($Reaction =~ m/^bio/) {  
4442                          $ObjectiveCoef = "1.0";                          $ObjectiveCoef = "1.0";
4443                  }                  }
4444                  my $LowerBound = -10000;                  my $LowerBound = -10000;
4445                  my $UpperBound = 10000;                  my $UpperBound = 10000;
4446                  my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateDataFromEquation($row->{"EQUATION"}->[0]);                  my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateDataFromEquation($rxnObj->equation());
4447                  my $Name = $Reaction;                  my $Name = $rxnObj->name();
                 if (defined($row->{"NAME"}->[0])) {  
                         $Name = $row->{"NAME"}->[0];  
4448                          $Name =~ s/[<>;&]//g;                          $Name =~ s/[<>;&]//g;
                 }  
4449                  my $Reversibility = "true";                  my $Reversibility = "true";
4450                  if (defined($mdlrow->{"DIRECTIONALITY"}->[0])) {                  if (defined($mdlrow->{"DIRECTIONALITY"}->[0])) {
4451                          if ($mdlrow->{"DIRECTIONALITY"}->[0] ne "<=>") {                          if ($mdlrow->{"DIRECTIONALITY"}->[0] ne "<=>") {
# Line 4232  Line 4458 
4458                                  $Reactants = $Temp;                                  $Reactants = $Temp;
4459                          }                          }
4460                  }                  }
4461                  print SBMLOUTPUT '<reaction id="'.$Reaction.'" name="'.$Name.'" reversible="'.$Reversibility.'">'."\n";                  print SBMLOUTPUT '<reaction id="'.$rxnObj->id().'" name="'.$Name.'" reversible="'.$Reversibility.'">'."\n";
4462                  print SBMLOUTPUT "<notes>\n";                  print SBMLOUTPUT "<notes>\n";
4463                  my $ECData = "";                  my $ECData = "";
4464                  if (defined($row->{"ENZYME"}->[0])) {                  if ($rxnObj->id() !~ m/^bio/) {
4465                          $ECData = $row->{"ENZYME"}->[0];                          if (defined($rxnObj->enzyme())) {
4466                                    my @ecList = split(/\|/,$rxnObj->enzyme());
4467                                    if (defined($ecList[1])) {
4468                                            $ECData = $ecList[1];
4469                                    }
4470                            }
4471                  }                  }
4472                  my $SubsystemData = "";                  my $SubsystemData = "";
4473                  if (defined($mdlrow->{"SUBSYSTEM"}->[0])) {                  if (defined($mdlrow->{"SUBSYSTEM"}->[0])) {
# Line 4301  Line 4532 
4532    
4533          my @ExchangeList = keys(%{$ExchangeHash});          my @ExchangeList = keys(%{$ExchangeHash});
4534          foreach my $ExCompound (@ExchangeList) {          foreach my $ExCompound (@ExchangeList) {
4535                  my $ExCompoundName = $ExCompound;                  my $cpdObj = $cpdMgr->get_objects({id=>$ExCompound})->[0];
4536                  my $Row = $cpdTbl->get_row_by_key($ExCompound,"DATABASE");                  if (!defined($cpdObj)) {
4537                  if (defined($Row) && defined($Row->{"NAME"}->[0])) {                          next;
                         $ExCompoundName = $Row->{"NAME"}->[0];  
                         $ExCompoundName =~ s/[<>;&]//g;  
4538                  }                  }
4539                    my $ExCompoundName = $cpdObj->name();
4540                    $ExCompoundName =~ s/[<>;&]//g;
4541                  $ObjectiveCoef = "0.0";                  $ObjectiveCoef = "0.0";
4542                  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";
4543                  print SBMLOUTPUT "\t".'<notes>'."\n";                  print SBMLOUTPUT "\t".'<notes>'."\n";
# Line 4387  Line 4618 
4618    
4619  =head3 patch_model  =head3 patch_model
4620  Definition:  Definition:
4621          FIGMODELTable:patch results = FIGMODELmodel->patch_model(FIGMODELTable:patch table);          FIGMODELmodel->patch_model([] -or- {} of patch arguments);
4622  Description:  Description:
4623  =cut  =cut
4624  sub patch_model {  sub patch_model {
4625          my ($self,$tbl) = @_;          my ($self,$arguments) = @_;
4626            if (($self->owner() eq "master" || $self->owner() eq "chenry") && $self->id() =~ m/Seed/) {
4627          #Instantiating table                  print "Model qualifies.";
4628          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();
4629          #Getting genome annotations                  my $reactionTable = $self->reaction_table();
4630          my $features = $self->figmodel()->database()->get_genome_feature_table($self->genome());                  my $row = $reactionTable->get_row_by_key("rxn05294","LOAD");
4631          #Gettubg reaction table                  if (defined($row)) {
4632          my $reactions = $self->reaction_table();                          $reactionTable->delete_row($reactionTable->row_index($row));
         #Checking for patched roles  
         for (my $i=0; $i < $tbl->size(); $i++) {  
                 my $row = $tbl->get_row($i);  
                 my @genes = $features->get_rows_by_key($row->{ROLE}->[0],"ROLES");  
                 if (@genes > 0) {  
                         for (my $j=0; $j < @{$row->{REACTIONS}};$j++) {  
                                 my $resultrxn = $results->get_row_by_key($row->{REACTIONS}->[$j],"Reactions");  
                                 if (!defined($resultrxn)) {  
                                         $resultrxn = $results->add_row({"Reactions"=>[$row->{REACTIONS}->[$j]],"Roles"=>[$row->{ROLE}->[0]]});  
                                 }  
                                 my $rxnrow = $reactions->get_row_by_key($row->{REACTIONS}->[$j],"LOAD");  
                                 if (defined($rxnrow) && !defined($resultrxn->{"Old genes"})) {  
                                         $resultrxn->{"Old genes"} = $rxnrow->{"ASSOCIATED PEG"};  
                                         if ($resultrxn->{"Old genes"}->[0] !~ m/GAP|BOF|UNIVERSAL|SPONTANEOUS/) {  
                                                 push(@{$resultrxn->{"Genes"}},@{$resultrxn->{"Old genes"}});  
                                         }  
                                 }  
                                 delete $resultrxn->{"Current gene set"};  
                                 if (defined($resultrxn->{"Genes"})) {  
                                         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"}});  
4633                          }                          }
4634                    $row = $reactionTable->get_row_by_key("rxn05295","LOAD");
4635                    if (defined($row)) {
4636                            $reactionTable->delete_row($reactionTable->row_index($row));
4637                  }                  }
4638                    $row = $reactionTable->get_row_by_key("rxn05296","LOAD");
4639                    if (defined($row)) {
4640                            $reactionTable->delete_row($reactionTable->row_index($row));
4641          }          }
4642                    $reactionTable->save();
4643          #Ensuring that the old model is preserved                  my $growth = $self->calculate_growth();
4644          $self->ArchiveModel();                  if ($growth =~ m/^NOGROWTH/ ||  $growth < 1e-5) {
4645          #Modifing the reaction list                          $self->error_message("patch_model:model did not grow, gapfilling!");
4646          for (my $i=0; $i < $results->size();$i++) {                          $self->GapFillModel(1);
4647                  my $row = $results->get_row($i);                          $growth = $self->calculate_growth();
4648                  my $rxnrow = $reactions->get_row_by_key($row->{"Reactions"}->[0],"LOAD");                          if ($growth =~ m/^NOGROWTH/ ||  $growth < 1e-5) {
4649                  if (defined($rxnrow)) {                                  $self->error_message("patch_model:model still did not grow!");
4650                          $rxnrow->{"ASSOCIATED PEG"} = $row->{"Genes"};                                  print "Patch failed.\n";
4651                  } else {                                  return;
                         $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"]});  
4652                  }                  }
4653          }          }
4654          $reactions->save();                  $self->PrintSBMLFile();
         $results->save();  
         $self->update_model_stats();  
4655          $self->PrintModelLPFile();          $self->PrintModelLPFile();
4656                    $self->PrintModelSimpleReactionTable();
4657          $self->run_default_model_predictions();          $self->run_default_model_predictions();
4658          #Returning results                  $self->update_model_stats();
4659          return $results;                  print "Patch complete.\n";
4660            }
4661  }  }
4662    
4663  =head3 translate_genes  =head3 translate_genes

Legend:
Removed from v.1.21  
changed lines
  Added in v.1.27

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3