View of /FigKernelScripts/check_model.pl

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

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

my $debug = 0;

unless (scalar (@ARGV) >= 1)
{
    print STDERR "usage: check_model <genome_id>\n";
    exit(-1);
}
chomp $genome_id;
my $fig;
if($orgdir)
{
    $fig = new FIGV($orgdir);
}
else
{
    $fig = new FIGV;
}
Scenario::set_fig($fig,$genome_id);
print STDERR "Running check_model on 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,0)};
my $filebase = $fig->model_directory($genome_id)."/Analysis/";

my %scenarios_used;
my %all_substrates;
foreach (@scenarios)
{
    map {$all_substrates{$_} = 1} @{$_->get_substrates_all()};
    $scenarios_used{$_->get_scenario_name()} = 0;
}

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

my %cofactor_cpds;
open(INPUTCO,$filebase."cofactors.txt")
    or die("Failed to open ".$filebase."cofactors.txt. Run find_model_cofactors $genome_id\n");
while(<INPUTCO>)
{
    chomp;
    my @line = split "\t", $_;
    chomp @line;
    $cofactor_cpds{$line[0]} = 1;
}
close(INPUTCO);

my %biomass_cpds;
open(BIOMASS,$filebase."biomass.txt")
    or die ("Failed to open ".$filebase."biomass.txt Run find_model_information $genome_id biomass");
while(<BIOMASS>)
{
    chomp;
    my @temp = split "\t", $_;
    chomp @temp;
    $biomass_cpds{$temp[0]} = 1;
}
close(BIOMASS);
print STDERR "\tStage 1 - Complete\n" if($debug);

print STDERR "\tStage 2 - Run Analysis Tools (fueled by caffine)\n" if($debug);
#Total list of what compounds we reach.
my @list_generated = keys %input_cpds;

my %possible_problems;
my %good_compounds;

# Here we are running all the scenarios we have inputs for
#  to analyze what biomass we can reach. We are also identifying 
#  any Scenarios that are not being utilized.

my %current_inputs;
foreach (keys %input_cpds)
{
    $current_inputs{$_} = 1;
}
my %current_cofactors;
my %previous_inputs;
foreach (keys %cofactor_cpds)
{
    $current_cofactors{$_} = 1;
}

while(scalar(keys %current_inputs))
{
    my ($temp_in,$temp_co) = generate_new_inputs(\%current_inputs,\%current_cofactors);
    #Check to make sure we have a different compound
    my $break = 1;
    foreach my $key (keys %$temp_in)
    {
	if(!defined $current_inputs{$key} && !defined $previous_inputs{$key})
	{
	    $break = 0;
	}
    }
    if($break) { last;} #break the while loop, we are looping

    undef %previous_inputs;
    map {$previous_inputs{$_} = 1;} sort keys %current_inputs;

    #undef %current_inputs;
    map {$current_inputs{$_} = 1;} sort keys %$temp_in;

    #undef %current_cofactors;
    map {$current_cofactors{$_} = 1;} sort keys %$temp_co;

    push @list_generated, keys %$temp_in;
    push @list_generated, keys %$temp_co;
    
}

foreach my $found (@list_generated)
{
    if(defined $biomass_cpds{$found})
    {
	$biomass_cpds{$found} = 2;
    }
}

# Analyze input compounds
# and also identify dead-end Scenarios
my %dead_end_scenarios;

foreach my $scenario (@scenarios)
{
    my $count = 0;
    my %used_inputs;
    foreach my $substrate (@{$scenario->get_substrates_all()})
    {
	if($input_cpds{$substrate} == 1)
	{
	    $used_inputs{$substrate} = 1;
	    $count++;
	}
    }

    if($count !=0 && $count != scalar(@{$scenario->get_substrates_all()}))
    {
       	map {$possible_problems{$_} = 1;} keys %used_inputs;
    }
    else
    {
	map {$good_compounds{$_} = 0;} keys %used_inputs;
    }

    #dead-end analysis
    my $bad_end = 1;
    foreach my $product (@{$scenario->get_products_all()})
    {
	if(defined $biomass_cpds{$product} || $all_substrates{$product} == 1)
	{
	    $bad_end = 0;
	    last;
	}
    }
    #Iff we reach this point, then this scenario produces no biomass
    # or compounds used by another scenario. It's all waste, could be a error.
    if($bad_end)
    {
	$dead_end_scenarios{$scenario->get_scenario_name()} = 1;
    }
}

open(RESULTS,">".$filebase."info_biomass.txt");
open(OUT, ">".$filebase."biomass_valid.txt");
foreach (sort keys %biomass_cpds)
{
    if($biomass_cpds{$_} != 2)
    {
	print RESULTS "Biomass $_ failed to be generated using inputs.txt\n";
    }
    else
    {
	print OUT "$_\n";
    }
}
close(RESULTS);
close(OUT);

open(OUT,">$filebase\scenarios_valid.txt");
open(SCEN_RESULTS,">".$filebase."info_scenarios.txt");
foreach my $scen_id (sort keys %scenarios_used)
{
    if($scenarios_used{$scen_id} == 0)
    {
	print SCEN_RESULTS "Scenario $scen_id unable to run with current inputs\n";
    }
    else
    {
	print OUT "$scen_id\n";
    }
}
print SCEN_RESULTS "\n\n#################\n\n\n";
foreach (sort keys %dead_end_scenarios)
{
    print SCEN_RESULTS "Scenario $_ is possibly a dead end(no outputs used by other scenarios or biomass.)\n";
}
close(SCEN_RESULTS);
close(OUT);

open(RESULTS_INPUT,">".$filebase."info_inputs.txt");
foreach (sort keys %possible_problems)
{
    if($possible_problems{$_} == 1 && !defined $good_compounds{$_})
    {
	print RESULTS_INPUT "Input $_ may be a false-positive. Perhaps a Scenario should generate this.\n";
    }
}
close(RESULTS_INPUT);

print STDERR "\tStage 2 - Complete! Results in info_inputs.txt, info_biomass.txt and info_scenarios.txt\n" if($debug);

sub generate_new_inputs
{
    my ($inputs,$cofactors) = @_;
    my %new_inputs;
    my %new_cofactors;
    foreach my $scenario (@scenarios)
    {
	my $count = 0;
	foreach my $substrate (@{$scenario->get_substrates_all()})
	{
	    if($inputs->{$substrate})
	    {
		$count++;
	    }
	    elsif( $cofactors->{$substrate})
	    {
		$count++;
	    }
	}

	if($count != scalar(@{$scenario->get_substrates_all()}))
	{
	    next;
	}
	else
	{
	    print STDERR "Scenario ".$scenario->get_scenario_name()." ran.\n" if($debug);
	    
	    ### here we want to mark what inputs we used 
	    foreach (@{$scenario->get_substrates()})
	    {
		$inputs->{$_} = 2;
	    }

	    $scenarios_used{$scenario->get_scenario_name()} = 1;
	    foreach my $product (@{$scenario->get_products()})
	    {
		$new_inputs{$product} = 1;
	    }
	    foreach my $product (@{$scenario->get_product_cofactors()})
	    {
		$new_cofactors{$product} = 1;
	    }
	}
    }
    #Add any inputs that aren't used to the return again
    foreach (keys %$inputs)
    {
	if($inputs->{$_} !=2)
	{
	    $new_inputs{$_} = 1;
	}
    }

    return (\%new_inputs,\%new_cofactors);
}