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

View of /FigKernelScripts/make_MPS_model.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (download) (as text) (annotate)
Wed Nov 12 18:29:31 2008 UTC (11 years ago) by dejongh
Branch: MAIN
CVS Tags: mgrast_dev_08112011, rast_rel_2009_05_18, mgrast_dev_08022011, rast_rel_2014_0912, myrast_rel40, mgrast_dev_05262011, rast_rel_2008_12_18, mgrast_dev_04082011, rast_rel_2010_0928, mgrast_version_3_2, mgrast_dev_12152011, mgrast_dev_06072011, rast_rel_2009_0925, rast_rel_2010_0526, rast_rel_2014_0729, mgrast_dev_02212011, rast_rel_2010_1206, mgrast_release_3_0, mgrast_dev_03252011, rast_rel_2010_0118, rast_rel_2009_02_05, rast_rel_2011_0119, mgrast_release_3_0_4, mgrast_release_3_0_2, mgrast_release_3_0_3, mgrast_release_3_0_1, mgrast_dev_03312011, mgrast_release_3_1_2, mgrast_release_3_1_1, mgrast_release_3_1_0, mgrast_dev_04132011, mgrast_dev_04012011, rast_rel_2009_07_09, rast_rel_2010_0827, myrast_33, rast_rel_2011_0928, mgrast_dev_04052011, mgrast_dev_02222011, rast_rel_2009_03_26, mgrast_dev_10262011, rast_rel_2008_11_24, HEAD
Changes since 1.6: +1 -1 lines
switch to scenario_directory

# _*_ Perl _*_
#
#
#
use strict;
use FIGV;
use Scenario;

my ($genome_id,$orgdir) = @ARGV;

my $debug = 0;

unless (scalar (@ARGV) >= 1)
{
    print STDERR "usage: make_MPS_model <genome_id>\n";
    exit(-1);
}
my $fig;
if($orgdir)
{
    $fig = new FIGV($orgdir);
}
else
{
    $fig = new FIGV;
}
chomp $genome_id;
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->scenario_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."secondary_compounds.txt")
    or die("Failed to open ".$filebase."secondary_compounds.txt");
while(<COFACT>)
{
    chomp;
    my @line = split "\t",$_;
    chomp @line;
    $input_cpds{$line[0]} = 1;
}
close(COFACT);

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);

#$biomass_cpds{"C00002"} = 1;

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;
foreach my $scenario (@scenarios)
{
    if(!defined $valid_scenarios{$scenario->get_scenario_name})
    {
	next;
    }
    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;
    }
}
#the ATP Synthase hack reaction
$core_metabolites{"C00009"} = 1;
$core_metabolites{"C00001"} = 1;
$core_metabolites{"C00008"} = 1;
$core_metabolites{"C00002"} = 1;
$reaction_reversible{"ATPSYN"} = 0;
$reaction_to_substrates{"ATPSYN"} = {};
$reaction_to_products{"ATPSYN"} = {};
$reaction_to_substrates{"ATPSYN"}->{"C00008"} = 1;
$reaction_to_substrates{"ATPSYN"}->{"C00009"} = 1;
$reaction_to_products{"ATPSYN"}->{"C00001"} = 1;
$reaction_to_products{"ATPSYN"}->{"C00002"} = 1;
#end hack  

open(OUT,">".$filebase.$genome_id.".mps");

print OUT "NAME          test\n";
print OUT "ROWS\n N  OBJ\n G  BIOMAS\n";

my @all_cpds;

#write the compound list.
foreach my $id (sort keys %input_cpds)
{
    print OUT " G  $id\n"; #regular cellular
    push @all_cpds, $id;
    $id = substr($id,1);
    print OUT " G  E$id\n"; #extra cellular
    push @all_cpds, "E".$id;
}
foreach my $id (sort keys %output_cpds)
{
    if(!$input_cpds{$id})
    {
	print OUT " G  $id\n"; #regular cellular
	push @all_cpds, $id;
	$id = substr($id,1);
	print OUT " G  E$id\n"; #extra cellular
	push @all_cpds, "E".$id;
    }
}
foreach my $id (sort keys %biomass_cpds)
{
    if(!$input_cpds{$id} && !$output_cpds{$id})
    {
	print OUT " G  $id\n"; #regular cellular
	push @all_cpds, $id;
	$id = substr($id,1);
	print OUT " G  E$id\n"; #extra cellular
	push @all_cpds, "E".$id;
    }
}

