# _*_ Perl _*_ # # # use strict; use FIGV; use Scenario; my ($genome_id,$orgdir) = @ARGV; my $debug = 0; unless (scalar (@ARGV) >= 1) { print STDERR "usage: make_SBML_model \n"; exit(-1); } chomp $genome_id; my $fig; if($orgdir) { $fig = new FIGV($orgdir); } else { $fig = new FIGV; } Scenario::set_fig($fig,$genome_id); print STDERR "Creating SBML model for genome $genome_id\n" if($debug); print STDERR "\tStage 1 - Loading Scenarios, model inputs and biomass\n" if($debug); my @scenarios = @{Scenario->get_genome_scenarios($genome_id,1)}; my $filebase = $fig->model_directory($genome_id)."/Analysis/"; my %input_cpds; open(INPUTS,$filebase."inputs.txt") or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs"); while() { chomp; my @line = split "\t", $_; chomp @line; $input_cpds{$line[0]} = 1; } close(INPUTS); #we are going to load up the cofactor compounds into inputs as well open(COFACT,$filebase."cofactors.txt") or die("Failed to open ".$filebase."cofactors.txt"); while() { chomp; my @line = split "\t",$_; chomp @line; $input_cpds{$line[0]} = 1; } close(COFACT); my %compound_id_to_name; open(IN,"$FIG_Config::global/Models/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt"); while($_ = ){ chomp($_); my($id,$name) = split("\t",$_); $compound_id_to_name{$id} = $name; } close(IN); my %biomass_cpds; open(BIOMASS,$filebase."biomass_valid.txt") or die ("Failed to open ".$filebase."biomass_valid.txt\n"); while() { chomp; my @temp = split "\t", $_; chomp @temp; $biomass_cpds{$temp[0]} = 1; } close(BIOMASS); open(IN,$filebase."scenarios_valid.txt"); my @lines = ; chomp @lines; close(IN); my %valid_scenarios; foreach (@lines) { $valid_scenarios{$_} = 1; } my %output_cpds; open(OUTS,$filebase."outputs.txt"); while() { chomp; my @temp = split "\t", $_; chomp @temp; $output_cpds{$temp[0]} = 1; } close(OUTS); #load core_metabolites, scenario reaction information my %core_metabolites; my %reaction_to_substrates; my %reaction_to_products; my %reaction_reversible; my %reaction_notes; foreach my $scenario (@scenarios) { if(!defined $valid_scenarios{$scenario->get_scenario_name}) { next; } $scenario->load_reaction_pegs($fig); my $subsystem = $scenario->get_subsystem_name(); foreach my $reaction (@{$scenario->get_reaction_ids()}) { $reaction_to_substrates{$reaction} = {}; $reaction_to_products{$reaction} = {}; my %substrates = %{$scenario->get_reaction_substrates($reaction)}; my %products = %{$scenario->get_reaction_products($reaction)}; $reaction_reversible{$reaction} = $scenario->get_reaction_reversibility($reaction); map {$core_metabolites{$_} = 1; $reaction_to_substrates{$reaction}->{$_} = $substrates{$_}; } keys %substrates; map {$core_metabolites{$_} = 1; $reaction_to_products{$reaction}->{$_} = $products{$_}; } keys %products; if(!defined $reaction_notes{$reaction}) { $reaction_notes{$reaction} = {}; } foreach ( @{$scenario->get_reaction_pegs($reaction) }) { $reaction_notes{$reaction}->{"SEED_PEG:$_"} = 1; my $roles = $fig->protein_subsystem_to_roles($_, $subsystem); foreach my $role (@$roles) { $reaction_notes{$reaction}->{"Functional Role:$role"} = 1; } } $reaction_notes{$reaction}->{"Subsystem/Scenario:".$scenario->get_scenario_name()} = 1; } } print STDERR "\tStage 1 - Complete\n" if($debug); print STDERR "\tStage 2 - Write SBML file\n" if($debug); open(OUT,">".$filebase.$genome_id."_FBA_Model.xml"); write_header(); write_species_list(); write_reactions_list(); print OUT "\n\n"; close(OUT); print STDERR "\tStage 2 - Complete\n" if($debug); sub write_header { my $print_id = $genome_id; $print_id =~ s/\.//g; print OUT " \n"; } sub write_species_list { print OUT "\n"; #input compounds foreach my $id (sort keys %input_cpds) { my $name = $compound_id_to_name{$id}; $name =~ s/\\"//g; make_species($id,$name,1,1); } #biomass compounds foreach my $id (sort keys(%biomass_cpds)) { my $name = $compound_id_to_name{$id}; $name =~ s/\\"//g; if(!$input_cpds{$id}) { make_species($id,$name,1,1); } } #output compounds foreach my $id (sort keys(%output_cpds)) { my $name = $compound_id_to_name{$id}; $name =~ s/\\"//g; if(!$input_cpds{$id} && !$biomass_cpds{$id}) { make_species($id,$name,1,1); } } #core_metabolites are all the compounds used in any scenario. foreach my $id (sort keys(%core_metabolites)) { my $name = $compound_id_to_name{$id}; $name =~ s/\\"//g; if(!$input_cpds{$id} && !$biomass_cpds{$id} && !$output_cpds{$id}) { make_species($id,$name,1,1); } } print OUT "\t\n"; print OUT "\n"; } sub write_reactions_list { print OUT "\n"; #make uptake/transport reactions for the input compounds (provide cpds for the system) foreach my $input (keys %input_cpds) { make_reaction($input."_transport_in",0,{$input."_e" => "1.000000"},{$input => "1.000000"},"0.000000",()); } #Make sink/transport reactions for the output compounds foreach my $output (keys %output_cpds) { make_reaction($output."_transport_out",0,{$output => "1.000000"},{$output."_e" => "1.000000"},"0.000000"); } #make sinks/transport reactions for all cpds not already outputed or biomassed foreach my $id (keys %core_metabolites) { if(!$output_cpds{$id} && !$biomass_cpds{$id}) { make_reaction($id."_sink",0,{$id => "1.000000"},{$id."_e" => "1.000000"},"0.000000"); } } #Now we write out the Scenario reactions foreach my $rid (keys %reaction_to_substrates) { make_reaction($rid,$reaction_reversible{$rid},$reaction_to_substrates{$rid},$reaction_to_products{$rid}, "0.000000",$reaction_notes{$rid}); } #now its time for the biomass reaction make_reaction("BIOMASS",0,\%biomass_cpds,{"biomass_out" => "1.000000"},"-1.000000"); print OUT "\n"; } sub make_species { my ($id,$name,$internal,$external) = @_; if($internal) { print OUT "\t\n"; } if($external) { print OUT "\t\n"; } } sub make_reaction { my($reaction_id,$reversible,$substrates,$products,$objective,$notes) = @_; if($reversible) { print OUT "\n"; } else { print OUT "\n"; } if(defined $notes) { print OUT "\t\n"; foreach my $note (keys %$notes) { print OUT "\t\t$note\n"; } print OUT "\t\n"; } print OUT "\t\n" if(keys %$substrates); foreach (keys %$substrates) { print OUT "\t\t{$_}\"/>\n"; } print OUT "\t\n" if(keys %$substrates); print OUT "\t\n" if(keys %$products); foreach (keys %$products) { print OUT "\t\t{$_}\"/>\n"; } print OUT "\t\n" if(keys %$products); print OUT "\t\n"; print OUT "\t\t\n". "\t\t\t\n". "\t\t\t\t LOWER_BOUND \n". "\t\t\t\t UPPER_BOUND \n". "\t\t\t\t OBJECTIVE_COEFFICIENT \n". "\t\t\t\t FLUX_VALUE \n". "\t\t\t\t REDUCED_COST \n". "\t\t\t\n\t\t\n"; print OUT "\t\t\n"; my $low_bound = "0.000000"; if($reversible) { $low_bound = "-999999.000000"; } print OUT "\t\t\t\n" . "\t\t\t\n". "\t\t\t\n". "\t\t\t\n". "\t\t\t\n"; print OUT "\t\t\n". "\t". "\n\n"; }