[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.2 - (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 :     Scenario::set_fig($fig);
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 : formsma 1.2 my $filebase = $fig->model_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 :     open(COFACT,$filebase."cofactors.txt")
48 :     or die("Failed to open ".$filebase."cofactors.txt");
49 :     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 :     print STDERR "\tStage 1 - Complete\n" if($debug);
151 :    
152 :     print STDERR "\tStage 2 - Write SBML file\n" if($debug);
153 :     open(OUT,">".$filebase.$genome_id."_FBA_Model.xml");
154 :     write_header();
155 :     write_species_list();
156 :     write_reactions_list();
157 :     print OUT "</model>\n</sbml>\n";
158 :     close(OUT);
159 :    
160 :     print STDERR "\tStage 2 - Complete\n" if($debug);
161 :    
162 :    
163 :     sub write_header
164 :     {
165 :    
166 :     my $print_id = $genome_id;
167 :     $print_id =~ s/\.//g;
168 :     print OUT "<?xml version=\"1.0\" encoding=\"UTF-8\"?>
169 :     <sbml xmlns=\"http://www.sbml.org/sbml/level2\" level=\"2\" version=\"1\" xmlns:html=\"http://www.w3.org/1999/xhtml\">
170 :     <model id=\"Genome_$print_id\" name=\"Genome_$print_id\">
171 :     <listOfUnitDefinitions>
172 :     <unitDefinition id=\"mmol_per_gDW_per_hr\">
173 :     <listOfUnits>
174 :     <unit kind=\"mole\" scale=\"-3\"/>
175 :     <unit kind=\"gram\" exponent=\"-1\"/>
176 :     <unit kind=\"second\" multiplier=\".00027777\" exponent=\"-1\"/>
177 :     </listOfUnits>
178 :     </unitDefinition>
179 :     </listOfUnitDefinitions>
180 :     <listOfCompartments>
181 :     <compartment id=\"Extra_organism\"/>
182 :     <compartment id=\"Cytosol\" outside=\"Extra_organism\"/>
183 :     </listOfCompartments>\n";
184 :     }
185 :    
186 :     sub write_species_list
187 :     {
188 :     print OUT "<listOfSpecies>\n";
189 :    
190 :     #input compounds
191 :     foreach my $id (sort keys %input_cpds)
192 :     {
193 :     my $name = $compound_id_to_name{$id};
194 :     $name =~ s/\\"//g;
195 :     make_species($id,$name,1,1);
196 :     }
197 :     #biomass compounds
198 :     foreach my $id (sort keys(%biomass_cpds))
199 :     {
200 :     my $name = $compound_id_to_name{$id};
201 :     $name =~ s/\\"//g;
202 :     if(!$input_cpds{$id})
203 :     {
204 :     make_species($id,$name,1,1);
205 :     }
206 :     }
207 :     #output compounds
208 :     foreach my $id (sort keys(%output_cpds))
209 :     {
210 :     my $name = $compound_id_to_name{$id};
211 :     $name =~ s/\\"//g;
212 :     if(!$input_cpds{$id} && !$biomass_cpds{$id})
213 :     {
214 :     make_species($id,$name,1,1);
215 :     }
216 :     }
217 :    
218 :     #core_metabolites are all the compounds used in any scenario.
219 :     foreach my $id (sort keys(%core_metabolites))
220 :     {
221 :     my $name = $compound_id_to_name{$id};
222 :     $name =~ s/\\"//g;
223 :     if(!$input_cpds{$id} && !$biomass_cpds{$id} && !$output_cpds{$id})
224 :     {
225 :     make_species($id,$name,1,1);
226 :     }
227 :     }
228 :    
229 :     print OUT "\t<species id=\"biomass_out\" name=\"BIOMASS\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n";
230 :    
231 :     print OUT "</listOfSpecies>\n";
232 :     }
233 :    
234 :     sub write_reactions_list
235 :     {
236 :     print OUT "<listOfReactions>\n";
237 :    
238 :     #make uptake/transport reactions for the input compounds (provide cpds for the system)
239 :     foreach my $input (keys %input_cpds)
240 :     {
241 :     make_reaction($input."_transport_in",0,{$input."_e" => "1.000000"},{$input => "1.000000"},"0.000000",());
242 :     }
243 :    
244 :     #Make sink/transport reactions for the output compounds
245 :     foreach my $output (keys %output_cpds)
246 :     {
247 :     make_reaction($output."_transport_out",0,{$output => "1.000000"},{$output."_e" => "1.000000"},"0.000000");
248 :     }
249 :    
250 :     #make sinks/transport reactions for all cpds not already outputed or biomassed
251 :     foreach my $id (keys %core_metabolites)
252 :     {
253 :     if(!$output_cpds{$id} && !$biomass_cpds{$id})
254 :     {
255 :     make_reaction($id."_sink",0,{$id => "1.000000"},{$id."_e" => "1.000000"},"0.000000");
256 :     }
257 :     }
258 :    
259 :     #Now we write out the Scenario reactions
260 :     foreach my $rid (keys %reaction_to_substrates)
261 :     {
262 :     make_reaction($rid,$reaction_reversible{$rid},$reaction_to_substrates{$rid},$reaction_to_products{$rid},
263 :     "0.000000",$reaction_notes{$rid});
264 :     }
265 :    
266 :     #now its time for the biomass reaction
267 :     make_reaction("BIOMASS",0,\%biomass_cpds,{"biomass_out" => "1.000000"},"-1.000000");
268 :    
269 :     print OUT "</listOfReactions>\n";
270 :     }
271 :    
272 :     sub make_species
273 :     {
274 :     my ($id,$name,$internal,$external) = @_;
275 :     if($internal)
276 :     {
277 :     print OUT "\t<species id=\"$id\" name=\"$name\" compartment=\"Cytosol\" charge=\"0\" boundaryCondition=\"false\"/>\n";
278 :     }
279 :     if($external)
280 :     {
281 :     print OUT "\t<species id=\"$id\_e\" name=\"$name\" compartment=\"Extra_organism\" charge=\"0\" boundaryCondition=\"false\"/>\n";
282 :     }
283 :     }
284 :    
285 :     sub make_reaction
286 :     {
287 :     my($reaction_id,$reversible,$substrates,$products,$objective,$notes) = @_;
288 :    
289 :    
290 :     if($reversible)
291 :     {
292 :     print OUT "<reaction id=\"$reaction_id\" reversible=\"true\">\n";
293 :     }
294 :     else
295 :     {
296 :     print OUT "<reaction id=\"$reaction_id\" reversible=\"false\">\n";
297 :     }
298 :     if(defined $notes)
299 :     {
300 :     print OUT "\t<notes>\n";
301 :     foreach my $note (keys %$notes)
302 :     {
303 :     print OUT "\t\t<html:p>$note</html:p>\n";
304 :     }
305 :     print OUT "\t</notes>\n";
306 :     }
307 :     print OUT "\t<listOfReactants>\n" if(keys %$substrates);
308 :     foreach (keys %$substrates)
309 :     {
310 :     print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$substrates->{$_}\"/>\n";
311 :     }
312 :     print OUT "\t</listOfReactants>\n" if(keys %$substrates);
313 :    
314 :     print OUT "\t<listOfProducts>\n" if(keys %$products);
315 :     foreach (keys %$products)
316 :     {
317 :     print OUT "\t\t<speciesReference species=\"$_\" stoichiometry=\"$products->{$_}\"/>\n";
318 :     }
319 :     print OUT "\t</listOfProducts>\n" if(keys %$products);
320 :    
321 :     print OUT "\t<kineticLaw>\n";
322 :    
323 :     print OUT "\t\t<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n".
324 :     "\t\t\t<apply>\n".
325 :     "\t\t\t\t<ci> LOWER_BOUND </ci>\n".
326 :     "\t\t\t\t<ci> UPPER_BOUND </ci>\n".
327 :     "\t\t\t\t<ci> OBJECTIVE_COEFFICIENT </ci>\n".
328 :     "\t\t\t\t<ci> FLUX_VALUE </ci>\n".
329 :     "\t\t\t\t<ci> REDUCED_COST </ci>\n".
330 :     "\t\t\t</apply>\n\t\t</math>\n";
331 :    
332 :    
333 :     print OUT "\t\t<listOfParameters>\n";
334 :     my $low_bound = "0.000000";
335 :     if($reversible)
336 :     {
337 :     $low_bound = "-999999.000000";
338 :     }
339 :     print OUT "\t\t\t<parameter id=\"LOWER_BOUND\" value=\"$low_bound\" units=\"mmol_per_gDW_per_hr\"/>\n" .
340 :     "\t\t\t<parameter id=\"UPPER_BOUND\" value=\"999999.000000\" units=\"mmol_per_gDW_per_hr\"/>\n".
341 :     "\t\t\t<parameter id=\"OBJECTIVE_COEFFICIENT\" value=\"$objective\"/>\n".
342 :     "\t\t\t<parameter id=\"FLUX_VALUE\" value=\"0.000000\" units=\"mmol_per_gDW_per_hr\"/>\n".
343 :     "\t\t\t<parameter id=\"REDUCED_COST\" value=\"0.000000\"/>\n";
344 :     print OUT "\t\t</listOfParameters>\n".
345 :     "\t</kineticLaw>".
346 :     "\n</reaction>\n";
347 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3