foreach my $id (sort keys %core_metabolites)
{
    if(!$input_cpds{$id} && !$output_cpds{$id} && !$biomass_cpds{$id})
    {
	print OUT " G  $id\n"; #regular cellular
	push @all_cpds, $id;
	$id = substr($id,1);
	print OUT " G  E$id\n"; #extra cellular
	push @all_cpds, "E".$id;
    }
}

#special compounds
print OUT " G  bioout\n";
push @all_cpds, "bioout";

print OUT "COLUMNS\n";
#for each reaction print ' 4 spaces   $rid 4 spaces   $cid  14 spaces        $stoich' 

my @all_rids;

#scenario reactions
foreach my $rid (sort keys %reaction_to_substrates)
{
    
    foreach (sort keys %{$reaction_to_substrates{$rid}})
    {
	if(!defined $reaction_to_products{$rid}->{$_})
	{
	    print OUT "    $rid    $_             -$reaction_to_substrates{$rid}->{$_}\n";
	}
    }
    foreach (sort keys %{$reaction_to_products{$rid}})
    {
	if(!defined $reaction_to_substrates{$rid}->{$_})
	{
	    print OUT "    $rid    $_              $reaction_to_products{$rid}->{$_}\n";
	}
    }
    push @all_rids, $rid;
    if($reaction_reversible{$rid})
    {
	foreach (sort keys %{$reaction_to_substrates{$rid}})
	{
	    print OUT "    $rid\_r  $_              $reaction_to_substrates{$rid}->{$_}\n";
	}
	foreach (sort keys %{$reaction_to_products{$rid}})
	{
	    print OUT "    $rid\_r  $_             -$reaction_to_products{$rid}->{$_}\n";
	}
	push @all_rids, $rid."_r";
    }

}

#transport reactions
foreach my $cid (keys %input_cpds)
{
    my $trans_id = "T".substr($cid,1);
    my $extern_id = "E".substr($cid,1);
    print OUT "    $trans_id    $cid              1\n";
    print OUT "    $trans_id    $extern_id             -1\n";
    #lets enable reverse transports to fix some sink rxn issues
    print OUT "    $trans_id\_r  $cid             -1\n";
    print OUT "    $trans_id\_r  $extern_id              1\n";
}
foreach my $cid (keys %output_cpds)
{
    if(!$input_cpds{$cid})
    {
	my $trans_id = "T".substr($cid,1);
	my $extern_id = "E".substr($cid,1);
	print OUT "    $trans_id    $cid             -1\n";
	print OUT "    $trans_id    $extern_id              1\n";
    }
}

#biomass reaction
foreach my $cid (keys %biomass_cpds)
{
    print OUT "    BIOMAS    $cid             -1\n";
}
print OUT "    BIOMAS    bioout              1\n";
#uptake reactions
foreach my $cid (keys %input_cpds)
{
    my $temp = substr($cid,1);
    print OUT "    U$temp    E$temp              1\n";
}
#sink reactions
foreach my $cid (keys %output_cpds)
{
    my $temp = substr($cid,1);
    print OUT "    S$temp    E$temp             -1\n";
}

#sink reactions for all?
foreach my $cid (@all_cpds)
{
    if($cid =~/E\d{5}/ || $cid eq "bioout")
    {
	next;
    } 
    my $temp  = substr($cid,1);
    if(!$input_cpds{$cid} && !$output_cpds{$cid})
    {
	print OUT "    T$temp    E$temp              1\n";
	print OUT "    T$temp    C$temp             -1\n";
    }
    if(!$output_cpds{$cid})
    {
	print OUT "    S$temp    E$temp             -1\n";
    }
}



#OBJ function and biomass
print OUT "    FOBJ      OBJ                 1\n";
print OUT "    FOBJ      bioout             -1\n";
print OUT "RHS\nRANGES\n";
foreach my $cid (@all_cpds)
{
    print OUT "    CBOUND    $cid              0\n";
}
print OUT "BOUNDS\n";
print OUT " UP RBOUND    BIOMAS            1000\n";
foreach my $rid (@all_rids)
{
    if(!($rid =~/_r$/))
    {
	print OUT " UP RBOUND    $rid            1000\n";
    }
    else
    {
	print OUT " UP RBOUND    $rid          1000\n";
    }
}
foreach my $cid (keys %input_cpds)
{
    $cid = substr($cid,1);
    print OUT " UP RBOUND    T$cid            1000\n";
    print OUT " UP RBOUND    U$cid            1000\n";
}
foreach my $cid (keys %output_cpds)
{
    $cid = substr($cid,1);
    print OUT " UP RBOUND    T$cid            1000\n";
    print OUT " UP RBOUND    S$cid            1000\n";
}


print OUT "ENDATA\n";	
close(OUT);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3