use strict; package FIGMODELmodel; =head1 FIGMODELmodel object =head2 Introduction Module for manipulating model objects. =head2 Core Object Methods =head3 new Definition: FIGMODELmodel = FIGMODELmodel->new(); Description: This is the constructor for the FIGMODELmodel object. =cut sub new { my ($class,$figmodel,$id) = @_; #Error checking first if (!defined($figmodel)) { print STDERR "FIGMODELmodel->new(undef,".$id."):figmodel must be defined to create a model object!\n"; return undef; } my $self = {_figmodel => $figmodel}; bless $self; if (!defined($id)) { $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,undef):id must be defined to create a model object"); return undef; } #Checking that the id exists my $modelHandle = $self->figmodel()->database()->get_object_manager("model"); if (!defined($modelHandle)) { $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not load MODELS table. Check database!"); return undef; } #If the id is a number, we get the model row by index if ($id =~ m/^\d+$/) { my $objects = $modelHandle->get_objects(); $self->{_data} = $objects->[$id]; $self->figmodel()->{_models}->{$id} = $self; } else { my $objects = $modelHandle->get_objects({id => $id}); if (!defined($objects->[0]) && $id =~ m/(.+)(V[^V]+)/) { $objects = $modelHandle->get_objects({id => $1}); if (!defined($objects->[0])) { $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!"); return undef; } $self->{_selectedversion} = $2; } elsif (!defined($objects->[0])) { $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find model ".$id." in database!"); return undef; } $self->{_data} = $objects->[0]; } if (!defined($self->{_data})) { $self->figmodel()->error_message("FIGMODELmodel->new(figmodel,".$id."):could not find specified id in database!"); return undef; } $self->figmodel()->{_models}->{$self->id()} = $self; return $self; } =head3 config Definition: ref::key value = FIGMODELmodel->config(string::key); Description: Trying to avoid using calls that assume configuration data is stored in a particular manner. Call this function to get file paths etc. =cut sub config { my ($self,$key) = @_; return $self->figmodel()->config($key); } =head3 fail Definition: -1 = FIGMODELmodel->fail(); Description: Standard return for failed functions. =cut sub fail { my ($self) = @_; return $self->figmodel()->fail(); } =head3 success Definition: 1 = FIGMODELmodel->success(); Description: Standard return for successful functions. =cut sub success { my ($self) = @_; return $self->figmodel()->success(); } =head3 figmodel Definition: FIGMODEL = FIGMODELmodel->figmodel(); Description: Returns a FIGMODEL object =cut sub figmodel { my ($self) = @_; return $self->{"_figmodel"}; } =head3 fig Definition: FIGMODEL = FIGMODELmodel->fig(); Description: Returns a FIG object =cut sub fig { my ($self) = @_; if (!defined($self->{_fig}) && $self->source() !~ /^MGRAST/) { if ($self->source() =~ /^RAST/) { $self->{"_fig"} = $self->figmodel()->fig();#$self->genome()); } else { $self->{"_fig"} = $self->figmodel()->fig(); } } return $self->{"_fig"}; } =head3 aquireModelLock Definition: FIGMODELmodel->aquireModelLock(); Description: Locks the database for alterations relating to the current model object =cut sub aquireModelLock { my ($self) = @_; $self->figmodel()->database()->genericLock($self->id()); } =head3 releaseModelLock Definition: FIGMODELmodel->releaseModelLock(); Description: Unlocks the database for alterations relating to the current model object =cut sub releaseModelLock { my ($self) = @_; $self->figmodel()->database()->genericUnlock($self->id()); } =head3 mgdata Definition: FIGMODEL = FIGMODELmodel->mgdata(); Description: Returns a mgrast database object =cut sub mgdata { my ($self) = @_; if (!defined($self->{_mgdata}) && $self->source() =~ /^MGRAST/) { require MGRAST; $self->{_mgdata} = $self->figmodel()->mgrast()->Job->get_objects( { 'genome_id' => $self->genome() } ) } return $self->{_mgdata}; } =head3 id Definition: string = FIGMODELmodel->id(); Description: Returns model id =cut sub id { my ($self) = @_; return $self->{_data}->id(); } =head3 owner Definition: string = FIGMODELmodel->owner(); Description: Returns the username for the model owner =cut sub owner { my ($self) = @_; return $self->{_data}->owner(); } =head3 name Definition: string = FIGMODELmodel->name(); Description: Returns the name of the organism or metagenome sample being modeled =cut sub name { my ($self) = @_; if (!defined($self->{_name})) { my $source = $self->source(); if ($source =~ /^MGRAST/) { $self->{_name} = "NA"; if (defined($self->mgdata())) { $self->{_name} = $self->mgdata()->genome_name; } } elsif ($source !~ /^RAST/) { $self->{_name} = $self->fig()->orgname_of_orgid($self->genome()); } else { $self->{_name} = $self->figmodel()->get_genome_stats($self->genome())->{NAME}->[0]; } } return $self->{_name}; } =head3 get_reaction_class Definition: string = FIGMODELmodel->get_reaction_class(string::reaction ID); Description: Returns reaction class =cut sub get_reaction_class { my ($self,$reaction,$nohtml,$brief_flux) = @_; if (!-e $self->directory()."ReactionClassification-".$self->id().".tbl") { if (!defined($self->{_reaction_classes})) { $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id()."-Complete.tbl",";","|",0,["REACTION"]); if (!defined($self->{_reaction_classes})) { return undef; } } my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION"); if (defined($ClassRow) && defined($ClassRow->{CLASS})) { my $class; my $min = $ClassRow->{MIN}->[0]; my $max = $ClassRow->{MAX}->[0]; if ($ClassRow->{CLASS}->[0] eq "Positive") { $class = "Essential =>"; $brief_flux ? $class.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $class.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[0] eq "Negative") { $class = "Essential <="; $brief_flux ? $class.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $class.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[0] eq "Positive variable") { $class = "Active =>"; $brief_flux ? $class.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $class.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[0] eq "Negative variable") { $class = "Active <="; $brief_flux ? $class.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $class.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[0] eq "Variable") { $class = "Active <=>"; $brief_flux ? $class.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $class.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[0] eq "Blocked") { $class = "Inactive"; } elsif ($ClassRow->{CLASS}->[0] eq "Dead") { $class = "Disconnected"; } if (!defined($nohtml) || $nohtml ne "1") { $class = "".$class.""; } return $class; } return undef; } if (!defined($self->{_reaction_classes})) { $self->{_reaction_classes} = $self->figmodel()->database()->load_table($self->directory()."ReactionClassification-".$self->id().".tbl",";","|",0,["REACTION"]); if (!defined($self->{_reaction_classes})) { return undef; } } my $ClassRow = $self->{_reaction_classes}->get_row_by_key($reaction,"REACTION"); my $classstring = ""; if (defined($ClassRow) && defined($ClassRow->{CLASS})) { for (my $i=0; $i < @{$ClassRow->{CLASS}};$i++) { if (length($classstring) > 0) { $classstring .= "
"; } my $NewClass; my $min = $ClassRow->{MIN}->[$i]; my $max = $ClassRow->{MAX}->[$i]; if ($ClassRow->{CLASS}->[$i] eq "Positive") { $NewClass = $ClassRow->{MEDIA}->[$i].":Essential =>"; $brief_flux ? $NewClass.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $NewClass.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[$i] eq "Negative") { $NewClass = $ClassRow->{MEDIA}->[$i].":Essential <="; $brief_flux ? $NewClass.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $NewClass.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[$i] eq "Positive variable") { $NewClass = $ClassRow->{MEDIA}->[$i].":Active =>"; $brief_flux ? $NewClass.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $NewClass.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[$i] eq "Negative variable") { $NewClass = $ClassRow->{MEDIA}->[$i].":Active <="; $brief_flux ? $NewClass.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $NewClass.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[$i] eq "Variable") { $NewClass = $ClassRow->{MEDIA}->[$i].":Active <=>"; $brief_flux ? $NewClass.="
[Flux: ".sprintf("%.3g",$min)." to ".sprintf("%.3g",$max)."]
" : $NewClass.="
[Flux: ".$min." to ".$max."]
"; } elsif ($ClassRow->{CLASS}->[$i] eq "Blocked") { $NewClass = $ClassRow->{MEDIA}->[$i].":Inactive"; } elsif ($ClassRow->{CLASS}->[$i] eq "Dead") { $NewClass = $ClassRow->{MEDIA}->[$i].":Disconnected"; } if (!defined($nohtml) || $nohtml ne "1") { $NewClass = "".$NewClass.""; } $classstring .= $NewClass; } } return $classstring; } =head3 get_biomass Definition: string = FIGMODELmodel->get_biomass(); Description: Returns data for the biomass reaction =cut sub get_biomass { my ($self) = @_; if (!defined($self->{_biomass})) { my $rxntbl = $self->reaction_table(); if (defined($rxntbl)) { for (my $i=0; $i < $rxntbl->size(); $i++) { if ($rxntbl->get_row($i)->{"LOAD"}->[0] =~ m/bio\d\d\d\d\d/) { $self->{_biomass} = $rxntbl->get_row($i)->{"LOAD"}->[0]; last; } } } } return $self->get_reaction_data($self->{_biomass}); } =head3 get_reaction_data Definition: string = FIGMODELmodel->get_reaction_data(string::reaction ID); Description: Returns model reaction data =cut sub get_reaction_data { my ($self,$reaction) = @_; if (!defined($self->reaction_table())) { return undef; } if ($reaction =~ m/^\d+$/) { $self->reaction_table()->get_row($reaction); } return $self->reaction_table()->get_row_by_key($reaction,"LOAD"); } =head3 get_reaction_id Definition: string = FIGMODELmodel->get_reaction_id(string::reaction ID); Description: Returns model reaction id =cut sub get_reaction_id { my ($self,$reaction) = @_; my $Data = $self->get_reaction_data($reaction); if (!defined($Data) || !defined($Data->{LOAD}->[0])) { return undef; } return $Data->{LOAD}->[0]; } =head3 load_model_table Definition: FIGMODELTable = FIGMODELmodel->load_model_table(string:table name,0/1:refresh the table)); Description: Returns the table specified by the input filename. Table will be stored in a file in the model directory. =cut sub load_model_table { my ($self,$name,$refresh) = @_; if (!defined($refresh)) { $refresh = 1; } if ($refresh == 1) { delete $self->{"_".$name}; } if (!defined($self->{"_".$name})) { my $tbldef = $self->figmodel()->config("$name"); if (!defined($tbldef)) { return undef; } my $filename = $self->directory().$name."-".$self->id().$self->selected_version().".tbl"; $self->{"_".$name} = $self->figmodel()->database()->load_table($filename,"\t","|",$tbldef->{headingline}->[0],$tbldef->{hashcolumns}); if (!defined($self->{"_".$name})) { if (defined($tbldef->{prefix})) { $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},"\t","|",join(@{$tbldef->{prefix}},"\n")); } else { $self->{"_".$name} = FIGMODELTable->new($tbldef->{columns},$filename,$tbldef->{hashcolumns},"\t","|"); } } } return $self->{"_".$name}; } =head3 create_table_prototype Definition: FIGMODELTable::table = FIGMODELmodel->create_table_prototype(string::table); Description: Returns a empty FIGMODELTable with all the metadata associated with the input table name =cut sub create_table_prototype { my ($self,$TableName) = @_; #Checking if the table definition exists in the FIGMODELconfig file if (!defined($self->figmodel()->config($TableName))) { $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:Definition not found for ".$TableName); return undef; } #Checking that this is a database table if (!defined($self->config($TableName)->{tabletype}) || $self->config($TableName)->{tabletype}->[0] ne "ModelTable") { $self->figmodel()->error_message("FIGMODELdatabase:create_table_prototype:".$TableName." is not a model table!"); return undef; } if (!defined($self->config($TableName)->{delimiter})) { $self->config($TableName)->{delimiter}->[0] = ";"; } if (!defined($self->config($TableName)->{itemdelimiter})) { $self->config($TableName)->{itemdelimiter}->[0] = "|"; } my $prefix; if (defined($self->config($TableName)->{prefix})) { $prefix = join("\n",@{$self->config($TableName)->{prefix}}); } 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); return $tbl; } =head3 get_reaction_number Definition: int = FIGMODELmodel->get_reaction_number(); Description: Returns the number of reactions in the model =cut sub get_reaction_number { my ($self) = @_; if (!defined($self->reaction_table())) { return 0; } return $self->reaction_table()->size(); } =head3 reaction_table Definition: FIGMODELTable = FIGMODELmodel->reaction_table(); Description: Returns FIGMODELTable with the reaction list for the model =cut sub reaction_table { my ($self) = @_; if (!defined($self->{_reaction_data})) { $self->{_reaction_data} = $self->figmodel()->database()->GetDBModel($self->id()); my $classTbl = $self->reaction_class_table(); for (my $i=0; $i < $classTbl->size(); $i++) { my $row = $classTbl->get_row($i); my $rxnRow = $self->{_reaction_data}->get_row_by_key($row->{"REACTION"}->[0],"LOAD"); for (my $j=0; $j < @{$row->{MEDIA}};$j++) { my $class = "Active <=>"; if ($row->{CLASS}->[$j] eq "Positive") { $class = "Essential =>"; } elsif ($row->{CLASS}->[$j] eq "Negative") { $class = "Essential <="; } elsif ($row->{CLASS}->[$j] eq "Blocked") { $class = "Inactive"; } elsif ($row->{CLASS}->[$j] eq "Positive variable") { $class = "Active =>"; } elsif ($row->{CLASS}->[$j] eq "Negative variable") { $class = "Active <="; } elsif ($row->{CLASS}->[$j] eq "Variable") { $class = "Active <=>"; } elsif ($row->{CLASS}->[$j] eq "Dead") { $class = "Dead"; } push(@{$rxnRow->{PREDICTIONS}},$row->{MEDIA}->[$j].":".$class); } } } return $self->{_reaction_data}; } =head3 feature_table Definition: FIGMODELTable = FIGMODELmodel->feature_table(); Description: Returns FIGMODELTable with the feature list for the model =cut sub feature_table { my ($self) = @_; if (!defined($self->{_feature_data})) { #Getting the genome feature list my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome()); if (!defined($FeatureTable)) { print STDERR "FIGMODELmodel:feature_table:Could not get features for genome ".$self->genome()." in database!"; return undef; } #Getting the reaction table for the model my $rxnTable = $self->reaction_table(); if (!defined($rxnTable)) { print STDERR "FIGMODELmodel:feature_table:Could not get reaction table for model ".$self->id()." in database!"; return undef; } #Cloning the feature table $self->{_feature_data} = $FeatureTable->clone_table_def(); $self->{_feature_data}->add_headings(($self->id()."REACTIONS",$self->id()."PREDICTIONS")); for (my $i=0; $i < $rxnTable->size(); $i++) { my $Row = $rxnTable->get_row($i); if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) { foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) { my $temp = $GeneSet; $temp =~ s/\+/|/g; $temp =~ s/\sAND\s/|/gi; $temp =~ s/\sOR\s/|/gi; $temp =~ s/[\(\)\s]//g; my @GeneList = split(/\|/,$temp); foreach my $Gene (@GeneList) { my $FeatureRow = $self->{_feature_data}->get_row_by_key("fig|".$self->genome().".".$Gene,"ID"); if (!defined($FeatureRow)) { $FeatureRow = $FeatureTable->get_row_by_key("fig|".$self->genome().".".$Gene,"ID"); if (defined($FeatureRow)) { $self->{_feature_data}->add_row($FeatureRow); } } if (defined($FeatureRow)) { $self->{_feature_data}->add_data($FeatureRow,$self->id()."REACTIONS",$Row->{"LOAD"}->[0],1); } } } } } #Loading predictions my @files = glob($self->directory()."EssentialGenes-".$self->id()."-*"); 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; } for (my $i=0; $i < $self->{_feature_data}->size(); $i++) { my $Row = $self->{_feature_data}->get_row($i); if ($Row->{ID}->[0] =~ m/(peg\.\d+)/) { my $gene = $1; if (defined($hash->{$gene})) { push(@{$Row->{$self->id()."PREDICTIONS"}},$media.":essential"); } else { push(@{$Row->{$self->id()."PREDICTIONS"}},$media.":nonessential"); } } } } } } return $self->{_feature_data}; } =head3 reaction_class_table Definition: FIGMODELTable = FIGMODELmodel->reaction_class_table(); Description: Returns FIGMODELTable with the reaction class data, and creates the table file if it does not exist =cut sub reaction_class_table { my ($self) = @_; 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}; } =head3 compound_class_table Definition: FIGMODELTable = FIGMODELmodel->compound_class_table(); Description: Returns FIGMODELTable with the compound class data, and creates the table file if it does not exist =cut sub compound_class_table { my ($self) = @_; 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}; } =head3 get_essential_genes Definition: [string::peg ID] = FIGMODELmodel->get_essential_genes(string::media condition); Description: Returns an reference to an array of the predicted essential genes during growth in the input media condition =cut sub get_essential_genes { my ($self,$media) = @_; if (!defined($media)) { $media = "Complete"; } if (!defined($self->{_essential_genes}->{$media})) { $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}; } =head3 compound_table Definition: FIGMODELTable = FIGMODELmodel->compound_table(); Description: Returns FIGMODELTable with the compound list for the model =cut sub compound_table { my ($self) = @_; if (!defined($self->{_compound_table})) { $self->{_compound_table} = $self->figmodel()->database()->GetDBModelCompounds($self->id()); } return $self->{_compound_table}; } =head3 get_compound_data Definition: {string:key=>[string]:values} = FIGMODELmodel->get_compound_data(string::compound ID); Description: Returns model compound data =cut sub get_compound_data { my ($self,$compound) = @_; if (!defined($self->compound_table())) { return undef; } if ($compound =~ m/^\d+$/) { return $self->compound_table()->get_row($compound); } return $self->compound_table()->get_row_by_key($compound,"DATABASE"); } =head3 get_feature_data Definition: {string:key=>[string]:values} = FIGMODELmodel->get_feature_data(string::feature ID); Description: Returns model feature data =cut sub get_feature_data { my ($self,$feature) = @_; if (!defined($self->feature_table())) { return undef; } if ($feature =~ m/^\d+$/) { return $self->feature_table()->get_row($feature); } if ($feature =~ m/(peg\.\d+)/) { $feature = $1; } return $self->feature_table()->get_row_by_key("fig|".$self->genome().".".$feature,"ID"); } =head3 genome Definition: string = FIGMODELmodel->genome(); Description: Returns model genome =cut sub genome { my ($self) = @_; return $self->{_data}->genome(); } =head3 source Definition: string = FIGMODELmodel->source(); Description: Returns model source =cut sub source { my ($self) = @_; return $self->{_data}->source(); } =head3 rights Definition: 1/0 = FIGMODELmodel->rights(string::username); Description: Returns 1 if the input user can view the model, and zero otherwise =cut sub rights { my ($self,$username) = @_; if ($self->public()) { return 1; } if (!defined($username) || $username eq "NONE") { return 0; } if (!defined($self->{_userrights}->{$username})) { $self->{_userrights}->{$self->{_data}->owner()} = 1; my @users = split(/\|/,$self->{_data}->users()); for (my $i=0; $i < @users;$i++) { $self->{_userrights}->{$users[$i]} = 1; } } return $self->{_userrights}->{$username}; } =head3 public Definition: 1/0 = FIGMODELmodel->public(); Description: Returns 1 if the model is public, and zero otherwise =cut sub public { my ($self) = @_; if ($self->{_data}->users() eq "all") { return 1; } return 0; } =head3 directory Definition: string = FIGMODELmodel->directory(); Description: Returns model directory =cut sub directory { my ($self) = @_; if (!defined($self->{_directory})) { my $userdirectory = ""; if ($self->owner() ne "master") { $userdirectory = $self->owner()."/"; } my $source = $self->source(); if ($source =~ /^MGRAST/) { $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/"; } elsif ($source =~ /^RAST/) { $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/"; } elsif ($source =~ /^SEED/) { $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->genome()."/"; } elsif ($source =~ /^PM/) { if (length($userdirectory) == 0) { $self->{_directory} = $self->figmodel()->config("imported model directory")->[0].$self->id()."/"; } else { $self->{_directory} = $self->figmodel()->config("organism directory")->[0].$userdirectory.$self->id()."/"; } } } return $self->{_directory}; } =head3 filename Definition: string = FIGMODELmodel->filename(); Description: Returns model filename =cut sub filename { my ($self) = @_; return $self->directory().$self->id().$self->selected_version().".txt"; } =head3 version Definition: string = FIGMODELmodel->version(); Description: Returns the version of the model =cut sub version { my ($self) = @_; if (!defined($self->{_version})) { if (!defined($self->{_selectedversion})) { $self->{_version} = "V".$self->{_data}->version().".".$self->{_data}->autocompleteVersion(); } else { $self->{_version} = $self->{_selectedversion}; } } return $self->{_version}; } =head3 selected_version Definition: string = FIGMODELmodel->selected_version(); Description: Returns the selected version of the model =cut sub selected_version { my ($self) = @_; if (!defined($self->{_selectedversion})) { return ""; } return $self->{_selectedversion}; } =head3 modification_time Definition: string = FIGMODELmodel->modification_time(); Description: Returns the selected version of the model =cut sub modification_time { my ($self) = @_; return $self->{_data}->modificationDate(); } =head3 gene_reactions Definition: string = FIGMODELmodel->gene_reactions(); Description: Returns the number of reactions added by the gap filling =cut sub gene_reactions { my ($self) = @_; return ($self->{_data}->reactions() - $self->{_data}->autoCompleteReactions() - $self->{_data}->spontaneousReactions() - $self->{_data}->gapFillReactions()); } =head3 total_compounds Definition: string = FIGMODELmodel->total_compounds(); Description: Returns the number of compounds in the model =cut sub total_compounds { my ($self) = @_; return $self->{_data}->compounds(); } =head3 gapfilling_reactions Definition: string = FIGMODELmodel->gapfilling_reactions(); Description: Returns the number of reactions added by the gap filling =cut sub gapfilling_reactions { my ($self) = @_; return ($self->{_data}->autoCompleteReactions()+$self->{_data}->gapFillReactions()); } =head3 total_reactions Definition: string = FIGMODELmodel->total_reactions(); Description: Returns the total number of reactions in the model =cut sub total_reactions { my ($self) = @_; return $self->{_data}->reactions(); } =head3 model_genes Definition: string = FIGMODELmodel->model_genes(); Description: Returns the number of genes mapped to one or more reactions in the model =cut sub model_genes { my ($self) = @_; return $self->{_data}->associatedGenes(); } =head3 class Definition: string = FIGMODELmodel->class(); Description: Returns the class of the model: gram positive, gram negative, other =cut sub class { my ($self) = @_; return $self->{_data}->cellwalltype(); } =head3 taxonomy Definition: string = FIGMODELmodel->taxonomy(); Description: Returns model taxonomy or biome if this is an metagenome model =cut sub taxonomy { my ($self) = @_; if (!defined($self->{_taxonomy})) { my $source = $self->source(); if ($source =~ /^MGRAST/) { $self->{_taxonomy} = "NA"; my $mgmodeldata = $self->figmodel()->database()->mg_model_data($self->genome()); if (defined($mgmodeldata)) { $self->{_taxonomy} = $mgmodeldata->biome(); } } else { $self->{_taxonomy} = $self->fig()->taxonomy_of($self->genome()); } } return $self->{_taxonomy}; } =head3 genome_size Definition: string = FIGMODELmodel->genome_size(); Description: Returns size of the modeled genome in KB =cut sub genome_size { my ($self) = @_; if (!defined($self->{_genome_size})) { my $source = $self->source(); if ($source =~ /^MGRAST/) { $self->{_genome_size} = "NA"; if (defined($self->mgdata())) { $self->{_genome_size} = $self->mgdata()->size; } } else { $self->{_genome_size} = $self->fig()->genome_szdna($self->genome()); } } return $self->{_genome_size}; } =head3 genome_genes Definition: string = FIGMODELmodel->genome_genes(); Description: Returns the number of genes in the modeled genome =cut sub genome_genes { my ($self) = @_; if (!defined($self->{_genome_genes})) { my $source = $self->source(); if ($source =~ /^MGRAST/) { $self->{_genome_genes} = "NA"; if (defined($self->mgdata())) { $self->{_genome_genes} = $self->mgdata()->genome_contig_count; } } else { $self->{_genome_genes} = $self->figmodel()->get_genome_stats($self->genome())->{"TOTAL GENES"}->[0]; } } return $self->{_genome_genes}; } =head3 run_default_model_predictions Definition: 0/1::status = FIGMODELmodel->run_default_model_predictions(string::media ID); Description: =cut sub run_default_model_predictions { my ($self,$Media) = @_; #Assuming complete media if none is provided if (!defined($Media)) { $Media = "Complete"; } #Predicting essentiality my $result = $self->figmodel()->RunFBASimulation($self->id(),"SINGLEKO",undef,undef,[$self->id()],[$Media]); #Checking that the table is defined and the output file exists if (defined($result) && defined($result->get_row(0)->{"ESSENTIALGENES"})) { $self->figmodel()->database()->print_array_to_file($self->directory()."EssentialGenes-".$self->id()."-".$Media.".tbl",[join("\n",@{$result->get_row(0)->{"ESSENTIALGENES"}})]); } else { $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not identify essential reactions for model ".$self->id().$self->selected_version()."."); return $self->figmodel()->fail(); } #Classifying reactions and compounds my $tbl = $self->classify_model_reactions($Media); if (!defined($tbl)) { $self->figmodel()->error_message("FIGMODELmodel:run_default_model_predictions:could not classify reactions for model ".$self->id().$self->selected_version()."."); return $self->figmodel()->fail(); } $tbl->save(); return $self->figmodel()->success(); } =head3 update_stats_for_gap_filling Definition: {string => [string]} = FIGMODELmodel->update_stats_for_gap_filling(int::gapfill time); Description: =cut sub update_stats_for_gap_filling { my ($self,$gapfilltime) = @_; $self->{_data}->autoCompleteTime($gapfilltime); $self->{_data}->autocompleteDate(time()); $self->{_data}->modificationDate(time()); my $version = $self->{_data}->autocompleteVersion(); $self->{_data}->autocompleteVersion($version+1); } =head3 update_stats_for_build Definition: {string => [string]} = FIGMODELmodel->update_stats_for_build(); Description: =cut sub update_stats_for_build { my ($self) = @_; $self->{_data}->builtDate(time()); $self->{_data}->modificationDate(time()); my $version = $self->{_data}->version(); $self->{_data}->version($version+1); } =head3 update_model_stats Definition: FIGMODELmodel->update_model_stats(); Description: =cut sub update_model_stats { my ($self) = @_; #Getting reaction table my $rxntbl = $self->reaction_table(); if (!defined($rxntbl)) { $self->figmodel()->error_message("FIGMODELmodel:update_model_stats:Could not load reaction list for ".$self->id()); return undef; } my $cpdtbl = $self->compound_table(); #Calculating all necessary stats my %GeneHash; my %NonpegHash; my %CompoundHash; my $spontaneousReactions = 0; my $gapFillReactions = 0; my $biologReactions = 0; my $transporters = 0; my $autoCompleteReactions = 0; my $associatedSubsystemGenes = 0; for (my $i=0; $i < $rxntbl->size(); $i++) { my $Row = $rxntbl->get_row($i); if (defined($Row) && defined($Row->{"ASSOCIATED PEG"})) { my $ReactionRow = $self->figmodel()->get_reaction($Row->{"LOAD"}->[0]); if (defined($ReactionRow->{"EQUATION"}->[0])) { #Checking for extracellular metabolites which indicate that this is a transporter if ($ReactionRow->{"EQUATION"}->[0] =~ m/\[e\]/) { $transporters++; } } #Identifying spontaneous/biolog/gapfilling/gene associated reactions if ($Row->{"ASSOCIATED PEG"}->[0] =~ m/BIOLOG/i) { $biologReactions++; } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GROW/i) { $gapFillReactions++; } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/SPONTANEOUS/i) { $spontaneousReactions++; } elsif ($Row->{"ASSOCIATED PEG"}->[0] =~ m/GAP/ || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNIVERSAL/i || $Row->{"ASSOCIATED PEG"}->[0] =~ m/UNKNOWN/i) { $autoCompleteReactions++; } else { foreach my $GeneSet (@{$Row->{"ASSOCIATED PEG"}}) { $_ = $GeneSet; my @GeneList = /(peg\.\d+)/g; foreach my $Gene (@GeneList) { if ($Gene =~ m/(peg\.\d+)/) { $GeneHash{$1} = 1; } else { $NonpegHash{$Gene} = 1; } } } } } } my @genes = keys(%GeneHash); my @othergenes = keys(%NonpegHash); #Setting the reaction count $self->{_data}->reactions($rxntbl->size()); #Setting the metabolite count $self->{_data}->compounds($rxntbl->size()); #Setting the gene count my $geneCount = @genes + @othergenes; $self->{_data}->associatedGenes($geneCount); #Setting remaining stats $self->{_data}->spontaneousReactions($spontaneousReactions); $self->{_data}->gapFillReactions($gapFillReactions); $self->{_data}->biologReactions($biologReactions); $self->{_data}->transporters($transporters); $self->{_data}->autoCompleteReactions($autoCompleteReactions); $self->{_data}->associatedSubsystemGenes($associatedSubsystemGenes); #Setting the model class my $class = ""; for (my $i=0; $i < @{$self->figmodel()->config("class list")}; $i++) { if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i]))) { if (defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->id()})) { $class = $self->figmodel()->config("class list")->[$i]; last; } if ($class eq "" && defined($self->figmodel()->config($self->figmodel()->config("class list")->[$i])->{$self->genome()})) { $class = $self->figmodel()->config("class list")->[$i]; } } } if ($class eq "") { $class = $self->figmodel()->get_genome_stats($self->genome())->{CLASS}->[0]; } if ($class eq "") { $class = "unknown"; } $self->{_data}->cellwalltype($class); } =head3 GapFillModel Definition: (success/fail) = FIGMODELmodel->GapFillModel(); Description: This function performs an optimization to identify the minimal set of reactions that must be added to a model in order for biomass to be produced by the biomass reaction in the model. Before running the gap filling, the existing model is backup in the same directory with the current version numbers appended. If the model has been gap filled previously, the previous gap filling reactions are removed prior to running the gap filling again. =cut sub GapFillModel { my ($self,$donotclear) = @_; #Setting status of model to gap filling my $OrganismID = $self->genome(); $self->set_status(1,"Auto completion running"); my $UniqueFilename = $self->figmodel()->filename(); my $StartTime = time(); #Archiving the existing model $self->ArchiveModel(); #Removing any gapfilling reactions that may be currently present in the model if (!defined($donotclear) || $donotclear != 1) { my $ModelTable = $self->reaction_table(); for (my $i=0; $i < $ModelTable->size(); $i++) { 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/) { $ModelTable->delete_row($i); $i--; } } $ModelTable->save(); } #Calling the MFAToolkit to run the gap filling optimization my $MinimalMediaTable = $self->figmodel()->database()->GetDBTable("MINIMAL MEDIA TABLE"); my $Media = "Complete"; if (defined($MinimalMediaTable->get_row_by_key($self->genome(),"Organism"))) { $Media = $MinimalMediaTable->get_row_by_key($self->genome(),"Organism")->{"Minimal media"}->[0]; #Loading media, changing bounds, saving media as a test media my $MediaTable = FIGMODELTable::load_table($self->config("Media directory")->[0].$Media.".txt",";","",0,["VarName"]); for (my $i=0; $i < $MediaTable->size(); $i++) { if ($MediaTable->get_row($i)->{"Min"}->[0] < 0) { $MediaTable->get_row($i)->{"Min"}->[0] = -10000; } if ($MediaTable->get_row($i)->{"Max"}->[0] > 0) { $MediaTable->get_row($i)->{"Max"}->[0] = 10000; } } $MediaTable->save($self->config("Media directory")->[0].$UniqueFilename."TestMedia.txt"); 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)); 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)); } #Parse the solutions file for the model and get the reaction list from it my $SolutionData = $self->figmodel()->LoadProblemReport($UniqueFilename); #Clearing the mfatoolkit output and log file $self->figmodel()->clearing_output($UniqueFilename,"GapFill".$self->id().".log"); if (!defined($SolutionData) || $SolutionData->size() == 0) { $self->set_status(1,"Auto completion failed: auto completion solution not found"); $self->figmodel()->error_message("FIGMODEL:GapFillModel: Gap filling report not found!"); return $self->fail(); } $SolutionData->save("/home/chenry/Solution.txt"); #Looking for the last printed recursive MILP solution for (my $i=($SolutionData->size()-1); $i >=0; $i--) { if (defined($SolutionData->get_row($i)->{"Notes"}) && $SolutionData->get_row($i)->{"Notes"}->[0] =~ m/^Recursive/) { my $AllSolutions = substr($SolutionData->get_row($i)->{"Notes"}->[0],15); my @TempThree = split(/\|/,$AllSolutions); if (@TempThree > 0 && $TempThree[0] =~ m/.+:(.+)/) { my @TempFour = split(/,/,$1); my $DirectionList; my $ReactionList; for (my $j=0; $j < @TempFour; $j++) { my $ID = ""; my $Sign = "=>"; if ($TempFour[$j] =~ m/\-(rxn\d\d\d\d\d)/) { $ID = $1; $Sign = "<="; } elsif ($TempFour[$j] =~ m/\+(rxn\d\d\d\d\d)/) { $ID = $1; } if ($ID ne "") { if ($self->figmodel()->reversibility_of_reaction($ID) ne $Sign) { $Sign = "<=>"; } push(@{$DirectionList},$Sign); push(@{$ReactionList},$ID); } } $self->figmodel()->IntegrateGrowMatchSolution($self->id(),undef,$ReactionList,$DirectionList,"GAP FILLING",0,1); } last; } } #Updating model stats with gap filling results my $ElapsedTime = time() - $StartTime; $self->figmodel()->database()->ClearDBModel($self->id(),1); #Determining why each gap filling reaction was added $self->figmodel()->IdentifyDependancyOfGapFillingReactions($self->id(),$Media); if (!defined($donotclear) || $donotclear != 1) { $self->figmodel()->AddBiologTransporters($self->id()); if ($self->id() !~ m/MGRast/) { $self->update_stats_for_gap_filling($ElapsedTime); } } else { $self->update_model_stats(); } #Printing the updated SBML file $self->PrintSBMLFile(); $self->PrintModelLPFile($self->id()); $self->set_status(1,"Auto completion successfully finished"); $self->run_default_model_predictions(); return $self->success(); } =head3 GapGenModel Definition: FIGMODELmodel->GapGenModel(); Description: Runs the gap generation algorithm to correct a single false positive prediction. Results are loaded into a table. =cut sub GapGenModel { my ($self,$Media,$KOList,$NoKOList,$Experiment,$SolutionLimit) = @_; #Enforcing nonoptional arguments if (!defined($Media)) { return undef; } if (!defined($KOList)) { $KOList->[0] = "none"; } if (!defined($NoKOList)) { $NoKOList->[0] = "none"; } if (!defined($Experiment)) { $Experiment= "ReactionKO"; } if (!defined($SolutionLimit)) { $SolutionLimit = "10"; } #Translating the KO lists into arrays if (ref($KOList) ne "ARRAY") { my $temp = $KOList; $KOList = (); push(@{$KOList},split(/[,;]/,$temp)); } my $noKOHash; if (defined($NoKOList) && ref($NoKOList) ne "ARRAY") { my $temp = $NoKOList; $NoKOList = (); push(@{$NoKOList},split(/[,;]/,$temp)); foreach my $rxn (@{$NoKOList}) { $noKOHash->{$rxn} = 1; } } #Checking if solutions exist for the input parameters $self->aquireModelLock(); my $tbl = $self->load_model_table("GapGenSolutions"); my $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist"); my $solutions; if (defined($solutionRow)) { #Checking if any solutions conform to the no KO list foreach my $solution (@{$solutionRow->{Solutions}}) { my @reactions = split(/,/,$solution); my $include = 1; foreach my $rxn (@reactions) { if ($rxn =~ m/(rxn\d\d\d\d\d)/) { if (defined($noKOHash->{$1})) { $include = 0; } } } if ($include == 1) { push(@{$solutions},$solution); } } } else { $solutionRow = {Media => [$Media],Experiment => [$Experiment],KOlist => [join(",",@{$KOList})]}; $tbl->add_row($solutionRow); $self->figmodel()->database()->save_table($tbl); } $self->releaseModelLock(); #Returning solution list of solutions were found if (defined($solutions) && @{$solutions} > 0) { return $solutions; } #Getting unique filename my $Filename = $self->figmodel()->filename(); #Running the gap generation system($self->figmodel()->GenerateMFAToolkitCommandLineCall($Filename,$self->id().$self->selected_version(),$Media,["GapGeneration"],{"Recursive MILP solution limit" => $SolutionLimit ,"Reactions that should always be active" => join(";",@{$NoKOList}),"Reactions to knockout" => join(";",@{$KOList}),"Reactions that are always blocked" => "none"},"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log",undef,undef)); my $ProblemReport = $self->figmodel()->LoadProblemReport($Filename); if (!defined($ProblemReport)) { $self->figmodel()->error_message("FIGMODEL:GapGenerationAlgorithm;No problem report;".$Filename.";".$self->id().$self->selected_version().";".$Media.";".$KOList.";".$NoKOList); return undef; } #Clearing the output folder and log file $self->figmodel()->clearing_output($Filename,"Gapgeneration-".$self->id().$self->selected_version()."-".$Filename.".log"); #Saving the solution $self->aquireModelLock(); $tbl = $self->load_model_table("GapGenSolutions"); $solutionRow = $tbl->get_table_by_key($Experiment,"Experiment")->get_table_by_key($Media,"Media")->get_row_by_key(join(",",@{$KOList}),"KOlist"); for (my $j=0; $j < $ProblemReport->size(); $j++) { if ($ProblemReport->get_row($j)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^)]+)/) { my @SolutionList = split(/\|/,$1); for (my $k=0; $k < @SolutionList; $k++) { if ($SolutionList[$k] =~ m/(\d+):(.+)/) { push(@{$solutionRow->{Solutions}},$2); push(@{$solutions},$2); } } } } $self->figmodel()->database()->save_table($tbl); $self->releaseModelLock(); return $solutions; } =head3 datagapfill Definition: success()/fail() = FIGMODELmodel->datagapfill(); Description: Run gapfilling on the input run specifications =cut sub datagapfill { my ($self,$GapFillingRunSpecs,$TansferFileSuffix) = @_; my $UniqueFilename = $self->figmodel()->filename(); if (defined($GapFillingRunSpecs) && @{$GapFillingRunSpecs} > 0) { system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id().$self->selected_version(),"NoBounds",["DataGapFilling"],{"Reactions to knockout" => $self->config("permanently knocked out reactions")->[0],"Gap filling runs" => join(";",@{$GapFillingRunSpecs})},"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,undef)); #Checking that the solution exists if (!-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt") { $self->figmodel()->error_message("FIGMODEL:GapFillingAlgorithm: Could not find MFA output file!"); $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",["Experiment;Solution index;Solution cost;Solution reactions"]); return undef; } my $GapFillResultTable = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt",";","",0,undef); if (defined($TansferFileSuffix)) { system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/GapFillingSolutionTable.txt ".$self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt"); } #If the system is not configured to preserve all logfiles, then the mfatoolkit output folder is deleted $self->figmodel()->clearing_output($UniqueFilename,"GapFilling-".$self->id().$self->selected_version()."-".$UniqueFilename.".log"); return $GapFillResultTable; } if (defined($TansferFileSuffix)) { $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-".$TansferFileSuffix.".txt",["Experiment;Solution index;Solution cost;Solution reactions"]); } return undef; } =head3 TestSolutions Definition: $model->TestSolutions($ModelID,$NumProcessors,$ProcessorIndex,$GapFill); Description: Example: =cut sub TestSolutions { my ($self,$OriginalErrorFilename,$GapFillResultTable) = @_; #Getting the filename my $UniqueFilename = $self->figmodel()->filename(); #Reading in the original error matrix which has the headings for the original model simulation my $OriginalErrorData; if (!defined($OriginalErrorFilename) || !-e $self->directory().$OriginalErrorFilename) { my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All"); $OriginalErrorData = [$HeadingVector,$Errorvector]; } else { $OriginalErrorData = $self->figmodel()->database()->load_single_column_file($self->directory().$OriginalErrorFilename,""); } my $HeadingHash; my @HeadingArray = split(/;/,$OriginalErrorData->[0]); my @OrigErrorArray = split(/;/,$OriginalErrorData->[1]); for (my $i=0; $i < @HeadingArray; $i++) { my @SubArray = split(/:/,$HeadingArray[$i]); $HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]} = $i; } #Scanning through the gap filling solutions my $TempVersion = "V".$UniqueFilename; my $ErrorMatrixLines; for (my $i=0; $i < $GapFillResultTable->size(); $i++) { print "Starting problem solving ".$i."\n"; my $ErrorLine = $GapFillResultTable->get_row($i)->{"Experiment"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution index"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution cost"}->[0].";".$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0]; #Integrating solution into test model my $ReactionArray; my $DirectionArray; my @ReactionList = split(/,/,$GapFillResultTable->get_row($i)->{"Solution reactions"}->[0]); my %SolutionHash; for (my $k=0; $k < @ReactionList; $k++) { if ($ReactionList[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) { my $Reaction = $2; my $Sign = $1; if (defined($SolutionHash{$Reaction})) { $SolutionHash{$Reaction} = "<=>"; } elsif ($Sign eq "-") { $SolutionHash{$Reaction} = "<="; } elsif ($Sign eq "+") { $SolutionHash{$Reaction} = "=>"; } else { $SolutionHash{$Reaction} = $Sign; } } } @ReactionList = keys(%SolutionHash); for (my $k=0; $k < @ReactionList; $k++) { push(@{$ReactionArray},$ReactionList[$k]); push(@{$DirectionArray},$SolutionHash{$ReactionList[$k]}); } print "Integrating solution!\n"; $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"Gapfilling ".$GapFillResultTable->get_row($i)->{"Experiment"}->[0],1,1); $self->PrintModelLPFile(); #Running the model against all available experimental data print "Running test model!\n"; my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All"); @HeadingArray = split(/;/,$HeadingVector); my @ErrorArray = @OrigErrorArray; my @TempArray = split(/;/,$Errorvector); for (my $j=0; $j < @HeadingArray; $j++) { my @SubArray = split(/:/,$HeadingArray[$j]); $ErrorArray[$HeadingHash->{$SubArray[0].":".$SubArray[1].":".$SubArray[2]}] = $TempArray[$j]; } $ErrorLine .= ";".$FalsePostives."/".$FalseNegatives.";".join(";",@ErrorArray); push(@{$ErrorMatrixLines},$ErrorLine); print "Finishing problem solving ".$i."\n"; } #Clearing out the test model if (-e $self->directory().$self->id().$TempVersion.".txt") { unlink($self->directory().$self->id().$TempVersion.".txt"); unlink($self->directory()."SimulationOutput".$self->id().$TempVersion.".txt"); } return $ErrorMatrixLines; } =head3 status Definition: int::model status = FIGMODELmodel->status(); Description: Returns the current status of the SEED model associated with the input genome ID. model status = 1: model exists model status = 0: model is being built model status = -1: model does not exist model status = -2: model build failed =cut sub status { my ($self) = @_; return $self->{_data}->{status}->[0]; } =head3 message Definition: string::model message = FIGMODELmodel->message(); Description: Returns a message associated with the models current status =cut sub message { my ($self) = @_; return $self->{_data}->{message}->[0]; } =head3 set_status Definition: (success/fail) = FIGMODELmodel->set_status(int::new status,string::status message); Description: Changes the current status of the SEED model new status = 1: model exists new status = 0: model is being built new status = -1: model does not exist new status = -2: model build failed =cut sub set_status { my ($self,$NewStatus,$Message) = @_; #Getting the model row from the MODELS table $self->{_data}->{status}->[0] = $NewStatus; $self->{_data}->{message}->[0] = $Message; $self->figmodel()->database()->update_row("MODELS",$self->{_data},"id"); return $self->config("SUCCESS")->[0]; } =head3 CreateSingleGenomeReactionList Definition: FIGMODEL->CreateSingleGenomeReactionList(); Description: This function uses fig calls to obtain a list of genes and functions for a genome, and it uses a file mapping reactions and functional roles to produce a reaction list. Example: =cut sub CreateSingleGenomeReactionList { my ($self,$RunGapFilling) = @_; #Creating directory if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/") { system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/"); } elsif ($self->owner() eq "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->genome()."/") { system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->genome()."/"); } if ($self->owner() ne "master" && !-d $self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/") { system("mkdir ".$self->figmodel()->config("organism directory")->[0].$self->owner()."/".$self->genome()."/"); } #Getting genome stats my $genomestats = $self->figmodel()->get_genome_stats($self->genome()); my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome()); if (!defined($FeatureTable)) { $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome features could not be accessed!"); return $self->fail(); } #Checking that the number of genes exceeds the minimum size if ($FeatureTable->size() < $self->config("minimum genome size for modeling")->[0]) { $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList: ".$self->id()." genome rejected as too small for modeling!"); return $self->fail(); } #Setting up needed variables my $OriginalModelTable = undef; #Locking model status table 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) { $self->figmodel()->error_message("FIGMODEL:CreateSingleGenomeReactionList:Model is already being built. Canceling current build."); $self->figmodel()->database()->UnlockDBTable("MODELS"); return $self->fail(); }elsif ($Row->{status}->[0] == 1) { $OriginalModelTable = $self->reaction_table(); $self->ArchiveModel(); $Row->{status}->[0] = 0; $Row->{message}->[0] = "Rebuilding preliminary reconstruction"; } else { $Row->{status}->[0] = 0; $Row->{message}->[0] = "Preliminary reconstruction"; } #Updating the status table $self->figmodel()->database()->save_table($ModelTable); $self->figmodel()->database()->UnlockDBTable("MODELS"); #Sorting GenomeData by gene location on the chromosome $FeatureTable->sort_rows("MIN LOCATION"); my ($ComplexHash,$SuggestedMappings,$UnrecognizedReactions,$ReactionHash); my %LastGenePosition; my $GeneRoles; for (my $j=0; $j < $FeatureTable->size(); $j++) { my $CurrentRow = $FeatureTable->get_row($j); #"ID","ALIASES","MIN LOCATION","MAX LOCATION","ROLES","SUBSYSTEMS","SUBSYSTEM CLASS" if (defined($CurrentRow)) { my $GeneID = $CurrentRow->{"ID"}->[0]; if ($GeneID =~ m/(peg\.\d+)/) { $GeneID = $1; } foreach my $Role (@{$CurrentRow->{"ROLES"}}) { if ($self->figmodel()->role_is_valid($Role) != 0) { push(@{$GeneRoles->{$GeneID}},$Role); my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role); if ($ReactionHashArrayRef != 0) { foreach my $Mapping (@{$ReactionHashArrayRef}) { if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) { if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) { if ($Mapping->{"MASTER"}->[0] eq 1) { #Creating a complex if consecutive genes have been assigned to the same reaction $ComplexHash->{$Mapping->{"REACTION"}->[0]}->{$Mapping->{"COMPLEX"}->[0]}->{$Role}->{$GeneID} = 1; if (!defined($LastGenePosition{$Mapping->{"REACTION"}->[0]})) { $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j; push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID); } elsif (($j-$LastGenePosition{$Mapping->{"REACTION"}->[0]}) < 3 && $LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) { my $CurrentComplex = pop(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}}); push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$CurrentComplex."+".$GeneID); } elsif ($LastGenePosition{$Mapping->{"REACTION"}->[0]} != $j) { push(@{$ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"GENES"}},$GeneID); } $LastGenePosition{$Mapping->{"REACTION"}->[0]} = $j; #Adding a subsystem for the reaction if ($self->figmodel()->subsystem_is_valid($Mapping->{"SUBSYSTEM"}->[0]) == 1) { ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SUBSYSTEMS"},$Mapping->{"SUBSYSTEM"}->[0]); if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 1) { if ($Mapping->{"SOURCE"}->[0] =~ m/Hope\sFiles/) { $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 1; } elsif ($Mapping->{"SOURCE"}->[0] =~ m/SEED/) { $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 2; } elsif (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) { $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3; } } } #Handling confidence if (!defined($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}) || $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] > 2) { if ($Mapping->{"SOURCE"}->[0] =~ m/MATT/) { $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 3; } elsif ($Mapping->{"SOURCE"}->[0] =~ m/CHRIS/) { $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 4; } else { $ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"CONFIDENCE"}->[0] = 5; } } #Parsing sources ($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},my $NumMatches) = $self->figmodel()->add_elements_unique($ReactionHash->{$Mapping->{"REACTION"}->[0]}->{"SOURCE"},split(/\|/,$Mapping->{"SOURCE"}->[0])); } else { push(@{$SuggestedMappings},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role); } } else { push(@{$UnrecognizedReactions},$GeneID."\t".$Mapping->{"REACTION"}->[0]."\t".$Role); } } } } } } } } #Creating nonadjacent complexes my @ReactionList = keys(%{$ReactionHash}); foreach my $Reaction (@ReactionList) { #If multiple genes are assigned to the reaction, we check if they should should be in a complex if (@{$ReactionHash->{$Reaction}->{"GENES"}} > 0 && defined($ComplexHash->{$Reaction})) { my $GeneArray; foreach my $Complex (keys(%{$ComplexHash->{$Reaction}})) { my %ComplexComponents; foreach my $CurrentGeneSet (@{$ReactionHash->{$Reaction}->{"GENES"}}) { my @GeneList = split(/\+/,$CurrentGeneSet); my %RoleHash; foreach my $Gene (@GeneList) { foreach my $Role (@{$GeneRoles->{$Gene}}) { if (defined($ComplexHash->{$Reaction}->{$Complex}->{$Role})) { $RoleHash{$Role} = 1; } } } if (keys(%RoleHash) > 0) { if (!defined($ComplexComponents{join("|",sort(keys(%RoleHash)))})) { my @RoleList = keys(%RoleHash); my @ComplexList = keys(%ComplexComponents); foreach my $ComplexSet (@ComplexList) { my @RoleList = split(/\|/,$ComplexSet); my $Match = 0; foreach my $SingleRole (@RoleList) { if (defined($RoleHash{$SingleRole})) { $Match = 1; last; } } if ($Match == 1) { foreach my $SingleRole (@RoleList) { $RoleHash{$SingleRole} = 1 } push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},@{$ComplexComponents{$ComplexSet}}); delete $ComplexComponents{$ComplexSet}; } } } push(@{$ComplexComponents{join("|",sort(keys(%RoleHash)))}},$CurrentGeneSet); } } my @Position; my @Options; my $Count = 0; foreach my $RoleSet (keys(%ComplexComponents)) { push(@Position,0); push(@{$Options[$Count]},@{$ComplexComponents{$RoleSet}}); $Count++; } my $Done = 0; $Count = 0; my $NewRelationship; while($Done == 0) { #Creating complex with current indecies $NewRelationship->[$Count] = $Options[0]->[$Position[0]]; for (my $i=1; $i < @Position; $i++) { $NewRelationship->[$Count] .= "+".$Options[$i]->[$Position[$i]]; } $NewRelationship->[$Count] = join("+",$self->figmodel()->remove_duplicates(split(/\+/,$NewRelationship->[$Count]))); $Count++; #Iterating indecies my $CurrentIndex = 0; while($CurrentIndex >= 0) { if ($CurrentIndex >= @Position) { $CurrentIndex = -1000; } elsif ($Position[$CurrentIndex]+1 == @{$Options[$CurrentIndex]}) { $Position[$CurrentIndex] = -1; $CurrentIndex++; } else { $Position[$CurrentIndex]++; $CurrentIndex--; } } if ($CurrentIndex == -1000) { $Done = 1; } } push(@{$GeneArray},@{$NewRelationship}); } @{$ReactionHash->{$Reaction}->{"GENES"}} = $self->figmodel()->remove_duplicates(@{$GeneArray}); } } #Getting the reaction table my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS"); #Creating the model reaction table my $NewModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n"); @ReactionList = keys(%{$ReactionHash}); foreach my $ReactionID (@ReactionList) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID); my $Subsystem = "NONE"; if (defined($ReactionHash->{$ReactionID}->{"SUBSYSTEMS"})) { $Subsystem = join("|",@{$ReactionHash->{$ReactionID}->{"SUBSYSTEMS"}}); } my $Source = "NONE"; if (defined($ReactionHash->{$ReactionID}->{"SOURCE"})) { $Source = join("|",@{$ReactionHash->{$ReactionID}->{"SOURCE"}}); } $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"]}); } #Adding the spontaneous and universal reactions foreach my $ReactionID (@{$self->config("spontaneous reactions")}) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID); #Checking if the reaction is already in the model if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) { $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]}); } } foreach my $ReactionID (@{$self->config("universal reactions")}) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID); #Checking if the reaction is already in the model if (!defined($NewModelTable->get_row_by_key($ReactionID,"LOAD"))) { $NewModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]}); } } #Checking if a biomass reaction already exists my $BiomassReactionRow = $self->BuildSpecificBiomassReaction(); if (!defined($BiomassReactionRow)) { $self->set_status(-2,"Preliminary reconstruction failed: could not generate biomass reaction"); $self->figmodel()->error_message("FIGMODEL:CreateModelReactionList: Could not generate biomass function for ".$self->id()."!"); return $self->fail(); } my $ReactionList = $BiomassReactionRow->{"ESSENTIAL REACTIONS"}; push(@{$ReactionList},$BiomassReactionRow->{DATABASE}->[0]); #Adding biomass reactions to the model table foreach my $BOFReaction (@{$ReactionList}) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($BOFReaction); #Checking if the reaction is already in the model if (!defined($NewModelTable->get_row_by_key($BOFReaction,"LOAD"))) { if ($BOFReaction =~ m/bio/) { $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["BOF"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [1],"REFERENCE" => ["Biomass objective function"],"NOTES" => ["NONE"]}); } else { $NewModelTable->add_row({"LOAD" => [$BOFReaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["INITIAL GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Initial gap filling"],"NOTES" => ["NONE"]}); } } } #Completing any incomplete reactions sets my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS"); for (my $i=0; $i < $ReactionSetTable->size(); $i++) { if (defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) { foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) { if (!defined($NewModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction); $NewModelTable->add_row({"LOAD" => [$Reaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["REACTION SET GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Added due to presence of ".$ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0]],"NOTES" => ["NONE"]}); } } } } #Now we compare the model to the previous model to determine if any differences exist that aren't gap filling reactions if (defined($OriginalModelTable)) { my $PerfectMatch = 1; my $ReactionCount = 0; for (my $i=0; $i < $OriginalModelTable->size(); $i++) { #We only check that nongapfilling reactions exist in the new model if ($OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] !~ m/GAP/ || $OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}->[0] eq "INITIAL GAP FILLING") { $ReactionCount++; my $Row = $NewModelTable->get_row_by_key($OriginalModelTable->get_row($i)->{"LOAD"}->[0],"LOAD"); if (defined($Row)) { #We check that the reaction directionality is identical if ($Row->{"DIRECTIONALITY"}->[0] ne $OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]) { if (defined($OriginalModelTable->get_row($i)->{"NOTES"}->[0]) && $OriginalModelTable->get_row($i)->{"NOTES"}->[0] =~ m/Directionality\sswitched\sfrom\s([^\s])/) { if ($1 ne $Row->{"DIRECTIONALITY"}->[0]) { print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$1." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n"; $PerfectMatch = 0; last; } } else { print "Directionality mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$OriginalModelTable->get_row($i)->{"DIRECTIONALITY"}->[0]." vs ".$Row->{"DIRECTIONALITY"}->[0]."\n"; $PerfectMatch = 0; last; } } #We check that the genes assigned to the reaction are identical if ($PerfectMatch == 1 && @{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}} != @{$Row->{"ASSOCIATED PEG"}}) { print "Gene associatation mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}}." vs ".@{$Row->{"ASSOCIATED PEG"}}."\n"; $PerfectMatch = 0; last; } if ($PerfectMatch == 1) { my @GeneSetOne = sort(@{$OriginalModelTable->get_row($i)->{"ASSOCIATED PEG"}}); my @GeneSetTwo = sort(@{$Row->{"ASSOCIATED PEG"}}); for (my $j=0; $j < @GeneSetOne; $j++) { if ($GeneSetOne[$j] ne $GeneSetTwo[$j]) { print "Gene mismatch for reaction ".$OriginalModelTable->get_row($i)->{"LOAD"}->[0].": ".$GeneSetOne[$j]." vs ".$GeneSetTwo[$j]."\n"; $PerfectMatch = 0; $i = $OriginalModelTable->size(); last; } } } } else { print "Original model contains an extra reaction:".$OriginalModelTable->get_row($i)->{"LOAD"}->[0]."\n"; $PerfectMatch = 0; last; } } } if ($PerfectMatch == 1 && $ReactionCount == $NewModelTable->size()) { #Bailing out of function as the model has not changed $self->set_status(1,"rebuild canceled because model has not changed"); return $self->success(); } } #Saving the new model to file $self->set_status(1,"Preliminary reconstruction complete"); $self->figmodel()->database()->save_table($NewModelTable); $self->{_reaction_data} = $NewModelTable; #Clearing the previous model from the cache $self->figmodel()->database()->ClearDBModel($self->id(),1); #Updating the model stats table $self->update_stats_for_build(); $self->PrintSBMLFile(); #Adding model to gapfilling queue if (defined($RunGapFilling) && $RunGapFilling == 1) { $self->set_status(1,"Autocompletion queued"); $self->figmodel()->add_job_to_queue("gapfillmodel?".$self->id(),"QSUB","cplex",$self->owner(),"BACK"); } return $self->success(); } =head3 CreateMetaGenomeReactionList Definition: (success/fail) = FIGMODELmodel->CreateMetaGenomeReactionList(); Description: This is the code called to create or update the reaction list for a metgenome model =cut sub CreateMetaGenomeReactionList { my ($self) = @_; #Checking if the metagenome file exists if (!-e $self->config("raw MGRAST directory")->[0].$self->genome().".summary") { $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome()); } #Loading metagenome data my $MGRASTData = $self->figmodel()->database()->load_multiple_column_file($self->config("raw MGRAST directory")->[0].$self->genome().".summary","\t"); if (!defined($MGRASTData)) { $self->error_message("FIGMODEL:CreateMetaGenomeReactionList: could not find raw data file for metagenome ".$self->genome()); } #Setting up needed variables my $OriginalModelTable = undef; #Checking status if ($self->status() < 0) { $self->set_status(0,"Preliminary reconstruction"); } elsif ($self->status() == 0) { $self->error_message("FIGMODEL->CreateModelReactionList:Model is already being built. Canceling current build."); return $self->fail(); } else { $OriginalModelTable = $self->reaction_table(); $self->ArchiveModel(); $self->set_status(0,"Rebuilding preliminary reconstruction"); } #Getting the reaction table my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS"); #Creating model table my $ModelTable = FIGMODELTable->new(["LOAD","DIRECTIONALITY","COMPARTMENT","ASSOCIATED PEG","SUBSYSTEM","CONFIDENCE","REFERENCE","NOTES"],$self->directory().$self->id().".txt",["LOAD"],";","|","REACTIONS\n"); for (my $i=0; $i < @{$MGRASTData};$i++) { #MD5,PEG,number of sims,role,sim e-scores my $Role = $MGRASTData->[$i]->[3]; my $MD5 = $MGRASTData->[$i]->[0]; my $peg = $MGRASTData->[$i]->[1]; my $sims = $MGRASTData->[$i]->[4]; $sims =~ s/;/,/g; #Checking for subsystems my $GeneSubsystems = $self->figmodel()->subsystems_of_role($Role); #Checking if there are reactions associated with this role my $ReactionHashArrayRef = $self->figmodel()->reactions_of_role($Role); if ($ReactionHashArrayRef != 0) { foreach my $Mapping (@{$ReactionHashArrayRef}) { if (defined($Mapping->{"REACTION"}) && defined($Mapping->{"MASTER"}) && defined($Mapping->{"SUBSYSTEM"}) && defined($Mapping->{"SOURCE"})) { if ($Mapping->{"REACTION"}->[0] =~ m/rxn\d\d\d\d\d/) { if ($Mapping->{"MASTER"}->[0] eq 1) { #Checking if the reaction is already in the model my $ReactionRow = $ModelTable->get_row_by_key($Mapping->{"REACTION"}->[0],"LOAD"); if (!defined($ReactionRow)) { $ReactionRow = {"LOAD" => [$Mapping->{"REACTION"}->[0]],"DIRECTIONALITY" => [$self->figmodel()->reversibility_of_reaction($Mapping->{"REACTION"}->[0])],"COMPARTMENT" => ["c"]}; $ModelTable->add_row($ReactionRow); } push(@{$ReactionRow->{"ASSOCIATED PEG"}},substr($peg,4)); push(@{$ReactionRow->{"REFERENCE"}},$MD5.":".$Role); push(@{$ReactionRow->{"CONFIDENCE"}},$sims); if (defined($GeneSubsystems)) { push(@{$ReactionRow->{"SUBSYSTEM"}},@{$GeneSubsystems}); } } } } } } } #Adding the spontaneous and universal reactions foreach my $ReactionID (@{$self->config("spontaneous reactions")}) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID); #Checking if the reaction is already in the model if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) { $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["SPONTANEOUS"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["SPONTANEOUS"],"NOTES" => ["NONE"]}); } } foreach my $ReactionID (@{$self->config("universal reactions")}) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($ReactionID); #Checking if the reaction is already in the model if (!defined($ModelTable->get_row_by_key($ReactionID,"LOAD"))) { $ModelTable->add_row({"LOAD" => [$ReactionID],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["UNIVERSAL"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [4],"REFERENCE" => ["UNIVERSAL"],"NOTES" => ["NONE"]}); } } #Completing any incomplete reactions sets my $ReactionSetTable = $self->figmodel()->database()->GetDBTable("REACTION SETS"); for (my $i=0; $i < $ReactionSetTable->size(); $i++) { if (defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) { foreach my $Reaction (@{$ReactionSetTable->get_row($i)->{"Dependant reactions"}}) { if (!defined($ModelTable->get_row_by_key($ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0],"LOAD"))) { #Getting the thermodynamic reversibility from the database my $Directionality = $self->figmodel()->reversibility_of_reaction($Reaction); $ModelTable->add_row({"LOAD" => [$Reaction],"DIRECTIONALITY" => [$Directionality],"COMPARTMENT" => ["c"],"ASSOCIATED PEG" => ["REACTION SET GAP FILLING"],"SUBSYSTEM" => ["NONE"],"CONFIDENCE" => [5],"REFERENCE" => ["Added due to presence of ".$ReactionSetTable->get_row($i)->{"Trigger reaction"}->[0]],"NOTES" => ["NONE"]}); } } } } #Clearing the previous model from the cache $self->figmodel()->database()->ClearDBModel($self->id(),1); $ModelTable->save(); return $self->success(); } =head3 ArchiveModel Definition: (success/fail) = FIGMODELmodel->ArchiveModel(); Description: This function archives the specified model in the model directory with the current version numbers appended. This function is used to preserve old versions of models prior to overwriting so new versions may be compared with old versions. =cut sub ArchiveModel { my ($self) = @_; #Checking that the model file exists if (!(-e $self->filename())) { $self->figmodel()->error_message("FIGMODEL:ArchiveModel: Model file ".$self->filename()." not found!"); return $self->fail(); } #Copying the model file system("cp ".$self->filename()." ".$self->directory().$self->id().$self->version().".txt"); } =head3 PrintModelDataToFile Definition: (success/fail) = FIGMODELmodel->PrintModelDataToFile(); Description: This function uses the MFAToolkit to print out all of the compound and reaction data for the input model. Some of the data printed by the toolkit is calculated internally in the toolkit and not stored in any files, so this data can only be retrieved through this function. The LoadModel function for example would not give you this data. =cut sub PrintModelDataToFile { my($self) = @_; #Running the MFAToolkit on the model file my $OutputIndex = $self->figmodel()->filename(); my $Command = $self->config("MFAToolkit executable")->[0]." parameterfile ../Parameters/Printing.txt resetparameter output_folder ".$OutputIndex.'/ LoadCentralSystem "'.$self->filename().'"'; system($Command); #Copying the model file printed by the toolkit out of the output directory and into the model directory if (!-e $self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().".txt") { $self->figmodel()->error_message("New model file not created due to an error. Check that the input modelfile exists."); $self->figmodel()->cleardirectory($OutputIndex); return $self->fail(); } $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex."/".$self->id().$self->selected_version().'.txt" "'.$self->directory().$self->id().$self->selected_version().'Data.txt"'; system($Command); $Command = 'cp "'.$self->config("MFAToolkit output directory")->[0].$OutputIndex.'/ErrorLog0.txt" "'.$self->directory().'ModelErrors.txt"'; system($Command); $self->figmodel()->cleardirectory($OutputIndex); return $self->success(); } =head2 Analysis Functions =head3 run_microarray_analysis Definition: int::status = FIGMODEL->run_microarray_analysis(string::media,string::job id,string::gene calls); Description: Runs microarray analysis attempting to turn off genes that are inactive in the microarray =cut sub run_microarray_analysis { my ($self,$media,$label,$index,$genecall) = @_; $genecall =~ s/_/:/g; $genecall =~ s/\//;/g; my $uniqueFilename = $self->figmodel()->filename(); my $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($uniqueFilename,$self->id(),$media,["ProductionMFA","ShewenellaExperiment"],{"Microarray assertions" => $label.";".$index.";".$genecall,"MFASolver" => "CPLEX","Network output location" => "/scratch/"},"MicroarrayAnalysis-".$uniqueFilename.".txt",undef,$self->selected_version()); system($command); my $filename = $self->figmodel()->config("MFAToolkit output directory")->[0].$uniqueFilename."/MicroarrayOutput-".$index.".txt"; if (-e $filename) { my $output = $self->figmodel()->database()->load_single_column_file($filename); if (defined($output->[0])) { my @array = split(/;/,$output->[0]); $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt"); return ($array[0],$array[1],$array[8].":".$array[2],$array[9].":".$array[3],$array[10].":".$array[4],$array[11].":".$array[5],$array[12].":".$array[6],$array[13].":".$array[7]); } print STDERR $filename." is empty!"; } print STDERR $filename." not found!"; $self->figmodel()->clearing_output($uniqueFilename,"MicroarrayAnalysis-".$uniqueFilename.".txt"); return undef; } =head3 find_minimal_pathways Definition: int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective); Description: Runs microarray analysis attempting to turn off genes that are inactive in the microarray =cut sub find_minimal_pathways { my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_; #Setting default media if (!defined($media)) { $media = "Complete"; } #Setting default solution number if (!defined($solutionnum)) { $solutionnum = "5"; } #Setting additional exchange fluxes if (!defined($additionalexchange) || length($additionalexchange) == 0) { if ($self->id() eq "iAF1260") { $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0"; } else { $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0]; } } #Translating objective my $objectivestring; if ($objective eq "ALL") { #Getting the list of universal building blocks my $buildingblocks = $self->config("universal building blocks"); my @objectives = keys(%{$buildingblocks}); #Getting the nonuniversal building blocks my $otherbuildingblocks = $self->config("nonuniversal building blocks"); my @array = keys(%{$otherbuildingblocks}); if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) { my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0]; if (defined($equation)) { for (my $i=0; $i < @array; $i++) { if (CORE::index($equation,$array[$i]) > 0) { push(@objectives,$array[$i]); } } } } for (my $i=0; $i < @objectives; $i++) { $self->find_minimal_pathways($media,$objectives[$i]); } return; } elsif ($objective eq "ENERGY") { $objectivestring = "MAX;FLUX;rxn00062;c;1"; } elsif ($objective =~ m/cpd\d\d\d\d\d/) { if ($objective =~ m/\[(\w)\]/) { $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1"; $additionalexchange .= ";".$objective."[".$1."]:-100:0"; } else { $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1"; $additionalexchange .= ";".$objective."[c]:-100:0"; } } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) { my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective); for (my $i=0; $i < @{$Products};$i++) { my $temp = $Products->[$i]->{"DATABASE"}->[0]; if ($additionalexchange !~ m/$temp/) { #$additionalexchange .= ";".$temp."[c]:-100:0"; } } for (my $i=0; $i < @{$Reactants};$i++) { print $Reactants->[$i]->{"DATABASE"}->[0]." started\n"; $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange); print $Reactants->[$i]->{"DATABASE"}->[0]." done\n"; } return; } #Adding additional drains if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) { $additionalexchange .= ";cpd15666[c]:0:100"; } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) { $additionalexchange .= ";cpd12370[c]:0:100"; } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) { $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100"; } #Running MFAToolkit my $filename = $self->figmodel()->filename(); my $command; if (defined($AllReversible) && $AllReversible == 1) { $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"Make all reactions reversible in MFA"=>1, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "CPLEX","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version()); } else { $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"Make all reactions reversible in MFA"=>0, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "CPLEX","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version()); } system($command); #Loading problem report my $results = $self->figmodel()->LoadProblemReport($filename); #Clearing output $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt"); if (!defined($results)) { print STDERR $objective." pathway results not found!\n"; return; } #Parsing output my @Array; my $row = $results->get_row(1); if (defined($row->{"Notes"}->[0])) { $_ = $row->{"Notes"}->[0]; @Array = /\d+:([^\|]+)\|/g; } #Writing output to file $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]); } =head3 find_minimal_pathways Definition: int::status = FIGMODEL->find_minimal_pathways(string::media,string::objective); Description: Runs microarray analysis attempting to turn off genes that are inactive in the microarray =cut sub find_minimal_pathways_two { my ($self,$media,$objective,$solutionnum,$AllReversible,$additionalexchange) = @_; #Setting default media if (!defined($media)) { $media = "Complete"; } #Setting default solution number if (!defined($solutionnum)) { $solutionnum = "5"; } #Setting additional exchange fluxes if (!defined($additionalexchange) || length($additionalexchange) == 0) { if ($self->id() eq "iAF1260") { $additionalexchange = "cpd03422[c]:-100:100;cpd01997[c]:-100:100;cpd11416[c]:-100:0;cpd15378[c]:-100:0;cpd15486[c]:-100:0"; } else { $additionalexchange = $self->figmodel()->config("default exchange fluxes")->[0]; } } #Translating objective my $objectivestring; if ($objective eq "ALL") { #Getting the list of universal building blocks my $buildingblocks = $self->config("universal building blocks"); my @objectives = keys(%{$buildingblocks}); #Getting the nonuniversal building blocks my $otherbuildingblocks = $self->config("nonuniversal building blocks"); my @array = keys(%{$otherbuildingblocks}); if (defined($self->get_biomass()) && defined($self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0]))) { my $equation = $self->figmodel()->get_reaction($self->get_biomass()->{"LOAD"}->[0])->{"EQUATION"}->[0]; if (defined($equation)) { for (my $i=0; $i < @array; $i++) { if (CORE::index($equation,$array[$i]) > 0) { push(@objectives,$array[$i]); } } } } for (my $i=0; $i < @objectives; $i++) { $self->find_minimal_pathways($media,$objectives[$i]); } return; } elsif ($objective eq "ENERGY") { $objectivestring = "MAX;FLUX;rxn00062;c;1"; } elsif ($objective =~ m/cpd\d\d\d\d\d/) { if ($objective =~ m/\[(\w)\]/) { $objectivestring = "MIN;DRAIN_FLUX;".$objective.";".$1.";1"; $additionalexchange .= ";".$objective."[".$1."]:-100:0"; } else { $objectivestring = "MIN;DRAIN_FLUX;".$objective.";c;1"; $additionalexchange .= ";".$objective."[c]:-100:0"; } } elsif ($objective =~ m/(rxn\d\d\d\d\d)/) { my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($objective); for (my $i=0; $i < @{$Products};$i++) { my $temp = $Products->[$i]->{"DATABASE"}->[0]; if ($additionalexchange !~ m/$temp/) { #$additionalexchange .= ";".$temp."[c]:-100:0"; } } for (my $i=0; $i < @{$Reactants};$i++) { print $Reactants->[$i]->{"DATABASE"}->[0]." started\n"; $self->find_minimal_pathways($media,$Reactants->[$i]->{"DATABASE"}->[0],$additionalexchange); print $Reactants->[$i]->{"DATABASE"}->[0]." done\n"; } return; } #Adding additional drains if (($objective eq "cpd15665" || $objective eq "cpd15667" || $objective eq "cpd15668" || $objective eq "cpd15669") && $additionalexchange !~ m/cpd15666/) { $additionalexchange .= ";cpd15666[c]:0:100"; } elsif ($objective eq "cpd11493" && $additionalexchange !~ m/cpd12370/) { $additionalexchange .= ";cpd12370[c]:0:100"; } elsif ($objective eq "cpd00166" && $additionalexchange !~ m/cpd01997/) { $additionalexchange .= ";cpd01997[c]:0:100;cpd03422[c]:0:100"; } #Running MFAToolkit my $filename = $self->figmodel()->filename(); my $command; if (defined($AllReversible) && $AllReversible == 1) { $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"use simple variable and constraint names"=>1,"Make all reactions reversible in MFA"=>1, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "SCIP","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version()); } else { $command = $self->figmodel()->GenerateMFAToolkitCommandLineCall($filename,$self->id(),$media,["ProductionMFA"],{"use simple variable and constraint names"=>1,"Make all reactions reversible in MFA"=>0, "Recursive MILP solution limit" => $solutionnum,"Generate pathways to objective" => 1,"MFASolver" => "SCIP","objective" => $objectivestring,"exchange species" => $additionalexchange},"MinimalPathways-".$media."-".$self->id().$self->selected_version().".txt",undef,$self->selected_version()); } print $command."\n"; system($command); #Loading problem report my $results = $self->figmodel()->LoadProblemReport($filename); #Clearing output $self->figmodel()->clearing_output($filename,"MinimalPathways-".$media."-".$self->id()."-".$objective.".txt"); if (!defined($results)) { print STDERR $objective." pathway results not found!\n"; return; } #Parsing output my @Array; my $row = $results->get_row(1); if (defined($row->{"Notes"}->[0])) { $_ = $row->{"Notes"}->[0]; @Array = /\d+:([^\|]+)\|/g; } #Writing output to file $self->figmodel()->database()->print_array_to_file($self->directory()."MinimalPathways-".$media."-".$objective."-".$self->id()."-".$AllReversible."-".$self->selected_version().".txt",[join("|",@Array)]); } sub combine_minimal_pathways { my ($self) = @_; my $tbl; if (-e $self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl") { $tbl = FIGMODELTable::load_table($self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",";","|",0,["Objective","Media","Reversible"]); } else { $tbl = FIGMODELTable->new(["Objective","Media","Reactions","Reversible","Shortest path","Number of essentials","Essentials","Length"],$self->directory()."MinimalPathwayTable-".$self->id().$self->selected_version().".tbl",["Objective","Media","Reversible"],";","|"); } my @files = glob($self->directory()."MinimalPathways-*"); for (my $i=0; $i < @files;$i++) { if ($files[$i] =~ m/MinimalPathways\-(\S+)\-(cpd\d\d\d\d\d)\-(\w+)\-(\d)\-/ || $files[$i] =~ m/MinimalPathways\-(\S+)\-(ENERGY)\-(\w+)\-(\d)\-/) { my $reactions = $self->figmodel()->database()->load_single_column_file($files[$i],""); if (defined($reactions) && @{$reactions} > 0 && length($reactions->[0]) > 0) { my $newrow = {"Objective"=>[$2],"Media"=>[$1],"Reversible"=>[$4]}; my $row = $tbl->get_table_by_key($newrow->{"Objective"}->[0],"Objective")->get_table_by_key($newrow->{"Media"}->[0],"Media")->get_row_by_key($newrow->{"Reversible"}->[0],"Reversible"); if (!defined($row)) { $row = $tbl->add_row($newrow); } $row->{Reactions} = $self->figmodel()->database()->load_single_column_file($files[$i],""); delete($row->{"Shortest path"}); delete($row->{"Number of essentials"}); delete($row->{"Essentials"}); delete($row->{"Length"}); for (my $j=0; $j < @{$row->{Reactions}}; $j++) { my @array = split(/,/,$row->{Reactions}->[$j]); $row->{"Length"}->[$j] = @array; if (!defined($row->{"Shortest path"}->[0]) || $row->{"Length"}->[$j] < $row->{"Shortest path"}->[0]) { $row->{"Shortest path"}->[0] = $row->{"Length"}->[$j]; } $row->{"Number of essentials"}->[0] = 0; for (my $k=0; $k < @array;$k++) { if ($array[$k] =~ m/(rxn\d\d\d\d\d)/) { my $class = $self->get_reaction_class($1,1); my $temp = $row->{Media}->[0].":Essential"; if ($class =~ m/$temp/) { $row->{"Number of essentials"}->[$j]++; if (!defined($row->{"Essentials"}->[$j]) && length($row->{"Essentials"}->[$j]) > 0) { $row->{"Essentials"}->[$j] = $array[$k]; } else { $row->{"Essentials"}->[$j] .= ",".$array[$k]; } } } } } } } } $tbl->save(); } =head3 calculate_growth Definition: string::growth = FIGMODELmodel->calculate_growth(string:media); Description: Calculating growth in the input media =cut sub calculate_growth { my ($self,$Media) = @_; my $UniqueFilename = $self->figmodel()->filename(); 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())); my $ProblemReport = $self->figmodel()->LoadProblemReport($UniqueFilename); my $Result; if (defined($ProblemReport)) { my $Row = $ProblemReport->get_row(0); if (defined($Row) && defined($Row->{"Objective"}->[0])) { if ($Row->{"Objective"}->[0] < 0.00000001) { $Result = "NOGROWTH:".$Row->{"Individual metabolites with zero production"}->[0]; } else { $Result = $Row->{"Objective"}->[0]; } } } return $Result; } =head3 classify_model_reactions Definition: (FIGMODELTable:Reaction classes,FIGMODELTable:Compound classes) = FIGMODELmodel->classify_model_reactions(string:media); Description: This function uses the MFAToolkit to minimize and maximize the flux through every reaction in the input model during minimal growth on the input media. The results are returned in a hash of strings where the keys are the reaction IDs and the strings are structured as follows: "Class;Min flux;Max flux". Possible values for "Class" include: 1.) Positive: these reactions are essential in the forward direction. 2.) Negative: these reactions are essential in the reverse direction. 3.) Positive variable: these reactions are nonessential, but they only ever proceed in the forward direction. 4.) Negative variable: these reactions are nonessential, but they only ever proceed in the reverse direction. 5.) Variable: these reactions are nonessential and proceed in the forward or reverse direction. 6.) Blocked: these reactions never carry any flux at all in the media condition tested. 7.) Dead: these reactions are disconnected from the network. =cut sub classify_model_reactions { my ($self,$Media,$SaveChanges) = @_; #Getting unique file for printing model output my $UniqueFilename = $self->figmodel()->filename(); #Running the MFAToolkit system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),$Media,["ProductionMFA"],{"identify dead ends" => 1,"find tight bounds" => 1,"MFASolver" => "GLPK"},"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log",undef,$self->selected_version())); #Reading in the output bounds file my ($ReactionTB,$CompoundTB,$DeadCompounds,$DeadEndCompounds,$DeadReactions); if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt") { $ReactionTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsReactionData0.txt",";","|",1,["DATABASE ID"]); } if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt") { $CompoundTB = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/MFAOutput/TightBoundsCompoundData0.txt",";","|",1,["DATABASE ID"]); } if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt") { $DeadReactions = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadReactions.txt","")); } if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt") { $DeadCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadMetabolites.txt","")); } if (-e $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt") { $DeadEndCompounds = $self->figmodel()->put_array_in_hash($self->figmodel()->database()->load_single_column_file($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/DeadEndMetabolites.txt","")); } if (!defined($ReactionTB) && !defined($CompoundTB)) { 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"; return (undef,undef); } #Clearing output $self->figmodel()->clearing_output($UniqueFilename,"Classify-".$self->id().$self->selected_version()."-".$UniqueFilename.".log"); #Creating the table objects that will hold the results of the reaction classification my $rxnclasstable = $self->reaction_class_table(); my $cpdclasstable = $self->compound_class_table(); #Loading the compound table if (defined($CompoundTB)) { for (my $i=0; $i < $CompoundTB->size(); $i++) { my $Row = $CompoundTB->get_row($i); if (defined($Row->{"DATABASE ID"})) { #Getting the compound row my $CpdRow = $cpdclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0].$Row->{COMPARTMENT}->[0],"COMPOUND",1); #Setting row values my $Max = 0; my $Min = 0; my $Class = "Unknown"; if (defined($DeadCompounds) && defined($DeadCompounds->{$Row->{"DATABASE ID"}->[0]})) { $Class = "Dead"; } elsif (defined($DeadEndCompounds) && defined($DeadEndCompounds->{$Row->{"DATABASE ID"}->[0]})) { $Class = "Deadend"; } elsif (defined($Row->{"Min DRAIN_FLUX"}) && defined($Row->{"Max DRAIN_FLUX"}) && $Row->{"Min DRAIN_FLUX"}->[0] ne "1e+07") { $Max = $Row->{"Max DRAIN_FLUX"}->[0]; $Min = $Row->{"Min DRAIN_FLUX"}->[0]; if ($Row->{"Min DRAIN_FLUX"}->[0] > 0.00000001) { $Class = "Positive"; } elsif ($Row->{"Max DRAIN_FLUX"}->[0] < -0.00000001) { $Class = "Negative"; } elsif ($Row->{"Min DRAIN_FLUX"}->[0] < -0.00000001) { if ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) { $Class = "Variable"; } else { $Max = 0; $Class = "Negative variable"; } } elsif ($Row->{"Max DRAIN_FLUX"}->[0] > 0.00000001) { $Min = 0; $Class = "Positive variable"; } else { $Min = 0; $Max = 0; $Class = "Blocked"; } } my $index = 0; if (defined($CpdRow->{MEDIA})) { for (my $i=0; $i < @{$CpdRow->{MEDIA}};$i++) { $index++; if ($CpdRow->{MEDIA}->[$i] eq $Media) { $index = $i; last; } } } $CpdRow->{MIN}->[$index] = $Min; $CpdRow->{MAX}->[$index] = $Max; $CpdRow->{CLASS}->[$index] = $Class; $CpdRow->{MEDIA}->[$index] = $Media; } } if (!defined($SaveChanges) || $SaveChanges == 1) { $cpdclasstable->save(); } } if (defined($ReactionTB)) { for (my $i=0; $i < $ReactionTB->size(); $i++) { my $Row = $ReactionTB->get_row($i); if (defined($Row->{"DATABASE ID"})) { #Getting the compound row my $Compartment = "c"; if (defined($Row->{COMPARTMENT}->[0])) { $Compartment = $Row->{COMPARTMENT}->[0]; } my $RxnRow = $rxnclasstable->get_row_by_key($Row->{"DATABASE ID"}->[0],"REACTION",1); my $Max = 0; my $Min = 0; my $Class = "Unknown"; if (defined($DeadReactions) && defined($DeadReactions->{$Row->{"DATABASE ID"}->[0]})) { $Class = "Dead"; } elsif (defined($Row->{"Min FLUX"}) && defined($Row->{"Max FLUX"})) { $Max = $Row->{"Max FLUX"}->[0]; $Min = $Row->{"Min FLUX"}->[0]; if ($Row->{"Min FLUX"}->[0] > 0.00000001) { $Class = "Positive"; } elsif ($Row->{"Max FLUX"}->[0] < -0.00000001) { $Class = "Negative"; } elsif ($Row->{"Min FLUX"}->[0] < -0.00000001) { if ($Row->{"Max FLUX"}->[0] > 0.00000001) { $Class = "Variable"; } else { $Max = 0; $Class = "Negative variable"; } } elsif ($Row->{"Max FLUX"}->[0] > 0.00000001) { $Min = 0; $Class = "Positive variable"; } else { $Min = 0; $Max = 0; $Class = "Blocked"; } } my $index = 0; if (defined($RxnRow->{MEDIA})) { for (my $i=0; $i < @{$RxnRow->{MEDIA}};$i++) { $index++; if ($RxnRow->{MEDIA}->[$i] eq $Media) { $index = $i; last; } } } $RxnRow->{MIN}->[$index] = $Min; $RxnRow->{MAX}->[$index] = $Max; $RxnRow->{CLASS}->[$index] = $Class; $RxnRow->{MEDIA}->[$index] = $Media; } } if (!defined($SaveChanges) || $SaveChanges == 1) { $rxnclasstable->save(); } } return ($rxnclasstable,$cpdclasstable); } =head3 RunAllStudiesWithDataFast Definition: (integer::false positives,integer::false negatives,integer::correct negatives,integer::correct positives,string::error vector,string heading vector) = FIGMODELmodel->RunAllStudiesWithDataFast(string::experiment,0/1::print result); Description: Simulates every experimental condition currently available for the model. =cut sub RunAllStudiesWithDataFast { my ($self,$Experiment,$PrintResults) = @_; #Printing lp and key file for model if (!-e $self->directory()."FBA-".$self->id().$self->selected_version().".lp") { $self->PrintModelLPFile(); } my $UniqueFilename = $self->figmodel()->filename(); #Determing the simulations that need to be run my $ExperimentalDataTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment); #Creating the table of jobs to submit my $JobArray = $self->GetSimulationJobTable($ExperimentalDataTable,$Experiment,$UniqueFilename); #Printing the job file if (!-d $self->config("MFAToolkit output directory")->[0].$UniqueFilename."/") { system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/"); } $JobArray->save(); #Running simulations system($self->config("mfalite executable")->[0]." ".$self->config("Reaction database directory")->[0]."masterfiles/MediaTable.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Jobfile.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt"); #Parsing the results my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef); if (!defined($Results)) { $self->figmodel()->error_message("FIGMODELmodel:RunAllStudiesWithDataFast:Could not find simulation results: ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt"); return undef; } my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector,$SimulationResults) = $self->EvaluateSimulationResults($Results,$ExperimentalDataTable); #Printing results to file $self->figmodel()->database()->save_table($SimulationResults,undef,undef,undef,"False negatives\tFalse positives\tCorrect negatives\tCorrect positives\n".$FalseNegatives."\t".$FalsePostives."\t".$CorrectNegatives."\t".$CorrectPositives."\n"); $self->figmodel()->clearing_output($UniqueFilename); return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector); } =head3 GetSimulationJobTable Definition: my $JobTable = $model->GetSimulationJobTable($Experiment,$PrintResults,$Version); Description: =cut sub GetSimulationJobTable { my ($self,$SimulationTable,$Experiment,$Folder) = @_; #Determing the simulations that need to be run if (!defined($SimulationTable)) { $SimulationTable = $self->figmodel()->GetExperimentalDataTable($self->genome(),$Experiment); if (!defined($SimulationTable)) { return undef; } } #Creating the job table my $JobTable = $self->figmodel()->CreateJobTable($Folder); for (my $i=0; $i < $SimulationTable->size(); $i++) { if ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Gene\sKO/) { my $Row = $JobTable->get_row_by_key("Gene KO","LABEL",1); $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1); } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Media\sgrowth/) { my $Row = $JobTable->get_row_by_key("Growth phenotype","LABEL",1); $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1); } elsif ($SimulationTable->get_row($i)->{"Heading"}->[0] =~ m/Interval\sKO/) { my $Row = $JobTable->get_row_by_key($SimulationTable->get_row($i)->{"Heading"}->[0],"LABEL",1); $JobTable->add_data($Row,"MEDIA",$SimulationTable->get_row($i)->{"Media"}->[0],1); $JobTable->add_data($Row,"GENE KO",$SimulationTable->get_row($i)->{"Experiment type"}->[0],1); } } #Filling in model specific elements of the job table for (my $i=0; $i < $JobTable->size(); $i++) { if ($JobTable->get_row($i)->{"LABEL"}->[0] =~ m/Gene\sKO/) { $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "SINGLEKO"; $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 1; } else { $JobTable->get_row($i)->{"RUNTYPE"}->[0] = "GROWTH"; $JobTable->get_row($i)->{"SAVE NONESSENTIALS"}->[0] = 0; } $JobTable->get_row($i)->{"LP FILE"}->[0] = $self->directory()."FBA-".$self->id().$self->selected_version(); $JobTable->get_row($i)->{"MODEL"}->[0] = $self->directory().$self->id().$self->selected_version().".txt"; $JobTable->get_row($i)->{"SAVE FLUXES"}->[0] = 0; } return $JobTable; } =head3 EvaluateSimulationResults Definition: (integer::false positives,integer::false negatives,integer::correct negatives,integer::correct positives,string::error vector,string heading vector,FIGMODELtable::simulation results) = FIGMODELmodel->EvaluateSimulationResults(FIGMODELtable::raw simulation results,FIGMODELtable::experimental data); Description: Compares simulation results with experimental data to produce a table indicating where predictions are incorrect. =cut sub EvaluateSimulationResults { my ($self,$Results,$ExperimentalDataTable) = @_; #Comparing experimental results with simulation results my $SimulationResults = FIGMODELTable->new(["Run result","Experiment type","Media","Experiment ID","Reactions knocked out"],$self->directory()."SimulationOutput".$self->id().$self->selected_version().".txt",["Experiment ID","Media"],"\t",",",undef); my $FalsePostives = 0; my $FalseNegatives = 0; my $CorrectNegatives = 0; my $CorrectPositives = 0; my @Errorvector; my @HeadingVector; my $ReactionKOWithGeneHash; for (my $i=0; $i < $Results->size(); $i++) { if ($Results->get_row($i)->{"LABEL"}->[0] eq "Gene KO") { if (defined($Results->get_row($i)->{"REACTION KO WITH GENES"})) { for (my $j=0; $j < @{$Results->get_row($i)->{"REACTION KO WITH GENES"}}; $j++) { my @Temp = split(/:/,$Results->get_row($i)->{"REACTION KO WITH GENES"}->[$j]); if (defined($Temp[1]) && length($Temp[1]) > 0) { $ReactionKOWithGeneHash->{$Temp[0]} = $Temp[1]; } } } if ($Results->get_row($i)->{"OBJECTIVE"}->[0] == 0) { for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) { my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading"); if (defined($Row)) { my $KOReactions = "none"; if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) { $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]}; } push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions); my $Status = "Unknown"; if ($Row->{"Growth"}->[0] > 0) { $Status = "False negative"; $FalseNegatives++; push(@Errorvector,3); } else { $Status = "False positive"; $FalsePostives++; push(@Errorvector,2); } $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]}); } } } else { for (my $j=0; $j < @{$Results->get_row($i)->{"ESSENTIALGENES"}}; $j++) { #print $j."\t".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]."\n"; my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"ESSENTIALGENES"}->[$j],"Heading"); if (defined($Row)) { my $KOReactions = "none"; if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]})) { $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"ESSENTIALGENES"}->[$j]}; } push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions); my $Status = "Unknown"; if ($Row->{"Growth"}->[0] > 0) { $Status = "False negative"; $FalseNegatives++; push(@Errorvector,3); } else { $Status = "Correct negative"; $CorrectNegatives++; push(@Errorvector,1); } $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]}); } } for (my $j=0; $j < @{$Results->get_row($i)->{"NONESSENTIALGENES"}}; $j++) { my $Row = $ExperimentalDataTable->get_row_by_key("Gene KO:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j],"Heading"); if (defined($Row)) { my $KOReactions = "none"; if (defined($ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]})) { $KOReactions = $ReactionKOWithGeneHash->{$Results->get_row($i)->{"NONESSENTIALGENES"}->[$j]}; } push(@HeadingVector,$Row->{"Heading"}->[0].":".$KOReactions); my $Status = "Unknown"; if ($Row->{"Growth"}->[0] > 0) { $Status = "Correct positive"; $CorrectPositives++; push(@Errorvector,0); } else { $Status = "False positive"; $FalsePostives++; push(@Errorvector,2); } $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Gene KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => [$KOReactions]}); } } } } elsif ($Results->get_row($i)->{"LABEL"}->[0] eq "Growth phenotype") { my $Row = $ExperimentalDataTable->get_row_by_key("Media growth:".$Results->get_row($i)->{"MEDIA"}->[0].":".$Results->get_row($i)->{"MEDIA"}->[0],"Heading"); if (defined($Row)) { push(@HeadingVector,$Row->{"Heading"}->[0].":none"); my $Status = "Unknown"; if ($Row->{"Growth"}->[0] > 0) { if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) { $Status = "Correct positive"; $CorrectPositives++; push(@Errorvector,0); } else { $Status = "False negative"; $FalseNegatives++; push(@Errorvector,3); } } else { if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) { $Status = "False positive"; $FalsePostives++; push(@Errorvector,2); } else { $Status = "Correct negative"; $CorrectNegatives++; push(@Errorvector,1); } } $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Media growth"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Media"}->[0]],"Reactions knocked out" => ["none"]}); } } elsif ($Results->get_row($i)->{"LABEL"}->[0] =~ m/Interval\sKO/ && defined($Results->get_row($i)->{"KOGENES"}->[0])) { my $Row = $ExperimentalDataTable->get_row_by_key($Results->get_row($i)->{"LABEL"}->[0],"Heading"); if (defined($Row)) { my $Status = "Unknown"; if ($Row->{"Growth"}->[0] > 0) { if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) { $Status = "Correct positive"; $CorrectPositives++; push(@Errorvector,0); } else { $Status = "False negative"; $FalseNegatives++; push(@Errorvector,3); } } else { if ($Results->get_row($i)->{"OBJECTIVE"}->[0] > 0) { $Status = "False positive"; $FalsePostives++; push(@Errorvector,2); } else { $Status = "Correct negative"; $CorrectNegatives++; push(@Errorvector,1); } } $SimulationResults->add_row({"Run result" => [$Status],"Experiment type" => ["Interval KO"],"Media" => [$Row->{"Media"}->[0]],"Experiment ID" => [$Row->{"Experiment ID"}->[0]],"Reactions knocked out" => ["none"]}); } } } return ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,join(";",@Errorvector),join(";",@HeadingVector),$SimulationResults); } =head3 InspectSolution Definition: $model->InspectSolution(string::gene knocked out,string::media condition,[string]::list of reactions); Description: =cut sub InspectSolution { my ($self,$GeneKO,$Media,$ReactionList) = @_; #Getting a directory for the results my $UniqueFilename = $self->figmodel()->filename(); system("mkdir ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/"); my $TempVersion = "V".$UniqueFilename; #Setting gene ko to none if no genes are to be knocked out if ($GeneKO !~ m/^peg\./) { $GeneKO = "none"; } #Implementing the input solution in the test model my $ReactionArray; my $DirectionArray; my %SolutionHash; for (my $k=0; $k < @{$ReactionList}; $k++) { if ($ReactionList->[$k] =~ m/(.+)(rxn\d\d\d\d\d)/) { my $Reaction = $2; my $Sign = $1; if (defined($SolutionHash{$Reaction})) { $SolutionHash{$Reaction} = "<=>"; } elsif ($Sign eq "-") { $SolutionHash{$Reaction} = "<="; } elsif ($Sign eq "+") { $SolutionHash{$Reaction} = "=>"; } else { $SolutionHash{$Reaction} = $Sign; } } } my @TempList = keys(%SolutionHash); for (my $k=0; $k < @TempList; $k++) { push(@{$ReactionArray},$TempList[$k]); push(@{$DirectionArray},$SolutionHash{$TempList[$k]}); } print "Integrating solution!\n"; $self->figmodel()->IntegrateGrowMatchSolution($self->id().$self->selected_version(),$self->directory().$self->id().$TempVersion.".txt",$ReactionArray,$DirectionArray,"SolutionInspection",1,1); #Printing lp and key file for model $self->PrintModelLPFile(); #Running FBA on the test model my $JobTable = $self->figmodel()->CreateJobTable($UniqueFilename); $JobTable->add_row({"LABEL" => ["TEST"],"RUNTYPE" => ["GROWTH"],"LP FILE" => [$self->directory()."FBA-".$self->id().$TempVersion],"MODEL" => [$self->directory().$self->id().$TempVersion.".txt"],"MEDIA" => [$Media],"REACTION KO" => ["none|".join("|",@{$ReactionList})],"GENE KO" => [$GeneKO],"SAVE FLUXES" => [0],"SAVE NONESSENTIALS" => [0]}); $JobTable->save(); #Running simulations system($self->config("mfalite executable")->[0]." ".$self->config("Reaction database directory")->[0]."masterfiles/MediaTable.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Jobfile.txt ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt"); #Parsing the results my $Results = $self->figmodel()->database()->load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt",";","\\|",0,undef); if (!defined($Results)) { $self->figmodel()->error_message("FIGMODELmodel:InspectSolution:Could not load problem report ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/Output.txt"); return undef; } #Making sure that the model grew with all reactions present my $Found = 0; for (my $i=0; $i < $Results->size(); $i++) { if (defined($Results->get_row($i)->{"KOGENES"}->[0]) && defined($Results->get_row($i)->{"KOREACTIONS"}->[0]) && $Results->get_row($i)->{"KOREACTIONS"}->[0] eq "none" && $Results->get_row($i)->{"KOGENES"}->[0] eq $GeneKO && $Results->get_row($i)->{"OBJECTIVE"}->[0] > 0.00001) { $Found = 1; } } if ($Found == 0) { print "Solution no longer valid\n"; return undef; } #Making sure all of the reactions added are still necessary my $FinalReactionList; for (my $k=0; $k < $Results->size(); $k++) { if (defined($Results->get_row($k)->{"KOGENES"}->[0]) && $Results->get_row($k)->{"KOGENES"}->[0] eq $GeneKO) { if (defined($Results->get_row($k)->{"KOREACTIONS"}->[0]) && $Results->get_row($k)->{"KOREACTIONS"}->[0] =~ m/rxn\d\d\d\d\d/ && $Results->get_row($k)->{"OBJECTIVE"}->[0] < 0.000001) { push(@{$FinalReactionList},$Results->get_row($k)->{"KOREACTIONS"}->[0]); } } } #Deleting extra files created unlink($self->directory()."FBA-".$self->id().$TempVersion.".lp"); unlink($self->directory()."FBA-".$self->id().$TempVersion.".key"); unlink($self->directory().$self->id().$TempVersion.".txt"); #Deleting the test model and the MFA folder $self->figmodel()->clearing_output($UniqueFilename); return $FinalReactionList; } =head3 GapFillingAlgorithm Definition: FIGMODELmodel->GapFillingAlgorithm(); Description: This is a wrapper for running the gap filling algorithm on any model in the database. The algorithm performs a gap filling for any false negative prediction of the avialable experimental data. This function is threaded to improve efficiency: one thread does nothing but using the MFAToolkit to fill gaps for every false negative prediction. The other thread reads in the gap filling solutions, builds a test model for each solution, and runs the test model against all available experimental data. This function prints two important output files in the Model directory: 1.) GapFillingOutput.txt: this is a summary of the results of the gap filling analysis 2.) GapFillingErrorMatrix.txt: this lists the correct and incorrect predictions for each gapfilling solution implemented in a test model. =cut sub GapFillingAlgorithm { my ($self) = @_; #First the input model version and model filename should be simulated and the false negatives identified my ($FalsePostives,$FalseNegatives,$CorrectNegatives,$CorrectPositives,$Errorvector,$HeadingVector) = $self->RunAllStudiesWithDataFast("All"); #Getting the filename my $UniqueFilename = $self->figmodel()->filename(); #Printing the original performance vector $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-OPEM".".txt",[$HeadingVector,$Errorvector]); my $PreviousGapFilling; if (-e $self->directory().$self->id().$self->selected_version()."-GFS.txt") { #Backing up the old solution file system("cp ".$self->directory().$self->id().$self->selected_version()."-GFS.txt ".$self->directory().$self->id().$self->selected_version()."-OldGFS.txt"); unlink($self->directory().$self->id().$self->selected_version()."-GFS.txt"); } if (-e $self->directory().$self->id().$self->selected_version()."-OldGFS.txt") { #Reading in the solution file from the previous gap filling if it exists $PreviousGapFilling = $self->figmodel()->database()->load_table($self->directory().$self->id().$self->selected_version()."-OldGFS.txt",";",",",0,["Experiment"]); } #Now we use the simulation output to make the gap filling run data my @Errors = split(/;/,$Errorvector); my @Headings = split(/;/,$HeadingVector); my $GapFillingRunSpecs = ""; my $Count = 0; my $RescuedPreviousResults; my $RunCount = 0; my $SolutionExistedCount = 0; my $AcceptedSolutions = 0; my $RejectedSolutions = 0; my $NoExistingSolutions = 0; for (my $i=0; $i < @Errors; $i++) { if ($Errors[$i] == 3) { my @HeadingDataArray = split(/:/,$Headings[$i]); if ($HeadingDataArray[2] !~ m/^peg\./ || $HeadingDataArray[3] ne "none") { my $SolutionFound = 0; if (defined($PreviousGapFilling) && defined($PreviousGapFilling->get_row_by_key($HeadingDataArray[2],"Experiment"))) { my @Rows = $PreviousGapFilling->get_rows_by_key($HeadingDataArray[2],"Experiment"); for (my $j=0; $j < @Rows; $j++) { if ($HeadingDataArray[2] =~ m/^peg\./) { my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"}); if (defined($ReactionList)) { print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n"; $SolutionFound++; push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList})); $AcceptedSolutions++; } else { $RejectedSolutions++; } } else { my $ReactionList = $self->InspectSolution($HeadingDataArray[2],$HeadingDataArray[1],$Rows[$j]->{"Solution reactions"}); if (defined($ReactionList)) { print join(",",@{$Rows[$j]->{"Solution reactions"}})."\t".join(",",@{$ReactionList})."\n"; $SolutionFound++; push(@{$RescuedPreviousResults},$Rows[$j]->{"Experiment"}->[0].";".$Rows[$j]->{"Solution index"}->[0].";".$Rows[$j]->{"Solution cost"}->[0].";".join(",",@{$ReactionList})); $AcceptedSolutions++; } else { $RejectedSolutions++; } } } } else { $NoExistingSolutions++; } if ($SolutionFound == 0) { $RunCount++; if (length($GapFillingRunSpecs) > 0) { $GapFillingRunSpecs .= ";"; } $GapFillingRunSpecs .= $HeadingDataArray[2].":".$HeadingDataArray[1].":".$HeadingDataArray[3]; } else { $SolutionExistedCount++; } } $Count++; } } #Updating the growmatch progress table my $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1); $Row->{"INITIAL FP"}->[0] = $FalsePostives; $Row->{"INITIAL FN"}->[0] = $FalseNegatives; $Row->{"GF TIMING"}->[0] = time()."-"; $Row->{"FN WITH SOL"}->[0] = $FalseNegatives-$NoExistingSolutions; $Row->{"FN WITH ACCEPTED SOL"}->[0] = $SolutionExistedCount; $Row->{"TOTAL ACCEPTED GF SOL"}->[0] = $AcceptedSolutions; $Row->{"TOTAL REJECTED GF SOL"}->[0] = $RejectedSolutions; $Row->{"FN WITH NO SOL"}->[0] = $NoExistingSolutions+$RejectedSolutions; $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM"); #Running the gap filling once to correct all false negative errors my $SolutionsFound = 0; my $GapFillingArray; push(@{$GapFillingArray},split(/;/,$GapFillingRunSpecs)); my $GapFillingResults = $self->datagapfill($GapFillingArray,"GFS"); if (defined($GapFillingResults)) { $SolutionsFound = 1; } if (defined($RescuedPreviousResults) && @{$RescuedPreviousResults} > 0) { #Printing previous solutions to GFS file $self->figmodel()->database()->print_array_to_file($self->directory().$self->id().$self->selected_version()."-GFS.txt",$RescuedPreviousResults,1); $SolutionsFound = 1; } #Recording the finishing of the gapfilling $Row = $self->figmodel()->database()->get_row_by_key("GROWMATCH TABLE",$self->genome(),"ORGANISM",1); $Row->{"GF TIMING"}->[0] .= time(); $self->figmodel()->database()->update_row("GROWMATCH TABLE",$Row,"ORGANISM"); if ($SolutionsFound == 1) { #Scheduling solution testing $self->figmodel()->add_job_to_queue("testsolutions?".$self->id().$self->selected_version()."?-1?GF","QSUB","fast","master","BACK"); } else { $self->figmodel()->error_message("No false negative predictions found. Data gap filling not necessary!"); } return $self->success(); } =head3 SolutionReconciliation Definition: FIGMODELmodel->SolutionReconciliation(); Description: This is a wrapper for running the solution reconciliation algorithm on any model in the database. The algorithm performs a reconciliation of any gap filling solutions to identify the combination of solutions that results in the optimal model. This function prints out one output file in the Model directory: ReconciliationOutput.txt: this is a summary of the results of the reconciliation analysis =cut sub SolutionReconciliation { my ($self,$GapFill,$Stage) = @_; #Setting the output filenames my $OutputFilename; my $OutputFilenameTwo; if ($GapFill == 1) { $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GFReconciliation.txt"; $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GFSRS.txt"; } else { $OutputFilename = $self->directory().$self->id().$self->selected_version()."-GGReconciliation.txt"; $OutputFilenameTwo = $self->directory().$self->id().$self->selected_version()."-GGSRS.txt"; } #In stage one, we run the reconciliation and create a test file to check combined solution performance if (!defined($Stage) || $Stage == 1) { my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE"); my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1); $Row->{"GF RECONCILATION TIMING"}->[0] = time()."-"; $GrowMatchTable->save(); $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE"); #Getting a unique filename my $UniqueFilename = $self->figmodel()->filename(); #Copying over the necessary files if ($GapFill == 1) { if (!-e $self->directory().$self->id().$self->selected_version()."-GFEM.txt") { print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GFEM.txt file not found. Could not reconcile!"; return 0; } if (!-e $self->directory().$self->id().$self->selected_version()."-OPEM.txt") { print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-OPEM.txt file not found. Could not reconcile!"; return 0; } system("cp ".$self->directory().$self->id().$self->selected_version()."-GFEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GFEM.txt"); system("cp ".$self->directory().$self->id().$self->selected_version()."-OPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt"); #Backing up and deleting the existing reconciliation file if (-e $OutputFilename) { system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGFReconciliation.txt"); unlink($OutputFilename); } } else { if (!-e $self->directory().$self->id().$self->selected_version()."-GGEM.txt") { print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGEM.txt file not found. Could not reconcile!"; return 0; } if (!-e $self->directory().$self->id().$self->selected_version()."-GGOPEM.txt") { print STDERR "FIGMODEL:SolutionReconciliation:".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt file not found. Could not reconcile!"; return 0; } system("cp ".$self->directory().$self->id().$self->selected_version()."-GGEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-GGEM.txt"); system("cp ".$self->directory().$self->id().$self->selected_version()."-GGOPEM.txt ".$self->figmodel()->config("MFAToolkit input files")->[0].$UniqueFilename."-OPEM.txt"); #Backing up and deleting the existing reconciliation file if (-e $OutputFilename) { system("cp ".$OutputFilename." ".$self->directory().$self->id().$self->selected_version()."-OldGGReconciliation.txt"); unlink($OutputFilename); } } #Running the reconciliation system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),"NONE",["SolutionReconciliation"],{"Solution data for model optimization" => $UniqueFilename},"Reconciliation".$UniqueFilename.".log",undef,$self->selected_version())); $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE"); $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1); $Row->{"GF RECONCILATION TIMING"}->[0] .= time(); $GrowMatchTable->save(); $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE"); #Loading the problem report from the reconciliation run my $ReconciliatonOutput = $self->figmodel()->LoadProblemReport($UniqueFilename); print $UniqueFilename."\n"; #Clearing output files $self->figmodel()->clearing_output($UniqueFilename,"Reconciliation".$UniqueFilename.".log"); $ReconciliatonOutput->save("/home/chenry/Test.txt"); #Checking the a problem report was found and was loaded if (!defined($ReconciliatonOutput) || $ReconciliatonOutput->size() < 1 || !defined($ReconciliatonOutput->get_row(0)->{"Notes"}->[0])) { print STDERR "FIGMODEL:SolutionReconciliation: MFAToolkit output from SolutionReconciliation of ".$self->id()." not found!\n\n"; return 0; } #Processing the solutions my $SolutionCount = 0; my $ReactionSetHash; my $SingleReactionHash; my $ReactionDataHash; for (my $n=0; $n < $ReconciliatonOutput->size(); $n++) { if (defined($ReconciliatonOutput->get_row($n)->{"Notes"}->[0]) && $ReconciliatonOutput->get_row($n)->{"Notes"}->[0] =~ m/^Recursive\sMILP\s([^;]+)/) { #Breaking up the solution into reaction sets my @ReactionSets = split(/\|/,$1); #Creating reaction lists for each set my $SolutionHash; for (my $i=0; $i < @ReactionSets; $i++) { if (length($ReactionSets[$i]) > 0) { my @Alternatives = split(/:/,$ReactionSets[$i]); for (my $j=1; $j < @Alternatives; $j++) { if (length($Alternatives[$j]) > 0) { push(@{$SolutionHash->{$Alternatives[$j]}},$Alternatives[0]); } } if (@Alternatives == 1) { $SingleReactionHash->{$Alternatives[0]}->{$SolutionCount} = 1; if (!defined($SingleReactionHash->{$Alternatives[0]}->{"COUNT"})) { $SingleReactionHash->{$Alternatives[0]}->{"COUNT"} = 0; } $SingleReactionHash->{$Alternatives[0]}->{"COUNT"}++; } } } #Identifying reactions sets and storing the sets in the reactions set hash foreach my $Solution (keys(%{$SolutionHash})) { my $SetKey = join(",",sort(@{$SolutionHash->{$Solution}})); if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount})) { $ReactionSetHash->{$SetKey}->{$SetKey}->{$SolutionCount} = 1; if (!defined($ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"})) { $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"} = 0; } $ReactionSetHash->{$SetKey}->{$SetKey}->{"COUNT"}++; } $ReactionSetHash->{$SetKey}->{$Solution}->{$SolutionCount} = 1; if (!defined($ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"})) { $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"} = 0; } $ReactionSetHash->{$SetKey}->{$Solution}->{"COUNT"}++; } $SolutionCount++; } } #Handling the scenario where no solutions were found if ($SolutionCount == 0) { print STDERR "FIGMODEL:SolutionReconciliation: Reconciliation unsuccessful. No solution found.\n\n"; return 0; } #Printing results without solution performance figures. Also printing solution test file open (RECONCILIATION, ">$OutputFilename"); #Printing the file heading print RECONCILIATION "DATABASE;DEFINITION;REVERSIBLITY;DELTAG;DIRECTION;NUMBER OF SOLUTIONS"; for (my $i=0; $i < $SolutionCount; $i++) { print RECONCILIATION ";Solution ".$i; } print RECONCILIATION "\n"; #Printing the singlet reactions first my $Solutions; print RECONCILIATION "SINGLET REACTIONS\n"; my @SingletReactions = keys(%{$SingleReactionHash}); for (my $j=0; $j < $SolutionCount; $j++) { $Solutions->[$j]->{"BASE"} = $j; } for (my $i=0; $i < @SingletReactions; $i++) { my $ReactionData; if (defined($ReactionDataHash->{$SingletReactions[$i]})) { $ReactionData = $ReactionDataHash->{$SingletReactions[$i]}; } else { my $Direction = substr($SingletReactions[$i],0,1); if ($Direction eq "+") { $Direction = "=>"; } else { $Direction = "<="; } my $Reaction = substr($SingletReactions[$i],1); $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t"); $ReactionData->{"DIRECTIONS"}->[0] = $Direction; $ReactionData->{"REACTIONS"}->[0] = $Reaction; if (!defined($ReactionData->{"DEFINITION"}->[0])) { $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN"; } if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) { $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN"; } if (!defined($ReactionData->{"DELTAG"}->[0])) { $ReactionData->{"DELTAG"}->[0] = "UNKNOWN"; } $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData; } print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$SingleReactionHash->{$SingletReactions[$i]}->{"COUNT"}; for (my $j=0; $j < $SolutionCount; $j++) { print RECONCILIATION ";"; if (defined($SingleReactionHash->{$SingletReactions[$i]}->{$j})) { $Solutions->[$j]->{$SingletReactions[$i]} = 1; $Solutions->[$j]->{"BASE"} = $j; print RECONCILIATION "|".$j."|"; } } print RECONCILIATION "\n"; } #Printing the reaction sets with alternatives print RECONCILIATION "Reaction sets with alternatives\n"; my @ReactionSets = keys(%{$ReactionSetHash}); foreach my $ReactionSet (@ReactionSets) { my $NewSolutions; my $BaseReactions; my $AltList = [$ReactionSet]; push(@{$AltList},keys(%{$ReactionSetHash->{$ReactionSet}})); for (my $j=0; $j < @{$AltList}; $j++) { my $CurrentNewSolutions; my $Index; if ($j == 0) { print RECONCILIATION "NEW SET\n"; } elsif ($AltList->[$j] ne $ReactionSet) { print RECONCILIATION "ALTERNATIVE SET\n"; #For each base solution in which this set is represented, we copy the base solution to the new solution my $NewSolutionCount = 0; for (my $k=0; $k < $SolutionCount; $k++) { if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) { if (defined($Solutions)) { $Index->{$k} = @{$Solutions} + $NewSolutionCount; } else { $Index->{$k} = $NewSolutionCount; } if (defined($NewSolutions) && @{$NewSolutions} > 0) { $Index->{$k} += @{$NewSolutions}; } $CurrentNewSolutions->[$NewSolutionCount] = {}; foreach my $Reaction (keys(%{$Solutions->[$k]})) { $CurrentNewSolutions->[$NewSolutionCount]->{$Reaction} = $Solutions->[$k]->{$Reaction}; } $NewSolutionCount++; } } } if ($j == 0 || $AltList->[$j] ne $ReactionSet) { my @SingletReactions = split(/,/,$AltList->[$j]); for (my $i=0; $i < @SingletReactions; $i++) { #Adding base reactions to base solutions and set reactions the new solutions if ($j == 0) { push(@{$BaseReactions},$SingletReactions[$i]); } else { for (my $k=0; $k < @{$CurrentNewSolutions}; $k++) { $CurrentNewSolutions->[$k]->{$SingletReactions[$i]} = 1; } } #Getting reaction data and printing reaction in output file my $ReactionData; if (defined($ReactionDataHash->{$SingletReactions[$i]})) { $ReactionData = $ReactionDataHash->{$SingletReactions[$i]}; } else { my $Direction = substr($SingletReactions[$i],0,1); if ($Direction eq "+") { $Direction = "=>"; } else { $Direction = "<="; } my $Reaction = substr($SingletReactions[$i],1); $ReactionData = FIGMODELObject->load($self->figmodel()->config("reaction directory")->[0].$Reaction,"\t"); $ReactionData->{"DIRECTIONS"}->[0] = $Direction; $ReactionData->{"REACTIONS"}->[0] = $Reaction; if (!defined($ReactionData->{"DEFINITION"}->[0])) { $ReactionData->{"DEFINITION"}->[0] = "UNKNOWN"; } if (!defined($ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0])) { $ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0] = "UNKNOWN"; } if (!defined($ReactionData->{"DELTAG"}->[0])) { $ReactionData->{"DELTAG"}->[0] = "UNKNOWN"; } $ReactionDataHash->{$SingletReactions[$i]} = $ReactionData; } print RECONCILIATION $ReactionData->{"REACTIONS"}->[0].";".$ReactionData->{"DEFINITION"}->[0].";".$ReactionData->{"THERMODYNAMIC REVERSIBILITY"}->[0].";".$ReactionData->{"DELTAG"}->[0].";".$ReactionData->{"DIRECTIONS"}->[0].";".$ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{"COUNT"}; for (my $k=0; $k < $SolutionCount; $k++) { print RECONCILIATION ";"; if (defined($ReactionSetHash->{$ReactionSet}->{$AltList->[$j]}->{$k})) { if ($j == 0) { print RECONCILIATION "|".$k."|"; } else { print RECONCILIATION "|".$Index->{$k}."|"; } } } print RECONCILIATION "\n"; } #Adding the current new solutions to the new solutions array if (defined($CurrentNewSolutions) && @{$CurrentNewSolutions} > 0) { push(@{$NewSolutions},@{$CurrentNewSolutions}); } } } #Adding the base reactions to all existing solutions for (my $j=0; $j < @{$Solutions}; $j++) { if (defined($ReactionSetHash->{$ReactionSet}->{$ReactionSet}->{$Solutions->[$j]->{"BASE"}})) { foreach my $SingleReaction (@{$BaseReactions}) { $Solutions->[$j]->{$SingleReaction} = 1; } } } #Adding the new solutions to the set of existing solutions push(@{$Solutions},@{$NewSolutions}); } close(RECONCILIATION); #Now printing a file that defines all of the solutions in a format the testsolutions function understands open (RECONCILIATION, ">$OutputFilenameTwo"); print RECONCILIATION "Experiment;Solution index;Solution cost;Solution reactions\n"; for (my $i=0; $i < @{$Solutions}; $i++) { delete($Solutions->[$i]->{"BASE"}); print RECONCILIATION "SR".$i.";".$i.";10;".join(",",keys(%{$Solutions->[$i]}))."\n"; } close(RECONCILIATION); $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE"); $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1); $Row->{"GF RECON TESTING TIMING"}->[0] = time()."-"; $Row->{"GF RECON SOLUTIONS"}->[0] = @{$Solutions}; $GrowMatchTable->save(); $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE"); #Scheduling the solution testing if ($GapFill == 1) { system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GFSR:BACK:fast:QSUB\""); } else { system($self->figmodel()->config("scheduler executable")->[0]." \"add:testsolutions?".$self->id().$self->selected_version()."?-1?GGSR:BACK:fast:QSUB\""); } } else { #Reading in the solution testing results my $Data; if ($GapFill == 1) { $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GFSREM.txt",""); } else { $Data = $self->figmodel()->database()->load_single_column_file($self->directory().$self->id().$self->selected_version()."-GGSREM.txt",""); } #Reading in the preliminate reconciliation report my $OutputData = $self->figmodel()->database()->load_single_column_file($OutputFilename,""); #Replacing the file tags with actual performance data my $Count = 0; for (my $i=0; $i < @{$Data}; $i++) { if ($Data->[$i] =~ m/^SR(\d+);.+;(\d+\/\d+);/) { my $Index = $1; my $Performance = $Index."/".$2; for (my $j=0; $j < @{$OutputData}; $j++) { $OutputData->[$j] =~ s/\|$Index\|/$Performance/g; } } } $self->figmodel()->database()->print_array_to_file($OutputFilename,$OutputData); my $GrowMatchTable = $self->figmodel()->database()->LockDBTable("GROWMATCH TABLE"); my $Row = $GrowMatchTable->get_row_by_key($self->genome(),"ORGANISM",1); $Row->{"GF RECON TESTING TIMING"}->[0] .= time(); $GrowMatchTable->save(); $self->figmodel()->database()->UnlockDBTable("GROWMATCH TABLE"); } return 1; } =head3 BuildSpecificBiomassReaction Definition: FIGMODELmodel->BuildSpecificBiomassReaction(); Description: =cut sub BuildSpecificBiomassReaction { my ($self) = @_; 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 my $genomestats = $self->figmodel()->get_genome_stats($self->genome()); my $Class = $genomestats->{CLASS}->[0]; my $Name = $genomestats->{NAME}->[0]; #Checking for phoenix variants my $PhoenixVariantTable = $self->figmodel()->database()->GetDBTable("Phoenix variants table"); my $Phoenix = 0; my @Rows = $PhoenixVariantTable->get_rows_by_key($OrganismID,"GENOME"); my $VariantHash; for (my $i=0; $i < @Rows; $i++) { $Phoenix = 1; if (defined($Rows[$i]->{"SUBSYSTEM"}) && defined($Rows[$i]->{"VARIANT"})) { $VariantHash->{$Rows[$i]->{"SUBSYSTEM"}->[0]} = $Rows[$i]->{"VARIANT"}->[0]; } } #Collecting genome data my $RoleHash; my $FeatureTable = $self->figmodel()->GetGenomeFeatureTable($self->genome()); for (my $i=0; $i < $FeatureTable->size(); $i++) { for (my $j=0; $j < @{$FeatureTable->get_row($i)->{"ROLES"}}; $j++) { $RoleHash->{$FeatureTable->get_row($i)->{"ROLES"}->[$j]} = 1; my $Subsystems = $self->figmodel()->subsystems_of_role($FeatureTable->get_row($i)->{"ROLES"}->[$j]); if (defined($Subsystems)) { for (my $k=0; $k < @{$Subsystems}; $k++) { if ($Phoenix == 0) { $VariantHash->{$Subsystems->[$k]} = 1; } } } } } #Scanning through the template item by item and determinine which biomass components should be added my $ComponentTypes; my $EquationData; my $BiomassReactionTemplateTable = $self->figmodel()->database()->GetDBTable("BIOMASS TEMPLATE"); for (my $i=0; $i < $BiomassReactionTemplateTable->size(); $i++) { my $Row = $BiomassReactionTemplateTable->get_row($i); if (defined($Row->{"INCLUSION CRITERIA"}->[0]) && $Row->{"INCLUSION CRITERIA"}->[0] eq "UNIVERSAL") { push(@{$EquationData},$Row); $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1; } else { my $Criteria = $Row->{"INCLUSION CRITERIA"}->[0]; my $End = 0; while ($End == 0) { if ($Criteria =~ m/^(.+)(AND)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(AND)\{([^{^}]+)\}$/ || $Criteria =~ m/^(.+)(OR)\{([^{^}]+)\}(.+)$/ || $Criteria =~ m/^(OR)\{([^{^}]+)\}$/) { print $Criteria."\n"; my $Start = ""; my $End = ""; my $Condition = $1; my $Data = $2; if ($1 ne "AND" && $1 ne "OR") { $Start = $1; $End = $4; $Condition = $2; $Data = $3; } my $Result = "YES"; if ($Condition eq "OR") { $Result = "NO"; } my @Array = split(/\|/,$Data); for (my $j=0; $j < @Array; $j++) { if ($Array[$j] eq "YES" && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Array[$j] eq "NO" && $Condition eq "AND") { $Result = "NO"; last; } elsif ($Array[$j] =~ m/^COMPOUND:(.+)/) { my $Match = 0; for (my $k=0; $k < @{$EquationData}; $k++) { if ($EquationData->[$k]->{"ID"}->[0] eq $1) { $Match = 1; last; } } if ($Match == 1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Match != 1 && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^!COMPOUND:(.+)/) { my $Match = 0; for (my $k=0; $k < @{$EquationData}; $k++) { if ($EquationData->[$k]->{"ID"}->[0] eq $1) { $Match = 1; last; } } if ($Match != 1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Match == 1 && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^NAME:(.+)/) { my $Comparison = $1; if ((!defined($Comparison) || !defined($Name) || $Name =~ m/$Comparison/) && $Condition eq "OR") { $Result = "YES"; last; } elsif (defined($Comparison) && defined($Name) && $Name !~ m/$Comparison/ && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^!NAME:(.+)/) { my $Comparison = $1; if ((!defined($Comparison) || !defined($Name) || $Name !~ m/$Comparison/) && $Condition eq "OR") { $Result = "YES"; last; } elsif (defined($Comparison) && defined($Name) && $Name =~ m/$Comparison/ && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^SUBSYSTEM:(.+)/) { my @SubsystemArray = split(/`/,$1); if (@SubsystemArray == 1) { if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "AND") { $Result = "NO"; last; } } else { my $Match = 0; for (my $k=1; $k < @SubsystemArray; $k++) { if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) { $Match = 1; last; } } if ($Match == 1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Match != 1 && $Condition eq "AND") { $Result = "NO"; last; } } } elsif ($Array[$j] =~ m/^!SUBSYSTEM:(.+)/) { my @SubsystemArray = split(/`/,$1); if (@SubsystemArray == 1) { if ((!defined($VariantHash->{$SubsystemArray[0]}) || $VariantHash->{$SubsystemArray[0]} eq -1) && $Condition eq "OR") { $Result = "YES"; last; } elsif (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} ne -1 && $Condition eq "AND") { $Result = "NO"; last; } } else { my $Match = 0; for (my $k=1; $k < @SubsystemArray; $k++) { if (defined($VariantHash->{$SubsystemArray[0]}) && $VariantHash->{$SubsystemArray[0]} eq $SubsystemArray[$k]) { $Match = 1; last; } } if ($Match != 1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Match == 1 && $Condition eq "AND") { $Result = "NO"; last; } } } elsif ($Array[$j] =~ m/^ROLE:(.+)/) { if (defined($RoleHash->{$1}) && $Condition eq "OR") { $Result = "YES"; last; } elsif (!defined($RoleHash->{$1}) && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^!ROLE:(.+)/) { if (!defined($RoleHash->{$1}) && $Condition eq "OR") { $Result = "YES"; last; } elsif (defined($RoleHash->{$1}) && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^CLASS:(.+)/) { if ($Class eq $1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Class ne $1 && $Condition eq "AND") { $Result = "NO"; last; } } elsif ($Array[$j] =~ m/^!CLASS:(.+)/) { if ($Class ne $1 && $Condition eq "OR") { $Result = "YES"; last; } elsif ($Class eq $1 && $Condition eq "AND") { $Result = "NO"; last; } } } $Criteria = $Start.$Result.$End; print $Criteria."\n"; } else { $End = 1; last; } } if ($Criteria eq "YES") { push(@{$EquationData},$Row); $ComponentTypes->{$Row->{"CLASS"}->[0]}->{$Row->{"ID"}->[0]} = 1; } } } #Building biomass equation my %Reactants; my %Products; foreach my $EquationRow (@{$EquationData}) { #First determine what the coefficient should be my $Coefficient; if (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/^[0-9\.]+$/) { $Coefficient = $EquationRow->{"COEFFICIENT"}->[0]; } elsif (defined($EquationRow->{"COEFFICIENT"}) && $EquationRow->{"COEFFICIENT"}->[0] =~ m/cpd\d\d\d\d\d/) { $Coefficient = 0; 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]}; } else { $Products{$EquationRow->{"ID"}->[0]} = $Coefficient; } } } my $Equation = ""; my @ReactantList = sort(keys(%Reactants)); for (my $i=0; $i < @ReactantList; $i++) { if (length($Equation) > 0) { $Equation .= " + "; } $Equation .= $self->figmodel()->format_coefficient($Reactants{$ReactantList[$i]})." ".$ReactantList[$i]; } $Equation .= " => "; my $First = 1; @ReactantList = sort(keys(%Products)); for (my $i=0; $i < @ReactantList; $i++) { if ($First == 0) { $Equation .= " + "; } $First = 0; $Equation .= $self->figmodel()->format_coefficient($Products{$ReactantList[$i]})." ".$ReactantList[$i]; } #Adding the biomass equation to the biomass table my $NewRow = $self->figmodel()->add_biomass_reaction($Equation,$self->id(),"Template:".$self->genome()); $biomassrxn = $NewRow; } return $biomassrxn; } =head3 PrintSBMLFile Definition: FIGMODELmodel->PrintSBMLFile(); Description: Printing file with model data in SBML format =cut sub PrintSBMLFile { my($self) = @_; #Opening the SBML file for printing my $Filename = $self->config("SBML files")->[0].$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"; } if (!open (SBMLOUTPUT, ">$Filename")) { return; } #Loading and parsing the model data my $ModelTable = $self->reaction_table(); if (!defined($ModelTable) || !defined($ModelTable->{"array"})) { print "Failed to load ".$self->id()."\n"; return; } #Adding intracellular metabolites that also need exchange fluxes to the exchange hash my $ExchangeHash = {"cpd11416" => "c"}; my %CompartmentsPresent; $CompartmentsPresent{"c"} = 1; my %CompoundList; my @ReactionList; my $ReactionTable = $self->figmodel()->database()->GetDBTable("REACTIONS"); for (my $i=0; $i < $ModelTable->size(); $i++) { my $Reaction = $ModelTable->get_row($i)->{"LOAD"}->[0]; if (defined($ReactionTable->get_row_by_key($Reaction,"DATABASE")) && defined($ReactionTable->get_row_by_key($Reaction,"DATABASE")->{"EQUATION"}->[0])) { push(@ReactionList,$Reaction); $_ = $ReactionTable->get_row_by_key($Reaction,"DATABASE")->{"EQUATION"}->[0]; my @MatchArray = /(cpd\d\d\d\d\d)/g; for (my $j=0; $j < @MatchArray; $j++) { $CompoundList{$MatchArray[$j]}->{"c"} = 1; } $_ = $ReactionTable->get_row_by_key($Reaction,"DATABASE")->{"EQUATION"}->[0]; @MatchArray = /(cpd\d\d\d\d\d\[\D\])/g; for (my $j=0; $j < @MatchArray; $j++) { if ($MatchArray[$j] =~ m/(cpd\d\d\d\d\d)\[(\D)\]/) { $CompartmentsPresent{lc($2)} = 1; $CompoundList{$1}->{lc($2)} = 1; } } } } #Printing header to SBML file my $ModelName = $self->id(); $ModelName =~ s/\./_/; print SBMLOUTPUT ''."\n"; print SBMLOUTPUT '' . "\n"; if (defined($self->figmodel()->database()->GetDBTable("MODEL STATS")->{$self->id()}->[0]->{"Organism name"}->[0])) { print SBMLOUTPUT '{$self->id()}->[0]->{"Organism name"}->[0].' SEED model">'."\n"; } else { print SBMLOUTPUT ''."\n"; } #Printing the unit data print SBMLOUTPUT "\n"; print SBMLOUTPUT "\t\n"; print SBMLOUTPUT "\t\t\n"; print SBMLOUTPUT "\t\t\t\n"; print SBMLOUTPUT "\t\t\t\n"; print SBMLOUTPUT "\t\t\t\n"; print SBMLOUTPUT "\t\t\n"; print SBMLOUTPUT "\t\n"; print SBMLOUTPUT "\n"; #Printing compartments for SBML file print SBMLOUTPUT ''."\n"; foreach my $Compartment (keys(%CompartmentsPresent)) { if (defined($self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Name"}->[0])) { my @OutsideList = split(/\//,$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Compartment}->[0]->{"Outside"}->[0]); my $Printed = 0; foreach my $Outside (@OutsideList) { if (defined($CompartmentsPresent{$Outside}) && defined($self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Outside}->[0]->{"Name"}->[0])) { print SBMLOUTPUT '{$Compartment}->[0]->{"Name"}->[0].'" outside="'.$self->figmodel()->database()->GetDBTable("COMPARTMENTS")->{$Outside}->[0]->{"Name"}->[0].'"/>'."\n"; $Printed = 1; last; } } if ($Printed eq 0) { print SBMLOUTPUT '{$Compartment}->[0]->{"Name"}->[0].'"/>'."\n"; } } } print SBMLOUTPUT ''."\n"; #Printing the list of metabolites involved in the model print SBMLOUTPUT ''."\n"; foreach my $Compound (keys(%CompoundList)) { my $Name = $Compound; my $Formula = ""; if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0])) { $Formula = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0]; } if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0])) { $Name = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0]; $Name =~ s/\s/_/; $Name .= "_".$Formula; } $Name =~ s/[<>;&\*]//; my $Charge = 0; if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0])) { $Charge = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0]; } foreach my $Compartment (keys(%{$CompoundList{$Compound}})) { if ($Compartment eq "e") { $ExchangeHash->{$Compound} = "e"; } print SBMLOUTPUT '{$Compartment}->[0]->{"Name"}->[0].'" charge="'.$Charge.'" boundaryCondition="false"/>'."\n"; } } #Printing the boundary species foreach my $Compound (keys(%{$ExchangeHash})) { my $Name = $Compound; my $Formula = ""; if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0])) { $Formula = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"FORMULA"}->[0]; } if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0])) { $Name = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"NAME"}->[0]; $Name =~ s/\s/_/; $Name .= "_".$Formula; } $Name =~ s/[<>;&\*]//; my $Charge = 0; if (defined($self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0])) { $Charge = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->{$Compound}->[0]->{"CHARGE"}->[0]; } print SBMLOUTPUT ''."\n"; } print SBMLOUTPUT ''."\n"; #Printing the list of reactions involved in the model my $ObjectiveCoef; print SBMLOUTPUT ''."\n"; foreach my $Reaction (@ReactionList) { $ObjectiveCoef = "0.0"; if (defined($self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"EQUATION"}->[0])) { if ($Reaction =~ m/^bio/) { $ObjectiveCoef = "1.0"; } my $LowerBound = -10000; my $UpperBound = 10000; my ($Reactants,$Products) = $self->figmodel()->GetReactionSubstrateData($Reaction); my $Name = $Reaction; if (defined($self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"NAME"}->[0])) { $Name = $self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"NAME"}->[0]; $Name =~ s/[<>;&]//g; } my $Reversibility = "true"; if (defined($ModelTable->{$Reaction}->[0]->{"DIRECTIONALITY"}->[0])) { if ($ModelTable->{$Reaction}->[0]->{"DIRECTIONALITY"}->[0] ne "<=>") { $LowerBound = 0; $Reversibility = "false"; } if ($ModelTable->{$Reaction}->[0]->{"DIRECTIONALITY"}->[0] eq "<=") { my $Temp = $Products; $Products = $Reactants; $Reactants = $Temp; } } print SBMLOUTPUT ''."\n"; print SBMLOUTPUT "\n"; my $ECData = ""; if (defined($self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"ENZYME"}->[0])) { $ECData = $self->figmodel()->database()->GetDBTable("REACTIONS")->{$Reaction}->[0]->{"ENZYME"}->[0]; } my $SubsystemData = ""; if (defined($ModelTable->{$Reaction}->[0]->{"SUBSYSTEM"}->[0])) { $SubsystemData = $ModelTable->{$Reaction}->[0]->{"SUBSYSTEM"}->[0]; } my $GeneAssociation = ""; my $ProteinAssociation = ""; if (defined($ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0])) { if (@{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}} == 1 && $ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0] !~ m/\+/) { $GeneAssociation = $ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[0]; } else { if (@{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}} > 1) { $GeneAssociation = "( "; } for (my $i=0; $i < @{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}}; $i++) { if ($i > 0) { $GeneAssociation .= " ) or ( "; } $GeneAssociation .= $ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}->[$i]; } if (@{$ModelTable->{$Reaction}->[0]->{"ASSOCIATED PEG"}} > 1) { $GeneAssociation .= " )"; } } $GeneAssociation =~ s/\+/ and /g; if ($GeneAssociation =~ m/\sor\s/ || $GeneAssociation =~ m/\sand\s/) { $GeneAssociation = "( ".$GeneAssociation." )"; } $ProteinAssociation = $GeneAssociation; if (defined($self->figmodel()->database()->GetDBTable("MODEL STATS")->{$self->id()}->[0]->{"Genome ID"}->[0])) { $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]); } } print SBMLOUTPUT "GENE_ASSOCIATION:".$GeneAssociation."\n"; print SBMLOUTPUT "PROTEIN_ASSOCIATION:".$ProteinAssociation."\n"; print SBMLOUTPUT "SUBSYSTEM: ".$SubsystemData."\n"; print SBMLOUTPUT "PROTEIN_CLASS: ".$ECData."\n"; print SBMLOUTPUT "\n"; print SBMLOUTPUT "\n"; foreach my $Reactant (@{$Reactants}) { print SBMLOUTPUT '[0]."_".$Reactant->{"COMPARTMENT"}->[0].'" stoichiometry="'.$Reactant->{"COEFFICIENT"}->[0].'"/>'."\n"; } print SBMLOUTPUT "\n"; print SBMLOUTPUT "\n"; foreach my $Product (@{$Products}) { print SBMLOUTPUT '[0]."_".$Product->{"COMPARTMENT"}->[0].'" stoichiometry="'.$Product->{"COEFFICIENT"}->[0].'"/>'."\n"; } print SBMLOUTPUT "\n"; print SBMLOUTPUT "\n"; print SBMLOUTPUT "\t\n"; print SBMLOUTPUT "\t\t\t FLUX_VALUE \n"; print SBMLOUTPUT "\t\n"; print SBMLOUTPUT "\t\n"; print SBMLOUTPUT "\t\t\n"; print SBMLOUTPUT "\t\t\n"; print SBMLOUTPUT "\t\t\n"; print SBMLOUTPUT "\t\t\n"; print SBMLOUTPUT "\t\n"; print SBMLOUTPUT "\n"; print SBMLOUTPUT ''."\n"; } } my @ExchangeList = keys(%{$ExchangeHash}); foreach my $ExCompound (@ExchangeList) { my $ExCompoundName = $ExCompound; my $Row = $self->figmodel()->database()->GetDBTable("COMPOUNDS")->get_row_by_key($ExCompound,"DATABASE"); if (defined($Row) && defined($Row->{"NAME"}->[0])) { $ExCompoundName = $Row->{"NAME"}->[0]; $ExCompoundName =~ s/[<>;&]//g; } $ObjectiveCoef = "0.0"; print SBMLOUTPUT ''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t\t".'GENE_ASSOCIATION: '."\n"; print SBMLOUTPUT "\t\t".'PROTEIN_ASSOCIATION: '."\n"; print SBMLOUTPUT "\t\t".'SUBSYSTEM: S_'."\n"; print SBMLOUTPUT "\t\t".'PROTEIN_CLASS: '."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t\t".''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t\t".''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT "\t\t".''."\n"; print SBMLOUTPUT "\t\t\t\t".' FLUX_VALUE '."\n"; print SBMLOUTPUT "\t\t".''."\n"; print SBMLOUTPUT "\t\t".''."\n"; print SBMLOUTPUT "\t\t\t".''."\n"; print SBMLOUTPUT "\t\t\t".''."\n"; print SBMLOUTPUT "\t\t\t".''."\n"; print SBMLOUTPUT "\t\t\t".''."\n"; print SBMLOUTPUT "\t\t".''."\n"; print SBMLOUTPUT "\t".''."\n"; print SBMLOUTPUT ''."\n"; } #Closing out the file print SBMLOUTPUT ''."\n"; print SBMLOUTPUT ''."\n"; print SBMLOUTPUT "\n"; close(SBMLOUTPUT); } =head3 PrintModelSimpleReactionTable Definition: success()/fail() FIGMODELmodel->PrintModelSimpleReactionTable(); Description: Prints the table of model data =cut sub PrintModelSimpleReactionTable { my ($self) = @_; my $rxntbl = $self->reaction_table(); my $tbl = $self->create_table_prototype("ModelSimpleReactionTable"); for (my $i=0; $i < $rxntbl->size(); $i++) { my $row = $rxntbl->get_row($i); $row->{DATABASE} = $row->{LOAD}; $tbl->add_row($row); } $tbl->save(); system("cp ".$tbl->filename()." ".$self->figmodel()->config("Model table download directory")->[0].$self->figmodel()->config("ModelSimpleReactionTable")->{filename_prefix}->[0]."-".$self->id().".txt"); } =head3 PrintModelLPFile Definition: success()/fail() FIGMODELmodel->PrintModelLPFile(); Description: Prints the lp file needed to run the model using the mpifba program =cut sub PrintModelLPFile { my ($self) = @_; #Printing lp and key file for model my $UniqueFilename = $self->figmodel()->filename(); #Printing the standard FBA file system($self->figmodel()->GenerateMFAToolkitCommandLineCall($UniqueFilename,$self->id(),"NoBounds",["ProdFullFBALP"],undef,$self->id().$self->selected_version()."-LPPrint.log",undef,$self->selected_version())); system("cp ".$self->config("MFAToolkit output directory")->[0].$UniqueFilename."/CurrentProblem.lp ".$self->directory()."FBA-".$self->id().$self->selected_version().".lp"); my $KeyTable = FIGMODELTable::load_table($self->config("MFAToolkit output directory")->[0].$UniqueFilename."/VariableKey.txt",";","|",0,undef); if (!defined($KeyTable)) { print STDERR "FIGMODEL:RunAllStudiesWithDataFast: ".$self->id()." LP file could not be printed.\n"; return 0; } $KeyTable->headings(["Variable type","Variable ID"]); $KeyTable->save($self->directory()."FBA-".$self->id().$self->selected_version().".key"); unlink($self->config("database message file directory")->[0].$self->id().$self->selected_version()."-LPPrint.log"); $self->figmodel()->clearing_output($UniqueFilename,"FBA-".$self->id().$self->selected_version().".lp"); } =head3 patch_model Definition: FIGMODELTable:patch results = FIGMODELmodel->patch_model(FIGMODELTable:patch table); Description: =cut sub patch_model { my ($self,$tbl) = @_; #Instantiating table my $results = FIGMODELTable->new(["Reactions","New genes","Old genes","Genes","Roles","Status"],$self->directory()."PatchResults-".$self->id().$self->selected_version().".tbl",["Reaction"],"\t",";",undef); #Getting genome annotations my $features = $self->figmodel()->database()->get_genome_feature_table($self->genome()); #Gettubg reaction table my $reactions = $self->reaction_table(); #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"}}); } } } #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"]}); } } $reactions->save(); $results->save(); $self->update_model_stats(); $self->PrintModelLPFile(); $self->run_default_model_predictions(); #Returning results return $results; } =head3 translate_genes Definition: FIGMODELmodel->translate_genes(); Description: =cut sub translate_genes { my ($self) = @_; #Loading gene translations if (!defined($self->{_gene_aliases})) { #Loading gene aliases from feature table my $tbl = $self->figmodel()->GetGenomeFeatureTable($self->genome()); if (defined($tbl)) { for (my $i=0; $i < $tbl->size(); $i++) { my $row = $tbl->get_row($i); if ($row->{ID}->[0] =~ m/(peg\.\d+)/) { my $geneID = $1; for (my $j=0; $j < @{$row->{ALIASES}}; $j++) { $self->{_gene_aliases}->{$row->{ALIASES}->[$j]} = $geneID; } } } } #Loading additional gene aliases from the database if (-e $self->figmodel()->config("Translation directory")->[0]."AdditionalAliases/".$self->genome().".txt") { my $AdditionalAliases = $self->figmodel()->database()->load_multiple_column_file($self->figmodel()->config("Translation directory")->[0]."AdditionalAliases/".$self->genome().".txt","\t"); for (my $i=0; $i < @{$AdditionalAliases}; $i++) { $self->{_gene_aliases}->{$AdditionalAliases->[$i]->[1]} = $AdditionalAliases->[$i]->[0]; } } } #Cycling through reactions and translating genes for (my $i=0; $i < $self->reaction_table()->size(); $i++) { my $row = $self->reaction_table()->get_row($i); if (defined($row->{"ASSOCIATED PEG"})) { for (my $j=0; $j < @{$row->{"ASSOCIATED PEG"}}; $j++) { my $Original = $row->{"ASSOCIATED PEG"}->[$j]; $Original =~ s/\sand\s/:/g; $Original =~ s/\sor\s/;/g; my @GeneNames = split(/[,\+\s\(\):;]/,$Original); foreach my $Gene (@GeneNames) { if (length($Gene) > 0 && defined($self->{_gene_aliases}->{$Gene})) { my $Replace = $self->{_gene_aliases}->{$Gene}; $Original =~ s/([^\w])$Gene([^\w])/$1$Replace$2/g; $Original =~ s/^$Gene([^\w])/$Replace$1/g; $Original =~ s/([^\w])$Gene$/$1$Replace/g; $Original =~ s/^$Gene$/$Replace/g; } } $Original =~ s/:/ and /g; $Original =~ s/;/ or /g; $row->{"ASSOCIATED PEG"}->[$j] = $Original; } } } #Archiving model and saving reaction table $self->ArchiveModel(); $self->reaction_table()->save(); } =head3 feature_web_data Definition: string:web output for feature/model connection = FIGMODELmodel->feature_web_data(FIGMODELfeature:feature); Description: =cut sub feature_web_data { my ($self,$feature) = @_; #First checking if the feature is in the model my $featureGenome = $feature->{GENOME}->[0]; if ($feature->{GENOME}->[0] =~ m/\>(\d+\.\d+)\genome()) { return "Not in model"; } my $data = $self->get_feature_data($feature->{ID}->[0]); if (!defined($data)) { return "Not in model"; } my $output; #Printing predictions if (defined($data->{$self->id()."PREDICTIONS"})) { #Parsing essentiality data my $essentialityData; if (defined($feature->{ESSENTIALITY})) { for (my $i=0; $i < @{$feature->{ESSENTIALITY}};$i++) { my @temp = split(/:/,$feature->{ESSENTIALITY}->[$i]); $essentialityData->{$temp[0]} = $temp[1]; } } my $predictionHash; for (my $i=0; $i < @{$data->{$self->id()."PREDICTIONS"}};$i++) { my @temp = split(/:/,$data->{$self->id()."PREDICTIONS"}->[$i]); if (defined($essentialityData->{$temp[0]})) { if ($temp[1] eq "essential" && $essentialityData->{$temp[0]} eq "essential") { push(@{$predictionHash->{"Correct negative"}},$temp[0]); } elsif ($temp[1] eq "nonessential" && $essentialityData->{$temp[0]} eq "essential") { push(@{$predictionHash->{"False positive"}},$temp[0]); } elsif ($temp[1] eq "essential" && $essentialityData->{$temp[0]} eq "nonessential") { push(@{$predictionHash->{"False negative"}},$temp[0]); } elsif ($temp[1] eq "nonessential" && $essentialityData->{$temp[0]} eq "nonessential") { push(@{$predictionHash->{"Correct positive"}},$temp[0]); } } else { push(@{$predictionHash->{$temp[1]}},$temp[0]); } } foreach my $key (keys(%{$predictionHash})) { my $string = $key; $string = ucfirst($string); push(@{$output},'{$key}}).'">'.$string.''); } } #Printing reactions if (defined($data->{$self->id()."REACTIONS"})) { for (my $i=0; $i < @{$data->{$self->id()."REACTIONS"}};$i++) { my $rxnData = $self->get_reaction_data($data->{$self->id()."REACTIONS"}->[$i]); my $reactionString = $self->figmodel()->web()->create_reaction_link($data->{$self->id()."REACTIONS"}->[$i],join(" or ",@{$rxnData->{"ASSOCIATED PEG"}}),$self->id()); if (defined($rxnData->{PREDICTIONS})) { my $predictionHash; for (my $i=0; $i < @{$rxnData->{PREDICTIONS}};$i++) { my @temp = split(/:/,$rxnData->{PREDICTIONS}->[$i]); push(@{$predictionHash->{$temp[1]}},$temp[0]); } $reactionString .= "("; foreach my $key (keys(%{$predictionHash})) { if ($key eq "Essential =>") { $reactionString .= '{$key}}).'">E=>,'; } elsif ($key eq "Essential <=") { $reactionString .= '{$key}}).'">E<=,'; } elsif ($key eq "Active =>") { $reactionString .= '{$key}}).'">A=>,'; } elsif ($key eq "Active <=") { $reactionString .= '{$key}}).'">A<=,'; } elsif ($key eq "Active <=>") { $reactionString .= '{$key}}).'">A,'; } elsif ($key eq "Inactive") { $reactionString .= '{$key}}).'">I,'; } elsif ($key eq "Dead") { $reactionString .= 'D,'; } } $reactionString =~ s/,$/)/; } push(@{$output},$reactionString); } } #Returning output return join("
",@{$output}); } =head3 remove_obsolete_reactions Definition: void FIGMODELmodel->remove_obsolete_reactions(); Description: =cut sub remove_obsolete_reactions { my ($self) = @_; (my $dummy,my $translation) = $self->figmodel()->put_two_column_array_in_hash($self->figmodel()->database()->load_multiple_column_file($self->figmodel()->config("Translation directory")->[0]."ObsoleteRxnIDs.txt","\t")); my $rxnTbl = $self->reaction_table(); if (defined($rxnTbl)) { for (my $i=0; $i < $rxnTbl->size(); $i++) { my $row = $rxnTbl->get_row($i); if (defined($translation->{$row->{LOAD}->[0]}) || defined($translation->{$row->{LOAD}->[0]."r"})) { my $direction = $row->{DIRECTION}->[0]; my $newRxn; if (defined($translation->{$row->{LOAD}->[0]."r"})) { $newRxn = $translation->{$row->{LOAD}->[0]."r"}; if ($direction eq "<=") { $direction = "=>"; } elsif ($direction eq "=>") { $direction = "<="; } } else { $newRxn = $translation->{$row->{LOAD}->[0]}; } #Checking if the new reaction is already in the model my $newRow = $rxnTbl->get_row_by_key($newRxn,"LOAD"); if (defined($newRow)) { #Handling direction if ($newRow->{DIRECTION}->[0] ne $direction) { $newRow->{DIRECTION}->[0] = "<=>"; } push(@{$row->{"ASSOCIATED PEG"}},@{$rxnTbl->get_row($i)->{"ASSOCIATED PEG"}}); } else { $rxnTbl->get_row($i)->{LOAD}->[0] = $newRxn; $rxnTbl->get_row($i)->{DIRECTION}->[0] = $direction; } } } $rxnTbl->save(); } } 1;