[Bio] / FigKernelPackages / UnvSubsys.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/UnvSubsys.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 package UnvSubsys;
2 :    
3 : overbeek 1.2 use Subsystem;
4 : overbeek 1.1 use Carp;
5 :     use FIG;
6 :    
7 :     use Data::Dumper;
8 :     use strict;
9 :    
10 :     sub new
11 :     {
12 : overbeek 1.8 my($class, $ssa, $fig, $active_subsetR, $focus, $show_clusters, $aliases, $peg_colors) = @_;
13 : overbeek 1.1
14 :     $ssa =~ s/ /_/g;
15 :    
16 :    
17 : parrello 1.5 ### { Roles =>Roles,
18 : overbeek 1.2 ### RoleIndex => ToRoleIndexHash,
19 :     ### RoleSubsets => ColSubsets,
20 :     ### Genomes => Genomes,
21 :     ### GenomeIndex => ToGenomeIndexHash,
22 :     ### PegHash => PegHash,
23 :     ### Colors => ColorHash,
24 :     ### Aliases => AliasHash,
25 :     ### Curator => Curator,
26 :     ### Notes => Notes,
27 :     ### Reactions => ReactionHash
28 :     ### }
29 : parrello 1.5 ###
30 :     ### Roles = pointer to a list of [Role,Abbrev,[ReactionURLs]]
31 :     ###
32 :     ### ToRoleIndexHash = a pointer to a hash: key=Role Value=RoleIndex
33 :     ###
34 :     ### ColSubsets = pointer to a list of [SubsetName,[RoleIndexesFrom0]]
35 :     ###
36 : overbeek 1.6 ### RowSubSets = pointer to a list of [SubsetName,[GenomeIndexesFrom0]]
37 :     ###
38 : parrello 1.5 ### Genomes is a pointer to a list of [Genome,Variant]
39 :     ###
40 :     ### ToGenomeIndexHash = a pointer to a hash: key=Genome value=GenomeIndex
41 :     ###
42 :     ### PegHash = a pointer to a hash of hashes such that $peg_hash->{$genome_index}->{$role_index} = a
43 :     ### pointer to a list of PEGs
44 :     ###
45 :     ### ColorHash is a hash: key=PEG value=color
46 :     ###
47 :     ### AliasHash is a hash: key=PEG value=aliases
48 :     ###
49 : overbeek 1.2 ### ReactionHash is a hash: key=Role value=[reaction-ids]
50 : overbeek 1.1
51 : overbeek 1.6 if ((ref($fig) eq "FIG") || ((ref($fig) eq 'SFXlate') && ($fig = $fig->{'fig'})))
52 : overbeek 1.1 {
53 : parrello 1.5 my $subsystem = new Subsystem($ssa,$fig,0);
54 :     my $curator = $subsystem->get_curator;
55 :     my $notes = $subsystem->get_notes;
56 :     $notes =~ s/ /\n/g;
57 :     my @roles = $subsystem->get_roles;
58 :     my $reactions = $subsystem->get_reactions;
59 :     my @genomes = $subsystem->get_genomes;
60 : overbeek 1.1 my @col_subsets = $subsystem->get_subset_namesC;
61 :    
62 : parrello 1.5 my $role_info = [];
63 :     my $roleH = {};
64 : overbeek 1.1
65 : parrello 1.5 my($i,$j,$subset,$peg);
66 :     for ($i=0; ($i < @roles); $i++)
67 :     {
68 :     my $role = $roles[$i];
69 :     my $abbrev = $subsystem->get_role_abbr( $subsystem->get_role_index( $role ) );
70 :     my $react = $reactions ? join(",", map { &HTML::reaction_link($_) } @{$reactions->{$role}}) : [];
71 :     push(@$role_info,[$role,$abbrev,$react]);
72 :     $roleH->{$role} = $i;
73 :     }
74 :    
75 :     my $subset_info = [];
76 :     foreach $subset (@col_subsets)
77 :     {
78 :     if ($subset ne 'All')
79 :     {
80 :     push(@$subset_info,[$subset,[map { $roleH->{$_} } $subsystem->get_subsetC_roles($subset)]]);
81 :     }
82 :     }
83 :    
84 :     my $genomes_info = [];
85 :     my $genomeH = {};
86 :     for ($i=0; ($i < @genomes); $i++)
87 :     {
88 :     my $genome = $genomes[$i];
89 :     my $variant = $subsystem->get_variant_code( $subsystem->get_genome_index( $genome ) );
90 :     push(@$genomes_info,[$genome,$variant]);
91 :     $genomeH->{$genome} = $i;
92 :     }
93 :    
94 :     my $pegH = {};
95 :     for ($i=0; ($i < @genomes); $i++)
96 :     {
97 :     my $genome = $genomes[$i];
98 :     for ($j=0; ($j < @roles); $j++)
99 :     {
100 :     my $role = $roles[$j];
101 :     my @pegs = $subsystem->get_pegs_from_cell($genome,$role);
102 :     $pegH->{$i}->{$j} = [@pegs];
103 :     }
104 :     }
105 : overbeek 1.7
106 :     my $row_subsets = &row_subsets($fig,$genomeH,$genomes_info);
107 : overbeek 1.11 my $active_genomes = &active_genomes($fig,$row_subsets,$active_subsetR,$focus,$genomeH,$genomes_info);
108 :    
109 : overbeek 1.8 my $colorsH;
110 :     if ($peg_colors)
111 :     {
112 :     $colorsH = {};
113 :     foreach $_ (@$peg_colors)
114 :     {
115 :     my($peg,$color) = @$_;
116 :     $colorsH->{$peg} = $color;
117 :     }
118 :     }
119 :     elsif ($show_clusters)
120 :     {
121 :     $colorsH = &set_colors($fig,$pegH,$active_genomes);
122 :     }
123 :     else
124 :     {
125 :     $colorsH = {};
126 :     }
127 :    
128 : overbeek 1.7 my $aliasesH = $aliases ? &set_aliases($fig,$pegH,$active_genomes) : {};
129 : overbeek 1.6 my $reactions = $subsystem->get_reactions;
130 : parrello 1.5 my $self = { Roles => $role_info,
131 :     RoleIndex => $roleH,
132 :     RoleSubsets => $subset_info,
133 :     Genomes => $genomes_info,
134 :     GenomeIndex => $genomeH,
135 : overbeek 1.6 GenomeSubsets => $row_subsets,
136 : parrello 1.5 PegHash => $pegH,
137 :     Colors => $colorsH,
138 :     Aliases => $aliasesH,
139 :     Curator => $curator,
140 :     Notes => $notes,
141 :     Reactions => $reactions
142 :     };
143 :     bless($self, $class);
144 :     return $self;
145 : overbeek 1.1 }
146 :     else
147 :     {
148 : parrello 1.5 return undef;
149 : overbeek 1.1 }
150 :     }
151 :    
152 : overbeek 1.6 sub get_subset_namesR {
153 :     my($self) = @_;
154 :    
155 :     return map { $_->[0] } @{$self->{GenomeSubsets}};
156 :     }
157 :    
158 :     sub get_subsetR {
159 :     my($self,$set) = @_;
160 :     my($i);
161 :    
162 :     my $sets = $self->{GenomeSubsets};
163 :     for ($i=0; ($i < @$sets) && ($sets->[$i]->[0] ne $set); $i++) {}
164 :     if ($i < @$sets)
165 :     {
166 :     return map { $self->{Genomes}->[$_]->[0] } @{$sets->[$i]->[1]}
167 :     }
168 :     return undef;
169 :     }
170 :    
171 :     sub get_subsetsR {
172 :     my($self) = @_;
173 :    
174 :     my $sets = $self->{GenomeSubsets};
175 :     my @pairs = ();
176 :     my $pair;
177 :     foreach $pair (@$sets)
178 :     {
179 :     my($id,$members) = @$pair;
180 :     push(@pairs,[$id,[map { $self->{Genomes}->[$_]->[0] } @$members]]);
181 :     }
182 :     return @pairs;
183 :     }
184 :    
185 :     sub row_subsets {
186 :     my ($fig,$genomeH,$genomes_info) = @_;
187 :    
188 :     my $subsets = [];
189 : overbeek 1.11 my $taxonomic_groups = $fig->taxonomic_groups_of_complete(5);
190 : overbeek 1.9
191 : overbeek 1.6 my($pair,$id,$members);
192 :     foreach $pair (@$taxonomic_groups)
193 :     {
194 :     ($id,$members) = @$pair;
195 :     my @mem = grep { defined($_) } map { $genomeH->{$_} } @$members;
196 : overbeek 1.10 if (@mem > 0)
197 : overbeek 1.6 {
198 :     push(@$subsets,[$id,[@mem]]);
199 :     }
200 :     }
201 :     return $subsets;
202 :     }
203 :    
204 : overbeek 1.4 sub set_aliases {
205 : overbeek 1.7 my($fig,$pegH,$active_genomes) = @_;
206 : overbeek 1.4 my($genomeI,$roleI,$pegs,$peg,$roleH);
207 :    
208 :     my $aliasesH = {};
209 :    
210 : overbeek 1.7 foreach $genomeI (grep { $active_genomes->{$_} } keys(%$pegH))
211 : overbeek 1.4 {
212 : parrello 1.5 $roleH = $pegH->{$genomeI};
213 :     foreach $roleI (keys(%$roleH))
214 :     {
215 :     $pegs = $roleH->{$roleI};
216 :     foreach $peg (@$pegs)
217 :     {
218 :     if (! $aliasesH->{$peg})
219 :     {
220 :     $aliasesH->{$peg} = scalar &ext_id($fig,$peg);
221 :     }
222 :     }
223 :     }
224 : overbeek 1.4 }
225 :     return $aliasesH;
226 :     }
227 :    
228 :     sub set_colors {
229 : overbeek 1.7 my($fig,$pegH,$active_genomes) = @_;
230 : overbeek 1.4 my($genomeI,$roleI,$pegs,$peg,$roleH,$peg,%pegs_in_genome);
231 :    
232 :     my $colorsH = {};
233 : overbeek 1.7 foreach $genomeI (grep { $active_genomes->{$_} } keys(%$pegH))
234 : overbeek 1.4 {
235 : parrello 1.5 undef %pegs_in_genome;
236 :     $roleH = $pegH->{$genomeI};
237 :     foreach $roleI (keys(%$roleH))
238 :     {
239 :     $pegs = $roleH->{$roleI};
240 :     foreach $peg (@$pegs)
241 :     {
242 :     $pegs_in_genome{$peg} = 1;
243 :     }
244 :     }
245 :    
246 :     my @pegs = keys(%pegs_in_genome);
247 :     my($tuple,$peg,$color);
248 :     my $colors_for_one_genome = &set_colors_for_genome($fig,\@pegs);
249 :     while (($peg,$color) = each %$colors_for_one_genome)
250 :     {
251 :     $colorsH->{$peg} = $colors_for_one_genome->{$peg};
252 :     }
253 : overbeek 1.4 }
254 :     return $colorsH;
255 :     }
256 :    
257 :     sub set_colors_for_genome {
258 :     my($fig,$pegs) = @_;
259 :     my($peg,@clusters,$cluster,@colors,$color,%seen,%conn,$x,$peg1,@pegs,$i);
260 :    
261 :     my $color_of = {};
262 :     foreach $peg (@$pegs) { $color_of->{$peg} = '#FFFFFF' }
263 :    
264 :     @pegs = keys(%$color_of); # Use of keys makes @pegs entries unique
265 :    
266 :     foreach $peg (@pegs)
267 :     {
268 : parrello 1.5 $conn{$peg} = [grep { $color_of->{$_} && ($_ ne $peg) } $fig->close_genes($peg,5000)];
269 : overbeek 1.4 }
270 :    
271 :     @clusters = ();
272 :     while ($peg = shift @pegs)
273 :     {
274 : parrello 1.5 if (! $seen{$peg})
275 :     {
276 :     $cluster = [$peg];
277 :     $seen{$peg} = 1;
278 :     for ($i=0; ($i < @$cluster); $i++)
279 :     {
280 :     my @tmp = grep { ! $seen{$_} } @{$conn{$cluster->[$i]}};
281 :     if (@tmp > 0)
282 :     {
283 :     foreach my $peg1 (@tmp) { $seen{$peg1} = 1 }
284 :     push(@$cluster,@tmp);
285 :     }
286 :     }
287 :     push(@clusters,$cluster);
288 :     }
289 : overbeek 1.4 }
290 :    
291 :     @colors = &cool_colors();
292 :    
293 :     @clusters = grep { @$_ > 1 } sort { @$a <=> @$b } @clusters;
294 :    
295 :     if (@clusters > @colors) { splice(@clusters,0,(@clusters - @colors)) } # make sure we have enough colors
296 :    
297 :     my($cluster);
298 :     foreach $cluster (@clusters)
299 :     {
300 : parrello 1.5 $color = shift @colors;
301 :     foreach $peg (@$cluster)
302 :     {
303 :     $color_of->{$peg} = $color;
304 :     }
305 : overbeek 1.4 }
306 :     return $color_of;
307 :     }
308 :    
309 :     sub cool_colors {
310 :     # just an array of "websafe" colors or whatever colors we want to use. Feel free to remove bad colors (hence the lines not being equal length!)
311 :     return (
312 :     '#C0C0C0', '#FF40C0', '#FF8040', '#FF0080', '#FFC040', '#40C0FF', '#40FFC0', '#C08080', '#C0FF00', '#00FF80', '#00C040',
313 :     "#6B8E23", "#483D8B", "#2E8B57", "#008000", "#006400", "#800000", "#00FF00", "#7FFFD4",
314 :     "#87CEEB", "#A9A9A9", "#90EE90", "#D2B48C", "#8DBC8F", "#D2691E", "#87CEFA", "#E9967A", "#FFE4C4", "#FFB6C1",
315 :     "#E0FFFF", "#FFA07A", "#DB7093", "#9370DB", "#008B8B", "#FFDEAD", "#DA70D6", "#DCDCDC", "#FF00FF", "#6A5ACD",
316 :     "#00FA9A", "#228B22", "#1E90FF", "#FA8072", "#CD853F", "#DC143C", "#FF6347", "#98FB98", "#4682B4",
317 :     "#D3D3D3", "#7B68EE", "#2F4F4F", "#FF7F50", "#FF69B4", "#BC8F8F", "#A0522D", "#DEB887", "#00DED1",
318 :     "#6495ED", "#800080", "#FFD700", "#F5DEB3", "#66CDAA", "#FF4500", "#4B0082", "#CD5C5C",
319 :     "#EE82EE", "#7CFC00", "#FFFF00", "#191970", "#FFFFE0", "#DDA0DD", "#00BFFF", "#DAA520", "#008080",
320 :     "#00FF7F", "#9400D3", "#BA55D3", "#D8BFD8", "#8B4513", "#3CB371", "#00008B", "#5F9EA0",
321 :     "#4169E1", "#20B2AA", "#8A2BE2", "#ADFF2F", "#556B2F",
322 :     "#F0FFFF", "#B0E0E6", "#FF1493", "#B8860B", "#FF0000", "#F08080", "#7FFF00", "#8B0000",
323 :     "#40E0D0", "#0000CD", "#48D1CC", "#8B008B", "#696969", "#AFEEEE", "#FF8C00", "#EEE8AA", "#A52A2A",
324 :     "#FFE4B5", "#B0C4DE", "#FAF0E6", "#9ACD32", "#B22222", "#FAFAD2", "#808080", "#0000FF",
325 :     "#000080", "#32CD32", "#FFFACD", "#9932CC", "#FFA500", "#F0E68C", "#E6E6FA", "#F4A460", "#C71585",
326 :     "#BDB76B", "#00FFFF", "#FFDAB9", "#ADD8E6", "#778899",
327 :     );
328 :     }
329 :    
330 :     sub ext_id {
331 :     my($fig,$peg) = @_;
332 :    
333 :     my @tmp;
334 :     my @aliases = $fig->feature_aliases($peg);
335 :     if ((@tmp = grep { $_ =~ /^uni\|/ } @aliases) > 0)
336 :     {
337 : parrello 1.5 @aliases = @tmp;
338 : overbeek 1.4 }
339 :     elsif ((@tmp = grep { $_ =~ /^sp\|/ } @aliases) > 0)
340 :     {
341 : parrello 1.5 @aliases = @tmp;
342 : overbeek 1.4 }
343 :     elsif ((@tmp = grep { $_ =~ /^gi\|/ } @aliases) > 0)
344 :     {
345 : parrello 1.5 @aliases = @tmp;
346 : overbeek 1.4 }
347 :     elsif ((@tmp = grep { $_ =~ /^kegg\|/ } @aliases) > 0)
348 :     {
349 : parrello 1.5 @aliases = @tmp;
350 : overbeek 1.4 }
351 :     else
352 :     {
353 : parrello 1.5 @aliases = ();
354 : overbeek 1.4 }
355 :    
356 :     if (wantarray())
357 :     {
358 : parrello 1.5 return @aliases;
359 : overbeek 1.4 }
360 :     else
361 :     {
362 : parrello 1.5 return $aliases[0];
363 : overbeek 1.4 }
364 :     }
365 :    
366 :    
367 : overbeek 1.2 sub subsystem_curator {
368 :     my($self) = @_;
369 :    
370 :     my $curator = $self->{Curator};
371 :     $curator =~ s/master://;
372 :     return $curator;
373 :     }
374 :    
375 :     sub get_roles {
376 :     my($self) = @_;
377 :    
378 :     return map { $_->[0] } @{$self->{Roles}};
379 :     }
380 :    
381 : overbeek 1.3 sub get_genome_index {
382 :     my($self,$genome) = @_;
383 :    
384 :     return $self->{GenomeIndex}->{$genome};
385 :     }
386 :    
387 :     sub get_genomes {
388 :     my($self) = @_;
389 :    
390 :     return map { $_->[0] } @{$self->{Genomes}};
391 :     }
392 :    
393 :     sub get_variant_code {
394 :     my($self,$genome) = @_;
395 :    
396 :     if ($genome =~ /^\d+$/)
397 :     {
398 : parrello 1.5 return $self->{Genomes}->[$genome]->[1];
399 : overbeek 1.3 }
400 :     else
401 :     {
402 : parrello 1.5 my $genomeI = $self->{GenomeIndex}->{$genome};
403 :     return $self->{Genomes}->[$genomeI]->[1];
404 : overbeek 1.3 }
405 :     }
406 :    
407 :     sub get_pegs_from_cell {
408 :     my($self,$genome,$role) = @_;
409 :    
410 :     my $genomeI = $self->{GenomeIndex}->{$genome};
411 :     my $roleI = $self->{RoleIndex}->{$role};
412 :    
413 :     my $pegs = $self->{PegHash}->{$genomeI}->{$roleI};
414 :     return $pegs ? @$pegs : ();
415 :     }
416 :    
417 :     sub get_notes {
418 :     my($self) = @_;
419 :    
420 :     return $self->{Notes};
421 :     }
422 :    
423 : overbeek 1.2 sub get_role_index {
424 :     my($self,$role) = @_;
425 :    
426 :     return $self->{RoleIndex}->{$role};
427 :     }
428 :    
429 :     sub get_role_abbr {
430 :     my($self,$roleI) = @_;
431 :    
432 :     if ($roleI !~ /^\d+$/)
433 :     {
434 : parrello 1.5 $roleI = $self->{RoleIndex}->{$roleI};
435 : overbeek 1.2 }
436 :     my $roles = $self->{Roles};
437 :     return $roles->[$roleI]->[1];
438 :     }
439 :    
440 :     sub get_reactions {
441 :     my($self) = @_;
442 :    
443 :     return $self->{Reactions};
444 :     }
445 :    
446 :     sub get_subset_namesC {
447 :     my($self) = @_;
448 :    
449 :     return map { $_->[0] } @{$self->{RoleSubsets}};
450 :     }
451 :    
452 :     sub get_subsetC_roles {
453 :     my($self,$subset) = @_;
454 :     my($i,$j);
455 :    
456 :     my $subset_info = $self->{RoleSubsets};
457 :     for ($i=0; ($i < @$subset_info) && ($subset_info->[$i]->[0] ne $subset); $i++) {}
458 :     if ($i < @$subset_info)
459 :     {
460 : parrello 1.5 my @roles = ();
461 :     foreach $j (@{$subset_info->[$i]->[1]})
462 :     {
463 :     push(@roles,$self->{Roles}->[$j]->[0]);
464 :     }
465 :     return @roles;
466 : overbeek 1.2 }
467 :     return undef;
468 :     }
469 :    
470 : overbeek 1.4 sub get_color_of {
471 :     my($self,$peg) = @_;
472 :    
473 :     return $self->{Colors}->{$peg};
474 :     }
475 :    
476 : overbeek 1.7 sub active_genomes {
477 : overbeek 1.11 my($fig,$row_subsets,$active_subsetR,$focus,$genomeH,$genomes_info) = @_;
478 : overbeek 1.7 my($i,@bestL);
479 :    
480 :     my $active_genomes = {};
481 : overbeek 1.11
482 : overbeek 1.7 if ($active_subsetR)
483 :     {
484 :     for ($i=0; ($i < @$row_subsets) && ($row_subsets->[$i]->[0] ne $active_subsetR); $i++) {}
485 :     if ($i < @$row_subsets)
486 :     {
487 :     @bestL = @{$row_subsets->[$i]->[1]};
488 :     }
489 :     else
490 :     {
491 :     $active_subsetR = 'All';
492 : overbeek 1.10 @bestL = map { $genomeH->{$_} } keys(%$genomeH);
493 : overbeek 1.7 }
494 :     }
495 :     else
496 :     {
497 : overbeek 1.11 if ((! $focus) || (! $fig->is_complete($focus)))
498 : overbeek 1.7 {
499 :     $active_subsetR = 'All';
500 : overbeek 1.10 @bestL = map { $genomeH->{$_} } keys(%$genomeH);
501 : overbeek 1.7 }
502 :     else
503 :     {
504 :     my $bestN = undef;
505 :     @bestL = ();
506 :     my $tuple;
507 :     foreach $tuple (@$row_subsets)
508 :     {
509 :     my($id,$genomeIs) = @$tuple;
510 :     for ($i=0; ($i < @$genomeIs) && ($genomes_info->[$genomeIs->[$i]]->[0] ne $focus); $i++) {}
511 :     if ($i < @$genomeIs)
512 :     {
513 :     if ((! $bestN) || (@$genomeIs < @bestL))
514 :     {
515 :     $bestN = $id;
516 :     @bestL = @$genomeIs;
517 :     }
518 :     }
519 : overbeek 1.11 else
520 :     {
521 :     # print &Dumper($id);
522 :     }
523 : overbeek 1.7 }
524 :     $active_subsetR = $bestN;
525 :     }
526 :     }
527 :    
528 :     my %active_genomes = map { $_ => 1 } @bestL;
529 :     return \%active_genomes;
530 :     }
531 :    
532 : overbeek 1.1 1;
533 :    
534 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3