[Bio] / FigKernelScripts / update_PSEED_subsystems.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/update_PSEED_subsystems.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 ########################################################################
2 : olson 1.2 use Time::HiRes 'gettimeofday';
3 :     use YAML::Any qw(DumpFile LoadFile);
4 : overbeek 1.1 use Data::Dumper;
5 : olson 1.3
6 :     =head1 NAME
7 :    
8 :     update_PSEED_subsystems
9 :    
10 :     =head1 SYNOPSIS
11 :    
12 :     update_PSEED_subsystems [-restart restart-file] reconstruction-directory
13 :    
14 :     =head1 DESCRIPTION
15 :    
16 :     Updates the PSEED subsystems using the reconstruction information in
17 :     C<reconstruction-directory>.
18 :    
19 :     Basic steps.
20 :    
21 :     Use &build_needed() to pull all subsystem information from the Sapling. This
22 :     build a data structure that has all the subsystems and roles
23 :     from the annotator SEED. It also loads the metabolic reconstruction
24 :     data computed earlier by L<compute_PSEED_reconstructions.pl.
25 :    
26 :     =cut
27 :    
28 :    
29 : overbeek 1.1 use FIG; ### This references the PSEED
30 :     my $fig = new FIG;
31 :     use Subsystem;
32 : olson 1.2 use Getopt::Long;
33 : overbeek 1.1
34 : olson 1.2 my $restart_file;
35 :     my $dry_run;
36 :     my $rc = GetOptions("dry-run" => \$dry_run,
37 :     "restart=s" => \$restart_file);
38 :    
39 :     ($rc && @ARGV == 1) or die "Usage: update_PSEED_subsystems [-restart restart-file] reconstruction-dir\n";
40 :    
41 :     my $recon_dir = shift;
42 :     -d $recon_dir or die "Reconstruction directory $recon_dir does not exist\n";
43 :     print "go: dryrun=$dry_run restart=$restart_file\n";
44 : overbeek 1.1 use SeedEnv; ### The servers go (more-or-less) to the ASEED
45 :     my $annO = ANNOserver->new();
46 : olson 1.2 my $sapO = SAPserver->new();
47 : overbeek 1.1
48 :     use strict;
49 : olson 1.2
50 :     my($vcs,$bindings,$rolesH,$roleAbbrH);
51 :    
52 :     if (defined($restart_file) && -s $restart_file)
53 :     {
54 :     print STDERR "Loading restart file...\n";
55 :     my $r = LoadFile($restart_file);
56 :     print STDERR "...done\n";
57 :     ($vcs, $bindings, $rolesH, $roleAbbrH) = @$r;
58 :     }
59 :     else
60 :     {
61 :     ($vcs,$bindings,$rolesH,$roleAbbrH) = &build_needed($sapO, $recon_dir);
62 :     if (defined($restart_file))
63 :     {
64 :     DumpFile($restart_file, [$vcs, $bindings, $rolesH, $roleAbbrH]);
65 :     }
66 :     }
67 : overbeek 1.1
68 :     my @subsysL = sort keys(%$vcs);
69 :     foreach my $sub (@subsysL)
70 :     {
71 : olson 1.2 print "\nProcess $sub\n";
72 : overbeek 1.1 my $roles = [sort keys(%{$rolesH->{$sub}})];
73 : olson 1.2 &sync_sub_and_roles($annO,$fig,$sub,$roles,$roleAbbrH);
74 : overbeek 1.1
75 :     my $altered = 0;
76 :     if (my $vcsH = $vcs->{$sub})
77 :     {
78 :     my $subO = new Subsystem($sub,$fig);
79 : olson 1.2 my @genomes_in = sort { $a cmp $b } $subO->get_genomes;
80 :     my @genomes_needed = sort { $a cmp $b } keys(%$vcsH);
81 : overbeek 1.1 my $i = 0;
82 :     my $n = 0;
83 :     while (($i < @genomes_in) || ($n < @genomes_needed))
84 :     {
85 : olson 1.2 if (($i < @genomes_in) && (($n == @genomes_needed) || ($genomes_in[$i] lt $genomes_needed[$n])))
86 : overbeek 1.1 {
87 : olson 1.2 print " Remove $genomes_in[$i]\n";
88 :     $subO->remove_genome($genomes_in[$i]) if !$dry_run;
89 : overbeek 1.1 $i++;
90 :     $altered = 1;
91 :     }
92 : olson 1.2 elsif (($n < @genomes_needed) && ($i == @genomes_in) || ($genomes_in[$i] gt $genomes_needed[$n]))
93 : overbeek 1.1 {
94 : olson 1.2 my $idx = $subO->add_genome($genomes_needed[$n]) if !$dry_run;
95 :     print " Add $genomes_needed[$n]\n";
96 :     $subO->set_variant_code($idx,$vcsH->{$genomes_needed[$n]}) if !$dry_run;
97 : overbeek 1.1 $altered = 1;
98 :     foreach my $role (@$roles)
99 :     {
100 : olson 1.2 if (my @pegs = keys(%{$bindings->{$genomes_needed[$n]}->{$role}}))
101 : overbeek 1.1 {
102 : olson 1.2 print " $role => @pegs\n";
103 :     $subO->set_pegs_in_cell($genomes_needed[$n],$role,[sort { &FIG::by_fig_id($a,$b) } @pegs]) if !$dry_run;
104 : overbeek 1.1 }
105 :     }
106 :     $n++;
107 :     }
108 : olson 1.2 elsif (($i < @genomes_in) && ($n < @genomes_needed) && ($genomes_in[$i] eq $genomes_needed[$n]))
109 : overbeek 1.1 {
110 : olson 1.2 print " Update $genomes_in[$i]\n";
111 :     my $vc = $subO->get_variant_code_for_genome($genomes_in[$i]) if !$dry_run;
112 : overbeek 1.1 if ($vc ne $vcsH->{$genomes_in[$i]})
113 :     {
114 :     my $idx = $subO->get_genome_index($genomes_in[$i]);
115 : olson 1.2
116 :     print " change vc from $vc to $vcsH->{$genomes_needed[$n]}\n";
117 :     $subO->set_variant_code($idx,$vcsH->{$genomes_needed[$n]}) if !$dry_run;
118 : overbeek 1.1 $altered = 1;
119 :     }
120 :    
121 :     foreach my $role (@$roles)
122 :     {
123 :     my @pegs_in = sort { &SeedUtils::by_fig_id($a,$b) }
124 :     $subO->get_pegs_from_cell($genomes_in[$i],$role);
125 :     my @pegs_needed = sort { &SeedUtils::by_fig_id($a,$b) }
126 :     keys(%{$bindings->{$genomes_needed[$n]}->{$role}});
127 :     if (@pegs_in != @pegs_needed)
128 :     {
129 :     $altered = 1;
130 : olson 1.2 print " $role => @pegs_needed\n";
131 :     $subO->set_pegs_in_cell($genomes_in[$i],$role,\@pegs_needed) if !$dry_run;
132 : overbeek 1.1 }
133 :     else
134 :     {
135 :     for (my $i=0; ($i < @pegs_needed) && ($pegs_needed[$i] eq $pegs_in[$i]); $i++) {}
136 :     if ($i < @pegs_needed)
137 :     {
138 :     $altered = 1;
139 : olson 1.2 print " $role => @pegs_needed\n";
140 :     $subO->set_pegs_in_cell($genomes_needed[$n],$role,\@pegs_needed) if !$dry_run;
141 : overbeek 1.1 }
142 :     }
143 :    
144 :     }
145 :     $n++;
146 :     $i++;
147 :     }
148 :     }
149 :    
150 :     if ($altered)
151 :     {
152 : olson 1.2 if (!$dry_run)
153 :     {
154 :     $subO->write_subsystem;
155 :     }
156 : overbeek 1.1 }
157 :     }
158 :     }
159 :    
160 :     sub build_needed {
161 : olson 1.2 my($sapO, $recon_dir) = @_;
162 : overbeek 1.1
163 : olson 1.2 my @files = <$recon_dir/*/*>;
164 :     my @genome_files = sort { $a->[0] cmp $b->[0] } map { m,/(\d+\.\d+)$, ; [$1, $_] } @files;
165 :    
166 :     #
167 :     # Load the subsystem role/abbreviation information from sapling.
168 :     #
169 :    
170 :     my $roleInfo = $sapO->query(-limit => 999999999,
171 :     -objects => 'Subsystem Includes Role',
172 :     -fields => [qw(Subsystem(id) Includes(abbreviation) Role(id))]);
173 :    
174 :     my $role_map = {};
175 :    
176 :     map { $role_map->{$_->[0], $_->[2]} = $_->[1] } @$roleInfo;
177 :    
178 : overbeek 1.1 my $vcs = {};
179 :     my $bindings = {};
180 :     my $rolesH = {};
181 : olson 1.2 foreach my $genome_file (@genome_files)
182 : overbeek 1.1 {
183 : olson 1.2 my($genome, $recon_file) = @$genome_file;
184 :    
185 :     my $meta_recon = LoadFile($recon_file);
186 :    
187 : overbeek 1.1 foreach my $tuple (@$meta_recon)
188 :     {
189 :     my($subvc,$role,$peg) = @$tuple;
190 :     $bindings->{$genome}->{$role}->{$peg} = 1;
191 :     my($sub,$vc) = ($subvc =~ /^(\S.*\S)\s*:(\S+)$/);
192 :     defined($sub) || die "sub=$sub vc=$vc subvc=$subvc";
193 :     $vcs->{$sub}->{$genome} = $vc;
194 :     $rolesH->{$sub}->{$role} = 1;
195 :     }
196 :     }
197 : olson 1.2 return ($vcs,$bindings,$rolesH, $role_map);
198 : overbeek 1.1 }
199 :    
200 :     sub sync_sub_and_roles {
201 : olson 1.2 my($annO,$fig,$sub,$roles,$roleAbbrH) = @_;
202 : overbeek 1.1
203 :     my $subO = new Subsystem($sub,$fig);
204 :     my $must_fix = 0;
205 :     if ($subO)
206 :     {
207 :     my @roles_in = sort $subO->get_roles;
208 :     my @roles_needed = sort @$roles;
209 :     if (@roles_in != @roles_needed)
210 :     {
211 :     $must_fix = 1;
212 :     }
213 :     else
214 :     {
215 :     my $i;
216 :     for ($i=0; ($i < @roles_in) && ($roles_in[$i] eq $roles_needed[$i]); $i++) {}
217 :     if ($i < @roles_in) { $must_fix = 1 }
218 :     }
219 :     }
220 :     else
221 :     {
222 :     $must_fix = 1;
223 : olson 1.2 $subO = new Subsystem($sub,$fig,1) if !$dry_run; ## create new subsystem
224 : overbeek 1.1 }
225 :    
226 :     if ($must_fix)
227 :     {
228 : olson 1.2 my $roles_abbrs = [map { [$_, $roleAbbrH->{$sub, $_}] } @$roles];
229 :     print "Set $sub @$roles\n";
230 :     $subO->set_roles($roles_abbrs) if !$dry_run;
231 :     if (!$dry_run)
232 :     {
233 :     $subO->write_subsystem;
234 :     }
235 : overbeek 1.1 }
236 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3