[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.1 - (download) (as text) (annotate)
Tue Jul 31 16:07:23 2007 UTC (12 years, 10 months ago) by formsma
Branch: MAIN
First version of script to automatically dectect essential reactions in a generated model - 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(<FLUX>)
    {
	if($_ =~ /^Objective/)
	{
	    my @temp = split " ",$_;
	    chomp @temp;
	    if($temp[3] > 0)
	    {
		$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