[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.5 - (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 : formsma 1.2 my $filebase = $fig->model_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 :     }
170 :     }
171 :    
172 :     foreach my $id (sort keys %core_metabolites)
173 :     {
174 :     if(!$input_cpds{$id} && !$output_cpds{$id} && !$biomass_cpds{$id})
175 :     {
176 :     print OUT " G $id\n"; #regular cellular
177 :     push @all_cpds, $id;
178 :     $id = substr($id,1);
179 :     print OUT " G E$id\n"; #extra cellular
180 :     push @all_cpds, "E".$id;
181 :     }
182 :     }
183 :    
184 :     #special compounds
185 :     print OUT " G bioout\n";
186 :     push @all_cpds, "bioout";
187 :    
188 :     print OUT "COLUMNS\n";
189 :     #for each reaction print ' 4 spaces $rid 4 spaces $cid 14 spaces $stoich'
190 :    
191 :     my @all_rids;
192 :    
193 :     #scenario reactions
194 :     foreach my $rid (sort keys %reaction_to_substrates)
195 :     {
196 : formsma 1.5
197 : formsma 1.1 foreach (sort keys %{$reaction_to_substrates{$rid}})
198 :     {
199 : formsma 1.5 if(!defined $reaction_to_products{$rid}->{$_})
200 :     {
201 :     print OUT " $rid $_ -$reaction_to_substrates{$rid}->{$_}\n";
202 :     }
203 : formsma 1.1 }
204 :     foreach (sort keys %{$reaction_to_products{$rid}})
205 :     {
206 : formsma 1.5 if(!defined $reaction_to_substrates{$rid}->{$_})
207 :     {
208 :     print OUT " $rid $_ $reaction_to_products{$rid}->{$_}\n";
209 :     }
210 : formsma 1.1 }
211 :     push @all_rids, $rid;
212 :     if($reaction_reversible{$rid})
213 :     {
214 :     foreach (sort keys %{$reaction_to_substrates{$rid}})
215 :     {
216 :     print OUT " $rid\_r $_ $reaction_to_substrates{$rid}->{$_}\n";
217 :     }
218 :     foreach (sort keys %{$reaction_to_products{$rid}})
219 :     {
220 :     print OUT " $rid\_r $_ -$reaction_to_products{$rid}->{$_}\n";
221 :     }
222 :     push @all_rids, $rid."_r";
223 :     }
224 :    
225 :     }
226 :    
227 :     #transport reactions
228 :     foreach my $cid (keys %input_cpds)
229 :     {
230 :     my $trans_id = "T".substr($cid,1);
231 :     my $extern_id = "E".substr($cid,1);
232 :     print OUT " $trans_id $cid 1\n";
233 :     print OUT " $trans_id $extern_id -1\n";
234 :     #lets enable reverse transports to fix some sink rxn issues
235 :     print OUT " $trans_id\_r $cid -1\n";
236 :     print OUT " $trans_id\_r $extern_id 1\n";
237 :     }
238 :     foreach my $cid (keys %output_cpds)
239 :     {
240 :     if(!$input_cpds{$cid})
241 :     {
242 :     my $trans_id = "T".substr($cid,1);
243 :     my $extern_id = "E".substr($cid,1);
244 :     print OUT " $trans_id $cid -1\n";
245 :     print OUT " $trans_id $extern_id 1\n";
246 :     }
247 :     }
248 :    
249 :     #biomass reaction
250 :     foreach my $cid (keys %biomass_cpds)
251 :     {
252 :     print OUT " BIOMAS $cid -1\n";
253 :     }
254 :     print OUT " BIOMAS bioout 1\n";
255 :     #uptake reactions
256 :     foreach my $cid (keys %input_cpds)
257 :     {
258 :     my $temp = substr($cid,1);
259 :     print OUT " U$temp E$temp 1\n";
260 :     }
261 :     #sink reactions
262 :     foreach my $cid (keys %output_cpds)
263 :     {
264 :     my $temp = substr($cid,1);
265 :     print OUT " S$temp E$temp -1\n";
266 :     }
267 :    
268 :     #sink reactions for all?
269 :     foreach my $cid (@all_cpds)
270 :     {
271 :     if($cid =~/E\d{5}/ || $cid eq "bioout")
272 :     {
273 :     next;
274 :     }
275 :     my $temp = substr($cid,1);
276 :     if(!$input_cpds{$cid} && !$output_cpds{$cid} && !$biomass_cpds{$cid})
277 :     {
278 :     print OUT " T$temp E$temp 1\n";
279 :     print OUT " T$temp C$temp -1\n";
280 :     }
281 :     elsif(!$output_cpds{$cid} && !$biomass_cpds{$cid})
282 :     {
283 :     print OUT " S$temp E$temp -1\n";
284 :     }
285 :     }
286 :    
287 :    
288 :    
289 :     #OBJ function and biomass
290 :     print OUT " FOBJ OBJ 1\n";
291 :     print OUT " FOBJ bioout -1\n";
292 :     print OUT "RHS\nRANGES\n";
293 :     foreach my $cid (@all_cpds)
294 :     {
295 :     print OUT " CBOUND $cid 0\n";
296 :     }
297 :     print OUT "BOUNDS\n";
298 :     print OUT " UP RBOUND BIOMAS 1000\n";
299 :     foreach my $rid (@all_rids)
300 :     {
301 :     if(!($rid =~/_r$/))
302 :     {
303 :     print OUT " UP RBOUND $rid 1000\n";
304 :     }
305 :     else
306 :     {
307 :     print OUT " UP RBOUND $rid 1000\n";
308 :     }
309 :     }
310 :     foreach my $cid (keys %input_cpds)
311 :     {
312 :     $cid = substr($cid,1);
313 :     print OUT " UP RBOUND T$cid 1000\n";
314 :     print OUT " UP RBOUND U$cid 1000\n";
315 :     }
316 :     foreach my $cid (keys %output_cpds)
317 :     {
318 :     $cid = substr($cid,1);
319 :     print OUT " UP RBOUND T$cid 1000\n";
320 :     print OUT " UP RBOUND S$cid 1000\n";
321 :     }
322 :     #sink reactions here eh.
323 :    
324 :    
325 :     print OUT "ENDATA\n";
326 :     close(OUT);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3