[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.2 - (download) (as text) (annotate)
Tue Jul 31 19:47:02 2007 UTC (12 years, 9 months ago) by formsma
Branch: MAIN
CVS Tags: rast_rel_2008_06_18, rast_rel_2008_06_16, rast_rel_2008_07_21, rast_2008_0924, rast_rel_2008_04_23, rast_rel_2008_09_30, mgrast_rel_2008_0924, mgrast_rel_2008_1110_v2, mgrast_rel_2008_0625, rast_rel_2008_10_09, rast_release_2008_09_29, mgrast_rel_2008_0806, mgrast_rel_2008_0923, mgrast_rel_2008_0919, mgrast_rel_2008_1110, rast_rel_2008_09_29, mgrast_rel_2008_0917, rast_rel_2008_10_29, rast_rel_2008_08_07
Changes since 1.1: +2 -6 lines
Fixed incorrect filehandle. Tested Good - Kevin

# _*_ 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->model_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