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

Annotation of /FigKernelScripts/make_MPS_model.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.7 - (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: make_MPS_model <genome_id>\n";
16 :     exit(-1);
17 :     }
18 : formsma 1.2 my $fig;
19 :     if($orgdir)
20 :     {
21 :     $fig = new FIGV($orgdir);
22 :     }
23 :     else
24 :     {
25 :     $fig = new FIGV;
26 :     }
27 : formsma 1.1 chomp $genome_id;
28 : formsma 1.3 Scenario::set_fig($fig,$genome_id);
29 : formsma 1.1 print STDERR "Creating SBML model for 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,1)};
32 : dejongh 1.7 my $filebase = $fig->scenario_directory($genome_id)."/Analysis/";
33 : formsma 1.1
34 :     my %input_cpds;
35 :     open(INPUTS,$filebase."inputs.txt")
36 :     or die ("Failed to open ".$filebase."inputs.txt. Run find_model_information $genome_id inputs");
37 :     while(<INPUTS>)
38 :     {
39 :     chomp;
40 :     my @line = split "\t", $_;
41 :     chomp @line;
42 :     $input_cpds{$line[0]} = 1;
43 :     }
44 :     close(INPUTS);
45 :     #we are going to load up the cofactor compounds into inputs as well
46 : formsma 1.4 open(COFACT,$filebase."secondary_compounds.txt")
47 :     or die("Failed to open ".$filebase."secondary_compounds.txt");
48 : formsma 1.1 while(<COFACT>)
49 :     {
50 :     chomp;
51 :     my @line = split "\t",$_;
52 :     chomp @line;
53 :     $input_cpds{$line[0]} = 1;
54 :     }
55 :     close(COFACT);
56 :    
57 :     my %biomass_cpds;
58 :     open(BIOMASS,$filebase."biomass_valid.txt")
59 :     or die ("Failed to open ".$filebase."biomass_valid.txt\n");
60 :     while(<BIOMASS>)
61 :     {
62 :     chomp;
63 :     my @temp = split "\t", $_;
64 :     chomp @temp;
65 :     $biomass_cpds{$temp[0]} = 1;
66 :     }
67 :     close(BIOMASS);
68 :    
69 : formsma 1.5 #$biomass_cpds{"C00002"} = 1;
70 :    
71 : formsma 1.1 open(IN,$filebase."scenarios_valid.txt");
72 :     my @lines = <IN>;
73 :     chomp @lines;
74 :     close(IN);
75 :     my %valid_scenarios;
76 :     foreach (@lines)
77 :     {
78 :     $valid_scenarios{$_} = 1;
79 :     }
80 :    
81 :     my %output_cpds;
82 :     open(OUTS,$filebase."outputs.txt");
83 :     while(<OUTS>)
84 :     {
85 :     chomp;
86 :     my @temp = split "\t", $_;
87 :     chomp @temp;
88 :     $output_cpds{$temp[0]} = 1;
89 :     }
90 :     close(OUTS);
91 :    
92 :    
93 :    
94 :     #load core_metabolites, scenario reaction information
95 :     my %core_metabolites;
96 :     my %reaction_to_substrates;
97 :     my %reaction_to_products;
98 :     my %reaction_reversible;
99 :     foreach my $scenario (@scenarios)
100 :     {
101 :     if(!defined $valid_scenarios{$scenario->get_scenario_name})
102 :     {
103 :     next;
104 :     }
105 :     foreach my $reaction (@{$scenario->get_reaction_ids()})
106 :     {
107 :     $reaction_to_substrates{$reaction} = {};
108 :     $reaction_to_products{$reaction} = {};
109 :    
110 :     my %substrates = %{$scenario->get_reaction_substrates($reaction)};
111 :     my %products = %{$scenario->get_reaction_products($reaction)};
112 :    
113 :     $reaction_reversible{$reaction} = $scenario->get_reaction_reversibility($reaction);
114 :     map {$core_metabolites{$_} = 1;
115 :     $reaction_to_substrates{$reaction}->{$_} = $substrates{$_};
116 :     } keys %substrates;
117 :     map {$core_metabolites{$_} = 1;
118 :     $reaction_to_products{$reaction}->{$_} = $products{$_};
119 :     } keys %products;
120 :     }
121 :     }
122 : formsma 1.5 #the ATP Synthase hack reaction
123 :     $core_metabolites{"C00009"} = 1;
124 :     $core_metabolites{"C00001"} = 1;
125 :     $core_metabolites{"C00008"} = 1;
126 :     $core_metabolites{"C00002"} = 1;
127 :     $reaction_reversible{"ATPSYN"} = 0;
128 :     $reaction_to_substrates{"ATPSYN"} = {};
129 :     $reaction_to_products{"ATPSYN"} = {};
130 :     $reaction_to_substrates{"ATPSYN"}->{"C00008"} = 1;
131 :     $reaction_to_substrates{"ATPSYN"}->{"C00009"} = 1;
132 :     $reaction_to_products{"ATPSYN"}->{"C00001"} = 1;
133 :     $reaction_to_products{"ATPSYN"}->{"C00002"} = 1;
134 :     #end hack
135 :    
136 : formsma 1.1 open(OUT,">".$filebase.$genome_id.".mps");
137 :    
138 :     print OUT "NAME test\n";
139 :     print OUT "ROWS\n N OBJ\n G BIOMAS\n";
140 :    
141 :     my @all_cpds;
142 :    
143 :     #write the compound list.
144 :     foreach my $id (sort keys %input_cpds)
145 :     {
146 :     print OUT " G $id\n"; #regular cellular
147 :     push @all_cpds, $id;
148 :     $id = substr($id,1);
149 :     print OUT " G E$id\n"; #extra cellular
150 :     push @all_cpds, "E".$id;
151 :     }
152 :     foreach my $id (sort keys %output_cpds)
153 :     {
154 :     if(!$input_cpds{$id})
155 :     {
156 :     print OUT " G $id\n"; #regular cellular
157 :     push @all_cpds, $id;
158 :     $id = substr($id,1);
159 :     print OUT " G E$id\n"; #extra cellular
160 :     push @all_cpds, "E".$id;
161 :     }
162 :     }
163 :     foreach my $id (sort keys %biomass_cpds)
164 :     {
165 :     if(!$input_cpds{$id} && !$output_cpds{$id})
166 :     {
167 :     print OUT " G $id\n"; #regular cellular
168 :     push @all_cpds, $id;
169 : formsma 1.6 $id = substr($id,1);
170 :     print OUT " G E$id\n"; #extra cellular
171 :     push @all_cpds, "E".$id;
172 : formsma 1.1 }
173 :     }
174 :    
175 :     foreach my $id (sort keys %core_metabolites)
176 :     {
177 :     if(!$input_cpds{$id} && !$output_cpds{$id} && !$biomass_cpds{$id})
178 :     {
179 :     print OUT " G $id\n"; #regular cellular
180 :     push @all_cpds, $id;
181 :     $id = substr($id,1);
182 :     print OUT " G E$id\n"; #extra cellular
183 :     push @all_cpds, "E".$id;
184 :     }
185 :     }
186 :    
187 :     #special compounds
188 :     print OUT " G bioout\n";
189 :     push @all_cpds, "bioout";
190 :    
191 :     print OUT "COLUMNS\n";
192 :     #for each reaction print ' 4 spaces $rid 4 spaces $cid 14 spaces $stoich'
193 :    
194 :     my @all_rids;
195 :    
196 :     #scenario reactions
197 :     foreach my $rid (sort keys %reaction_to_substrates)
198 :     {
199 : formsma 1.5
200 : formsma 1.1 foreach (sort keys %{$reaction_to_substrates{$rid}})
201 :     {
202 : formsma 1.5 if(!defined $reaction_to_products{$rid}->{$_})
203 :     {
204 :     print OUT " $rid $_ -$reaction_to_substrates{$rid}->{$_}\n";
205 :     }
206 : formsma 1.1 }
207 :     foreach (sort keys %{$reaction_to_products{$rid}})
208 :     {
209 : formsma 1.5 if(!defined $reaction_to_substrates{$rid}->{$_})
210 :     {
211 :     print OUT " $rid $_ $reaction_to_products{$rid}->{$_}\n";
212 :     }
213 : formsma 1.1 }
214 :     push @all_rids, $rid;
215 :     if($reaction_reversible{$rid})
216 :     {
217 :     foreach (sort keys %{$reaction_to_substrates{$rid}})
218 :     {
219 :     print OUT " $rid\_r $_ $reaction_to_substrates{$rid}->{$_}\n";
220 :     }
221 :     foreach (sort keys %{$reaction_to_products{$rid}})
222 :     {
223 :     print OUT " $rid\_r $_ -$reaction_to_products{$rid}->{$_}\n";
224 :     }
225 :     push @all_rids, $rid."_r";
226 :     }
227 :    
228 :     }
229 :    
230 :     #transport reactions
231 :     foreach my $cid (keys %input_cpds)
232 :     {
233 :     my $trans_id = "T".substr($cid,1);
234 :     my $extern_id = "E".substr($cid,1);
235 :     print OUT " $trans_id $cid 1\n";
236 :     print OUT " $trans_id $extern_id -1\n";
237 :     #lets enable reverse transports to fix some sink rxn issues
238 :     print OUT " $trans_id\_r $cid -1\n";
239 :     print OUT " $trans_id\_r $extern_id 1\n";
240 :     }
241 :     foreach my $cid (keys %output_cpds)
242 :     {
243 :     if(!$input_cpds{$cid})
244 :     {
245 :     my $trans_id = "T".substr($cid,1);
246 :     my $extern_id = "E".substr($cid,1);
247 :     print OUT " $trans_id $cid -1\n";
248 :     print OUT " $trans_id $extern_id 1\n";
249 :     }
250 :     }
251 :    
252 :     #biomass reaction
253 :     foreach my $cid (keys %biomass_cpds)
254 :     {
255 :     print OUT " BIOMAS $cid -1\n";
256 :     }
257 :     print OUT " BIOMAS bioout 1\n";
258 :     #uptake reactions
259 :     foreach my $cid (keys %input_cpds)
260 :     {
261 :     my $temp = substr($cid,1);
262 :     print OUT " U$temp E$temp 1\n";
263 :     }
264 :     #sink reactions
265 :     foreach my $cid (keys %output_cpds)
266 :     {
267 :     my $temp = substr($cid,1);
268 :     print OUT " S$temp E$temp -1\n";
269 :     }
270 :    
271 :     #sink reactions for all?
272 :     foreach my $cid (@all_cpds)
273 :     {
274 :     if($cid =~/E\d{5}/ || $cid eq "bioout")
275 :     {
276 :     next;
277 :     }
278 :     my $temp = substr($cid,1);
279 : formsma 1.6 if(!$input_cpds{$cid} && !$output_cpds{$cid})
280 : formsma 1.1 {
281 :     print OUT " T$temp E$temp 1\n";
282 :     print OUT " T$temp C$temp -1\n";
283 :     }
284 : formsma 1.6 if(!$output_cpds{$cid})
285 : formsma 1.1 {
286 :     print OUT " S$temp E$temp -1\n";
287 :     }
288 :     }
289 :    
290 :    
291 :    
292 :     #OBJ function and biomass
293 :     print OUT " FOBJ OBJ 1\n";
294 :     print OUT " FOBJ bioout -1\n";
295 :     print OUT "RHS\nRANGES\n";
296 :     foreach my $cid (@all_cpds)
297 :     {
298 :     print OUT " CBOUND $cid 0\n";
299 :     }
300 :     print OUT "BOUNDS\n";
301 :     print OUT " UP RBOUND BIOMAS 1000\n";
302 :     foreach my $rid (@all_rids)
303 :     {
304 :     if(!($rid =~/_r$/))
305 :     {
306 :     print OUT " UP RBOUND $rid 1000\n";
307 :     }
308 :     else
309 :     {
310 :     print OUT " UP RBOUND $rid 1000\n";
311 :     }
312 :     }
313 :     foreach my $cid (keys %input_cpds)
314 :     {
315 :     $cid = substr($cid,1);
316 :     print OUT " UP RBOUND T$cid 1000\n";
317 :     print OUT " UP RBOUND U$cid 1000\n";
318 :     }
319 :     foreach my $cid (keys %output_cpds)
320 :     {
321 :     $cid = substr($cid,1);
322 :     print OUT " UP RBOUND T$cid 1000\n";
323 :     print OUT " UP RBOUND S$cid 1000\n";
324 :     }
325 :    
326 :    
327 :     print OUT "ENDATA\n";
328 :     close(OUT);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3