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

Annotation of /FigKernelScripts/check_model.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (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 : formsma 1.4 my %scenarios_id;
36 : formsma 1.1 my %all_substrates;
37 :     foreach (@scenarios)
38 :     {
39 :     map {$all_substrates{$_} = 1} @{$_->get_substrates_all()};
40 :     $scenarios_used{$_->get_scenario_name()} = 0;
41 : formsma 1.4 $scenarios_id{$_->get_id()} = 0;
42 : formsma 1.1 }
43 :    
44 :     my %input_cpds;
45 :     open(INPUTS,$filebase."inputs.txt")
46 :     or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs");
47 :     while(<INPUTS>)
48 :     {
49 :     chomp;
50 :     my @line = split "\t", $_;
51 :     chomp @line;
52 :     $input_cpds{$line[0]} = 1;
53 :     }
54 :     close(INPUTS);
55 :    
56 : formsma 1.4 #ATP Snythase HACK, lets add ATP as a input
57 : mkubal 1.5
58 :     #cofactors just used for testing connectivity, not actually provided to model
59 : formsma 1.4 $input_cpds{"C00002"} = 1;
60 : mkubal 1.5 $input_cpds{"C00028"} = 1;
61 :     $input_cpds{"C00030"} = 1;
62 :     $input_cpds{"C00035"} = 1;
63 :     $input_cpds{"C00040"} = 1;
64 :     $input_cpds{"C00044"} = 1;
65 :     $input_cpds{"C00126"} = 1;
66 :     $input_cpds{"C00138"} = 1;
67 :     $input_cpds{"C00139"} = 1;
68 :     $input_cpds{"C00239"} = 1;
69 :     $input_cpds{"C00342"} = 1;
70 :     $input_cpds{"C00343"} = 1;
71 :     $input_cpds{"C00361"} = 1;
72 :     $input_cpds{"C01070"} = 1;
73 :    
74 :    
75 :     #inputs not produced or imported but required for connecting all scenarios
76 :     =pod
77 :     $input_cpds{"C00141"} = 1;
78 :     $input_cpds{"C00099"} = 1;
79 :     $input_cpds{"C00049"} = 1;
80 :     $input_cpds{"C00040"} = 1;
81 :     $input_cpds{"C00054"} = 1;
82 :     $input_cpds{"C00059"} = 1;
83 :     $input_cpds{"C00116"} = 1;
84 :     $input_cpds{"C00163"} = 1;
85 :     $input_cpds{"C00185"} = 1;
86 :     $input_cpds{"C00206"} = 1;
87 :     $input_cpds{"C00208"} = 1;
88 :     $input_cpds{"C00234"} = 1;
89 :     $input_cpds{"C00235"} = 1;
90 :     $input_cpds{"C00239"} = 1;
91 :     $input_cpds{"C00250"} = 1;
92 :     $input_cpds{"C00266"} = 1;
93 :     $input_cpds{"C00281"} = 1;
94 :     $input_cpds{"C00288"} = 1;
95 :     $input_cpds{"C00301"} = 1;
96 :     $input_cpds{"C00314"} = 1;
97 :     $input_cpds{"C00416"} = 1;
98 :     $input_cpds{"C00534"} = 1;
99 :     $input_cpds{"C00670"} = 1;
100 :     $input_cpds{"C00672"} = 1;
101 :     $input_cpds{"C00794"} = 1;
102 :     $input_cpds{"C00898"} = 1;
103 :     $input_cpds{"C00988"} = 1;
104 :     $input_cpds{"C01070"} = 1;
105 :     $input_cpds{"C01096"} = 1;
106 :     $input_cpds{"C01137"} = 1;
107 :     $input_cpds{"C01157"} = 1;
108 :     $input_cpds{"C01233"} = 1;
109 :     $input_cpds{"C02191"} = 1;
110 :     $input_cpds{"C02232"} = 1;
111 :     $input_cpds{"C02315"} = 1;
112 :     $input_cpds{"C02336"} = 1;
113 :     $input_cpds{"C02350"} = 1;
114 :     $input_cpds{"C03150"} = 1;
115 :     $input_cpds{"C03543"} = 1;
116 :     $input_cpds{"C04121"} = 1;
117 :     $input_cpds{"C04146"} = 1;
118 :     $input_cpds{"C04489"} = 1;
119 :     $input_cpds{"C04688"} = 1;
120 :     $input_cpds{"C04824"} = 1;
121 :     $input_cpds{"C05223"} = 1;
122 :     $input_cpds{"C05761"} = 1;
123 :     $input_cpds{"C05898"} = 1;
124 :     $input_cpds{"C06143"} = 1;
125 :     $input_cpds{"C06187"} = 1;
126 :     $input_cpds{"C06188"} = 1;
127 :     $input_cpds{"C11457"} = 1;
128 :     $input_cpds{"C14818"} = 1;
129 :    
130 :    
131 :     $input_cpds{"C00002"} = 1;
132 :     $input_cpds{"C00001"} = 1;
133 :     $input_cpds{"C00003"} = 1;
134 :     $input_cpds{"C00007"} = 1;
135 :     $input_cpds{"C00009"} = 1;
136 :     $input_cpds{"C00023"} = 1;
137 :     $input_cpds{"C00024"} = 1;
138 :     $input_cpds{"C00034"} = 1;
139 :     $input_cpds{"C00035"} = 1;
140 :     $input_cpds{"C00044"} = 1;
141 :     $input_cpds{"C00076"} = 1;
142 :     $input_cpds{"C00080"} = 1;
143 :     $input_cpds{"C00082"} = 1;
144 :     $input_cpds{"C00087"} = 1;
145 :     $input_cpds{"C00114"} = 1;
146 :     $input_cpds{"C00115"} = 1;
147 :     $input_cpds{"C00120"} = 1;
148 :     $input_cpds{"C00126"} = 1;
149 :     $input_cpds{"C00138"} = 1;
150 :     $input_cpds{"C00139"} = 1;
151 :     $input_cpds{"C00153"} = 1;
152 :     $input_cpds{"C00229"} = 1;
153 :     $input_cpds{"C00238"} = 1;
154 :     $input_cpds{"C00255"} = 1;
155 :     $input_cpds{"C00283"} = 1;
156 :     $input_cpds{"C00305"} = 1;
157 :     $input_cpds{"C00342"} = 1;
158 :     $input_cpds{"C00343"} = 1;
159 :     $input_cpds{"C00361"} = 1;
160 :     $input_cpds{"C00378"} = 1;
161 :     $input_cpds{"C00504"} = 1;
162 :     $input_cpds{"C00568"} = 1;
163 :     $input_cpds{"C00864"} = 1;
164 :     $input_cpds{"C05776"} = 1;
165 :    
166 :     $input_cpds{"C00141"} = 1;
167 :     $input_cpds{"C00099"} = 1;
168 :     $input_cpds{"C00049"} = 1;
169 :     $input_cpds{"C00040"} = 1;
170 :     $input_cpds{"C00054"} = 1;
171 :     $input_cpds{"C00059"} = 1;
172 :     $input_cpds{"C00116"} = 1;
173 :     $input_cpds{"C00163"} = 1;
174 :     $input_cpds{"C00185"} = 1;
175 :     $input_cpds{"C00206"} = 1;
176 :     $input_cpds{"C00208"} = 1;
177 :     $input_cpds{"C00234"} = 1;
178 :     $input_cpds{"C00235"} = 1;
179 :     $input_cpds{"C00239"} = 1;
180 :     $input_cpds{"C00250"} = 1;
181 :     $input_cpds{"C00266"} = 1;
182 :     $input_cpds{"C00281"} = 1;
183 :     $input_cpds{"C00288"} = 1;
184 :     $input_cpds{"C00301"} = 1;
185 :     $input_cpds{"C00314"} = 1;
186 :     $input_cpds{"C00416"} = 1;
187 :     $input_cpds{"C00534"} = 1;
188 :     $input_cpds{"C00670"} = 1;
189 :     $input_cpds{"C00672"} = 1;
190 :     $input_cpds{"C00794"} = 1;
191 :     $input_cpds{"C00898"} = 1;
192 :     $input_cpds{"C00988"} = 1;
193 :     $input_cpds{"C01070"} = 1;
194 :     $input_cpds{"C01096"} = 1;
195 :     $input_cpds{"C01137"} = 1;
196 :     $input_cpds{"C01157"} = 1;
197 :     $input_cpds{"C01233"} = 1;
198 :     $input_cpds{"C02191"} = 1;
199 :     $input_cpds{"C02232"} = 1;
200 :     $input_cpds{"C02315"} = 1;
201 :     $input_cpds{"C02336"} = 1;
202 :     $input_cpds{"C02350"} = 1;
203 :     $input_cpds{"C03150"} = 1;
204 :     $input_cpds{"C03543"} = 1;
205 :     $input_cpds{"C04121"} = 1;
206 :     $input_cpds{"C04146"} = 1;
207 :     $input_cpds{"C04489"} = 1;
208 :     $input_cpds{"C04688"} = 1;
209 :     $input_cpds{"C04824"} = 1;
210 :     $input_cpds{"C05223"} = 1;
211 :     $input_cpds{"C05761"} = 1;
212 :     $input_cpds{"C05898"} = 1;
213 :     $input_cpds{"C06143"} = 1;
214 :     $input_cpds{"C06187"} = 1;
215 :     $input_cpds{"C06188"} = 1;
216 :     $input_cpds{"C11457"} = 1;
217 :     $input_cpds{"C14818"} = 1;
218 :    
219 :     =cut
220 : formsma 1.4
221 : formsma 1.1 my %cofactor_cpds;
222 : formsma 1.4 open(INPUTCO,$filebase."secondary_compounds.txt")
223 :     or die("Failed to open ".$filebase."secondary_compounds.txt. Run find_model_cofactors $genome_id\n");
224 : formsma 1.1 while(<INPUTCO>)
225 :     {
226 :     chomp;
227 :     my @line = split "\t", $_;
228 :     chomp @line;
229 :     $cofactor_cpds{$line[0]} = 1;
230 :     }
231 :     close(INPUTCO);
232 :    
233 :     my %biomass_cpds;
234 :     open(BIOMASS,$filebase."biomass.txt")
235 :     or die ("Failed to open ".$filebase."biomass.txt Run find_model_information $genome_id biomass");
236 :     while(<BIOMASS>)
237 :     {
238 :     chomp;
239 :     my @temp = split "\t", $_;
240 :     chomp @temp;
241 :     $biomass_cpds{$temp[0]} = 1;
242 :     }
243 :     close(BIOMASS);
244 : formsma 1.4
245 :     my %compound_id_to_name;
246 :     open(IN,"$FIG_Config::global/Models/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt");
247 :     while($_ = <IN>){
248 :     chomp($_);
249 :     my($id,$name) = split("\t",$_);
250 :     $compound_id_to_name{$id} = $name;
251 :     }
252 :     close(IN);
253 :    
254 : formsma 1.1 print STDERR "\tStage 1 - Complete\n" if($debug);
255 :    
256 :     print STDERR "\tStage 2 - Run Analysis Tools (fueled by caffine)\n" if($debug);
257 :     #Total list of what compounds we reach.
258 :     my @list_generated = keys %input_cpds;
259 : formsma 1.4 push @list_generated , keys %cofactor_cpds;
260 : formsma 1.1 my %possible_problems;
261 :     my %good_compounds;
262 :    
263 :     # Here we are running all the scenarios we have inputs for
264 :     # to analyze what biomass we can reach. We are also identifying
265 :     # any Scenarios that are not being utilized.
266 :    
267 :     my %current_inputs;
268 :     foreach (keys %input_cpds)
269 :     {
270 :     $current_inputs{$_} = 1;
271 :     }
272 :     my %current_cofactors;
273 :     foreach (keys %cofactor_cpds)
274 :     {
275 :     $current_cofactors{$_} = 1;
276 :     }
277 :    
278 :     while(scalar(keys %current_inputs))
279 :     {
280 : formsma 1.4 my ($temp_in,$temp_co,$break) = generate_new_inputs(\%current_inputs,\%current_cofactors);
281 : formsma 1.1
282 : formsma 1.4 map {$current_inputs{$_} = 1;} keys %$temp_in;
283 :     map {$current_cofactors{$_} = 1;} keys %$temp_co;
284 : formsma 1.1
285 :     push @list_generated, keys %$temp_in;
286 :     push @list_generated, keys %$temp_co;
287 : formsma 1.4
288 :     print STDERR "Rechecking with inputs\n" if($debug);
289 :     if($break)
290 :     {
291 :     last;
292 :     }
293 : formsma 1.1 }
294 : formsma 1.4 print STDERR "Finished checking scenarios\n" if($debug);
295 : formsma 1.1 foreach my $found (@list_generated)
296 :     {
297 :     if(defined $biomass_cpds{$found})
298 :     {
299 :     $biomass_cpds{$found} = 2;
300 :     }
301 :     }
302 :    
303 :     # Analyze input compounds
304 :     # and also identify dead-end Scenarios
305 :     my %dead_end_scenarios;
306 :    
307 :     foreach my $scenario (@scenarios)
308 :     {
309 :     my $count = 0;
310 :     my %used_inputs;
311 :     foreach my $substrate (@{$scenario->get_substrates_all()})
312 :     {
313 :     if($input_cpds{$substrate} == 1)
314 :     {
315 :     $used_inputs{$substrate} = 1;
316 :     $count++;
317 :     }
318 :     }
319 :    
320 :     if($count !=0 && $count != scalar(@{$scenario->get_substrates_all()}))
321 :     {
322 :     map {$possible_problems{$_} = 1;} keys %used_inputs;
323 :     }
324 :     else
325 :     {
326 :     map {$good_compounds{$_} = 0;} keys %used_inputs;
327 :     }
328 :    
329 :     #dead-end analysis
330 :     my $bad_end = 1;
331 :     foreach my $product (@{$scenario->get_products_all()})
332 :     {
333 :     if(defined $biomass_cpds{$product} || $all_substrates{$product} == 1)
334 :     {
335 :     $bad_end = 0;
336 :     last;
337 :     }
338 :     }
339 :     #Iff we reach this point, then this scenario produces no biomass
340 :     # or compounds used by another scenario. It's all waste, could be a error.
341 :     if($bad_end)
342 :     {
343 :     $dead_end_scenarios{$scenario->get_scenario_name()} = 1;
344 :     }
345 :     }
346 :    
347 :     open(RESULTS,">".$filebase."info_biomass.txt");
348 :     open(OUT, ">".$filebase."biomass_valid.txt");
349 :     foreach (sort keys %biomass_cpds)
350 :     {
351 :     if($biomass_cpds{$_} != 2)
352 :     {
353 : formsma 1.4 print RESULTS "Biomass $_ $compound_id_to_name{$_} failed to be generated using inputs.txt\n";
354 : formsma 1.1 }
355 :     else
356 :     {
357 : formsma 1.4 print OUT "$_\t$compound_id_to_name{$_}\n";
358 : formsma 1.1 }
359 :     }
360 :     close(RESULTS);
361 :     close(OUT);
362 :    
363 :     open(OUT,">$filebase\scenarios_valid.txt");
364 :     open(SCEN_RESULTS,">".$filebase."info_scenarios.txt");
365 :     foreach my $scen_id (sort keys %scenarios_used)
366 :     {
367 :     if($scenarios_used{$scen_id} == 0)
368 :     {
369 :     print SCEN_RESULTS "Scenario $scen_id unable to run with current inputs\n";
370 :     }
371 :     else
372 :     {
373 :     print OUT "$scen_id\n";
374 :     }
375 :     }
376 :     print SCEN_RESULTS "\n\n#################\n\n\n";
377 :     foreach (sort keys %dead_end_scenarios)
378 :     {
379 :     print SCEN_RESULTS "Scenario $_ is possibly a dead end(no outputs used by other scenarios or biomass.)\n";
380 :     }
381 :     close(SCEN_RESULTS);
382 :     close(OUT);
383 :    
384 : formsma 1.4 #open(RESULTS_INPUT,">".$filebase."info_inputs.txt");
385 :     #foreach (sort keys %possible_problems)
386 :     #{
387 :     # if($possible_problems{$_} == 1 && !defined $good_compounds{$_})
388 :     # {
389 :     # print RESULTS_INPUT "Input $_ may be a false-positive. Perhaps a Scenario should generate this.\n";
390 :     # }
391 :     #}
392 :     #close(RESULTS_INPUT);
393 : formsma 1.1
394 :     print STDERR "\tStage 2 - Complete! Results in info_inputs.txt, info_biomass.txt and info_scenarios.txt\n" if($debug);
395 :    
396 :     sub generate_new_inputs
397 :     {
398 :     my ($inputs,$cofactors) = @_;
399 :     my %new_inputs;
400 :     my %new_cofactors;
401 : formsma 1.4 my $break = 1;
402 : formsma 1.1 foreach my $scenario (@scenarios)
403 :     {
404 : formsma 1.4 if($scenarios_id{$scenario->get_id()} == 1)
405 :     {
406 :     next;
407 :     }
408 : formsma 1.1 my $count = 0;
409 :     foreach my $substrate (@{$scenario->get_substrates_all()})
410 :     {
411 :     if($inputs->{$substrate})
412 :     {
413 :     $count++;
414 :     }
415 : formsma 1.4 elsif($cofactors->{$substrate})
416 : formsma 1.1 {
417 :     $count++;
418 :     }
419 :     }
420 :    
421 :     if($count != scalar(@{$scenario->get_substrates_all()}))
422 :     {
423 :     next;
424 :     }
425 :     else
426 :     {
427 :     print STDERR "Scenario ".$scenario->get_scenario_name()." ran.\n" if($debug);
428 : formsma 1.4 $break = 0;
429 :     $scenarios_id{$scenario->get_id()} = 1;
430 : formsma 1.1 $scenarios_used{$scenario->get_scenario_name()} = 1;
431 :     foreach my $product (@{$scenario->get_products()})
432 :     {
433 :     $new_inputs{$product} = 1;
434 :     }
435 :     foreach my $product (@{$scenario->get_product_cofactors()})
436 :     {
437 :     $new_cofactors{$product} = 1;
438 :     }
439 :     }
440 :     }
441 : formsma 1.4 return (\%new_inputs,\%new_cofactors,$break);
442 : formsma 1.1 }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3