[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.1 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3