[Bio] / FigKernelScripts / make_SBML_model.pl Repository:
ViewVC logotype

View of /FigKernelScripts/make_SBML_model.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (download) (as text) (annotate)
Mon Jul 23 20:33:49 2007 UTC (12 years, 4 months ago) by formsma
Branch: MAIN
Changes since 1.2: +1 -1 lines
Fixed errors on filehandleing with RASt - Kevin

# _*_ 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 <genome_id>\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(<INPUTS>)
{
    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(<COFACT>)
{
    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($_ = <IN>){
    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(<BIOMASS>)
{
    chomp;
    my @temp = split "\t", $_;
    chomp @temp;
    $biomass_cpds{$temp[0]} = 1;
}
close(BIOMASS);

open(IN,$filebase."scenarios_valid.txt");
my @lines = <IN>;
chomp @lines;
close(IN);
my %valid_scenarios;
foreach (@lines)
{
    $valid_scenarios{$_} = 1;
}

my %output_cpds;
open(OUTS,$filebase."outputs.txt");
while(<OUTS>)
{
    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 "</model>\n</sbml>\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 "<?xml version=\"1.0\" encoding=\"UTF-8\"?>
       <sbml xmlns=\"http://www.sbml.org/sbml/level2\" level=\"2\" version=\"1\" xmlns:html=\"http://www.w3.org/1999/xhtml\">
       <model id=\"Genome_$print_id\" name=\"Genome_$print_id\"> 
       <listOfUnitDefinitions>
	    <unitDefinition id=\"mmol_per_gDW_per_hr\">
		<listOfUnits>
			<unit kind=\"mole\" scale=\"-3\"/>
			<unit kind=\"gram\" exponent=\"-1\"/>
			<unit kind=\"second\" multiplier=\".00027777\" exponent=\"-1\"/>
		</listOfUnits>
	    </unitDefinition>
       </listOfUnitDefinitions>
       <listOfCompartments>
	     <compartment id=\"Extra_organism\"/>
   	     <compartment id=\"Cytosol\" outside=\"Extra_organism\"/>
       </listOfCompartments>\n";
}

sub write_species_list
{
    print OUT "<listOfSpecies>\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<species id=\"biomass_out\" name=\"BIOMASS\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n";
    
    print OUT "</listOfSpecies>\n";
}

sub write_reactions_list
{
    print OUT "<listOfReactions>\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 "</listOfReactions>\n";
}

sub make_species
{
    my ($id,$name,$internal,$external) = @_;
    if($internal)
    {
	print OUT "\t<species id=\"$id\" name=\"$name\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n";
    }
    if($external)
    {
	print OUT "\t<species id=\"$id\_e\" name=\"$name\" compartment=\"Extra_organism\" charge=\"0\" boundaryCondition=\"false\"/>\n";
    }
}

sub make_reaction
{
    my($reaction_id,$reversible,$substrates,$products,$objective,$notes) = @_;
    

    if($reversible)
    {
	print OUT "<reaction id=\"$reaction_id\" reversible=\"true\">\n";
    }
    else
    {
	print OUT "<reaction id=\"$reaction_id\" reversible=\"false\">\n";
    }
    if(defined $notes)
    {
	print OUT "\t<notes>\n";
	foreach my $note (keys %$notes)
	{
	    print OUT "\t\t<html:p>$note</html:p>\n";
	}
	print OUT "\t</notes>\n";
    }
    print OUT "\t<listOfReactants>\n" if(keys %$substrates);
    foreach (keys %$substrates)
    {
    	print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$substrates->{$_}\"/>\n";
    }
    print OUT "\t</listOfReactants>\n" if(keys %$substrates);
    
    print OUT "\t<listOfProducts>\n" if(keys %$products);
    foreach (keys %$products)
    {
	print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$products->{$_}\"/>\n";
    }
    print OUT "\t</listOfProducts>\n" if(keys %$products);
    
    print OUT "\t<kineticLaw>\n";
    
    print OUT "\t\t<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n".
	"\t\t\t<apply>\n".
	"\t\t\t\t<ci> LOWER_BOUND </ci>\n".
	"\t\t\t\t<ci> UPPER_BOUND </ci>\n".
	"\t\t\t\t<ci> OBJECTIVE_COEFFICIENT </ci>\n".
	"\t\t\t\t<ci> FLUX_VALUE </ci>\n".
	"\t\t\t\t<ci> REDUCED_COST </ci>\n".
	"\t\t\t</apply>\n\t\t</math>\n";
    
    
    print OUT "\t\t<listOfParameters>\n";
    my $low_bound = "0.000000";
    if($reversible)
    {
	$low_bound = "-999999.000000";
    }
    print OUT "\t\t\t<parameter id=\"LOWER_BOUND\" value=\"$low_bound\" units=\"mmol_per_gDW_per_hr\"/>\n" .
	"\t\t\t<parameter id=\"UPPER_BOUND\" value=\"999999.000000\" units=\"mmol_per_gDW_per_hr\"/>\n".
	"\t\t\t<parameter id=\"OBJECTIVE_COEFFICIENT\" value=\"$objective\"/>\n".
	"\t\t\t<parameter id=\"FLUX_VALUE\" value=\"0.000000\" units=\"mmol_per_gDW_per_hr\"/>\n".
	"\t\t\t<parameter id=\"REDUCED_COST\" value=\"0.000000\"/>\n";
    print OUT "\t\t</listOfParameters>\n".
	"\t</kineticLaw>".
	"\n</reaction>\n";
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3