[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.3 - (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 :     open(COFACT,$filebase."cofactors.txt")
47 :     or die("Failed to open ".$filebase."cofactors.txt");
48 :     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 :     open(IN,$filebase."scenarios_valid.txt");
70 :     my @lines = <IN>;
71 :     chomp @lines;
72 :     close(IN);
73 :     my %valid_scenarios;
74 :     foreach (@lines)
75 :     {
76 :     $valid_scenarios{$_} = 1;
77 :     }
78 :    
79 :     my %output_cpds;
80 :     open(OUTS,$filebase."outputs.txt");
81 :     while(<OUTS>)
82 :     {
83 :     chomp;
84 :     my @temp = split "\t", $_;
85 :     chomp @temp;
86 :     $output_cpds{$temp[0]} = 1;
87 :     }
88 :     close(OUTS);
89 :    
90 :    
91 :    
92 :     #load core_metabolites, scenario reaction information
93 :     my %core_metabolites;
94 :     my %reaction_to_substrates;
95 :     my %reaction_to_products;
96 :     my %reaction_reversible;
97 :     foreach my $scenario (@scenarios)
98 :     {
99 :     if(!defined $valid_scenarios{$scenario->get_scenario_name})
100 :     {
101 :     next;
102 :     }
103 :     foreach my $reaction (@{$scenario->get_reaction_ids()})
104 :     {
105 :     $reaction_to_substrates{$reaction} = {};
106 :     $reaction_to_products{$reaction} = {};
107 :    
108 :     my %substrates = %{$scenario->get_reaction_substrates($reaction)};
109 :     my %products = %{$scenario->get_reaction_products($reaction)};
110 :    
111 :     $reaction_reversible{$reaction} = $scenario->get_reaction_reversibility($reaction);
112 :     map {$core_metabolites{$_} = 1;
113 :     $reaction_to_substrates{$reaction}->{$_} = $substrates{$_};
114 :     } keys %substrates;
115 :     map {$core_metabolites{$_} = 1;
116 :     $reaction_to_products{$reaction}->{$_} = $products{$_};
117 :     } keys %products;
118 :     }
119 :     }
120 :     open(OUT,">".$filebase.$genome_id.".mps");
121 :    
122 :     print OUT "NAME test\n";
123 :     print OUT "ROWS\n N OBJ\n G BIOMAS\n";
124 :    
125 :     my @all_cpds;
126 :    
127 :     #write the compound list.
128 :     foreach my $id (sort keys %input_cpds)
129 :     {
130 :     print OUT " G $id\n"; #regular cellular
131 :     push @all_cpds, $id;
132 :     $id = substr($id,1);
133 :     print OUT " G E$id\n"; #extra cellular
134 :     push @all_cpds, "E".$id;
135 :     }
136 :     foreach my $id (sort keys %output_cpds)
137 :     {
138 :     if(!$input_cpds{$id})
139 :     {
140 :     print OUT " G $id\n"; #regular cellular
141 :     push @all_cpds, $id;
142 :     $id = substr($id,1);
143 :     print OUT " G E$id\n"; #extra cellular
144 :     push @all_cpds, "E".$id;
145 :     }
146 :     }
147 :     foreach my $id (sort keys %biomass_cpds)
148 :     {
149 :     if(!$input_cpds{$id} && !$output_cpds{$id})
150 :     {
151 :     print OUT " G $id\n"; #regular cellular
152 :     push @all_cpds, $id;
153 :     }
154 :     }
155 :    
156 :     foreach my $id (sort keys %core_metabolites)
157 :     {
158 :     if(!$input_cpds{$id} && !$output_cpds{$id} && !$biomass_cpds{$id})
159 :     {
160 :     print OUT " G $id\n"; #regular cellular
161 :     push @all_cpds, $id;
162 :     $id = substr($id,1);
163 :     print OUT " G E$id\n"; #extra cellular
164 :     push @all_cpds, "E".$id;
165 :     }
166 :     }
167 :    
168 :     #special compounds
169 :     print OUT " G bioout\n";
170 :     push @all_cpds, "bioout";
171 :    
172 :     print OUT "COLUMNS\n";
173 :     #for each reaction print ' 4 spaces $rid 4 spaces $cid 14 spaces $stoich'
174 :    
175 :     my @all_rids;
176 :    
177 :     #scenario reactions
178 :     foreach my $rid (sort keys %reaction_to_substrates)
179 :     {
180 :     foreach (sort keys %{$reaction_to_substrates{$rid}})
181 :     {
182 :     print OUT " $rid $_ -$reaction_to_substrates{$rid}->{$_}\n";
183 :     }
184 :     foreach (sort keys %{$reaction_to_products{$rid}})
185 :     {
186 :     print OUT " $rid $_ $reaction_to_products{$rid}->{$_}\n";
187 :     }
188 :     push @all_rids, $rid;
189 :     if($reaction_reversible{$rid})
190 :     {
191 :     foreach (sort keys %{$reaction_to_substrates{$rid}})
192 :     {
193 :     print OUT " $rid\_r $_ $reaction_to_substrates{$rid}->{$_}\n";
194 :     }
195 :     foreach (sort keys %{$reaction_to_products{$rid}})
196 :     {
197 :     print OUT " $rid\_r $_ -$reaction_to_products{$rid}->{$_}\n";
198 :     }
199 :     push @all_rids, $rid."_r";
200 :     }
201 :    
202 :     }
203 :    
204 :     #transport reactions
205 :     foreach my $cid (keys %input_cpds)
206 :     {
207 :     my $trans_id = "T".substr($cid,1);
208 :     my $extern_id = "E".substr($cid,1);
209 :     print OUT " $trans_id $cid 1\n";
210 :     print OUT " $trans_id $extern_id -1\n";
211 :     #lets enable reverse transports to fix some sink rxn issues
212 :     print OUT " $trans_id\_r $cid -1\n";
213 :     print OUT " $trans_id\_r $extern_id 1\n";
214 :     }
215 :     foreach my $cid (keys %output_cpds)
216 :     {
217 :     if(!$input_cpds{$cid})
218 :     {
219 :     my $trans_id = "T".substr($cid,1);
220 :     my $extern_id = "E".substr($cid,1);
221 :     print OUT " $trans_id $cid -1\n";
222 :     print OUT " $trans_id $extern_id 1\n";
223 :     }
224 :     }
225 :    
226 :     #biomass reaction
227 :     foreach my $cid (keys %biomass_cpds)
228 :     {
229 :     print OUT " BIOMAS $cid -1\n";
230 :     }
231 :     print OUT " BIOMAS bioout 1\n";
232 :     #uptake reactions
233 :     foreach my $cid (keys %input_cpds)
234 :     {
235 :     my $temp = substr($cid,1);
236 :     print OUT " U$temp E$temp 1\n";
237 :     }
238 :     #sink reactions
239 :     foreach my $cid (keys %output_cpds)
240 :     {
241 :     my $temp = substr($cid,1);
242 :     print OUT " S$temp E$temp -1\n";
243 :     }
244 :    
245 :     #sink reactions for all?
246 :     foreach my $cid (@all_cpds)
247 :     {
248 :     if($cid =~/E\d{5}/ || $cid eq "bioout")
249 :     {
250 :     next;
251 :     }
252 :     my $temp = substr($cid,1);
253 :     if(!$input_cpds{$cid} && !$output_cpds{$cid} && !$biomass_cpds{$cid})
254 :     {
255 :     print OUT " T$temp E$temp 1\n";
256 :     print OUT " T$temp C$temp -1\n";
257 :     }
258 :     elsif(!$output_cpds{$cid} && !$biomass_cpds{$cid})
259 :     {
260 :     print OUT " S$temp E$temp -1\n";
261 :     }
262 :     }
263 :    
264 :    
265 :    
266 :     #OBJ function and biomass
267 :     print OUT " FOBJ OBJ 1\n";
268 :     print OUT " FOBJ bioout -1\n";
269 :     print OUT "RHS\nRANGES\n";
270 :     foreach my $cid (@all_cpds)
271 :     {
272 :     print OUT " CBOUND $cid 0\n";
273 :     }
274 :     print OUT "BOUNDS\n";
275 :     print OUT " UP RBOUND BIOMAS 1000\n";
276 :     foreach my $rid (@all_rids)
277 :     {
278 :     if(!($rid =~/_r$/))
279 :     {
280 :     print OUT " UP RBOUND $rid 1000\n";
281 :     }
282 :     else
283 :     {
284 :     print OUT " UP RBOUND $rid 1000\n";
285 :     }
286 :     }
287 :     foreach my $cid (keys %input_cpds)
288 :     {
289 :     $cid = substr($cid,1);
290 :     print OUT " UP RBOUND T$cid 1000\n";
291 :     print OUT " UP RBOUND U$cid 1000\n";
292 :     }
293 :     foreach my $cid (keys %output_cpds)
294 :     {
295 :     $cid = substr($cid,1);
296 :     print OUT " UP RBOUND T$cid 1000\n";
297 :     print OUT " UP RBOUND S$cid 1000\n";
298 :     }
299 :     #sink reactions here eh.
300 :    
301 :    
302 :     print OUT "ENDATA\n";
303 :     close(OUT);

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3