[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.32, Thu Apr 2 01:19:02 2009 UTC revision 1.33, Mon Nov 16 17:56:21 2009 UTC
# Line 1  Line 1 
   
1  ## _*_ Perl _*_ ##  ## _*_ Perl _*_ ##
2  #  #
3  # model.pm  # model.pm
# Line 41  Line 40 
40  use FIG;  use FIG;
41  use Subsystem;  use Subsystem;
42  use File::Path;  use File::Path;
43    use Data::Dumper;
44    
45  my $fig = eval("FIG->new()"); ## So that this will compile properly on non-FIG systems!  my $fig = eval("FIG->new()"); ## So that this will compile properly on non-FIG systems!
46    
# Line 150  Line 150 
150      }      }
151      else      else
152      {      {
153          my %reactions_for_genome = get_reactions_for_genome_in_subsystem($genome,$ss_name);          my $reactions_for_genome = get_reactions_for_genome_in_subsystem($genome,$ss_name);
154          map { $ss_reactions{$_} = 1 } keys %reactions_for_genome if keys %reactions_for_genome;          map { $ss_reactions{$_} = 1 } keys %$reactions_for_genome if keys %$reactions_for_genome;
155      }      }
156    
157      map { $sc_inputs{$_} = 1 } @{$scenario_data->{inputs}};      map { $sc_inputs{$_} = 1 } @{$scenario_data->{inputs}};
# Line 2695  Line 2695 
2695    
2696  sub load_superset_file  sub load_superset_file
2697  {  {
2698        # only need to load once during execution context
2699        return \%ss_to_superset if defined %ss_to_superset;
2700    
2701      %superset_to_ss = ();      %superset_to_ss = ();
2702      %ss_to_superset = ();      %ss_to_superset = ();
2703    
# Line 2709  Line 2712 
2712          $ss_to_superset{$line[1]} = $line[0];          $ss_to_superset{$line[1]} = $line[0];
2713          push(@{$superset_to_ss{$line[0]}},$line[1]);          push(@{$superset_to_ss{$line[0]}},$line[1]);
2714      }      }
2715      close(FILE);      close FILE;
2716    
2717      return \%ss_to_superset      return \%ss_to_superset
2718  }  }
# Line 2777  Line 2780 
2780      }      }
2781      else      else
2782      {      {
2783          my %reactions_for_genome = get_reactions_for_genome_in_subsystem($genome,$ss_name);          my $reactions_for_genome = get_reactions_for_genome_in_subsystem($genome,$ss_name);
2784          map { $ss_reactions{$_} = 1 } keys %reactions_for_genome if keys %reactions_for_genome;          map { $ss_reactions{$_} = 1 } keys %$reactions_for_genome if keys %$reactions_for_genome;
2785      }      }
2786    
2787      # first find paths in the All directory that should be valid for the genome      # first find paths in the All directory that should be valid for the genome
# Line 3063  Line 3066 
3066    
3067      foreach my $name (@ss_scenarios)      foreach my $name (@ss_scenarios)
3068      {      {
3069            print STDERR "\t Running $name\n";
3070          run_scenario($genome,$superset,$subsystem,$name,$find_first);          run_scenario($genome,$superset,$subsystem,$name,$find_first);
3071      }      }
3072  }  }
# Line 3339  Line 3343 
3343          {          {
3344              $length = 4;              $length = 4;
3345          }          }
         print STDERR "Expanding $path \n" if $debug;  
3346          my @parts = split "/", $path;          my @parts = split "/", $path;
3347          shift @parts; # get ride of the first blank entry from /$genome          shift @parts; # get ride of the first blank entry from /$genome
3348          $length =$length - scalar @parts;          $length =$length - scalar @parts;
3349          my $genome = shift @parts;          my $genome = shift @parts;
3350          $path = join "/" , @parts;          $path = join "/" , @parts;
         print STDERR "Length : $length Path: $path\n" if $debug;  
