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

Annotation of /FigKernelScripts/make_model_spreadsheets.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : mkubal 1.1 # _*_ Perl _*_
2 :     #
3 :     #
4 :     #
5 :     use strict;
6 :     use FIGV;
7 :     use Scenario;
8 :    
9 :     my ($genome_id,$orgdir) = @ARGV;
10 :    
11 :     my $debug = 0;
12 :    
13 :    
14 :     unless (scalar (@ARGV) >= 1)
15 :     {
16 :     print STDERR "usage: make_model_spreadsheets <genome_id>\n";
17 :     exit(-1);
18 :     }
19 :     chomp $genome_id;
20 :     my $fig;
21 :     if($orgdir)
22 :     {
23 :     $fig = new FIGV($orgdir);
24 :     }
25 :     else
26 :     {
27 :     $fig = new FIGV;
28 :     }
29 :     Scenario::set_fig($fig,$genome_id);
30 :     print STDERR "Creating model spreadsheets forgenome $genome_id\n" if($debug);
31 :     print STDERR "\tStage 1 - Loading Scenarios, model inputs and biomass\n" if($debug);
32 :     my @scenarios = @{Scenario->get_genome_scenarios($genome_id,1)};
33 : dejongh 1.2 my $filebase = $fig->scenario_directory($genome_id)."/Analysis/";
34 : mkubal 1.1
35 :     #library of transported compounds in other models
36 :     open(LIB,$FIG_Config::global."/Models/transports_in.txt.original");
37 :     my %lib_input_cpds;
38 :     while(<LIB>)
39 :     {
40 :     chomp;
41 :     my @line = split "\t", $_;
42 :     chomp @line;
43 :     $lib_input_cpds{$line[0]} = 1;
44 :     }
45 :     close(LIB);
46 :    
47 :     my %evidence_of_transport_peg;
48 :     my %evidence_of_transport_annotation;
49 :     open(EVIDENCE,$FIG_Config::global."/Models/transport_evidence.txt");
50 :     #open(EVIDENCE,$FIG_Config::global."/Models/transport_evidence.txt");
51 :     while(<EVIDENCE>)
52 :     {
53 :     chomp;
54 :     my ($cid,$annotation) = split "\t", $_;
55 :     if($evidence_of_transport_annotation{$cid}){
56 :     my $ref = $evidence_of_transport_annotation{$cid};
57 :     push(@$ref,$annotation);
58 :     }
59 :     else{
60 :     $evidence_of_transport_annotation{$cid} = [$annotation];
61 :     }
62 :    
63 :     my ($peg_index_data,undef) = $fig->search_index($annotation);
64 :     foreach my $peg (map { $_->[0] } @$peg_index_data)
65 :     {
66 :     if($peg =~/$genome_id/){
67 :     if($evidence_of_transport_peg{$cid}){
68 :     my $ref = $evidence_of_transport_peg{$cid};
69 :     push(@$ref,$peg);
70 :     }
71 :     else{
72 :     $evidence_of_transport_peg{$cid} = [$peg];
73 :     }
74 :     }
75 :     }
76 :     }
77 :    
78 :     my %input_cpds;
79 :     open(INPUTS,$filebase."inputs.txt")
80 :     or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs");
81 :     while(<INPUTS>)
82 :     {
83 :     chomp;
84 :     my @line = split "\t", $_;
85 :     chomp @line;
86 :     $input_cpds{$line[0]} = 1;
87 :     }
88 :     close(INPUTS);
89 :     #we are going to load up the cofactor compounds into inputs as well
90 :     my %cofactor_cpds;
91 :     open(COFACT,$filebase."secondary_compounds.txt")
92 :     or die("Failed to open ".$filebase."secondary_compounds.txt");
93 :     while(<COFACT>)
94 :     {
95 :     chomp;
96 :     my @line = split "\t",$_;
97 :     chomp @line;
98 :     $input_cpds{$line[0]} = 1;
99 :     $cofactor_cpds{$line[0]} = 1;
100 :     }
101 :     close(COFACT);
102 :    
103 :    
104 :     my %compound_id_to_name;
105 :     open(IN,"$FIG_Config::global/Models/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt");
106 :     while($_ = <IN>){
107 :     chomp($_);
108 :     my($id,$name) = split("\t",$_);
109 :     $compound_id_to_name{$id} = $name;
110 :     }
111 :     close(IN);
112 :    
113 :     my %biomass_cpds;
114 :     open(BIOMASS,$filebase."biomass_valid.txt")
115 :     or die ("Failed to open ".$filebase."biomass_valid.txt\n");
116 :     while(<BIOMASS>)
117 :     {
118 :     chomp;
119 :     my @temp = split "\t", $_;
120 :     chomp @temp;
121 :     $biomass_cpds{$temp[0]} = 1;
122 :     }
123 :     close(BIOMASS);
124 :    
125 :     open(IN,$filebase."scenarios_valid.txt");
126 :     my @lines = <IN>;
127 :     chomp @lines;
128 :     close(IN);
129 :     my %valid_scenarios;
130 :     foreach (@lines)
131 :     {
132 :     $valid_scenarios{$_} = 1;
133 :     }
134 :    
135 :     my %output_cpds;
136 :     open(OUTS,$filebase."outputs.txt");
137 :     while(<OUTS>)
138 :     {
139 :     chomp;
140 :     my @temp = split "\t", $_;
141 :     chomp @temp;
142 :     $output_cpds{$temp[0]} = 1;
143 :     }
144 :     close(OUTS);
145 :    
146 :     #load core_metabolites, scenario reaction information
147 :     my %core_metabolites;
148 :     my %reaction_to_substrates;
149 :     my %reaction_to_products;
150 :     my %reaction_reversible;
151 :     my %reaction_notes;
152 :     foreach my $scenario (@scenarios)
153 :     {
154 :     if(!defined $valid_scenarios{$scenario->get_scenario_name})
155 :     {
156 :     next;
157 :     }
158 :     $scenario->load_reaction_pegs($fig);
159 :     my $subsystem = $scenario->get_subsystem_name();
160 :     foreach my $reaction (@{$scenario->get_reaction_ids()})
161 :     {
162 :     $reaction_to_substrates{$reaction} = {};
163 :     $reaction_to_products{$reaction} = {};
164 :    
165 :     my %substrates = %{$scenario->get_reaction_substrates($reaction)};
166 :     my %products = %{$scenario->get_reaction_products($reaction)};
167 :    
168 :     $reaction_reversible{$reaction} = $scenario->get_reaction_reversibility($reaction);
169 :     map {$core_metabolites{$_} = 1;
170 :     $reaction_to_substrates{$reaction}->{$_} = $substrates{$_};
171 :     } keys %substrates;
172 :     map {$core_metabolites{$_} = 1;
173 :     $reaction_to_products{$reaction}->{$_} = $products{$_};
174 :     } keys %products;
175 :     if(!defined $reaction_notes{$reaction})
176 :     {
177 :     $reaction_notes{$reaction} = {};
178 :     }
179 :     foreach ( @{$scenario->get_reaction_pegs($reaction) })
180 :     {
181 :     $reaction_notes{$reaction}->{"SEED_PEG:$_"} = 1;
182 :     my $roles = $fig->protein_subsystem_to_roles($_, $subsystem);
183 :     foreach my $role (@$roles)
184 :     {
185 :     $reaction_notes{$reaction}->{"Functional Role:$role"} = 1;
186 :     }
187 :     }
188 :     $reaction_notes{$reaction}->{"Subsystem/Scenario:".$scenario->get_scenario_name()} = 1;
189 :    
190 :     }
191 :     }
192 :    
193 :     #the ATP Synthase hack reaction
194 :     $core_metabolites{"C00009"} = 1;
195 :     $core_metabolites{"C00001"} = 1;
196 :     $core_metabolites{"C00008"} = 1;
197 :     $core_metabolites{"C00002"} = 1;
198 :     $reaction_reversible{"ATPSYN"} = 0;
199 :     $reaction_to_substrates{"ATPSYN"} = {};
200 :     $reaction_to_products{"ATPSYN"} = {};
201 :     $reaction_to_substrates{"ATPSYN"}->{"C00008"} = 1;
202 :     $reaction_to_substrates{"ATPSYN"}->{"C00009"} = 1;
203 :     $reaction_to_products{"ATPSYN"}->{"C00001"} = 1;
204 :     $reaction_to_products{"ATPSYN"}->{"C00002"} = 1;
205 :     $reaction_notes{"ATPSYN"} = {};
206 :     $reaction_notes{"ATPSYN"}->{"Manual Reaction:This fixes problems because we lack electron transport reactions."} = 1;
207 :     #end hack
208 :    
209 :     print STDERR "\tStage 1 - Complete\n" if($debug);
210 :    
211 :     print STDERR "\tStage 2 - Write Spreadsheets file\n" if($debug);
212 :     open(OUT,">".$filebase."compounds_spreadsheet.txt");
213 :     write_species_list();
214 :     close(OUT);
215 :    
216 :     open(OUT,">".$filebase."reactions_spreadsheet.txt");
217 :     write_reactions_list();
218 :     close(OUT);
219 :    
220 :     open(OUT,">".$filebase."transports_spreadsheet.txt");
221 :     write_transports_list();
222 :     close(OUT);
223 :    
224 :     sub write_species_list
225 :     {
226 :     print OUT "KEGG ID\tName\n";
227 :    
228 :     #input compounds
229 :     foreach my $id (sort keys %input_cpds)
230 :     {
231 :     my $name = $compound_id_to_name{$id};
232 :     $name =~ s/\\"//g;
233 :     print OUT "$id\t$name\n";
234 :     }
235 :     #biomass compounds
236 :     foreach my $id (sort keys(%biomass_cpds))
237 :     {
238 :     my $name = $compound_id_to_name{$id};
239 :     $name =~ s/\\"//g;
240 :     if(!$input_cpds{$id})
241 :     {
242 :     print OUT "$id\t$name\n";
243 :     }
244 :     }
245 :     #output compounds
246 :     foreach my $id (sort keys(%output_cpds))
247 :     {
248 :     my $name = $compound_id_to_name{$id};
249 :     $name =~ s/\\"//g;
250 :     if(!$input_cpds{$id} && !$biomass_cpds{$id})
251 :     {
252 :     print OUT "$id\t$name\n";
253 :     }
254 :     }
255 :    
256 :     #core_metabolites are all the compounds used in any scenario.
257 :     foreach my $id (sort keys(%core_metabolites))
258 :     {
259 :     my $name = $compound_id_to_name{$id};
260 :     $name =~ s/\\"//g;
261 :     if(!$input_cpds{$id} && !$biomass_cpds{$id} && !$output_cpds{$id})
262 :     {
263 :     print OUT "$id\t$name\n";
264 :     }
265 :     }
266 :     }
267 :    
268 :     sub write_reactions_list
269 :     {
270 :     print OUT "KEGG RID\tGene\tRole\tSubsystem/Scenario(s)\n";
271 :    
272 :     #Now we write out the Scenario reactions
273 :     foreach my $rid (keys %reaction_to_substrates)
274 :     {
275 :     make_reaction($rid,$reaction_reversible{$rid},$reaction_to_substrates{$rid},$reaction_to_products{$rid},
276 :     "0.000000",$reaction_notes{$rid});
277 :     }
278 :    
279 :     my @biomass_list;
280 :     open(IN,"$FIG_Config::global/Models/ec_biomass.txt");
281 :     while($_ = <IN>){
282 :     chomp($_);
283 :     my($id,$name,$stoich) = split("\t",$_);
284 :     $stoich = $stoich * (.000001);
285 :     push(@biomass_list,"$stoich $id");
286 :     }
287 :    
288 :     print OUT "\n";
289 :     my $biomass_string = join(" + ",@biomass_list);
290 :     print OUT "BIOMAS\t$biomass_string\n";
291 :     close(IN);
292 :    
293 :     }
294 :    
295 :     sub write_transports_list
296 :     {
297 :     print OUT "Direction\tCompound ID\tName\tClass\tGene(s)\tRole\n";
298 :    
299 :     foreach my $input (keys %input_cpds)
300 :     {
301 :     if(!$cofactor_cpds{$input}){
302 :     my $name = $compound_id_to_name{$input};
303 :     if($evidence_of_transport_peg{$input}){
304 :     my $peg_ref = $evidence_of_transport_peg{$input};
305 :     my $genes = join(",",@$peg_ref);
306 :     my $annotation_ref = $evidence_of_transport_annotation{$input};
307 :     my $annotations = join(",",@$annotation_ref);
308 :     print OUT "IN\t$input\t$name\tevidence_of_input\t$genes\t$annotations\n";
309 :     }
310 :     elsif($lib_input_cpds{$input}){
311 :     print OUT "IN\t$input\t$name\tcommon_primary_input\n";
312 :     }
313 :     else{
314 :     print OUT "IN\t$input\t$name\tabstract_input\n";
315 :     }
316 :     }
317 :     }
318 :    
319 :     foreach my $input (keys %cofactor_cpds)
320 :     {
321 :     my $name = $compound_id_to_name{$input};
322 :     if($evidence_of_transport_peg{$input}){
323 :     my $peg_ref = $evidence_of_transport_peg{$input};
324 :     my $genes = join(",",@$peg_ref);
325 :     my $annotation_ref = $evidence_of_transport_annotation{$input};
326 :     my $annotations = join(",",@$annotation_ref);
327 :     print OUT "IN\t$input\t$name\tevidence_of_input\t$genes\t$annotations\n";
328 :     }
329 :    
330 :     else{
331 :     print OUT "IN\t$input\t$name\tcommon_secondary_input\n";
332 :     }
333 :     }
334 :    
335 :     #Make sink/transport reactions for the output compounds
336 :     foreach my $output (keys %output_cpds)
337 :     {
338 :     my $name = $compound_id_to_name{$output};
339 :     print OUT "OUT\t$output\t$name\tcommon_output\n";
340 :     }
341 :    
342 :     #make sinks/transport reactions for all cpds not already outputed or biomassed
343 :     foreach my $id (keys %core_metabolites)
344 :     {
345 :     my $name = $compound_id_to_name{$id};
346 :     if(!$output_cpds{$id} && !$biomass_cpds{$id})
347 :     {
348 :     print OUT "OUT\t$id\t$name\tabstract_output\n";
349 :     }
350 :     }
351 :     }
352 :    
353 :     sub make_reaction
354 :     {
355 :     my($reaction_id,$reversible,$substrates,$products,$objective,$notes) = @_;
356 :    
357 :    
358 :     if($reversible)
359 :     {
360 :     print OUT "$reaction_id\ttrue\t";
361 :     }
362 :     else
363 :     {
364 :     print OUT "$reaction_id\tfalse\t";
365 :     }
366 :     if(defined $notes)
367 :     {
368 :     my @temp_array = ();
369 :     my $temp_string;
370 :     foreach my $note (keys %$notes)
371 :     {
372 :     if($note =~/SEED_PEG:(.*)/){push(@temp_array,$1);}
373 :     }
374 :     $temp_string = join(",",@temp_array);
375 :     print OUT "$temp_string\t";
376 :    
377 :     @temp_array = ();
378 :     foreach my $note (keys %$notes)
379 :     {
380 :     if($note =~/Functional Role:(.*)/){push(@temp_array,$1);}
381 :     }
382 :     $temp_string = join(",",@temp_array);
383 :     print OUT "$temp_string\t";
384 :    
385 :     @temp_array = ();
386 :     foreach my $note (keys %$notes)
387 :     {
388 :     if($note =~/Subsystem\/Scenario:(.*)/){push(@temp_array,$1);}
389 :     }
390 :     $temp_string = join(",",@temp_array);
391 :     print OUT "$temp_string";
392 :    
393 :     }
394 :     print OUT "\n"
395 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3