[Bio] / FigKernelPackages / model.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/model.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.20, Mon Nov 3 18:38:44 2008 UTC revision 1.21, Wed Nov 5 15:24:12 2008 UTC
# Line 38  Line 38 
38  use strict;  use strict;
39  use Scenario;  use Scenario;
40  use FIG;  use FIG;
41    use FIGV;
42  use Subsystem;  use Subsystem;
43  use File::Path;  use File::Path;
44    
# Line 63  Line 64 
64    
65  #Flip this bit to enable debugging  #Flip this bit to enable debugging
66  my $debug = int($ENV{HOPE_DEBUG});  my $debug = int($ENV{HOPE_DEBUG});
67  #$debug = 0;  $debug = 0;
68    
69  sub set_fig  sub set_fig
70  {  {
# Line 3869  Line 3870 
3870      return \%peg_to_scenario;      return \%peg_to_scenario;
3871  }  }
3872    
3873    sub analyze_network_for_genome
3874    {
3875        my ($genome_id, $transports, $biomass) = @_;
3876    
3877        my ($node_type, $node_connections);
3878        my $debug = 0;
3879        my $unreachable = 0;
3880        chomp $genome_id;
3881        my $fig = new FIGV;
3882    
3883        my @transports = split ("_", $transports);
3884        my @biomass = split ("_", $biomass);
3885    
3886        my %compound_id_to_name;
3887        open(IN,"$FIG_Config::global/Models/All/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt");
3888        while(<IN>)
3889        {
3890            chomp($_);
3891            my($id,$name) = split("\t");
3892            $compound_id_to_name{$id} = $name;
3893        }
3894        close(IN);
3895    
3896        Scenario::set_fig($fig,$genome_id);
3897        print STDERR "Running check_model_hope on genome $genome_id\n" if($debug);
3898        print STDERR "\tStage 1 - Loading Scenarios, transports and biomass\n" if($debug);
3899        my @scenarios = @{Scenario->get_genome_scenarios($genome_id,0)};
3900    
3901        my %scenarios_used;
3902        my %scenarios_id;
3903    
3904        foreach (@scenarios)
3905        {
3906            $scenarios_used{$_->get_scenario_name()} = 0;
3907            $scenarios_id{$_->get_id()} = 0;
3908        }
3909    
3910        my %biomass_lib_cpds;
3911        foreach my $biomass_cpd (@biomass)
3912        {
3913            $biomass_lib_cpds{$biomass_cpd} = 1;
3914        }
3915    
3916        my (%transportable_cpds, %transportable_in_biomass);
3917        foreach my $transport_cpd (@transports)
3918        {
3919            if (defined $biomass_lib_cpds{$transport_cpd})
3920            {
3921                $transportable_in_biomass{$transport_cpd} = 1;
3922            }
3923            else
3924            {
3925                $transportable_cpds{$transport_cpd} = 1;
3926            }
3927        }
3928    
3929        print STDERR "\tStage 1 - Complete\n" if($debug);
3930        print STDERR "\tStage 2 - Run Analysis Tools\n" if($debug);
3931    
3932    #   Total list of what compounds we reach.
3933        my %list_generated;
3934    
3935    #   Here we are running all the scenarios we have inputs for
3936    #   to analyze what biomass we can reach. We are also identifying
3937    #   any Scenarios that are not being utilized.
3938    
3939        my %current_inputs;
3940        foreach (keys %transportable_cpds)
3941        {
3942            # keep track of the paths that lead to each compound
3943            $current_inputs{$_}->{"[]transport of ".$_."/".$compound_id_to_name{$_}."[$_]"} = 1;
3944        }
3945    
3946        while(1)
3947        {
3948            my ($temp_in,$break) = generate_new_inputs(\%current_inputs);
3949    
3950            foreach my $cpd (keys %$temp_in)
3951            {
3952                map { $current_inputs{$cpd}->{$_} = 1 }  keys %{$temp_in->{$cpd}};
3953            }
3954    
3955            map { $list_generated{$_} = 1 } keys %$temp_in;
3956    
3957            print STDERR "Rechecking with inputs\n" if($debug);
3958            if($break)
3959            {
3960                last;
3961            }
3962        }
3963        print STDERR "Finished checking organism scenarios\n" if($debug);
3964    
3965        my %biomass_reached;
3966    
3967        foreach my $found (keys %list_generated)
3968        {
3969            if(defined $biomass_lib_cpds{$found})
3970            {
3971                $biomass_reached{$found} = $current_inputs{$found};
3972            }
3973        }
3974    
3975        foreach (keys %transportable_cpds)
3976        {
3977            push (@{$node_type->{'transport'}}, $_."/".$compound_id_to_name{$_});
3978        }
3979    
3980        print STDERR "\tStage 2 - Complete!\n" if($debug);
3981        print STDERR "\tStage 3 - Loading All Scenarios\n" if($debug);
3982    
3983        @scenarios = @{Scenario->get_genome_scenarios("All",0)};
3984    
3985        my %old_scenarios_used = %scenarios_used;
3986        my %old_scenarios_id = %scenarios_id;
3987    
3988        my %old_biomass_reached;
3989        foreach my $cid (keys %biomass_reached)
3990        {
3991            my %inner_hash_copy = %{$biomass_reached{$cid}};
3992            $old_biomass_reached{$cid} = \%inner_hash_copy
3993        }
3994    
3995        my %old_current_inputs;
3996        foreach my $cid (keys %current_inputs)
3997        {
3998            my %inner_hash_copy = %{$current_inputs{$cid}};
3999            $old_current_inputs{$cid} = \%inner_hash_copy;
4000        }
4001    
4002    
4003        %list_generated = {};
4004    
4005        foreach (@scenarios)
4006        {
4007            $scenarios_used{$_->get_scenario_name()} = 0 if ! exists $scenarios_used{$_->get_scenario_name()};
4008            $scenarios_id{$_->get_id()} = 0 if ! exists $scenarios_id{$_->get_id()};
4009        }
4010    
4011        while(1)
4012        {
4013            my ($temp_in,$break) = generate_new_inputs(\%current_inputs);
4014    
4015            foreach my $cpd (keys %$temp_in)
4016            {
4017                map { $current_inputs{$cpd}->{$_} = 1 } keys %{$temp_in->{$cpd}};
4018            }
4019    
4020            map { $list_generated{$_} = 1 } keys %$temp_in;
4021    
4022            print STDERR "Rechecking with inputs\n" if($debug);
4023            if($break)
4024            {
4025                last;
4026            }
4027        }
4028        print STDERR "Finished checking all scenarios\n" if($debug);
4029    
4030        my %biomass_reached;
4031    
4032        foreach my $found (keys %list_generated)
4033        {
4034            if(defined $biomass_lib_cpds{$found})
4035            {
4036                $biomass_reached{$found} = $current_inputs{$found};
4037            }
4038        }
4039    
4040        my (%cpds_to_print);
4041        my %org_cpds_to_check;
4042        map { $org_cpds_to_check{$_} = 1 } keys %old_biomass_reached;
4043        foreach my $cpd (sort keys %biomass_reached)
4044        {
4045            next if exists $old_biomass_reached{$cpd};
4046    
4047            $cpds_to_print{$cpd} = 1;
4048            my $done = 0;
4049    
4050            my @cpd_list = ($cpd);
4051            while (! $done)
4052            {
4053                my @new_substrates = collect_substrates(\%current_inputs, @cpd_list);
4054                @cpd_list = ();
4055                foreach my $substrate (@new_substrates)
4056                {
4057                if (exists $old_current_inputs{$substrate})
4058                {
4059                    $org_cpds_to_check{$substrate} = 1;
4060                    next;
4061                }
4062                push (@cpd_list, $substrate) if ! exists $cpds_to_print{$substrate};
4063                $cpds_to_print{$substrate} = 1;
4064                }
4065                $done = 1 if scalar @new_substrates == 0;
4066            }
4067        }
4068    
4069        foreach my $cpd (keys %cpds_to_print)
4070        {
4071            foreach my $path (keys %{$current_inputs{$cpd}})
4072            {
4073                my ($substrates,$scenario,$products) = $path =~ /\[(.*)\](.*)\[(.*)\]/;
4074    
4075                foreach my $substrate (split ",", $substrates)
4076                {
4077                    foreach my $product (split ",", $products)
4078                    {
4079                        $node_connections->{$substrate."/".$compound_id_to_name{$substrate}}->$scenario = 1;
4080                        $node_connections->{$scenario}->$product."/".$compound_id_to_name{$product} = 1;
4081                    }
4082                }
4083            }
4084        }
4085    
4086        %cpds_to_print = ();
4087    
4088        foreach my $cpd (sort keys %org_cpds_to_check)
4089        {
4090            $cpds_to_print{$cpd} = 1;
4091            my $done = 0;
4092    
4093            my @cpd_list = ($cpd);
4094            while (! $done)
4095            {
4096                my @new_substrates = collect_substrates(\%old_current_inputs, @cpd_list);
4097                @cpd_list = ();
4098                foreach my $substrate (@new_substrates)
4099                {
4100                    push (@cpd_list, $substrate) if ! exists $cpds_to_print{$substrate};
4101                    $cpds_to_print{$substrate} = 1;
4102                }
4103                $done = 1 if scalar @new_substrates == 0;
4104            }
4105        }
4106    
4107        foreach my $cpd (keys %cpds_to_print)
4108        {
4109            foreach my $path (keys %{$old_current_inputs{$cpd}})
4110            {
4111                my ($substrates,$scenario,$products) = $path =~ /\[(.*)\](.*)\[(.*)\]/;
4112    
4113                foreach my $substrate (split ",", $substrates)
4114                {
4115                    foreach my $product (split ",", $products)
4116                    {
4117                        $node_connections->{$substrate."/".$compound_id_to_name{$substrate}}->$scenario = 1;
4118                        $node_connections->{$scenario}->$product."/".$compound_id_to_name{$product} = 1;
4119                    }
4120                }
4121            }
4122        }
4123    
4124        foreach (sort keys %biomass_lib_cpds)
4125        {
4126            if(! exists $biomass_reached{$_} && ! exists $old_biomass_reached{$_})
4127            {
4128                print STDERR "Biomass $_/".$compound_id_to_name{$_}." not reachable using current transports" if ($debug);
4129                $unreachable++;
4130                if (exists $transportable_in_biomass{$_})
4131                {
4132                    print STDERR " (transportable)" if ($debug);
4133                    push (@{$node_type->{'biomass not reached (transportable)'}}, $_."/".$compound_id_to_name{$_});
4134                }
4135                else
4136                {
4137                    push (@{$node_type->{'biomass not reached'}}, $_."/".$compound_id_to_name{$_});
4138                }
4139                print STDERR "\n" if ($debug);
4140            }
4141            elsif ($old_biomass_reached{$_})
4142            {
4143                push (@{$node_type->{'biomass reached'}}, $_."/".$compound_id_to_name{$_});
4144            }
4145            else
4146            {
4147                # do the paths all include a non-organism scenario as the last step?
4148                my $last_not_org = 1;
4149    
4150                foreach my $path_set (keys %{$biomass_reached{$_}})
4151                {
4152                    my @paths = split ";", $path_set;
4153                    my $last_path = $paths[-1];
4154                    my ($substrates,$scenario,$products) = $last_path =~ /\[(.*)\](.*)\[(.*)\]/;
4155    
4156                    if (exists $old_scenarios_used{$scenario})
4157                    {
4158                        $last_not_org = 0;
4159                        last;
4160                    }
4161                }
4162    
4163                if ($last_not_org == 1)
4164                {
4165                    push (@{$node_type->{'biomass reached all last_not_org'}}, $_."/".$compound_id_to_name{$_});
4166                }
4167                else
4168                {
4169                    push (@{$node_type->{'biomass reached all'}}, $_."/".$compound_id_to_name{$_});
4170                }
4171            }
4172        }
4173    
4174        foreach my $scen_name (sort keys %scenarios_used)
4175        {
4176            if($scenarios_used{$scen_name} == 1)
4177            {
4178                if (exists ($old_scenarios_used{$scen_name}))
4179                {
4180                    push (@{$node_type->{'scenario present'}}, $scen_name);
4181                }
4182                else
4183                {
4184                    push (@{$node_type->{'scenario present all'}}, $scen_name);
4185                }
4186            }
4187        }
4188    
4189        foreach my $cpd (keys %current_inputs)
4190        {
4191            # compounds present with All scenarios that aren't biomass
4192            if (! exists($old_current_inputs{$cpd}) && ! exists ($biomass_reached{$cpd}))
4193            {
4194                push (@{$node_type->{'compound present all'}}, $cpd."/".$compound_id_to_name{$cpd});
4195            }
4196        }
4197    
4198        print STDERR "\tStage 3 - Complete!\n" if($debug);
4199    
4200        return ($node_type, $node_connections);
4201    
4202        sub generate_new_inputs
4203        {
4204            my ($inputs) = @_;
4205            my %new_inputs;
4206            my $break = 1;
4207    
4208            my %used_scenario_paths_this_time;
4209    
4210            foreach my $scenario (@scenarios)
4211            {
4212                if($scenarios_id{$scenario->get_id()} == 1) # || $scenarios_used{$scenario->get_scenario_name()} != 0)
4213                {
4214                    next;
4215                }
4216    
4217                my $sids = join ",", sort @{$scenario->get_substrates};
4218                my $pids = join ",", sort @{$scenario->get_products};
4219                my $cid_scenario_name = "[".$sids."]".$scenario->get_scenario_name()."[".$pids."]";
4220    
4221                my @matched;
4222    
4223                foreach my $substrate (@{$scenario->get_substrates()})
4224                {
4225                    if(exists $inputs->{$substrate})
4226                    {
4227                        push @matched, $substrate;
4228                    }
4229                }
4230    
4231                if (scalar @matched != scalar(@{$scenario->get_substrates()}))
4232                {
4233                    next;
4234                }
4235                else
4236                {
4237                    print STDERR "Scenario ".$cid_scenario_name." ran.\n" if($debug);
4238                    $break = 0;
4239                    $scenarios_id{$scenario->get_id()} = 1;
4240                    $scenarios_used{$scenario->get_scenario_name()} = 1;
4241    
4242                    # local determination that we've used an equivalent path already in this
4243                    # time through the subroutine, so no need to process further
4244                    next if (exists $used_scenario_paths_this_time{$cid_scenario_name});
4245                    $used_scenario_paths_this_time{$cid_scenario_name} = 1;
4246    
4247                    foreach my $product (@{$scenario->get_products()})
4248                    {
4249                        $new_inputs{$product}->{$cid_scenario_name} = 1;
4250                    }
4251                }
4252            }
4253            return (\%new_inputs,$break);
4254        }
4255    
4256        sub collect_substrates
4257        {
4258            my $input_hash = shift @_;
4259            my @cpds = @_;
4260            my %substrate_hash;
4261    
4262            foreach my $cpd (@cpds)
4263            {
4264                foreach my $path (keys %{$input_hash->{$cpd}})
4265                {
4266                    my ($substrates,$scenario,$products) = $path =~ /\[(.*)\](.*)\[(.*)\]/;
4267    
4268                    foreach my $substrate (split ",", $substrates)
4269                    {
4270                        $substrate_hash{$substrate} = 1;
4271                    }
4272                }
4273            }
4274            return keys %substrate_hash;
4275    
4276    
4277        }
4278    }
4279    
4280  sub get_reactions_for_genome_subsystem  sub get_reactions_for_genome_subsystem
4281  {  {
4282      my ($genome_id, $ss_name) = @_;      my ($genome_id, $ss_name) = @_;

Legend:
Removed from v.1.20  
changed lines
  Added in v.1.21

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3