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

View of /FigKernelScripts/find_essential_reactions.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (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.2: +1 -1 lines
switch to scenario_directory

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

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

my $debug = 0;

#We want to basically generate a bunch of MPS files
# that are missing one reaction and determine what ones are nessisary for flux.
#


unless (scalar (@ARGV) >= 1)
{
    print STDERR "usage: find_essential_reactions <genome_id>\n";
    exit(-1);
}
my $fig; 
if($orgdir)
{
    $fig = new FIGV($orgdir);
}
else
{
    $fig = new FIGV;
}
chomp $genome_id;

my $filebase = $fig->scenario_directory($genome_id)."/Analysis/";

my %reactions_tested;
my %flux_values;
while(1)
{
#open file for reading and writing.
    open(NEW,">".$filebase."test.mps");
    open(MPS,"<".$filebase."$genome_id.mps");
    my $current_reaction = "";
    while(<MPS>)
    {
	chomp;
	if(!($_ =~ /(R\d{5}[\ |_r])/) || defined $reactions_tested{$1}
	   || ($1 ne $current_reaction && $current_reaction ne ""))
	{
	    print NEW $_."\n";
	}
	elsif($current_reaction eq "" && defined $1)
	{
	    print "Mached $1 as reaction\n";
	    $current_reaction = $1;
	}
    }
    close(MPS);
    close(NEW);
    if($current_reaction eq "") {last;}
    print "Testing reaction $current_reaction\n";
    $reactions_tested{$current_reaction} = 1;
    system(("/vol/biotools/bin/glpsol","--max",$filebase."test.mps","-o",$filebase."test.output"));
    
    #determine if we have flux!
    my $flux;
    open(OUT,$filebase."test.output");
    while(<OUT>)
    {
	if($_ =~ /^Objective/)
	{
	    my @temp = split " ",$_;
	    chomp @temp;
	    $flux_values{$current_reaction} = $temp[3];
	}
    }
    close(OUT);
}
open(OUT,">".$filebase."essential_reactions.txt");
foreach (sort keys %flux_values)
{
    if($flux_values{$_} == 0)
    {
	print OUT "$_ needed for biomass\n";
    }
}
close(OUT);
system(("rm",$filebase."test.mps"));

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3