[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.1 - (download) (as text) (annotate)
Tue Jul 17 16:17:00 2007 UTC (12 years, 10 months ago) by formsma
Branch: MAIN
Initial revision. Creates a metabolic model for given genome in MPS format - Kevin

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

my ($genome_id) = @ARGV;

my $debug = 0;
my $fig = new FIG;

unless (scalar (@ARGV) == 1)
{
    print STDERR "usage: make_MPS_model <genome_id>\n";
    exit(-1);
}
chomp $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_Config::var/Models/$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 %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;
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;
    }
}
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;
    }
}

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}})
    {
	print OUT "    $rid    $_             -$reaction_to_substrates{$rid}->{$_}\n";
    }
    foreach (sort keys %{$reaction_to_products{$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} && !$biomass_cpds{$cid})
    {
	print OUT "    T$temp    E$temp              1\n";
	print OUT "    T$temp    C$temp             -1\n";
    }
    elsif(!$output_cpds{$cid} && !$biomass_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";
}
#sink reactions here eh.


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3