3351          my @temp = @{expand_recursive($genome, $path, $length)};          my @temp = @{expand_recursive($genome, $path, $length)};
3352          push @final_paths , @temp if scalar @temp > 0;          push @final_paths , @temp if scalar @temp > 0;
3353      }      }
# Line 3366  Line 3368 
3368      {      {
3369          #read this path, and pull out all the sub-directories.          #read this path, and pull out all the sub-directories.
3370          my $model_dir = get_scenario_directory($genome) . "/$path/";          my $model_dir = get_scenario_directory($genome) . "/$path/";
3371            return [] if (! -d $model_dir);
         print STDERR "reading directory $model_dir\n" if $debug;  
3372          opendir (DIR, $model_dir) or die("$model_dir");          opendir (DIR, $model_dir) or die("$model_dir");
         print STDERR "reading directory $model_dir\n" if $debug;  
3373          @sub_dirs  = readdir DIR;          @sub_dirs  = readdir DIR;
3374          close DIR;          close DIR;
         print STDERR "Found: " if $debug;  
         print STDERR @sub_dirs , "," if $debug;  
3375      }      }
3376      else      else
3377      {      {
# Line 3385  Line 3383 
3383      foreach my $sub_path (@sub_dirs)      foreach my $sub_path (@sub_dirs)
3384      {      {
3385          next if $sub_path =~ /^\.$/ || $sub_path =~ /^\.\.$/;          next if $sub_path =~ /^\.$/ || $sub_path =~ /^\.\.$/;
         print STDERR "Calling on $sub_path , $count \n" if $debug;  
3386          push @to_return , @{expand_recursive($genome, "$path/$sub_path",$count)};          push @to_return , @{expand_recursive($genome, "$path/$sub_path",$count)};
3387      }      }
     print STDERR "Returning" ,@to_return if $debug;  
3388      return \@to_return;      return \@to_return;
3389    
3390  }  }
# Line 3798  Line 3794 
3794      {      {
3795          my $subsystem = $fig->get_subsystem($subsystem_name);          my $subsystem = $fig->get_subsystem($subsystem_name);
3796          next if(!defined $subsystem);          next if(!defined $subsystem);
3797          my %reactions_for_ss = get_reactions_for_genome_in_subsystem($genome_id,$subsystem_name);          my $reactions_for_ss = get_reactions_for_genome_in_subsystem($genome_id,$subsystem_name);
3798          next if(! keys %reactions_for_ss);          next if(! keys %$reactions_for_ss);
3799          foreach my $reaction (keys %reactions_for_ss)          foreach my $reaction (keys %$reactions_for_ss)
3800          {          {
3801              if(defined $reaction_to_pegs{$reaction})              if(defined $reaction_to_pegs{$reaction})
3802              {              {
3803                  push @{$reaction_to_pegs{$reaction}} , @{$reactions_for_ss{$reaction}};                  push @{$reaction_to_pegs{$reaction}} , @{$reactions_for_ss->{$reaction}};
3804              }              }
3805              else              else
3806              {              {
3807                  $reaction_to_pegs{$reaction} = $reactions_for_ss{$reaction};                  $reaction_to_pegs{$reaction} = $reactions_for_ss->{$reaction};
3808              }              }
3809          }          }
3810      }      }
# Line 3854  Line 3850 
3850  {  {
3851      my ($genome_id, $transports, $biomass) = @_;      my ($genome_id, $transports, $biomass) = @_;
3852    
3853      my ($node_type, $node_connections);      my ($node_type, $node_connections, %compound_ids, $path_info, $reaction_presence);
     my $debug = 0;  
3854      my $unreachable = 0;      my $unreachable = 0;
3855      chomp $genome_id;      chomp $genome_id;
3856    
3857      my @transports = split ("_", $transports);      Scenario->set_fig($fig);
     my @biomass = split ("_", $biomass);  
   
     my %compound_id_to_name;  
     open(IN,"$FIG_Config::global/Models/All/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt");  
     while(<IN>)  
     {  
         chomp($_);  
         my($id,$name) = split("\t");  
         $compound_id_to_name{$id} = $name;  
     }  
     close(IN);  
   
     Scenario::set_fig($fig,$genome_id);  
3858      print STDERR "Running check_model_hope on genome $genome_id\n" if($debug);      print STDERR "Running check_model_hope on genome $genome_id\n" if($debug);
3859      print STDERR "\tStage 1 - Loading Scenarios, transports and biomass\n" if($debug);      print STDERR "\tStage 1 - Loading Scenarios, transports and biomass\n" if($debug);
3860      my @scenarios = @{Scenario->get_genome_scenarios($genome_id,0)};      my @scenarios = @{Scenario->get_genome_scenarios($genome_id,0)};
# Line 3882  Line 3864 
3864    
3865      foreach (@scenarios)      foreach (@scenarios)
3866      {      {
3867            print STDERR "\t\tLoading ", $_->get_scenario_name(), "\n" if $debug;
3868          $scenarios_used{$_->get_scenario_name()} = 0;          $scenarios_used{$_->get_scenario_name()} = 0;
3869          $scenarios_id{$_->get_id()} = 0;          $scenarios_id{$_->get_id()} = 0;
3870      }      }
3871    
3872      my %biomass_lib_cpds;      my %biomass_lib_cpds;
3873      foreach my $biomass_cpd (@biomass)      foreach my $biomass_cpd (@$biomass)
3874      {      {
3875          $biomass_lib_cpds{$biomass_cpd} = 1;          $biomass_lib_cpds{$biomass_cpd} = 1;
3876      }      }
3877    
3878      my (%transportable_cpds, %transportable_in_biomass);      my %transportable_cpds = map { $_ => 1 } @$transports;
     foreach my $transport_cpd (@transports)  
     {  
         if (defined $biomass_lib_cpds{$transport_cpd})  
         {  
             $transportable_in_biomass{$transport_cpd} = 1;  
         }  
         else  
         {  
             $transportable_cpds{$transport_cpd} = 1;  
         }  
     }  
3879    
3880      print STDERR "\tStage 1 - Complete\n" if($debug);      print STDERR "\tStage 1 - Complete\n" if($debug);
3881      print STDERR "\tStage 2 - Run Analysis Tools\n" if($debug);      print STDERR "\tStage 2 - Run Analysis Tools\n" if($debug);
# Line 3915  Line 3887 
3887  #   to analyze what biomass we can reach. We are also identifying  #   to analyze what biomass we can reach. We are also identifying
3888  #   any Scenarios that are not being utilized.  #   any Scenarios that are not being utilized.
3889    
3890        # current_inputs contains the compounds that we have "reached" so far in
3891        # traversing the network; we preload transported compounds, as well as
3892        # some cofactors that are "main" inputs to scenarios
3893      my %current_inputs;      my %current_inputs;
3894        $current_inputs{"C00342"}->{"[]cofactor (thioredoxin)[C00342]"} = 1;
3895    
3896    
3897      foreach (keys %transportable_cpds)      foreach (keys %transportable_cpds)
3898      {      {
3899          # keep track of the paths that lead to each compound          # keep track of the paths that lead to each compound
3900          $current_inputs{$_}->{"[]transport of ".$_."/".$compound_id_to_name{$_}."[$_]"} = 1;          my @names = $fig->names_of_compound($_);
3901            my $name = @names > 0 ? $names[0] : $_;
3902            $current_inputs{$_}->{"[]transport of ${name}[$_]"} = 1;
3903      }      }
3904    
3905        my $iter_num = 1;
3906    
3907      while(1)      while(1)
3908      {      {
3909          my ($temp_in,$break) = generate_new_inputs(\%current_inputs);          my ($temp_in,$break) = generate_new_inputs(\%current_inputs);
# Line 3933  Line 3915 
3915    
3916          map { $list_generated{$_} = 1 } keys %$temp_in;          map { $list_generated{$_} = 1 } keys %$temp_in;
3917    
         print STDERR "Rechecking with inputs\n" if($debug);  
3918          if($break)          if($break)
3919          {          {
3920              last;              last;
3921          }          }
3922    
3923            print STDERR "Rechecking organism scenarios with inputs\n" if($debug);
3924            $iter_num++;
3925      }      }
3926      print STDERR "Finished checking organism scenarios\n" if($debug);      print STDERR "Finished checking organism scenarios\n" if($debug);
3927    
# Line 3953  Line 3937 
3937    
3938      foreach (keys %transportable_cpds)      foreach (keys %transportable_cpds)
3939      {      {
3940          push (@{$node_type->{'transport'}}, $_."/".$compound_id_to_name{$_});          push (@{$node_type->{'transport'}}, $_);
3941            $compound_ids{$_} = 1;
3942      }      }
3943    
3944      print STDERR "\tStage 2 - Complete!\n" if($debug);      print STDERR "\tStage 2 - Complete!\n" if($debug);
# Line 3989  Line 3974 
3974    
3975      while(1)      while(1)
3976      {      {
3977            print STDERR "Rechecking all scenarios with inputs\n" if($debug);
3978            $iter_num++;
3979    
3980          my ($temp_in,$break) = generate_new_inputs(\%current_inputs);          my ($temp_in,$break) = generate_new_inputs(\%current_inputs);
3981    
3982          foreach my $cpd (keys %$temp_in)          foreach my $cpd (keys %$temp_in)
# Line 3998  Line 3986 
3986    
3987          map { $list_generated{$_} = 1 } keys %$temp_in;          map { $list_generated{$_} = 1 } keys %$temp_in;
3988    
         print STDERR "Rechecking with inputs\n" if($debug);  
3989          if($break)          if($break)
3990          {          {
3991              last;              last;
# Line 4055  Line 4042 
4042              {              {
4043                  foreach my $product (split ",", $products)                  foreach my $product (split ",", $products)
4044                  {                  {
4045                      $node_connections->{$substrate."/".$compound_id_to_name{$substrate}}->$scenario = 1;                      print STDERR "cpd: $cpd  scenario: $scenario  path: $path  substrate: $substrate  product: $product\n" if ($debug);
4046                      $node_connections->{$scenario}->$product."/".$compound_id_to_name{$product} = 1;                      $node_connections->{$substrate}->{$scenario} = 1;
4047                        $compound_ids{$substrate} = 1;
4048                        $node_connections->{$scenario}->{$product} = 1;
4049                        $compound_ids{$product} = 1;
4050                  }                  }
4051              }              }
4052          }          }
# Line 4093  Line 4083 
4083              {              {
4084                  foreach my $product (split ",", $products)                  foreach my $product (split ",", $products)
4085                  {                  {
4086                      $node_connections->{$substrate."/".$compound_id_to_name{$substrate}}->$scenario = 1;                      $node_connections->{$substrate}->{$scenario} = 1;
4087                      $node_connections->{$scenario}->$product."/".$compound_id_to_name{$product} = 1;                      $compound_ids{$substrate} = 1;
4088                        $node_connections->{$scenario}->{$product} = 1;
4089                        $compound_ids{$product} = 1;
4090                  }                  }
4091              }              }
4092          }          }
# Line 4102  Line 4094 
4094    
4095      foreach (sort keys %biomass_lib_cpds)      foreach (sort keys %biomass_lib_cpds)
4096      {      {
4097          if(! exists $biomass_reached{$_} && ! exists $old_biomass_reached{$_})          if ($old_biomass_reached{$_})
4098          {          {
4099              print STDERR "Biomass $_/".$compound_id_to_name{$_}." not reachable using current transports" if ($debug);              push (@{$node_type->{'biomass reached'}}, $_);
4100              $unreachable++;              $compound_ids{$_} = 1;
             if (exists $transportable_in_biomass{$_})  
             {  
                 print STDERR " (transportable)" if ($debug);  
                 push (@{$node_type->{'biomass not reached (transportable)'}}, $_."/".$compound_id_to_name{$_});  
4101              }              }
4102              else          elsif ($biomass_reached{$_})
4103              {              {
4104                  push (@{$node_type->{'biomass not reached'}}, $_."/".$compound_id_to_name{$_});              push (@{$node_type->{'biomass reached all'}}, $_);
4105                $compound_ids{$_} = 1;
4106              }              }
             print STDERR "\n" if ($debug);  
4107          }          }
         elsif ($old_biomass_reached{$_})  
         {  
             push (@{$node_type->{'biomass reached'}}, $_."/".$compound_id_to_name{$_});  
         }  
         else  
         {  
             # do the paths all include a non-organism scenario as the last step?  
             my $last_not_org = 1;  
4108    
4109              foreach my $path_set (keys %{$biomass_reached{$_}})      my %hope_ss_info_cache;
             {  
                 my @paths = split ";", $path_set;  
                 my $last_path = $paths[-1];  
                 my ($substrates,$scenario,$products) = $last_path =~ /\[(.*)\](.*)\[(.*)\]/;  
4110    
4111                  if (exists $old_scenarios_used{$scenario})      foreach my $ss_and_scen (sort keys %scenarios_used)
4112                  {                  {
4113                      $last_not_org = 0;          if($scenarios_used{$ss_and_scen} == 1)
4114                      last;          {
4115                  }              print STDERR "$ss_and_scen was used\n";
4116                my ($ss_name, $scen_name) = split "/", $ss_and_scen;
4117    
4118                if (! exists $hope_ss_info_cache{$ss_name}) {
4119                    $hope_ss_info_cache{$ss_name}->{"kegg_info"} = $fig->get_scenario_info($scen_name);
4120                    $hope_ss_info_cache{$ss_name}->{"org_reactions"} = get_reactions_for_genome_in_subsystem_fast($genome_id, $ss_name);
4121              }              }
4122    
4123              if ($last_not_org == 1)              my $org_ss_reactions = $hope_ss_info_cache{$ss_name}->{"org_reactions"};
4124                my @scenarios; # path info for the scenarios
4125    
4126                if (exists ($old_scenarios_used{$ss_and_scen}))
4127              {              {
4128                  push (@{$node_type->{'biomass reached all last_not_org'}}, $_."/".$compound_id_to_name{$_});                  push (@{$node_type->{'scenario present'}}, $ss_and_scen);
4129                    @scenarios = @{Scenario->get_genome_scenarios_for_scenario($genome_id, $ss_name, $scen_name, 0)};
4130              }              }
4131              else              else
4132              {              {
4133                  push (@{$node_type->{'biomass reached all'}}, $_."/".$compound_id_to_name{$_});                  # determine how many reactions the organism is missing
4134              }                  @scenarios = @{Scenario->get_genome_scenarios_for_scenario("All", $ss_name, $scen_name, 0)};
4135                    my @percentiles = ("0", "25", "50", "75");
4136                    my $percentile = 3; # assume the greatest number of missing reactions
4137    
4138                    foreach my $scenario (@scenarios)
4139                    {
4140                        my @reactions = @{$scenario->get_path_info()};
4141                        map {s/_\w//} @reactions;
4142                        my $num_missing_rxns = 0;
4143                        map {$num_missing_rxns++ unless exists $org_ss_reactions->{$_}} @reactions;
4144                        my $sc_percentile = int( ($num_missing_rxns / scalar @reactions) * 4 );
4145                        $percentile = $sc_percentile < $percentile ? $sc_percentile : $percentile;
4146          }          }
4147    
4148                    push (@{$node_type->{"scenario present all $percentiles[$percentile]"}}, $ss_and_scen);
4149      }      }
4150    
4151      foreach my $scen_name (sort keys %scenarios_used)              # get the path info for scenarios
4152      {              my @scenario_path_info;
4153          if($scenarios_used{$scen_name} == 1)  
4154                foreach my $scenario (@scenarios)
4155          {          {
4156              if (exists ($old_scenarios_used{$scen_name}))                  my $map_ids = $hope_ss_info_cache{$ss_name}->{"kegg_info"}->{"scenarios"}->{$scen_name}->{"map_ids"};
4157                    my $reactions = $scenario->get_path_info();
4158                    my @rxns_and_cpds;
4159    
4160                    foreach my $reaction (@$reactions)
4161              {              {
4162                  push (@{$node_type->{'scenario present'}}, $scen_name);                      my ($rxn_id, $dir) = split "_", $reaction;
4163                        my ($subside, $prodside) = $dir eq "R" ? (0,1) : (1,0);
4164                        my (@main_subs, @main_prods, @reaction_info);
4165                        my @subs_info = $fig->reaction2comp($rxn_id,$subside,$map_ids);
4166                        map { if ($_->[2]) { push @main_subs, $_->[0]; $compound_ids{$_->[0]} = 1 } } @subs_info;
4167                        if (@main_subs == 0) {
4168                            # no main info
4169                            map { push @main_subs, $_->[0]; $compound_ids{$_->[0]} = 1 } @subs_info;
4170              }              }
4171              else                      my @prods_info = $fig->reaction2comp($rxn_id,$prodside,$map_ids);
4172              {                      map { if ($_->[2]) { push @main_prods, $_->[0]; $compound_ids{$_->[0]} = 1 } } @prods_info;
4173                  push (@{$node_type->{'scenario present all'}}, $scen_name);                      if (@main_prods == 0) {
4174                            # no main info
4175                            map { push @main_prods, $_->[0]; $compound_ids{$_->[0]} = 1 } @prods_info;
4176                        }
4177    
4178                        push @reaction_info, @main_subs, $rxn_id, @main_prods;
4179                        push @rxns_and_cpds, (join " ", @reaction_info);
4180                        $reaction_presence->{$rxn_id} = exists $org_ss_reactions->{$rxn_id} ? "reaction present" : "reaction absent";
4181              }              }
4182    
4183                    push @scenario_path_info, (join ";", @rxns_and_cpds);
4184                }
4185    
4186                $path_info->{$ss_name."/".$scen_name} = join "&", @scenario_path_info;
4187          }          }
4188      }      }
4189    
# Line 4170  Line 4192 
4192          # compounds present with All scenarios that aren't biomass          # compounds present with All scenarios that aren't biomass
4193          if (! exists($old_current_inputs{$cpd}) && ! exists ($biomass_reached{$cpd}))          if (! exists($old_current_inputs{$cpd}) && ! exists ($biomass_reached{$cpd}))
4194          {          {
4195              push (@{$node_type->{'compound present all'}}, $cpd."/".$compound_id_to_name{$cpd});              push (@{$node_type->{'compound present all'}}, $cpd);
4196                $compound_ids{$cpd} = 1;
4197          }          }
4198      }      }
4199    
4200      print STDERR "\tStage 3 - Complete!\n" if($debug);      print STDERR "\tStage 3 - Complete!\n" if($debug);
4201    
4202      return ($node_type, $node_connections);      return ($node_type, $node_connections, \%compound_ids, $path_info, $reaction_presence);
4203    
4204      sub generate_new_inputs      sub generate_new_inputs
4205      {      {
# Line 4193  Line 4216 
4216                  next;                  next;
4217              }              }
4218    
4219              my $sids = join ",", sort @{$scenario->get_substrates};              my $sids = join ",", sort keys %{$scenario->get_substrates};
4220              my $pids = join ",", sort @{$scenario->get_products};              my $pids = join ",", sort keys %{$scenario->get_products};
4221              my $cid_scenario_name = "[".$sids."]".$scenario->get_scenario_name()."[".$pids."]";              my $cid_scenario_name = "[".$sids."]".$scenario->get_scenario_name()."[".$pids."]";
4222    
4223              my @matched;              my @matched;
4224    
4225              foreach my $substrate (@{$scenario->get_substrates()})              foreach my $substrate (keys %{$scenario->get_substrates()})
4226              {              {
4227                  if(exists $inputs->{$substrate})                  if(exists $inputs->{$substrate})
4228                  {                  {
# Line 4207  Line 4230 
4230                  }                  }
4231              }              }
4232    
4233              if (scalar @matched != scalar(@{$scenario->get_substrates()}))              if (scalar @matched != scalar(keys %{$scenario->get_substrates()}))
4234              {              {
4235                    print STDERR "($iter_num) Scenario ".$cid_scenario_name." did not run.  matched: @matched\n" if($debug);
4236                  next;                  next;
4237              }              }
4238              else              else
4239              {              {
4240                  print STDERR "Scenario ".$cid_scenario_name." ran.\n" if($debug);                  print STDERR "($iter_num) Scenario ".$cid_scenario_name." ran. matched: @matched\n" if($debug);
4241                  $break = 0;                  $break = 0;
4242                  $scenarios_id{$scenario->get_id()} = 1;                  $scenarios_id{$scenario->get_id()} = 1;
4243                  $scenarios_used{$scenario->get_scenario_name()} = 1;                  $scenarios_used{$scenario->get_scenario_name()} = 1;
# Line 4223  Line 4247 
4247                  next if (exists $used_scenario_paths_this_time{$cid_scenario_name});                  next if (exists $used_scenario_paths_this_time{$cid_scenario_name});
4248                  $used_scenario_paths_this_time{$cid_scenario_name} = 1;                  $used_scenario_paths_this_time{$cid_scenario_name} = 1;
4249    
4250                  foreach my $product (@{$scenario->get_products()})                  foreach my $product (keys %{$scenario->get_products()})
4251                  {                  {
4252                      $new_inputs{$product}->{$cid_scenario_name} = 1;                      $new_inputs{$product}->{$cid_scenario_name} = 1;
4253                  }                  }
# Line 4288  Line 4312 
4312          map { $genome_reactions{$_} = 1 } @included_reactions;          map { $genome_reactions{$_} = 1 } @included_reactions;
4313      }      }
4314    
4315      return %genome_reactions;      return \%genome_reactions;
4316    }
4317    
4318    sub get_reactions_for_genome_in_subsystem_fast
4319    {
4320            my ($genome_id, $ss_name) = @_;
4321            my $ss_dir = Subsystem::get_dir_from_name($ss_name);
4322            my (%hope_reactions, %genome_reactions);
4323    
4324            # this is taken directly from FIGRules::GetHopeReactions with a few modifications
4325            # this is called when a subsystem object is constructed normally,
4326            # but here all we need are the reactions so this is quite a bit faster.
4327            # returns the same reactions as get_reactions_for_genome_in_subsystem
4328            my $hopeFileName = "$ss_dir/hope_reactions";
4329            if (-f $hopeFileName)
4330            {
4331                    # Open the hope reaction file.
4332                    open(RXNFILE, "<$hopeFileName");
4333                    # Loop through the hope reaction file.
4334                    while (<RXNFILE>)
4335                    {
4336                            # Parse out the role and the reaction list.
4337                            if ($_ =~ /^(\S.*\S)\t(\S+)/)
4338                            {
4339                                    my ($role, $reactions) = ($1, $2);
4340                                    # FIGRules checks to see if the role is in this subsystem but we can't. Watch for this!
4341                                    push( @{$hope_reactions{$role}}, split(/,\s*/,$reactions) );
4342                            }
4343                    }
4344                    close(RXNFILE);
4345            }
4346    
4347            foreach my $role (keys %hope_reactions)
4348        {
4349            my @peg_list = $fig->seqs_with_role($role, "", $genome_id);
4350    
4351            if (scalar @peg_list > 0)
4352            {
4353                foreach my $reaction (@{$hope_reactions{$role}})
4354                {
4355                    push @{$genome_reactions{$reaction}}, @peg_list;
4356                }
4357            }
4358        }
4359    
4360            return \%genome_reactions;
4361  }  }
4362    
4363  sub get_scenario_directory  sub get_scenario_directory
# Line 4307  Line 4376 
4376  sub check_hope_info  sub check_hope_info
4377  {  {
4378      my ($ss) = @_;      my ($ss) = @_;
4379      my (%invalid_roles, %invalid_reactions, $invalid_classification);      my (%invalid_roles, %invalid_reactions, $invalid_classification, %invalid_kegg_maps, %reaction_ec_mismatch);
4380      my $sub = $fig->get_subsystem($ss);      my $sub = $fig->get_subsystem($ss);
4381    
4382      if (! defined $sub)      if (! defined $sub)
4383      {      {
4384          return ( {}, {}, "Subsystem no longer exists" );          return ( {}, {}, "Subsystem no longer exists", {}, {} );
4385      }      }
4386    
4387      my $classification = $sub->get_classification->[0];      my $classification = $sub->get_classification->[0];
# Line 4338  Line 4407 
4407      {      {
4408          my @additional_reactions = $sub->get_hope_additional_reactions($name);          my @additional_reactions = $sub->get_hope_additional_reactions($name);
4409          map { push @{$reactions{$_}}, "ADDITIONAL_REACTIONS" } @additional_reactions;          map { push @{$reactions{$_}}, "ADDITIONAL_REACTIONS" } @additional_reactions;
4410    
4411            foreach my $id ($sub->get_hope_map_ids($name))
4412            {
4413                my $map_name = $fig->map_name("map".$id);
4414                $invalid_kegg_maps{$id} = 1 if $map_name eq "";
4415            }
4416      }      }
4417    
4418    
# Line 4347  Line 4422 
4422          {          {
4423              $invalid_reactions{$reaction} = $reactions{$reaction};              $invalid_reactions{$reaction} = $reactions{$reaction};
4424          }          }
4425            else
4426            {
4427                my %kegg_ecs;
4428                @kegg_ecs{$fig->catalyzed_by($reaction)} = ();
4429                foreach my $role ( @{$reactions{$reaction}})
4430                {
4431                    if ($role =~ /(\d+\.\d+\.\d+\.\d+)/ && ! exists($kegg_ecs{$1})) {
4432                        push @{$reaction_ec_mismatch{$reaction}}, $role;
4433                    }
4434                }
4435            }
4436      }      }
4437    
4438      return (\%invalid_roles, \%invalid_reactions, $invalid_classification);      return (\%invalid_roles, \%invalid_reactions, $invalid_classification, \%invalid_kegg_maps, \%reaction_ec_mismatch);
4439    
4440  }  }
4441    

Legend:
Removed from v.1.32  
changed lines
  Added in v.1.33

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3