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

Annotation of /FigKernelScripts/make_SBML_model.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (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
14 :     unless (scalar (@ARGV) >= 1)
15 : formsma 1.1 {
16 :     print STDERR "usage: make_SBML_model <genome_id>\n";
17 :     exit(-1);
18 :     }
19 :     chomp $genome_id;
20 : formsma 1.2 my $fig;
21 :     if($orgdir)
22 :     {
23 :     $fig = new FIGV($orgdir);
24 :     }
25 :     else
26 :     {
27 :     $fig = new FIGV;
28 :     }
29 : formsma 1.3 Scenario::set_fig($fig,$genome_id);
30 : formsma 1.1 print STDERR "Creating SBML model for genome $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.6 my $filebase = $fig->scenario_directory($genome_id)."/Analysis/";
34 : formsma 1.1
35 :     my %input_cpds;
36 :     open(INPUTS,$filebase."inputs.txt")
37 :     or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs");
38 :     while(<INPUTS>)
39 :     {
40 :     chomp;
41 :     my @line = split "\t", $_;
42 :     chomp @line;
43 :     $input_cpds{$line[0]} = 1;
44 :     }
45 :     close(INPUTS);
46 :     #we are going to load up the cofactor compounds into inputs as well
47 : formsma 1.4 open(COFACT,$filebase."secondary_compounds.txt")
48 :     or die("Failed to open ".$filebase."secondary_compounds.txt");
49 : formsma 1.1 while(<COFACT>)
50 :     {
51 :     chomp;
52 :     my @line = split "\t",$_;
53 :     chomp @line;
54 :     $input_cpds{$line[0]} = 1;
55 :     }
56 :     close(COFACT);
57 :    
58 :    
59 :     my %compound_id_to_name;
60 :     open(IN,"$FIG_Config::global/Models/compound_id_to_name.txt") or die("Failed to open compound_id_to_name.txt");
61 :     while($_ = <IN>){
62 :     chomp($_);
63 :     my($id,$name) = split("\t",$_);
64 :     $compound_id_to_name{$id} = $name;
65 :     }
66 :     close(IN);
67 :    
68 :     my %biomass_cpds;
69 :     open(BIOMASS,$filebase."biomass_valid.txt")
70 :     or die ("Failed to open ".$filebase."biomass_valid.txt\n");
71 :     while(<BIOMASS>)
72 :     {
73 :     chomp;
74 :     my @temp = split "\t", $_;
75 :     chomp @temp;
76 :     $biomass_cpds{$temp[0]} = 1;
77 :     }
78 :     close(BIOMASS);
79 :    
80 :     open(IN,$filebase."scenarios_valid.txt");
81 :     my @lines = <IN>;
82 :     chomp @lines;
83 :     close(IN);
84 :     my %valid_scenarios;
85 :     foreach (@lines)
86 :     {
87 :     $valid_scenarios{$_} = 1;
88 :     }
89 :    
90 :     my %output_cpds;
91 :     open(OUTS,$filebase."outputs.txt");
92 :     while(<OUTS>)
93 :     {
94 :     chomp;
95 :     my @temp = split "\t", $_;
96 :     chomp @temp;
97 :     $output_cpds{$temp[0]} = 1;
98 :     }
99 :     close(OUTS);
100 :    
101 :    
102 :    
103 :     #load core_metabolites, scenario reaction information
104 :     my %core_metabolites;
105 :     my %reaction_to_substrates;
106 :     my %reaction_to_products;
107 :     my %reaction_reversible;
108 :     my %reaction_notes;
109 :     foreach my $scenario (@scenarios)
110 :     {
111 :     if(!defined $valid_scenarios{$scenario->get_scenario_name})
112 :     {
113 :     next;
114 :     }
115 :     $scenario->load_reaction_pegs($fig);
116 :     my $subsystem = $scenario->get_subsystem_name();
117 :     foreach my $reaction (@{$scenario->get_reaction_ids()})
118 :     {
119 :     $reaction_to_substrates{$reaction} = {};
120 :     $reaction_to_products{$reaction} = {};
121 :    
122 :     my %substrates = %{$scenario->get_reaction_substrates($reaction)};
123 :     my %products = %{$scenario->get_reaction_products($reaction)};
124 :    
125 :     $reaction_reversible{$reaction} = $scenario->get_reaction_reversibility($reaction);
126 :     map {$core_metabolites{$_} = 1;
127 :     $reaction_to_substrates{$reaction}->{$_} = $substrates{$_};
128 :     } keys %substrates;
129 :     map {$core_metabolites{$_} = 1;
130 :     $reaction_to_products{$reaction}->{$_} = $products{$_};
131 :     } keys %products;
132 :     if(!defined $reaction_notes{$reaction})
133 :     {
134 :     $reaction_notes{$reaction} = {};
135 :     }
136 :     foreach ( @{$scenario->get_reaction_pegs($reaction) })
137 :     {
138 :     $reaction_notes{$reaction}->{"SEED_PEG:$_"} = 1;
139 :     my $roles = $fig->protein_subsystem_to_roles($_, $subsystem);
140 :     foreach my $role (@$roles)
141 :     {
142 :     $reaction_notes{$reaction}->{"Functional Role:$role"} = 1;
143 :     }
144 :     }
145 :     $reaction_notes{$reaction}->{"Subsystem/Scenario:".$scenario->get_scenario_name()} = 1;
146 :    
147 :     }
148 :     }
149 :    
150 : formsma 1.5 #the ATP Synthase hack reaction
151 :     $core_metabolites{"C00009"} = 1;
152 :     $core_metabolites{"C00001"} = 1;
153 :     $core_metabolites{"C00008"} = 1;
154 :     $core_metabolites{"C00002"} = 1;
155 :     $reaction_reversible{"ATPSYN"} = 0;
156 :     $reaction_to_substrates{"ATPSYN"} = {};
157 :     $reaction_to_products{"ATPSYN"} = {};
158 :     $reaction_to_substrates{"ATPSYN"}->{"C00008"} = 1;
159 :     $reaction_to_substrates{"ATPSYN"}->{"C00009"} = 1;
160 :     $reaction_to_products{"ATPSYN"}->{"C00001"} = 1;
161 :     $reaction_to_products{"ATPSYN"}->{"C00002"} = 1;
162 :     $reaction_notes{"ATPSYN"} = {};
163 :     $reaction_notes{"ATPSYN"}->{"Manual Reaction:This fixes problems because we lack electron transport reactions."} = 1;
164 :     #end hack
165 :    
166 : formsma 1.1 print STDERR "\tStage 1 - Complete\n" if($debug);
167 :    
168 :     print STDERR "\tStage 2 - Write SBML file\n" if($debug);
169 :     open(OUT,">".$filebase.$genome_id."_FBA_Model.xml");
170 :     write_header();
171 :     write_species_list();
172 :     write_reactions_list();
173 :     print OUT "</model>\n</sbml>\n";
174 :     close(OUT);
175 :    
176 :     print STDERR "\tStage 2 - Complete\n" if($debug);
177 :    
178 :    
179 :     sub write_header
180 :     {
181 :    
182 :     my $print_id = $genome_id;
183 :     $print_id =~ s/\.//g;
184 :     print OUT "<?xml version=\"1.0\" encoding=\"UTF-8\"?>
185 :     <sbml xmlns=\"http://www.sbml.org/sbml/level2\" level=\"2\" version=\"1\" xmlns:html=\"http://www.w3.org/1999/xhtml\">
186 :     <model id=\"Genome_$print_id\" name=\"Genome_$print_id\">
187 :     <listOfUnitDefinitions>
188 :     <unitDefinition id=\"mmol_per_gDW_per_hr\">
189 :     <listOfUnits>
190 :     <unit kind=\"mole\" scale=\"-3\"/>
191 :     <unit kind=\"gram\" exponent=\"-1\"/>
192 :     <unit kind=\"second\" multiplier=\".00027777\" exponent=\"-1\"/>
193 :     </listOfUnits>
194 :     </unitDefinition>
195 :     </listOfUnitDefinitions>
196 :     <listOfCompartments>
197 :     <compartment id=\"Extra_organism\"/>
198 :     <compartment id=\"Cytosol\" outside=\"Extra_organism\"/>
199 :     </listOfCompartments>\n";
200 :     }
201 :    
202 :     sub write_species_list
203 :     {
204 :     print OUT "<listOfSpecies>\n";
205 :    
206 :     #input compounds
207 :     foreach my $id (sort keys %input_cpds)
208 :     {
209 :     my $name = $compound_id_to_name{$id};
210 :     $name =~ s/\\"//g;
211 :     make_species($id,$name,1,1);
212 :     }
213 :     #biomass compounds
214 :     foreach my $id (sort keys(%biomass_cpds))
215 :     {
216 :     my $name = $compound_id_to_name{$id};
217 :     $name =~ s/\\"//g;
218 :     if(!$input_cpds{$id})
219 :     {
220 :     make_species($id,$name,1,1);
221 :     }
222 :     }
223 :     #output compounds
224 :     foreach my $id (sort keys(%output_cpds))
225 :     {
226 :     my $name = $compound_id_to_name{$id};
227 :     $name =~ s/\\"//g;
228 :     if(!$input_cpds{$id} && !$biomass_cpds{$id})
229 :     {
230 :     make_species($id,$name,1,1);
231 :     }
232 :     }
233 :    
234 :     #core_metabolites are all the compounds used in any scenario.
235 :     foreach my $id (sort keys(%core_metabolites))
236 :     {
237 :     my $name = $compound_id_to_name{$id};
238 :     $name =~ s/\\"//g;
239 :     if(!$input_cpds{$id} && !$biomass_cpds{$id} && !$output_cpds{$id})
240 :     {
241 :     make_species($id,$name,1,1);
242 :     }
243 :     }
244 :    
245 :     print OUT "\t<species id=\"biomass_out\" name=\"BIOMASS\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n";
246 :    
247 :     print OUT "</listOfSpecies>\n";
248 :     }
249 :    
250 :     sub write_reactions_list
251 :     {
252 :     print OUT "<listOfReactions>\n";
253 :    
254 :     #make uptake/transport reactions for the input compounds (provide cpds for the system)
255 :     foreach my $input (keys %input_cpds)
256 :     {
257 :     make_reaction($input."_transport_in",0,{$input."_e" => "1.000000"},{$input => "1.000000"},"0.000000",());
258 :     }
259 :    
260 :     #Make sink/transport reactions for the output compounds
261 :     foreach my $output (keys %output_cpds)
262 :     {
263 :     make_reaction($output."_transport_out",0,{$output => "1.000000"},{$output."_e" => "1.000000"},"0.000000");
264 :     }
265 :    
266 :     #make sinks/transport reactions for all cpds not already outputed or biomassed
267 :     foreach my $id (keys %core_metabolites)
268 :     {
269 :     if(!$output_cpds{$id} && !$biomass_cpds{$id})
270 :     {
271 :     make_reaction($id."_sink",0,{$id => "1.000000"},{$id."_e" => "1.000000"},"0.000000");
272 :     }
273 :     }
274 :    
275 :     #Now we write out the Scenario reactions
276 :     foreach my $rid (keys %reaction_to_substrates)
277 :     {
278 :     make_reaction($rid,$reaction_reversible{$rid},$reaction_to_substrates{$rid},$reaction_to_products{$rid},
279 :     "0.000000",$reaction_notes{$rid});
280 :     }
281 :    
282 :     #now its time for the biomass reaction
283 :     make_reaction("BIOMASS",0,\%biomass_cpds,{"biomass_out" => "1.000000"},"-1.000000");
284 :    
285 :     print OUT "</listOfReactions>\n";
286 :     }
287 :    
288 :     sub make_species
289 :     {
290 :     my ($id,$name,$internal,$external) = @_;
291 :     if($internal)
292 :     {
293 :     print OUT "\t<species id=\"$id\" name=\"$name\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n";
294 :     }
295 :     if($external)
296 :     {
297 :     print OUT "\t<species id=\"$id\_e\" name=\"$name\" compartment=\"Extra_organism\" charge=\"0\" boundaryCondition=\"false\"/>\n";
298 :     }
299 :     }
300 :    
301 :     sub make_reaction
302 :     {
303 :     my($reaction_id,$reversible,$substrates,$products,$objective,$notes) = @_;
304 :    
305 :    
306 :     if($reversible)
307 :     {
308 :     print OUT "<reaction id=\"$reaction_id\" reversible=\"true\">\n";
309 :     }
310 :     else
311 :     {
312 :     print OUT "<reaction id=\"$reaction_id\" reversible=\"false\">\n";
313 :     }
314 :     if(defined $notes)
315 :     {
316 :     print OUT "\t<notes>\n";
317 :     foreach my $note (keys %$notes)
318 :     {
319 :     print OUT "\t\t<html:p>$note</html:p>\n";
320 :     }
321 :     print OUT "\t</notes>\n";
322 :     }
323 :     print OUT "\t<listOfReactants>\n" if(keys %$substrates);
324 :     foreach (keys %$substrates)
325 :     {
326 :     print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$substrates->{$_}\"/>\n";
327 :     }
328 :     print OUT "\t</listOfReactants>\n" if(keys %$substrates);
329 :    
330 :     print OUT "\t<listOfProducts>\n" if(keys %$products);
331 :     foreach (keys %$products)
332 :     {
333 :     print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$products->{$_}\"/>\n";
334 :     }
335 :     print OUT "\t</listOfProducts>\n" if(keys %$products);
336 :    
337 :     print OUT "\t<kineticLaw>\n";
338 :    
339 :     print OUT "\t\t<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n".
340 :     "\t\t\t<apply>\n".
341 :     "\t\t\t\t<ci> LOWER_BOUND </ci>\n".
342 :     "\t\t\t\t<ci> UPPER_BOUND </ci>\n".
343 :     "\t\t\t\t<ci> OBJECTIVE_COEFFICIENT </ci>\n".
344 :     "\t\t\t\t<ci> FLUX_VALUE </ci>\n".
345 :     "\t\t\t\t<ci> REDUCED_COST </ci>\n".
346 :     "\t\t\t</apply>\n\t\t</math>\n";
347 :    
348 :    
349 :     print OUT "\t\t<listOfParameters>\n";
350 :     my $low_bound = "0.000000";
351 :     if($reversible)
352 :     {
353 :     $low_bound = "-999999.000000";
354 :     }
355 :     print OUT "\t\t\t<parameter id=\"LOWER_BOUND\" value=\"$low_bound\" units=\"mmol_per_gDW_per_hr\"/>\n" .
356 :     "\t\t\t<parameter id=\"UPPER_BOUND\" value=\"999999.000000\" units=\"mmol_per_gDW_per_hr\"/>\n".
357 :     "\t\t\t<parameter id=\"OBJECTIVE_COEFFICIENT\" value=\"$objective\"/>\n".
358 :     "\t\t\t<parameter id=\"FLUX_VALUE\" value=\"0.000000\" units=\"mmol_per_gDW_per_hr\"/>\n".
359 :     "\t\t\t<parameter id=\"REDUCED_COST\" value=\"0.000000\"/>\n";
360 :     print OUT "\t\t</listOfParameters>\n".
361 :     "\t</kineticLaw>".
362 :     "\n</reaction>\n";
363 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3