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

Annotation of /FigKernelScripts/wb.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 ########################################################################
2 :     # -*- perl -*-
3 :     ########################################################################
4 :     # Copyright (c) 2003-2008 University of Chicago and Fellowship
5 :     # for Interpretations of Genomes. All Rights Reserved.
6 :     #
7 :     # This file is part of the SEED Toolkit.
8 :     #
9 :     # The SEED Toolkit is free software. You can redistribute
10 :     # it and/or modify it under the terms of the SEED Toolkit
11 :     # Public License.
12 :     #
13 :     # You should have received a copy of the SEED Toolkit Public License
14 :     # along with this program; if not write to the University of Chicago
15 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
16 :     # Genomes at veronika@thefig.info or download a copy from
17 :     # http://www.theseed.org/LICENSE.TXT.
18 :     ########################################################################
19 :    
20 :     use Carp;
21 :     use Data::Dumper;
22 :     use Time::HiRes qw(gettimeofday);
23 :     use Time::Local;
24 :     use gjoseqlib;
25 :     use gjoalignment;
26 :     use representative_sequences;
27 :     use gjoalignandtree;
28 :    
29 :     use SeedEnv;
30 :     my $sapO = SAPserver->new;
31 :    
32 :     my($user);
33 :    
34 :     # usage: wb [-echo] [-time] [-ATdir ATdir] [command]
35 :    
36 :     $echo = 0;
37 :     $time_cmds = 0;
38 :     $ATdir = '/home/overbeek/Ross/GaryTrees/WB/ATdir';
39 :    
40 :     while ((@ARGV > 0) && ($ARGV[0] =~ /^-/))
41 :     {
42 :     $arg = shift @ARGV;
43 :     if ($arg =~ /^-time/i) { $time_cmds = 1 }
44 :     if ($arg =~ /^-echo/i) { $echo = 1 }
45 :     if ($arg =~ /^-ATdir/)
46 :     {
47 :     $ATdir = shift @ARGV;
48 :     }
49 :     }
50 :    
51 :     my($t1,$t2);
52 :     if (@ARGV > 0) { $req = join( " ", @ARGV ); }
53 :     while ( (defined($req) && $req) || ((@ARGV == 0) && ($req = &get_req)) )
54 :     {
55 :     if ($time_cmds)
56 :     {
57 :     $t1 = gettimeofday;
58 :     }
59 :    
60 :     if ($req =~ /^\s*h\s*$/ || $req =~ /^\s*help\s*$/)
61 :     {
62 :     &help;
63 :     }
64 :     elsif ($req =~ /^\s*show_subsys\s*$/)
65 :     {
66 :     &show_all_subsys($sapO);
67 :     }
68 :     elsif ($req =~ /^\s*show_FR_index\s*$/)
69 :     {
70 :     &show_FR_index($ATdir);
71 :     }
72 :     elsif ($req =~ /^\s*build_reps\s+(\d+)\s*$/)
73 :     {
74 :     &build_reps($sapO,$ATdir,$1);
75 :     }
76 :     elsif ($req =~ /^\s*summarize_FR\s+(\d+)\s*$/)
77 :     {
78 :     &summarize_FR($ATdir,$1);
79 :     }
80 :     elsif ($req =~ /^\s*get_FR\s+\'(\S[^\']*\S)\'\s+\'(.*)\'\s*$/)
81 :     {
82 :     my $subsys = $1;
83 :     my @roles = split(/\'\s+\'/,$2);
84 :     &get_FR($sapO,$ATdir,$subsys,\@roles);
85 :     }
86 :     elsif ($req =~ /^\s*align1\s+(\d+)\.(\d+)\s*$/)
87 :     {
88 :     my ($index, $subI) = ($1, $2);
89 :     &align1($sapO, $ATdir, $index, $subI);
90 :     }
91 :     elsif ($req =~ /^\s*trim1\s+(\d+)\.(\d+)\s*$/)
92 :     {
93 :     my ($index, $subI) = ($1, $2);
94 :     &trim1($sapO, $ATdir, $index, $subI);
95 :     }
96 :     elsif ($req =~ /^\s*process\s+\'(\S[^\']*\S)\'\s+\'(.*)\'\s*$/)
97 :     {
98 :     my $subsys = $1;
99 :     my @roles = split(/\'\s+\'/,$2);
100 :     my $index = &get_FR($sapO,$ATdir,$subsys,\@roles);
101 :     &build_reps($sapO,$ATdir,$index);
102 :     opendir(TMP,"$ATdir/$index") || die "could not open $ATdir/$index";
103 :     my @subIs = map { ($_ =~ /^\d+\.(\d+)$/) ? $1 : () } readdir(TMP);
104 :     closedir(TMP);
105 :     foreach my $subI (@subIs) {
106 :     &align1($sapO,$ATdir,$index,$subI);
107 :     &trim1($sapO,$ATdir,$index,$subI);
108 :     }
109 :     }
110 :     elsif ($req =~ /^\s*extend_FR\s+(\d+)\s*$/)
111 :     {
112 :     my $index = $1;
113 :     &extend_FR($sapO,$ATdir,$index);
114 :     }
115 :     elsif ($req =~ /^\s*show_roles_in_subsys\s+\'(\S.*\S)\'\s*$/)
116 :     {
117 :     my $subsys = $1;
118 :     my $ssH = $sapO->subsystem_roles( -ids => [$subsys] );
119 :     if (my $ss_entry = $ssH->{$subsys})
120 :     {
121 :     foreach my $role (@$ss_entry)
122 :     {
123 :     print "$role\n";
124 :     }
125 :     }
126 :     else
127 :     {
128 :     print "Invalid subsystem\n";
129 :     }
130 :     }
131 :     else
132 :     {
133 :     print "Invalid command\n";
134 :     }
135 :     print "\n";
136 :     $req = "";
137 :    
138 :     if ($time_cmds)
139 :     {
140 :     $t2 = gettimeofday;
141 :     print $t2-$t1," seconds to execute command\n\n";
142 :     }
143 :     }
144 :    
145 :     sub padded {
146 :     my($x,$n) = @_;
147 :    
148 :     if (length($x) < $n)
149 :     {
150 :     return $x . (" " x ($n - length($x)));
151 :     }
152 :     return $x;
153 :     }
154 :    
155 :     sub get_req {
156 :     my($x);
157 :    
158 :     print "?? ";
159 :     $x = <STDIN>;
160 :     while (defined($x) && ($x =~ /^h$/i) )
161 :     {
162 :     &help;
163 :     print "?? ";
164 :     $x = <STDIN>;
165 :     }
166 :    
167 :     if ((! defined($x)) || ($x =~ /^\s*[qQxX]/))
168 :     {
169 :     return "";
170 :     }
171 :     else
172 :     {
173 :     if ($echo)
174 :     {
175 :     print ">> $x\n";
176 :     }
177 :     return $x;
178 :     }
179 :     }
180 :    
181 :     sub show_all_subsys {
182 :     my($sapO) = @_;
183 :    
184 :     my $ssH = $sapO->all_subsystems();
185 :     my @ss = sort keys(%$ssH);
186 :     foreach my $subsys (@ss)
187 :     {
188 :     print "$subsys\n";
189 :     }
190 :     print "\n";
191 :     }
192 :    
193 :     sub build_fr_entry {
194 :     my($sapO,$ATdir,$subsys,$roles) = @_;
195 :    
196 :     my $ts = time;
197 :    
198 :     &SeedUtils::verify_dir($ATdir);
199 :     if (! -e "$ATdir/FR-index") { system "touch $ATdir/FR-index" }
200 :    
201 :     my @index = `cat $ATdir/FR-index`;
202 :     my $nxt = @index + 1;
203 :     my @peg_fasta_entries = &get_ss_entries($sapO,$subsys,$roles);
204 :     if (@peg_fasta_entries > 0)
205 :     {
206 :     open(INDEX,">>$ATdir/FR-index") || die "Something is amiss with $ATdir/FR-index";
207 :     print INDEX join("\t", $nxt, $ts, $subsys, @$roles),"\n";
208 :     close(INDEX);
209 :    
210 :     if (-e "$ATdir/$nxt") { system "rm -rf $ATdir/$nxt" }
211 :     mkdir("$ATdir/$nxt", 0777) || die "Could not make $ATdit/$nxt";
212 :     &gjoseqlib::print_alignment_as_fasta("$ATdir/$nxt/seqs",\@peg_fasta_entries);
213 :     return $nxt;
214 :     }
215 :     else
216 :     {
217 :     return 0;
218 :     }
219 :     }
220 :    
221 :     sub get_FR {
222 :     my($sapO,$ATdir,$subsys,$roles) = @_;
223 :    
224 :     my $index;
225 :     if ((@$roles > 0) && ($index = &build_fr_entry($sapO,$ATdir,$subsys,$roles)))
226 :     {
227 :     print "Built FR entry for index $index\n";
228 :     }
229 :     else
230 :     {
231 :     print "Failed to build entry for $subsys\n";
232 :     foreach my $role (@$roles)
233 :     {
234 :     print "\t$role\n";
235 :     }
236 :     }
237 :     return $index;
238 :     }
239 :    
240 :     sub get_ss_entries {
241 :     my($sapO,$subsys,$roles) = @_;
242 :    
243 :     my $roleH = $sapO->pegs_implementing_roles( -subsystem => $subsys, -roles => $roles );
244 :     my @pegs = map { @{$roleH->{$_}} } keys(%$roleH);
245 :     my $idH = $sapO->ids_to_sequences( -protein => 1, -ids => \@pegs );
246 :     return map { [$_,'',$idH->{$_}] } keys(%$idH);
247 :     }
248 :    
249 :     sub show_FR_index {
250 :     my($ATdir) = @_;
251 :    
252 :     if (! -s "$ATdir/FR-index")
253 :     {
254 :     print "Sorry, no index exists for the FR entries\n";
255 :     }
256 :     else
257 :     {
258 :     my @index = `cat $ATdir/FR-index`;
259 :     foreach $_ (@index)
260 :     {
261 :     print $_;
262 :     }
263 :     print "\n";
264 :     }
265 :     }
266 :    
267 :     sub build_reps {
268 :     my($sapO,$ATdir,$index) = @_;
269 :    
270 :     my $dir = "$ATdir/$index";
271 :     my $seqsF = "$dir/seqs";
272 :     -s $seqsF || die "Could not find $seqsF";
273 :    
274 :     my @seqs = &gjoseqlib::read_fasta($seqsF);
275 :     @seqs = &gjoalignandtree::order_sequences_by_length(\@seqs, { order => 'median_outward' } );
276 :     my %to_tuple = map { $_->[0] => $_ } @seqs;
277 :     printf "Building reps for %d sequences...\n", scalar@seqs;
278 :    
279 :     my($reps,$repH) = &representative_sequences::rep_seq(\@seqs,{ max_ref_sim => 0.8 } );
280 :     my @repIds = map {$_->[0] } @$reps;
281 :    
282 :     &gjoseqlib::print_alignment_as_fasta("$dir/reps-0.8",$reps);
283 :     open(GROUPS,">$dir/groups-0.8") || die "Could not make $dir/groups-0.8";
284 :     foreach my $repId (@repIds)
285 :     {
286 :     print GROUPS join("\t",($repId, @{$repH->{$repId}})),"\n";
287 :     }
288 :     close(GROUPS);
289 :     printf "Built %d reps with max similarity 0.8 \n", scalar@$reps;
290 :    
291 :     my($reps2,$repH2) = &representative_sequences::rep_seq($reps,{ max_ref_sim => 0.2 } );
292 :     &gjoseqlib::print_alignment_as_fasta("$dir/reps-0.2", $reps2);
293 :     open(GROUPS,">$dir/groups-0.2") || die "Could not make $dir/groups-0.2";
294 :     foreach my $repId (map { $_->[0] } @$reps2)
295 :     {
296 :     print GROUPS join("\t",($repId, @{$repH2->{$repId}})),"\n";
297 :     }
298 :     close(GROUPS);
299 :     printf "Built %d rep(s) with max similarity 0.2 \n", scalar@$reps2;
300 :    
301 :     if (@$reps2 > 1)
302 :     {
303 :     my $subI = 1;
304 :     my @groups = `cat $dir/groups-0.2`;
305 :     @groups = sort { split(/\t/, $b) <=> split(/\t/, $a) } @groups;
306 :     foreach my $group (@groups) {
307 :     my $subdir = "$dir/$index.$subI";
308 :     &SeedUtils::verify_dir($subdir);
309 :     chomp($group);
310 :     my @repIds = split(/\t/, $group);
311 :     open(GROUPS,">$subdir/groups-0.8") || die "Could not make $subdir/groups-0.8";
312 :     my @seqIds;
313 :     foreach my $repId (@repIds)
314 :     {
315 :     print GROUPS join("\t",($repId, @{$repH->{$repId}})),"\n";
316 :     @seqIds = (@seqIds, $repId, @{$repH->{$repId}});
317 :     }
318 :     close(GROUPS);
319 :     &gjoseqlib::print_alignment_as_fasta("$subdir/seqs", [ map { $to_tuple{$_} } @seqIds ] );
320 :     &gjoseqlib::print_alignment_as_fasta("$subdir/reps-0.8", [ map { $to_tuple{$_} } @repIds ] );
321 :     printf "Built sub dir $index.$subI containing %5d sequences with %3d reps\n", scalar@seqIds, scalar@repIds;
322 :     $subI++;
323 :     }
324 :     }
325 :     elsif (@$reps2 == 1)
326 :     {
327 :     my $subdir = "$dir/$index.1";
328 :     &SeedUtils::verify_dir($subdir);
329 :     system "cp $dir/seqs $subdir/";
330 :     system "cp $dir/reps-0.8 $subdir/";
331 :     system "cp $dir/groups-0.8 $subdir/";
332 :     }
333 :    
334 :     }
335 :    
336 :     sub summarize_FR {
337 :     my($ATdir,$index) = @_;
338 :    
339 :     my $dir = "$ATdir/$index";
340 :     my $seqsF = "$dir/seqs";
341 :    
342 :     if (! -d $dir)
343 :     {
344 :     print "Directory $dir does not exist\n\n";
345 :     }
346 :     elsif (! -s $seqsF)
347 :     {
348 :     print "FR directory is empty\n\n";
349 :     }
350 :     else
351 :     {
352 :     my $nseqs = &num_seqs_in_fasta($seqsF);
353 :     print "FR is associated with $nseqs sequences\n";
354 :     if (! -s "$dir/groups-0.8" || ! -s "$dir/reps-0.8" || ! -s "$dir/groups-0.2")
355 :     {
356 :     print "Reps have not been built for FR\n";
357 :     }
358 :     else
359 :     {
360 :     my @groups = `cat $dir/groups-0.8`;
361 :     printf "%3d reps with max similarity 0.8\n", scalar@groups;
362 :     my @groups2 = `cat $dir/groups-0.2`;
363 :     if (@groups2 == 1)
364 :     {
365 :     printf "Only 1 rep with max similarity 0.2\n";
366 :     }
367 :     else
368 :     {
369 :     printf "%3d reps with max similarity 0.2\n", scalar@groups2;
370 :     my $complete = 1;
371 :     foreach (1..@groups2)
372 :     {
373 :     if (! -s "$dir/$index.$_/seqs" || ! -s "$dir/$index.$_/groups-0.8")
374 :     {
375 :     $complete = 0;
376 :     last;
377 :     }
378 :     }
379 :     if (!$complete)
380 :     {
381 :     print "subgroups incomplete. consider rebuilding reps.\n";
382 :     }
383 :     else
384 :     {
385 :     foreach (1..@groups2)
386 :     {
387 :     my $subdir = "$dir/$index.$_";
388 :     my $nseqs = &num_seqs_in_fasta("$subdir/seqs");
389 :     my $nreps = &num_seqs_in_fasta("$subdir/reps-0.8");
390 :     printf "Subgroup $index.$_ contains %5d sequences with %3d reps\n", $nseqs, $nreps;
391 :     }
392 :     }
393 :     }
394 :     }
395 :     }
396 :     print "\n";
397 :     }
398 :    
399 :     sub align1 {
400 :     my ($sapO, $ATdir, $index, $subI) = @_;
401 :     my $dir = "$ATdir/$index/$index.$subI";
402 :     print "Aligning reps in $index.$subI...\n";
403 :     my @reps = &gjoseqlib::read_fasta("$dir/reps-0.8");
404 :     my $align = &gjoalignment::align_with_muscle(\@reps);
405 :     &gjoseqlib::print_alignment_as_fasta("$dir/align1", $align);
406 :     }
407 :    
408 :     sub trim1 {
409 :     my ($sapO, $ATdir, $index, $subI) = @_;
410 :     my $dir = "$ATdir/$index/$index.$subI";
411 :     print "Trimming align1 in $index.$subI...\n";
412 :     my @align = &gjoseqlib::read_fasta("$dir/align1");
413 :     my $trim = &gjoalignandtree::trim_align_to_median_ends(\@align, { begin => 1, end => 1 } );
414 :     &gjoseqlib::print_alignment_as_fasta("$dir/trim1", $trim);
415 :     }
416 :    
417 :     sub num_seqs_in_fasta {
418 :     my ($fasta) = @_;
419 :     my $nseqs = `grep ">" $fasta |wc -l`;
420 :     chomp($nseqs);
421 :     return $nseqs;
422 :     }
423 :    
424 :     sub extend_FR {
425 :     my ($sapO,$ATdir,$index) = @_;
426 :     my $dir = "$ATdir/$index";
427 :     if (! -s "$dir/seqs") {
428 :     print "Nothing to extend on: $dir/seqs does not exist.\n";
429 :     } else {
430 :     if (! -s "$dir/reps-0.8" || ! -s "$dir/groups-0.8" || ! -s "$dir/groups-0.2") {
431 :     &build_reps($sapO, $ATdir, $index);
432 :     }
433 :     my @groups2 = `cat $dir/groups-0.2`;
434 :     my $complete = 1;
435 :     foreach my $subI (1 .. @groups2) {
436 :     my $subdir = "$dir/$index.$subI";
437 :     if (! -s "$subdir/seqs" || ! -s "$subdir/groups-0.8" || ! -s "$subdir/reps-0.8") {
438 :     $complete = 0;
439 :     last;
440 :     }
441 :     }
442 :     if (!$complete) {
443 :     &build_reps($sapO, $ATdir, $index);
444 :     }
445 :     foreach my $subI (1 .. @groups2) {
446 :     my $subdir = "$dir/$index.$subI";
447 :     &align1($sapO, $ATdir, $index, $subI) unless -s "$subdir/align1";
448 :     &trim1 ($sapO, $ATdir, $index, $subI) unless -s "$subdir/trim1";
449 :     }
450 :     print "FR $index is now complete\n";
451 :     }
452 :     }
453 :    
454 :     sub help {
455 :     print <<END;
456 :     align1 FRindex.subIndex
457 :     build_reps FRindex
458 :     extend_FR FRindex
459 :     get_FR 'Subsystem' 'FR' ... 'FR'
460 :     process 'Subsystem' 'FR' ... 'FR'
461 :     show_roles_in_subsys 'Subsystem'
462 :     show_subsys
463 :     show_FR_index
464 :     summarize_FR FRindex
465 :     trim1 FRindex.subIndex
466 :    
467 :     END
468 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3