[Bio] / FigKernelScripts / check_model.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/check_model.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (view) (download) (as text)

1 : formsma 1.1 # _*_ Perl _*_
2 :     #
3 :     #
4 :     #
5 :     use strict;
6 : formsma 1.2 use FIGV;
7 : formsma 1.1 use Scenario;
8 :    
9 : formsma 1.2 my ($genome_id,$orgdir) = @ARGV;
10 : formsma 1.1
11 :     my $debug = 0;
12 :    
13 : formsma 1.2 unless (scalar (@ARGV) >= 1)
14 : formsma 1.1 {
15 :     print STDERR "usage: check_model <genome_id>\n";
16 :     exit(-1);
17 :     }
18 :     chomp $genome_id;
19 : formsma 1.2 my $fig;
20 :     if($orgdir)
21 :     {
22 :     $fig = new FIGV($orgdir);
23 :     }
24 :     else
25 :     {
26 :     $fig = new FIGV;
27 :     }
28 : formsma 1.3 Scenario::set_fig($fig,$genome_id);
29 : formsma 1.1 print STDERR "Running check_model on genome $genome_id\n" if($debug);
30 :     print STDERR "\tStage 1 - Loading Scenarios, model inputs and biomass\n" if($debug);
31 :     my @scenarios = @{Scenario->get_genome_scenarios($genome_id,0)};
32 : formsma 1.2 my $filebase = $fig->model_directory($genome_id)."/Analysis/";
33 : formsma 1.1
34 :     my %scenarios_used;
35 :     my %all_substrates;
36 :     foreach (@scenarios)
37 :     {
38 :     map {$all_substrates{$_} = 1} @{$_->get_substrates_all()};
39 :     $scenarios_used{$_->get_scenario_name()} = 0;
40 :     }
41 :    
42 :     my %input_cpds;
43 :     open(INPUTS,$filebase."inputs.txt")
44 :     or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs");
45 :     while(<INPUTS>)
46 :     {
47 :     chomp;
48 :     my @line = split "\t", $_;
49 :     chomp @line;
50 :     $input_cpds{$line[0]} = 1;
51 :     }
52 :     close(INPUTS);
53 :    
54 :     my %cofactor_cpds;
55 :     open(INPUTCO,$filebase."cofactors.txt")
56 :     or die("Failed to open ".$filebase."cofactors.txt. Run find_model_cofactors $genome_id\n");
57 :     while(<INPUTCO>)
58 :     {
59 :     chomp;
60 :     my @line = split "\t", $_;
61 :     chomp @line;
62 :     $cofactor_cpds{$line[0]} = 1;
63 :     }
64 :     close(INPUTCO);
65 :    
66 :     my %biomass_cpds;
67 :     open(BIOMASS,$filebase."biomass.txt")
68 :     or die ("Failed to open ".$filebase."biomass.txt Run find_model_information $genome_id biomass");
69 :     while(<BIOMASS>)
70 :     {
71 :     chomp;
72 :     my @temp = split "\t", $_;
73 :     chomp @temp;
74 :     $biomass_cpds{$temp[0]} = 1;
75 :     }
76 :     close(BIOMASS);
77 :     print STDERR "\tStage 1 - Complete\n" if($debug);
78 :    
79 :     print STDERR "\tStage 2 - Run Analysis Tools (fueled by caffine)\n" if($debug);
80 :     #Total list of what compounds we reach.
81 :     my @list_generated = keys %input_cpds;
82 :    
83 :     my %possible_problems;
84 :     my %good_compounds;
85 :    
86 :     # Here we are running all the scenarios we have inputs for
87 :     # to analyze what biomass we can reach. We are also identifying
88 :     # any Scenarios that are not being utilized.
89 :    
90 :     my %current_inputs;
91 :     foreach (keys %input_cpds)
92 :     {
93 :     $current_inputs{$_} = 1;
94 :     }
95 :     my %current_cofactors;
96 :     my %previous_inputs;
97 :     foreach (keys %cofactor_cpds)
98 :     {
99 :     $current_cofactors{$_} = 1;
100 :     }
101 :    
102 :     while(scalar(keys %current_inputs))
103 :     {
104 :     my ($temp_in,$temp_co) = generate_new_inputs(\%current_inputs,\%current_cofactors);
105 :     #Check to make sure we have a different compound
106 :     my $break = 1;
107 :     foreach my $key (keys %$temp_in)
108 :     {
109 :     if(!defined $current_inputs{$key} && !defined $previous_inputs{$key})
110 :     {
111 :     $break = 0;
112 :     }
113 :     }
114 :     if($break) { last;} #break the while loop, we are looping
115 :    
116 :     undef %previous_inputs;
117 :     map {$previous_inputs{$_} = 1;} sort keys %current_inputs;
118 :    
119 :     #undef %current_inputs;
120 :     map {$current_inputs{$_} = 1;} sort keys %$temp_in;
121 :    
122 :     #undef %current_cofactors;
123 :     map {$current_cofactors{$_} = 1;} sort keys %$temp_co;
124 :    
125 :     push @list_generated, keys %$temp_in;
126 :     push @list_generated, keys %$temp_co;
127 :    
128 :     }
129 :    
130 :     foreach my $found (@list_generated)
131 :     {
132 :     if(defined $biomass_cpds{$found})
133 :     {
134 :     $biomass_cpds{$found} = 2;
135 :     }
136 :     }
137 :    
138 :     # Analyze input compounds
139 :     # and also identify dead-end Scenarios
140 :     my %dead_end_scenarios;
141 :    
142 :     foreach my $scenario (@scenarios)
143 :     {
144 :     my $count = 0;
145 :     my %used_inputs;
146 :     foreach my $substrate (@{$scenario->get_substrates_all()})
147 :     {
148 :     if($input_cpds{$substrate} == 1)
149 :     {
150 :     $used_inputs{$substrate} = 1;
151 :     $count++;
152 :     }
153 :     }
154 :    
155 :     if($count !=0 && $count != scalar(@{$scenario->get_substrates_all()}))
156 :     {
157 :     map {$possible_problems{$_} = 1;} keys %used_inputs;
158 :     }
159 :     else
160 :     {
161 :     map {$good_compounds{$_} = 0;} keys %used_inputs;
162 :     }
163 :    
164 :     #dead-end analysis
165 :     my $bad_end = 1;
166 :     foreach my $product (@{$scenario->get_products_all()})
167 :     {
168 :     if(defined $biomass_cpds{$product} || $all_substrates{$product} == 1)
169 :     {
170 :     $bad_end = 0;
171 :     last;
172 :     }
173 :     }
174 :     #Iff we reach this point, then this scenario produces no biomass
175 :     # or compounds used by another scenario. It's all waste, could be a error.
176 :     if($bad_end)
177 :     {
178 :     $dead_end_scenarios{$scenario->get_scenario_name()} = 1;
179 :     }
180 :     }
181 :    
182 :     open(RESULTS,">".$filebase."info_biomass.txt");
183 :     open(OUT, ">".$filebase."biomass_valid.txt");
184 :     foreach (sort keys %biomass_cpds)
185 :     {
186 :     if($biomass_cpds{$_} != 2)
187 :     {
188 :     print RESULTS "Biomass $_ failed to be generated using inputs.txt\n";
189 :     }
190 :     else
191 :     {
192 :     print OUT "$_\n";
193 :     }
194 :     }
195 :     close(RESULTS);
196 :     close(OUT);
197 :    
198 :     open(OUT,">$filebase\scenarios_valid.txt");
199 :     open(SCEN_RESULTS,">".$filebase."info_scenarios.txt");
200 :     foreach my $scen_id (sort keys %scenarios_used)
201 :     {
202 :     if($scenarios_used{$scen_id} == 0)
203 :     {
204 :     print SCEN_RESULTS "Scenario $scen_id unable to run with current inputs\n";
205 :     }
206 :     else
207 :     {
208 :     print OUT "$scen_id\n";
209 :     }
210 :     }
211 :     print SCEN_RESULTS "\n\n#################\n\n\n";
212 :     foreach (sort keys %dead_end_scenarios)
213 :     {
214 :     print SCEN_RESULTS "Scenario $_ is possibly a dead end(no outputs used by other scenarios or biomass.)\n";
215 :     }
216 :     close(SCEN_RESULTS);
217 :     close(OUT);
218 :    
219 :     open(RESULTS_INPUT,">".$filebase."info_inputs.txt");
220 :     foreach (sort keys %possible_problems)
221 :     {
222 :     if($possible_problems{$_} == 1 && !defined $good_compounds{$_})
223 :     {
224 :     print RESULTS_INPUT "Input $_ may be a false-positive. Perhaps a Scenario should generate this.\n";
225 :     }
226 :     }
227 :     close(RESULTS_INPUT);
228 :    
229 :     print STDERR "\tStage 2 - Complete! Results in info_inputs.txt, info_biomass.txt and info_scenarios.txt\n" if($debug);
230 :    
231 :     sub generate_new_inputs
232 :     {
233 :     my ($inputs,$cofactors) = @_;
234 :     my %new_inputs;
235 :     my %new_cofactors;
236 :     foreach my $scenario (@scenarios)
237 :     {
238 :     my $count = 0;
239 :     foreach my $substrate (@{$scenario->get_substrates_all()})
240 :     {
241 :     if($inputs->{$substrate})
242 :     {
243 :     $count++;
244 :     }
245 :     elsif( $cofactors->{$substrate})
246 :     {
247 :     $count++;
248 :     }
249 :     }
250 :    
251 :     if($count != scalar(@{$scenario->get_substrates_all()}))
252 :     {
253 :     next;
254 :     }
255 :     else
256 :     {
257 :     print STDERR "Scenario ".$scenario->get_scenario_name()." ran.\n" if($debug);
258 :    
259 :     ### here we want to mark what inputs we used
260 :     foreach (@{$scenario->get_substrates()})
261 :     {
262 :     $inputs->{$_} = 2;
263 :     }
264 :    
265 :     $scenarios_used{$scenario->get_scenario_name()} = 1;
266 :     foreach my $product (@{$scenario->get_products()})
267 :     {
268 :     $new_inputs{$product} = 1;
269 :     }
270 :     foreach my $product (@{$scenario->get_product_cofactors()})
271 :     {
272 :     $new_cofactors{$product} = 1;
273 :     }
274 :     }
275 :     }
276 :     #Add any inputs that aren't used to the return again
277 :     foreach (keys %$inputs)
278 :     {
279 :     if($inputs->{$_} !=2)
280 :     {
281 :     $new_inputs{$_} = 1;
282 :     }
283 :     }
284 :    
285 :     return (\%new_inputs,\%new_cofactors);
286 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3