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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3