[Bio] / DeJonghStuff / load_kegg.pl Repository:
ViewVC logotype

Annotation of /DeJonghStuff/load_kegg.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (view) (download) (as text)

1 : dejongh 1.4 #
2 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
3 :     # for Interpretations of Genomes. All Rights Reserved.
4 :     #
5 :     # This file is part of the SEED Toolkit.
6 :     #
7 :     # The SEED Toolkit is free software. You can redistribute
8 :     # it and/or modify it under the terms of the SEED Toolkit
9 :     # Public License.
10 :     #
11 :     # You should have received a copy of the SEED Toolkit Public License
12 :     # along with this program; if not write to the University of Chicago
13 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
14 :     # Genomes at veronika@thefig.info or download a copy from
15 :     # http://www.theseed.org/LICENSE.TXT.
16 :     #
17 :    
18 : dejongh 1.1
19 :     # -*- perl -*-
20 :    
21 :     ###########################################
22 :     use strict;
23 :    
24 :     # KEY RELATIONAL TABLES:
25 :     #
26 :     # 1. compound(Cid,Priority,Name) % we have a prioritized set of names for a compound
27 :     # 2. comp_to_CAS(Cid,CASid) % connection to chemical abstract society [optional]
28 :     # 3. reaction(Rid,Reversible) % nonreversible go from substrates to products
29 :     # 4. reaction_to_compound(Rid,0/1[substrate or product],Cid,Stoich,main_compound [part of major transformation])
30 :     # 5. reaction_to_role(Rid,FunctionalRole)
31 :     #
32 :    
33 :     use FIG;
34 :     my $fig = new FIG;
35 :    
36 :     my $usage = "usage: load_kegg";
37 :    
38 :     use Tracer;
39 : dejongh 1.4 TSetup('1 *', 'WARN');
40 : dejongh 1.1 &load_ec_and_map_data;
41 :     &load_compounds;
42 :     &load_reactions;
43 :     &load_catalyzes;
44 :    
45 :     undef $fig;
46 :    
47 :     sub load_ec_and_map_data {
48 :    
49 :     Open(\*TMPIN, "<$FIG_Config::data/KEGG/enzyme");
50 :     Open(\*ECMAP,">$FIG_Config::temp/ec_map.table");
51 :    
52 :     Trace("Reading KEGG enzymes.") if T(2);
53 :     my($ec,%name,$map);
54 :     $/ = "\n///\n";
55 :     while (defined($_ = <TMPIN>))
56 :     {
57 :     if ($_ =~ /ENTRY\s+EC\s+(\d+\.\d+\.\d+\.\d+)/s)
58 :     {
59 :     $ec = $1;
60 :     while ($_ =~ /PATH:\s+(map\d+)\s+(\S[^\n]+\S)/sg)
61 :     {
62 : dejongh 1.2 print ECMAP "$ec\t$1\n";
63 : dejongh 1.1 $name{$1} = $2;
64 :     }
65 :     }
66 :     }
67 :     $/ = "\n";
68 :     close(TMPIN);
69 :     close(ECMAP);
70 :    
71 :     Trace("Writing map table.") if T(2);
72 :     Open(\*MAP, ">$FIG_Config::temp/map_name.table");
73 :    
74 :     foreach $map (keys(%name))
75 :     {
76 :     print MAP "$map\t$name{$map}\n";
77 :     }
78 :     close(MAP);
79 :    
80 :     $fig->reload_table('all', "ec_map",
81 :     "ec varchar(100), map varchar(100)",
82 :     { index_ec_map_ec => "ec", index_ec_map_map => "map" },
83 :     "$FIG_Config::temp/ec_map.table");
84 :     unlink("$FIG_Config::temp/ec_map.table");
85 :     $fig->reload_table('all', "map_name",
86 :     "map varchar(100) UNIQUE NOT NULL, mapname varchar(200), primary key ( map )",
87 :     { }, "$FIG_Config::temp/map_name.table");
88 :     unlink("$FIG_Config::temp/map_name.table");
89 :     }
90 :    
91 :     sub load_compounds {
92 :    
93 :     Trace("Loading compounds.") if T(2);
94 :    
95 :     Open(\*TMPIN, "<$FIG_Config::data/KEGG/compound");
96 :     Open(\*COMP, ">$FIG_Config::temp/comp_name.table");
97 :     Open(\*CAS, ">$FIG_Config::temp/comp_cas.table");
98 :    
99 :     my($cid,$name,$cas,$names,$tmp,$n,$entry);
100 :     $/ = "\n///\n";
101 :     while (defined($entry = <TMPIN>))
102 :     {
103 : dejongh 1.3 if ($entry =~ /ENTRY\s+(C\d+).*\nNAME\s+(\S[^\n]*)\n((\s+(\S[^\n]*\S)\n)*)/s)
104 : dejongh 1.1 {
105 :     $cid = $1;
106 :     $names = $2;
107 :    
108 :     if ($3)
109 :     {
110 :     $tmp = $3;
111 :     chop $tmp;
112 :     $tmp =~ s/^\s+/ /;
113 :     $names = $names . $tmp;
114 :     $names =~ s/\n\s+/ /g;
115 :     $names =~ s/- /-/g;
116 :     }
117 :    
118 :     $n = 1;
119 :     foreach $name (map { $_ =~ s/^\s+//; $_ =~ s/\s+$//; $_ } split(/;/,$names))
120 :     {
121 :     print COMP "$cid\t$n\t$name\n";
122 :     $n++;
123 :     if (length $name > 200) { print "$cid, $name\n" }
124 :     }
125 :     }
126 :    
127 :     if ($entry =~ /DBLINKS\s+CAS:\s+(\S+)/s)
128 :     {
129 :     print CAS "$cid\t$1\n";
130 :     }
131 :     }
132 :     $/ = "\n";
133 :     close(TMPIN);
134 :     close(COMP);
135 :     close(CAS);
136 :    
137 :     $fig->reload_table('all', "comp_name",
138 :     "cid varchar(7), pos integer, name varchar(200)",
139 :     { index_comp_name_cid => "cid",
140 :     index_comp_name_name => "name" },
141 :     "$FIG_Config::temp/comp_name.table");
142 :     unlink("$FIG_Config::temp/comp_name.table");
143 :    
144 :     $fig->reload_table('all', "comp_cas",
145 :     "cid varchar(7), cas varchar(100)",
146 :     { index_comp_cas_cid => "cid",
147 :     index_comp_cas_cas => "cas" },
148 :     "$FIG_Config::temp/comp_cas.table");
149 :     unlink("$FIG_Config::temp/comp_cas.table");
150 :     }
151 :    
152 :     sub load_reactions {
153 :    
154 : dejongh 1.4 my($react,$path,$sub,$prod,@sub,@prod,$subs,$prods,$dir);
155 : dejongh 1.1 my($cid,$n,$main,%reaction,$x);
156 :    
157 :     Trace("Loading reactions.") if T(2);
158 :     Open(\*REACTION, "<$FIG_Config::data/KEGG/reaction.lst");
159 : dejongh 1.4 Open(\*RMAPFORMULA, "<$FIG_Config::data/KEGG/reaction_mapformula.lst");
160 : dejongh 1.1 Open(\*R2C, ">$FIG_Config::temp/reaction_to_compound.table");
161 :     Open(\*REV, ">$FIG_Config::temp/rev.table");
162 :    
163 :     Trace("Reading reaction list file.") if T(2);
164 :     while (defined($_ = <REACTION>))
165 :     {
166 :     if ($_ =~ /(R\d+):\s+(\S.*\S)\s+<=>\s+(\S.*\S)/)
167 :     {
168 :     $react = $1;
169 :     $sub = $2;
170 :     $prod = $3;
171 :     @sub = split(/\s+\+\s+/,$sub);
172 :     @prod = split(/\s+\+\s+/,$prod);
173 :     @sub = map { $_ =~ /^(([\dmn\(]\S*)\s+)?([CG]\d+)/; $2 ? [$3,$2,0] : [$3,1,0] } @sub;
174 :     @prod = map { $_ =~ /^(([\dmn\(]\S*)\s+)?([CG]\d+)/; $2 ? [$3,$2,0] : [$3,1,0] } @prod;
175 :     $reaction{$react} = [[@sub],[@prod]];
176 :     }
177 :     else
178 :     {
179 :     Trace("Invalid reaction format: $_") if T(1);
180 :     }
181 :     }
182 :     close(REACTION);
183 :    
184 :     Trace("Reading main reaction file.") if T(2);
185 : dejongh 1.4
186 :     my %reversibility;
187 :    
188 :     while (defined($_ = <RMAPFORMULA>))
189 : dejongh 1.1 {
190 : dejongh 1.4 if ($_ =~ /^(R\d+):\s+(\d+):\s+(\S.*\S)\s(\<?=\>?)\s(\S.*\S)/)
191 : dejongh 1.1 {
192 :     $react = $1;
193 : dejongh 1.4 $path = $2;
194 :     $sub = $3;
195 :     $dir = $4;
196 :     $prod = $5;
197 : dejongh 1.1
198 :     if (exists($reaction{$react}))
199 :     {
200 :     $subs = $reaction{$react}->[0];
201 :     $prods = $reaction{$react}->[1];
202 : dejongh 1.4 my $rc = &mark_main($sub,$subs);
203 :    
204 :     if ($rc == 0)
205 :     {
206 :     $rc = &mark_main($sub,$prods) && &mark_main($prod,$subs);
207 :     }
208 :     else
209 :     {
210 :     $rc = &mark_main($prod,$prods);
211 :     }
212 :    
213 :     if ($rc == 0)
214 :     {
215 :     print "Can't handle $_\n";
216 :     }
217 :    
218 :     # entry for the reaction in the context of the pathway
219 :     $reversibility{$react.".rn".$path} = $dir;
220 :    
221 :     # since there can be multiple entries per reaction, with different
222 :     # reversibility info (in the context of different pathways),
223 :     # reversible trumps non-reversible for the general reaction entry
224 : dejongh 1.1 if (($dir eq "<=") || ($dir eq "=>"))
225 :     {
226 : dejongh 1.4 if (! defined($reversibility{$react}))
227 :     {
228 :     $reversibility{$react} = $dir;
229 :     }
230 :     elsif ($reversibility{$react} ne "<=>" && $reversibility{$react} ne $dir)
231 :     {
232 :     $reversibility{$react} = "<=>";
233 :     }
234 : dejongh 1.1 }
235 : dejongh 1.4 else
236 : dejongh 1.1 {
237 : dejongh 1.4 $reversibility{$react} = "<=>";
238 : dejongh 1.1 }
239 :     }
240 :     }
241 :     }
242 : dejongh 1.4 close(RMAPFORMULA);
243 :    
244 :     foreach $react (keys %reversibility)
245 :     {
246 :     print REV "$react\t$reversibility{$react}\n";
247 :     }
248 :    
249 : dejongh 1.1 close(REV);
250 :    
251 :     Trace("Connecting reactions to compounds.") if T(2);
252 :     foreach $react (sort keys(%reaction))
253 :     {
254 :     ($subs,$prods) = @{$reaction{$react}};
255 :     foreach $x (@$subs)
256 :     {
257 :     ($cid,$n,$main) = @$x;
258 :     print R2C "$react\t0\t$cid\t$n\t$main\n";
259 :     }
260 :    
261 :     foreach $x (@$prods)
262 :     {
263 :     ($cid,$n,$main) = @$x;
264 :     print R2C "$react\t1\t$cid\t$n\t$main\n";
265 :     }
266 :     }
267 :     close(R2C);
268 :    
269 :     $fig->reload_table('all', "reaction_to_compound",
270 :     "rid varchar(8), setn char(1), cid varchar(8), stoich char(6), main char(1)",
271 :     { index_reaction_to_compound_rid => "rid",
272 :     index_reaction_to_compound_cid => "cid" },
273 :     "$FIG_Config::temp/reaction_to_compound.table");
274 :     unlink("$FIG_Config::temp/reaction_to_compound.table");
275 :    
276 :     $fig->reload_table('all', "reversible",
277 : dejongh 1.4 "rid varchar(16) UNIQUE NOT NULL, reversible char(3), primary key(rid)",
278 : dejongh 1.1 { }, "$FIG_Config::temp/rev.table");
279 :     unlink("$FIG_Config::temp/rev.table");
280 :     Trace("Reactions processed.") if T(2);
281 :     }
282 :    
283 :     sub mark_main {
284 :     my($main,$set) = @_;
285 :     my($cid,$i);
286 :    
287 :     foreach $cid (split(/\s+\+\s+/,$main))
288 :     {
289 :     for ($i=0; ($i < @$set) && ($set->[$i]->[0] ne $cid); $i++) {}
290 :     if ($i == @$set)
291 :     {
292 : dejongh 1.4 # compound id was not found in the reaction
293 :     return 0;
294 : dejongh 1.1 }
295 :     else
296 :     {
297 :     $set->[$i]->[2] = 1;
298 :     }
299 :     }
300 : dejongh 1.4
301 :     return 1;
302 : dejongh 1.1 }
303 :    
304 :     sub load_catalyzes {
305 :    
306 :     my($entry);
307 :     Open(\*REAC, "<$FIG_Config::data/KEGG/reaction");
308 :     Open(\*REAC2ENZ, ">$FIG_Config::temp/reaction_to_enzyme.table");
309 :    
310 :     Trace("Reading KEGG reaction file.") if T(2);
311 :     my($rid,$ec,@ecs,$ecs);
312 :     $/ = "\n///\n";
313 :     while (defined($entry = <REAC>))
314 :     {
315 :     if ($entry =~ /ENTRY\s+(R\d+).*\nENZYME\s+(\S[^a-zA-Z\/]+)/s)
316 :     {
317 :     $rid = $1;
318 :     $ecs = $2;
319 : dejongh 1.4
320 : dejongh 1.1 foreach $ec (split(/\s+/,$ecs))
321 :     {
322 :     print REAC2ENZ "$rid\t$ec\n";
323 :     }
324 :     }
325 :     }
326 :     $/ = "\n";
327 :     close(REAC);
328 :     close(REAC2ENZ);
329 :    
330 :     $fig->reload_table('all', "reaction_to_enzyme",
331 :     "rid varchar(8), role varchar(100)",
332 :     { index_reaction_to_enzyme_rid => "rid",
333 :     index_reaction_to_enzyme_role => "role" },
334 :     "$FIG_Config::temp/reaction_to_enzyme.table");
335 :     unlink("$FIG_Config::temp/reaction_to_enzyme.table");
336 :     Trace("Enzyme reactions loaded.") if T(2);
337 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3