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

Annotation of /FigKernelScripts/rapid_propagation4.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.1 # -*- perl -*-
2 :     ########################################################################
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     ########################################################################
18 :    
19 :     $SIG{HUP} = 'IGNORE'; # ... Force running 'nohup' ...
20 :    
21 :     use strict;
22 :     use warnings;
23 :    
24 :     use FIGV;
25 :     use FIG_Config;
26 :     use GenomeMeta;
27 :     use SAPserver;
28 :     use ANNOserver;
29 :     use gjoseqlib;
30 :    
31 :     use Getopt::Long;
32 :     use File::Basename;
33 :     use Carp;
34 :    
35 :     $ENV{PATH} .= ":" . join(":", $FIG_Config::ext_bin, $FIG_Config::bin);
36 :    
37 :     #
38 :     # Make figfams.pm report details of family data used, for any FF based code.
39 :     #
40 :     $ENV{REPORT_FIGFAM_DETAILS} = 1;
41 :    
42 :     $0 =~ m/([^\/]+)$/;
43 :     my $self = $1;
44 :     my $usage = "$self [--keep --glimmerV=[2,3] --kmerDataset=ReleaseID --code=num --errdir=dir --tmpdir=dir --chkdir=dir --meta=metafile] FromDir NewDir [NewDir must not exist]";
45 :    
46 :    
47 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
48 :     # Explanation of the various work directories
49 :     # ===========================================
50 :     # $origD --- the original raw genome directory
51 :     # $newD --- the destination for the completed rapid propagation
52 :     # $procD --- the temp directory in which the intermediate computations are held
53 :     # $tmpD --- the genome directory held within $procD
54 :     # $errD --- the directory to which error outputs are written to.
55 :     # $chkD --- the directory "Checkpoint" directories are copied to.
56 :     #
57 :     # In RAST we are called with the following assignments:
58 :     # $origD = jobdir/raw/<genomeid>
59 :     # $newD = jobdir/rp/<genomeid>
60 :     # $procD = /scratch/tmprp.job<jobnumber>.<pid>
61 :     # $errD = jobir/rp.errors
62 :     # $chkD = jobdir/Restart
63 :     #-----------------------------------------------------------------------
64 :    
65 :     my ($tmpD, $restartD);
66 :     my ($meta_file, $meta);
67 :    
68 :     my $bin = $FIG_Config::bin;
69 :    
70 :     my $procD = "$FIG_Config::temp/tmp_rp.$$";
71 :     my $errD = "";
72 :     my $chkD = "";
73 :    
74 :     my $code = 11;
75 :     my $domain = "";
76 :    
77 :     my $keep = 0;
78 :     my $fix_fs = 0;
79 :     my $backfill_gaps = 0;
80 :     my $glimmer_version = 3;
81 :     my $num_neighbors = 30;
82 :    
83 :     my $NR;
84 :     my $kmerDataset = defined($FIG_Config::kmerDataset) ? $FIG_Config::kmerDataset : q();
85 :    
86 :     my $rc = GetOptions(
87 :     "keep!" => \$keep,
88 :     "fix_fs!" => \$fix_fs,
89 :     "backfill!" => \$backfill_gaps,
90 :     "glimmerV=s" => \$glimmer_version,
91 :     "code=s" => \$code,
92 :     "domain=s" => \$domain,
93 :     "tmpdir=s" => \$procD,
94 :     "errdir=s" => \$errD,
95 :     "chkdir=s" => \$chkD,
96 :     "meta=s" => \$meta_file,
97 :     "nr=s" => \$NR,
98 :     "kmerDataset=s" => \$kmerDataset,
99 :     );
100 :     if (!$rc) {
101 :     die "\n usage: $usage\n\n";
102 :     }
103 :     my $old = $keep ? " old" : "";
104 :    
105 :    
106 :     if (defined($code)) {
107 :     if ($code !~ /^\d+$/) {
108 :     die "Genetic code must be numeric\n";
109 :     }
110 :     }
111 :     else {
112 :     warn "Warning: No genetic code passed to $self. Defaulting to 11.\n";
113 :     $code = 11;
114 :     }
115 :    
116 :    
117 :     #...Set $glimmer_version if user does not want the default version=3
118 :     if ($glimmer_version == 2) {
119 :     $ENV{RAST_GLIMMER_VERSION} = $glimmer_version;
120 :     }
121 :     elsif ($glimmer_version == 3) {
122 :     $ENV{RAST_GLIMMER_VERSION} = $glimmer_version;
123 :     }
124 :     else {
125 :     die "GLIMMER version $glimmer_version not supported. Only versions 2 and 3 supported";
126 :     }
127 :    
128 :    
129 :     #...Handle the mandatory arguments...
130 :     my $trouble = 0;
131 :     my ($origD, $newD) = @ARGV;
132 :    
133 :     if (!-d $origD) {
134 :     $trouble = 1;
135 :     warn "$origD does not exist";
136 :     }
137 :    
138 :     if (-d $newD) {
139 :     $trouble = 1;
140 :     warn "$newD already exists";
141 :     }
142 :    
143 :     die "\n\n usage: $usage\n\n" if $trouble;
144 :    
145 :    
146 :     my $taxonID = basename($origD);
147 :     if ($taxonID !~ /^\d+\.\d+$/) {
148 :     die "FromDir must end in a valid genome identifier (e.g., \"83333.1\"\n";
149 :     }
150 :    
151 :     &FIG::verify_dir($procD);
152 :    
153 :    
154 :    
155 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
156 :     # ... Define temporary, error, and checkpoint directories.
157 :     #-----------------------------------------------------------------------
158 :     $tmpD = "$procD/$taxonID";
159 :     &FIG::verify_dir($tmpD);
160 :    
161 :     if (! $errD) {
162 :     $errD = $tmpD;
163 :     }
164 :    
165 :     if (! $chkD) {
166 :     $chkD = $errD;
167 :     }
168 :    
169 :     $restartD = "$chkD/$taxonID.restart";
170 :    
171 :    
172 :    
173 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
174 :     # ... Set metafile variables ...
175 :     #-----------------------------------------------------------------------
176 :     if ($meta_file) {
177 :     $meta = new GenomeMeta($taxonID, $meta_file);
178 :    
179 :     $meta->set_metadata("status.rp", "in_progress");
180 :     $meta->add_log_entry("rp", ["Begin processing", $origD, $procD, $taxonID]);
181 :    
182 :     $meta->set_metadata("rp.glimmer_version", $glimmer_version);
183 :    
184 :     if ($fix_fs) {
185 :     $meta->set_metadata("correction.frameshifts", 1);
186 :     }
187 :    
188 :     if ($backfill_gaps) {
189 :     $meta->set_metadata("correction.backfill_gaps", 1);
190 :     }
191 :    
192 :     $ENV{DEBUG} = $meta->get_metadata('env.debug') || 0;
193 :     $ENV{VERBOSE} = $meta->get_metadata('env.verbose') || 0;
194 :     }
195 :     else {
196 :     die "$self now requires a meta file";
197 :     }
198 :    
199 :     my $restart = 0;
200 :     if (-d $restartD) {
201 :     warn "NOTE: This is a restart run\n";
202 :     $meta->add_log_entry("rp", ["Restart of previous job"]);
203 :    
204 :     $restart = 1;
205 :     system "rm -r $tmpD";
206 :     system(qq(cp -R $restartD $tmpD))
207 :     && die qq(Could not copy $restartD ==> $tmpD);
208 :     }
209 :     else {
210 :     &my_run("/bin/cp -R $origD $procD");
211 :     &make_restart($chkD, $tmpD, $restartD, qq(orig));
212 :    
213 :     if (! $keep) {
214 :     &my_run("rm -fR $tmpD/Features/* $tmpD/assigned_functions");
215 :     }
216 :     }
217 :    
218 :    
219 :     if ($code) {
220 :     print STDERR "Using genetic code $code\n" if $ENV{VERBOSE};
221 :    
222 :     open( GENETIC_CODE, ">$tmpD/GENETIC_CODE" )
223 :     || die "Could not write-open $tmpD/GENETIC_CODE";
224 :     print GENETIC_CODE "$code\n";
225 :     close(GENETIC_CODE) || die "Could not close $tmpD/GENETIC_CODE";
226 :     }
227 :    
228 :    
229 :    
230 :     ########################################################################
231 :     #...Main Body of Code...
232 :     ########################################################################
233 :    
234 :     my $figV = FIGV->new($tmpD);
235 :     my $sapO = SAPserver->new();
236 :     my $annoO = ANNOserver->new();
237 :    
238 :     my $trans_table = FIG::genetic_code($code);
239 :    
240 :     my $contigs_file = qq($tmpD/contigs);
241 :     my $contig_lens = {};
242 :     %$contig_lens = map { $_->[0] => length($_->[2]) } read_fasta($contigs_file);
243 :    
244 :     my $min_contig_len = 2_000;
245 :     my $min_tot_dna = 100_000;
246 :     my $max_frac_short = 0.90;
247 :    
248 :     my $result;
249 :     my $initial_fasta_txt;
250 :     my $initial_calls;
251 :     if ((not &is_too_small( $contig_lens, $min_contig_len, $min_tot_dna)) &&
252 :     (not &is_metagenome($contig_lens, $min_contig_len, $max_frac_short)) &&
253 :     (not &is_raw_reads( $contig_lens, $min_contig_len)) &&
254 :     (0 == &run_find_genes($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code,
255 :     $chkD,$tmpD,$restartD,$errD,
256 :     q(find_genes_based_on_kmers.stage-0.stderr),
257 :     $kmerDataset))
258 :     ) {
259 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
260 :     #... Large genome ...
261 :     #-----------------------------------------------------------------------
262 :     if (!$ keep) {
263 :     &make_restart($chkD,$tmpD,$restartD, qq(find_genes_based_on_kmers.stage-0));
264 :     }
265 :     &find_rnas($keep, $meta,$figV,$sapO,$annoO, $domain,$taxonID, $chkD,$tmpD,$restartD,$errD, q(find_rnas.stderr));
266 :     &find_special_proteins($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD);
267 :     &post_process_rnas_and_glimmer($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code,
268 :     $chkD,$tmpD,$restartD,$errD,
269 :     q(postprocess_rna_and_glimmer.stderr),
270 :     $kmerDataset);
271 :    
272 :     &correct_frameshifts($keep,$fix_fs, $meta,$figV,$sapO,$annoO, $chkD,$tmpD,$restartD,$errD);
273 :     &find_and_backfill_missing_and_miscalled($keep,$backfill_gaps, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD);
274 :     }
275 :     elsif (0 == &run_mga($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $kmerDataset)) {
276 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
277 :     #... Small genome, plasmid, or fragment ...
278 :     #-----------------------------------------------------------------------
279 :     if (!$keep) {
280 :     if (($domain eq q(Bacteria)) || ($domain eq q(Archaea))) {
281 :     &find_rnas($keep, $meta,$figV,$sapO,$annoO, $domain,$taxonID, $chkD,$tmpD,$restartD,$errD, q(find_rnas.stderr));
282 :     }
283 :     &post_process_mga($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code,
284 :     $chkD,$tmpD,$restartD,$errD,
285 :     qq(postprocess_mga.stderr), $kmerDataset);
286 :    
287 :     &correct_frameshifts($keep,$fix_fs, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD);
288 :     &find_and_backfill_missing_and_miscalled($keep,$backfill_gaps, $meta,$figV,$sapO,$annoO, $taxonID,$code,
289 :     $chkD,$tmpD,$restartD,$errD);
290 :     }
291 :     }
292 :     else {
293 :     &kmer_approach();
294 :     }
295 :    
296 :     &cleanup($bin, $keep, $origD,$restartD, $tmpD,$errD,$newD);
297 :    
298 :     $meta->set_metadata("status.rp", "complete");
299 :     $meta->set_metadata("rp.running", "no");
300 :    
301 :     exit(0);
302 :    
303 :    
304 :    
305 :    
306 :     ########################################################################
307 :     #...Utility routines...
308 :     ########################################################################
309 :     sub date_stamp {
310 :     my ($sec,$min,$hour,$mday,$mon,$year,$wday, $yday,$isdst) = localtime(time);
311 :     return sprintf ("%4d-%02d-%02d %02d:%02d:%02d",
312 :     ($year+1900, $mon+1, $mday, $hour, $min, $sec)
313 :     );
314 :     }
315 :    
316 :     sub log {
317 :     my($dir,$message) = @_;
318 :    
319 :     my $whole_message = scalar(localtime(time)) . ": $message";
320 :    
321 :     open( LOG, ">>$dir/log");
322 :     print LOG "$whole_message\n";
323 :     close(LOG);
324 :    
325 :     $meta->add_log_entry("rp", $whole_message) if $meta;
326 :     }
327 :    
328 :     sub my_run {
329 :     my($cmd, $nofatal) = @_;
330 :     $cmd =~ s/\n/ /gs;
331 :    
332 :     use Carp qw( cluck );
333 :     if ($cmd =~ m{^([^/]+?)(\s+.*)$}) {
334 :     my $prog = $1;
335 :     my $rest = $2;
336 :    
337 :     my $changed;
338 :     for my $dir (split(/:/, $ENV{PATH})) {
339 :     my $path = "$dir/$prog";
340 :     if (-x $path) {
341 :     warn "using $path\n";
342 :     $cmd = "$path $rest";
343 :     ++$changed;
344 :     last;
345 :     }
346 :     }
347 :     warn "Cmd running bare command $prog\n" unless $changed;
348 :     }
349 :    
350 :     if ($ENV{FIG_VERBOSE}) {
351 :     print STDERR (&date_stamp(), q(: running ), $cmd, qq(\n\n));
352 :     }
353 :    
354 :     if ($meta) {
355 :     $meta->add_log_entry("rp", ['run_start', $cmd]);
356 :     }
357 :    
358 :     my $rc = system($cmd);
359 :     if ($meta) {
360 :     $meta->add_log_entry("rp", ['run_finish', $cmd, $rc]);
361 :     }
362 :    
363 :     if ($rc != 0) {
364 :     my $msg;
365 :     (undef, undef, $msg) = &FIG::interpret_error_code($rc);
366 :     if ($nofatal) {
367 :     cluck "WARNING, rc=$rc, reason=$msg: $cmd";
368 :     }
369 :     else {
370 :     confess "FAILED, rc=$rc, reason=$msg: $cmd";
371 :     }
372 :     }
373 :     }
374 :    
375 :    
376 :     sub is_too_small {
377 :     my ($contigs_lens, $min_contig_len, $min_tot_dna) = @_;
378 :    
379 :     my $tot_dna = 0;
380 :     foreach my $contig (keys %$contig_lens) {
381 :     if ($contig_lens->{$contig} >= $min_contig_len) {
382 :     $tot_dna += $contig_lens->{$contig};
383 :     }
384 :     }
385 :    
386 :     my $rc = ($tot_dna < $min_tot_dna) || 0;
387 :     print STDERR qq(too_small: tot_dna=$tot_dna, min_tot_dna=$min_tot_dna ==> rc=$rc\n\n)
388 :     if $ENV{VERBOSE};
389 :     return $rc
390 :     }
391 :    
392 :    
393 :     sub is_metagenome {
394 :     my ($contigs_lens, $min_contig_len, $max_frac_short) = @_;
395 :    
396 :     my $tot_dna = 0;
397 :     my $tot_short_dna = 0;
398 :     foreach my $contig (keys %$contig_lens) {
399 :     $tot_dna += $contig_lens->{$contig};
400 :     if ($contig_lens->{$contig} <= $min_contig_len) {
401 :     $tot_short_dna += $contig_lens->{$contig};
402 :     }
403 :     }
404 :    
405 :     my $rc = ($tot_short_dna >= $max_frac_short * $tot_dna) || 0;
406 :     print STDERR qq(is_metagenome: tot_dna=$tot_dna, tot_short_dna=$tot_short_dna, max_frac_short=$max_frac_short ==> rc=$rc\n\n)
407 :     if $ENV{VERBOSE};
408 :     return $rc;
409 :     }
410 :    
411 :    
412 :     sub is_raw_reads {
413 :     my ($contigs_lens, $min_contig_len) = @_;
414 :    
415 :     my $tot_long_dna = 0;
416 :     foreach my $contig (keys %$contig_lens) {
417 :     if ($contig_lens->{$contig} >= $min_contig_len) {
418 :     $tot_long_dna += $contig_lens->{$contig};
419 :     }
420 :     }
421 :    
422 :     my $rc = ($tot_long_dna == 0) || 0;
423 :     print STDERR qq(is_raw_reads: tot_long_dna=$tot_long_dna ==> rc=$rc\n\n)
424 :     if $ENV{VERBOSE};
425 :     return $rc;
426 :     }
427 :    
428 :    
429 :    
430 :    
431 :     sub make_restart {
432 :     my($chkD, $tmpD, $restartD, $checkpoint_extension) = @_;
433 :     print STDERR (qq(\nmake_restart:\n),
434 :     qq( extension=\'$checkpoint_extension\'\n),
435 :     qq( chkD=\'$chkD\'\n),
436 :     qq( tmpD=\'$tmpD\'\n),
437 :     qq( restartD=\'$restartD\'\n\n)
438 :     ) if ($ENV{DEBUG} && (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1)));
439 :    
440 :     if (-d $restartD) {
441 :     system(q(/bin/rm), q(-fR), $restartD)
442 :     && confess qq(Could not remove restartD=\'$restartD\');
443 :     }
444 :    
445 :     system(q(/bin/cp), q(-R), $tmpD, $restartD)
446 :     && die qq(Could not copydir tmpD=\'$tmpD\' to restartD=\'$restartD\');
447 :    
448 :     system(q(/bin/touch), $restartD)
449 :     && die qq(Could not touch timestamp for restartD=\'$restartD\');
450 :    
451 :     if ($chkD && (-d $chkD) && $checkpoint_extension) {
452 :     my $name = basename($tmpD);
453 :     my $chkD_extended = qq($chkD/$name.$checkpoint_extension);
454 :    
455 :     system(q(/bin/cp), q(-R), $tmpD, $chkD_extended)
456 :     && die qq(Could not copydir tmpD=\'$tmpD\' to chkD=\'$chkD_extended\');
457 :    
458 :     system(q(/bin/touch), $chkD_extended)
459 :     && die qq(Could not touch timestamp for chkD=\'$chkD_extended\');
460 :     }
461 :    
462 :     return 1;
463 :     }
464 :    
465 :    
466 :     sub find_rnas {
467 :     my ($keep, $meta,$figV,$sapO,$annoO, $domain,$taxonID, $chkD,$tmpD,$restartD,$errD, $errfile_name) = @_;
468 :     my ($msg, $tmp, $genus, $species);
469 :    
470 :     #... Skip if not a prokaryote ...
471 :     return 1 if ($domain && ($domain ne q(Bacteria)) && ($domain ne q(Archaea)));
472 :    
473 :     #...
474 :     if (-s "$tmpD/Features/peg/tbl") {
475 :     ($tmp, $genus, $species) = &find_nearest_neighbor($figV, $sapO, $tmpD, $num_neighbors);
476 :     $msg = qq(find_nearest_neighbor returns: domain=$domain, genus=$genus, species=$species);
477 :     warn($msg, qq(\n\n)) if $ENV{VERBOSE};
478 :     &log($tmpD, $msg) if $ENV{VERBOSE};
479 :     }
480 :    
481 :     $domain = substr($domain, 0, 1) || $tmp;
482 :     $genus = quotemeta($genus) || q(Unknown);
483 :     $species = quotemeta($species) || q(sp.);
484 :    
485 :     my $contigs_fh;
486 :     my $contigs_file = qq($tmpD/contigs);
487 :     open($contigs_fh, qq(<$contigs_file))
488 :     || die qq(Could not read-open \'$contigs_file\');
489 :    
490 :     my $results = $annoO->find_rnas(-input => $contigs_fh,
491 :     -domain => $domain,
492 :     -genus => $genus,
493 :     -species => $species,
494 :     );
495 :     close($contigs_fh);
496 :    
497 :     my ($rna_fasta, $rna_tuples) = @$results;
498 :    
499 :     if (not $rna_fasta) {
500 :     $msg = q(No RNAs found);
501 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
502 :     &log($tmpD, $msg);
503 :     }
504 :     else {
505 :     my $tmp_rna_fh;
506 :     my $tmp_rna_fasta_file = qq($FIG_Config::temp/tmp_rna.$$.fasta);
507 :     open($tmp_rna_fh, qq(>$tmp_rna_fasta_file))
508 :     || die qq(Could not write-open \'$tmp_rna_fasta_file\');
509 :     print $tmp_rna_fh $rna_fasta;
510 :     close($tmp_rna_fh);
511 :    
512 :     open($tmp_rna_fh, qq(<$tmp_rna_fasta_file))
513 :     || die qq(Could not read-open \'$tmp_rna_fasta_file\');
514 :    
515 :     my %seq_of = map { $_->[0] => $_->[2] } read_fasta($tmp_rna_fh);
516 :    
517 :     foreach my $feature (@$rna_tuples) {
518 :     my ($id, $contig, $beg, $end, $func) = @$feature;
519 :     my $loc = join(q(_), ($contig, $beg, $end));
520 :     my $fid = $figV->add_feature(q(master:rast), $taxonID, q(rna), $loc, q(), $seq_of{$id},
521 :     qq($tmpD/called_by), q(find_rnas));
522 :     if ($fid) {
523 :     $figV->assign_function($fid, q(master:rast), $func);
524 :     }
525 :     else {
526 :     die (qq(Could not add feature=$id\n), Dumper($rna_tuples, \%seq_of));
527 :     }
528 :     }
529 :     }
530 :    
531 :     return 1;
532 :     }
533 :    
534 :     sub find_special_proteins {
535 :     my ($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD) = @_;
536 :    
537 :     my @contigL = &gjoseqlib::read_fasta(qq($tmpD/contigs));
538 :    
539 :     my $selenoL = $annoO->find_special_proteins(-contigs => \@contigL,
540 :     -templates => q(selenoprotein),
541 :     -comment => q(selenoprotein),
542 :     );
543 :    
544 :     my $pyrroL = $annoO->find_special_proteins(-contigs => \@contigL,
545 :     -templates => q(pyrrolysine),
546 :     -comment => q(pyrrolysoprotein),
547 :     );
548 :     open(SPECIAL, qq(>$tmpD/special_pegs))
549 :     || die qq(Could not write-open special_pegs file '$tmpD/special_pegs');
550 :     foreach my $prot (@$selenoL, @$pyrroL) {
551 :     my $type = $prot->{comment};
552 :     my $func = $prot->{reference_def};
553 :    
554 :     my $fid = $figV->add_feature(q(master:rast), $taxonID, q(peg), $prot->{location}, q(),
555 :     $prot->{sequence}, qq($tmpD/called_by), q(find_special_proteins));
556 :     if ($fid) {
557 :     print SPECIAL qq($fid\t$type\n);
558 :     $figV->assign_function($fid, q(master:rast), $func, qq($tmpD/proposed_functions));
559 :     }
560 :     else {
561 :     die (qq(Could not add feature=$fid\n), Dumper($prot));
562 :     }
563 :     }
564 :     close(SPECIAL);
565 :    
566 :     return 1;
567 :     }
568 :    
569 :    
570 :     sub run_find_genes {
571 :     my ($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $errfile_name, $kmerDataset) = @_;
572 :     # use IPC::Run qw( run timeout );
573 :    
574 :     my $old = $keep ? q(old) : q();
575 :     my $kmerDataSwitch = $kmerDataset ? qq(-kmerDataset=$kmerDataset) : q();
576 :    
577 :     my $rc = &my_run(qq(find_genes_based_on_kmers $kmerDataSwitch $tmpD $tmpD/found $old \> $errD/$errfile_name 2\>\&1));
578 :    
579 :     if ($ENV{VERBOSE}) {
580 :     if ($rc == 0) {
581 :     print STDERR qq(run_find_genes succeeded:\trc=$rc);
582 :     }
583 :     else {
584 :     print STDERR qq(run_find_genes failed:\trc=$rc);
585 :     }
586 :     }
587 :    
588 :     return $rc;
589 :     }
590 :    
591 :    
592 :     sub post_process_rnas_and_glimmer {
593 :     my ($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $errfile_name, $kmerDataset) = @_;
594 :     my ($msg, $rc);
595 :     my $kmerDataSwitch = $kmerDataset ? qq(-kmerDataset=$kmerDataset) : q();
596 :    
597 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
598 :     #...Promote remaining ORFs to PEGs...
599 :     #-----------------------------------------------------------------------
600 :     if (! $meta->get_metadata('status.rp.postprocess_rnas_and_glimmer')) {
601 :     if ($keep) {
602 :     #...If just re-annotating an existing set of calls,
603 :     # write assigned_functions for the PEG IDs listed in the found file.
604 :    
605 :     system "rm -fR $tmpD/Features/orf";
606 :     open(IN,"<$tmpD/found") || die "could not open $tmpD/found";
607 :     open(OUT,">$tmpD/assigned_functions") || die "could not open $tmpD/assigned_functions";
608 :     my %seen;
609 :     while (defined($_ = <IN>)) {
610 :     chomp;
611 :     my($peg,$func) = split(/\t/,$_);
612 :     if (! $seen{$peg}) {
613 :     $seen{$peg} = 1;
614 :     print OUT "$peg\t$func\n";
615 :     }
616 :     }
617 :     close(IN);
618 :     close(OUT);
619 :     }
620 :     else {
621 :     &my_run("$bin/find_genes_based_on_kmers $kmerDataSwitch $tmpD $tmpD/found > $errD/find_genes_based_on_kmers.stage-1.stderr 2>&1");
622 :     &log($tmpD, 'Finished first pass of finding genes matching families found in close genomes');
623 :     &make_restart($chkD, $tmpD, $restartD, qq(find_genes_based_on_kmers.stage-1));
624 :    
625 :     &my_run("$bin/find_genes_based_on_kmers $kmerDataSwitch $tmpD $tmpD/found > $errD/find_genes_based_on_kmers.stage-2.stderr 2>&1");
626 :     &log($tmpD, 'Finished second pass of finding genes matching families found in close genomes');
627 :     &make_restart($chkD, $tmpD, $restartD, qq(find_genes_based_on_kmers.stage-2));
628 :    
629 :     #...Promote the remaining ORFs to PEGs (writes functions to the assigned_functions file)
630 :     print STDERR (qq(before promotion:\t), `wc $tmpD/Features/*/tbl`, qq(\n))
631 :     if ($ENV{DEBUG} && (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1)));
632 :     &my_run("$bin/promote_orfs_to_pegs $tmpD $tmpD/found > $errD/promote_orfs.stderr 2>&1");
633 :     print STDERR (qq(after promotion:\t), `wc $tmpD/Features/*/tbl`, qq(\n))
634 :     if ($ENV{DEBUG} && (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1)));
635 :     &log($tmpD, 'Finished promoting genes that could not be placed in any FIGfam');
636 :     }
637 :    
638 :     &make_restart($chkD, $tmpD, $restartD, qq(promote_orfs_to_pegs));
639 :     $meta->set_metadata('status.rp.promote_orfs_to_pegs',1);
640 :     }
641 :    
642 :     return 1;
643 :     }
644 :    
645 :     sub cleanup {
646 :     my ($bin, $keep, $origD, $restartD, $tmpD, $errD, $newD) = @_;
647 :    
648 :     if (!$keep) {
649 :     if (!-d qq($tmpD/Features/peg)) { mkdir(qq($tmpD/Features/peg)); }
650 :     if (!-e qq($tmpD/Features/peg/tbl)) { open(TMP, qq(>$tmpD/Features/peg/tbl)); }
651 :     if (!-e qq($tmpD/Features/peg/fasta)) { open(TMP, qq(>$tmpD/Features/peg/fasta)); }
652 :     close(TMP);
653 :     }
654 :    
655 :    
656 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
657 :     # make_genome_dir_for close will copy $tmpD to $newD, and rename the copied
658 :     # assigned_functions to proposed_functions.
659 :     #
660 :     # In the case where we are not keeping genecalls, we wish to recreate the
661 :     # assigned_functions from the original raw directory.
662 :     # We can use the tbl-based peg mapping that make_genome_dir_for_close also uses
663 :     # to do this; and in fact do this inline here, after make_genome_dir_for_close.
664 :     #-----------------------------------------------------------------------
665 :     &my_run("$bin/make_genome_dir_for_close $origD $tmpD $newD 2> $errD/make_genome_dir_for_close.stderr");
666 :     &make_restart($chkD, $newD, $restartD, qq(make_genome_dir_for_close));
667 :    
668 :     &my_run("$bin/renumber_features -print_map $newD > $errD/renumber_features.map 2> $errD/renumber_features.stderr");
669 :     &make_restart($chkD, $newD, $restartD, qq(renumber_features));
670 :    
671 :     # &my_run("cat $newD/proposed*functions | $bin/rapid_subsystem_inference $newD/Subsystems 2> $errD/rapid_subsystem_inference.stderr");
672 :     # &log($tmpD, 'Finished inferring subsystems');
673 :    
674 :     if (-s "$tmpD/Features/peg/fasta") {
675 :     &my_run("$FIG_Config::ext_bin/formatdb -p -i $tmpD/Features/peg/fasta", 1);
676 :     }
677 :    
678 :     system("/bin/cp -R $procD/* $errD/");
679 :     unless ($ENV{DEBUG}) {
680 :     system("/bin/rm -fR $procD");
681 :     }
682 :    
683 :     return 1;
684 :     }
685 :    
686 :    
687 :     sub run_mga {
688 :     my ($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $errfile_name, $kmerDataset) = @_;
689 :     my $msg;
690 :    
691 :     return (1) if ($code != 11);
692 :     my $trans_table = $figV->standard_genetic_code();
693 :    
694 :     my $contigs_file = qq($tmpD/contigs);
695 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
696 :     #...Run the fragments or plasmid sequences through the ORF caller (metagene)..
697 :     # Create training data
698 :     #-----------------------------------------------------------------------------
699 :     if (! $meta->get_metadata('status.rp.run_mga')) {
700 :     if ($keep) {
701 :     #...Do nothing --- Annotations will be handled by &post_process_mga();
702 :     }
703 :     else {
704 :     $msg = q(Starting to run MGA);
705 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
706 :     &log($tmpD, $msg);
707 :    
708 :     my $mga_pipe;
709 :     open($mga_pipe, qq($FIG_Config::mga_bin/mga $contigs_file -s |))
710 :     || die qq(Could not pipe-out open MGA);
711 :    
712 :     my @mga_hits;
713 :     my $line = <$mga_pipe>; print STDERR (q(MGA: ), $line) if ($ENV{DEBUG} && defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
714 :     while ($line && ($line =~ m/^\#\s+(\S+)/o)) {
715 :     my $contig_id = $1;
716 :     $line = <$mga_pipe>; print STDERR (q(MGA: ), $line) if ($ENV{DEBUG} && defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
717 :     $line = <$mga_pipe>; print STDERR (q(MGA: ), $line) if ($ENV{DEBUG} && defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
718 :     $line = <$mga_pipe>; print STDERR (q(MGA: ), $line) if ($ENV{DEBUG} && defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
719 :     while ($line && $line !~ /^\#/) {
720 :     chomp $line;
721 :     my ($id, $start, $end, $strand, $frame, $trunc) = split(/\t/, $line);
722 :     push @mga_hits, [$id, $contig_id, $start, $end, $strand, $frame, $trunc, $line];
723 :     $line = <$mga_pipe>; print STDERR (q(MGA: ), $line) if ($ENV{DEBUG} && defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
724 :     }
725 :     }
726 :     close($mga_pipe);
727 :    
728 :     $msg = q(MGA completed);
729 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
730 :     &log($tmpD, $msg);
731 :    
732 :     foreach my $call (@mga_hits) {
733 :     my ($id, $contig_id, $left, $right, $strand, $frame, $trunc, $line) = @$call;
734 :     print STDERR (join(qq(\t), ($id, $contig_id, $left, $right, $strand, $frame, $trunc, $line)), qq(\n))
735 :     if ($ENV{DEBUG} && defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
736 :    
737 :     my ($beg, $end);
738 :     if ($strand eq "+") {
739 :     ($beg, $end) = ($left, $right);
740 :     }
741 :     else {
742 :     ($beg, $end) = ($right, $left);
743 :     }
744 :    
745 :     my $sign = ($end <=> $beg);
746 :     print STDERR qq(RP4: beg=$beg,\tend=$end,\tsign=$sign\n)
747 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
748 :    
749 :     if ($trunc eq qq(01)) {
750 :     print STDERR qq(Fixed truncated START:\t($beg, $end) ==> ) if $ENV{VERBOSE};
751 :     $beg = $end - $sign * ( 3 * int( (1 + abs($end-$beg)) / 3) - 1);
752 :     print STDERR qq(($beg, $end):\t$line\n) if $ENV{VERBOSE};
753 :     }
754 :     elsif ($trunc eq qq(10)) {
755 :     warn qq(Fixed truncated STOP:\t($beg, $end) ==> ) if $ENV{VERBOSE};
756 :     $end = $beg + $sign * ( 3 * int( (1 + abs($end-$beg)) / 3) - 1);
757 :     print STDERR qq(($beg, $end):\t$line\n) if $ENV{VERBOSE};
758 :     }
759 :     elsif ($trunc eq qq(00)) {
760 :     warn qq(Fixed double-truncated partial:\t($beg, $end) ==> )
761 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
762 :     $beg = $beg + $sign * $frame;
763 :     $end = $beg + $sign * ( 3 * int( (1 + abs($end-$beg)) / 3) - 1);
764 :     print STDERR qq(($beg, $end):\t$line\n) if $ENV{VERBOSE};
765 :     }
766 :     print STDERR qq(===> beg=$beg,\tend=$end,\tsign=$sign\n\n)
767 :     if (defined($ENV{VERBOSE}) && ($ENV{VERBOSE} > 1));
768 :    
769 :     my $loc = join(q(_), ($contig_id, $beg, $end));
770 :     my $seq = $figV->dna_seq($taxonID, $loc);
771 :     my $pep = $figV->translate($seq, $trans_table, 1);
772 :     $pep =~ s/\*$//o;
773 :    
774 :     if (my $fid = $figV->add_feature(q(master:rast), $taxonID, q(orf), $loc, q(), $pep)) {
775 :     #...Suceeded, do nothing...
776 :     }
777 :     else {
778 :     die (q(Could not add ORF: ), $loc, qq(\t), $pep) unless $fid;
779 :     }
780 :     }
781 :     }
782 :     }
783 :    
784 :     &make_restart($chkD, $tmpD, $restartD, qq(run_mga));
785 :     $meta->set_metadata('status.rp.run_mga',1);
786 :    
787 :     return 0;
788 :     }
789 :    
790 :    
791 :     sub post_process_mga {
792 :     my ($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $errfile_name, $kmerDataset) = @_;
793 :     my ($msg);
794 :    
795 :     my $pegs_file;
796 :     if ($keep) {
797 :     #...do stuff
798 :     }
799 :     else {
800 :     #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
801 :     #...Run the training data from the metagene caller through the ORF caller (GLIMMER)..
802 :     #-----------------------------------------------------------------------------
803 :     my $tmp_orf_tbl = qq($tmpD/Features/orf/tbl);
804 :     if (! $meta->get_metadata('status.rp.annoserver.call_genes')) {
805 :     if (!-s $tmp_orf_tbl) {
806 :     $msg = q(No MGA ORFs found to train GLIMMER with);
807 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
808 :     &log($tmpD, $msg);
809 :     }
810 :     else {
811 :     $msg = q(Using MGA as training data for 'annoO->find_genes');
812 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
813 :     &log($tmpD, $msg);
814 :    
815 :     my %tmp_trainH = map { m/^(\S+)\s+(\S+)/o ? ($1 => $2) : () } &SeedUtils::file_read($tmp_orf_tbl);
816 :     # print STDERR Dumper(\%tmp_trainH);
817 :    
818 :     my ($contigs_fh, $retval_pair, $fasta_string, $prot_callsL);
819 :     my $contigs_file = qq($tmpD/contigs);
820 :     open($contigs_fh, qq(<$contigs_file))
821 :     || die qq(Could not read-open contigs file \'$contigs_file\');
822 :     if ($retval_pair = $annoO->call_genes(-input => $contigs_fh,
823 :     -trainingLocations => \%tmp_trainH,
824 :     -geneticCode => $code
825 :     )
826 :     ) {
827 :     ($fasta_string, $prot_callsL) = @$retval_pair;
828 :     my @prot_seqs = &gjoseqlib::read_fasta(\$fasta_string);
829 :     my %loc_of_prot = map { $_->[0] => join(q(_), ($_->[1], $_->[2], $_->[3])) } @$prot_callsL;
830 :     my %seq_of_prot = map { $_->[0] => $_->[2] } @prot_seqs;
831 :     # print STDERR Dumper($prot_callsL, \@prot_seqs, \%loc_of_prot, \%seq_of_prot);
832 :    
833 :     if (@prot_seqs > 0) {
834 :     $msg = q(Completed 'annoO->call_genes');
835 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
836 :     &log($tmpD, $msg);
837 :    
838 :     my $assignmentsH = &get_kmer_based_functions($annoO, 8, 2, 2, \@prot_seqs, $kmerDataset);
839 :    
840 :     if (&add_pegs($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $errfile_name,
841 :     \%loc_of_prot, \%seq_of_prot, $assignmentsH, $keep)
842 :     ) {
843 :     system("rm -fR $tmpD/Features/orf")
844 :     && die qq(Could not remove ORF directory \'$tmpD/Features/orf\');
845 :     }
846 :     else {
847 :     die qq(Could not add PEGs);
848 :     }
849 :     }
850 :     else {
851 :     $msg = q(FAILED 'annoO->find_genes');
852 :     warn ($msg, qq(\n)) if $ENV{VERBOSE};
853 :     &log($tmpD, $msg);
854 :     }
855 :     }
856 :     }
857 :     }
858 :     &make_restart($chkD, $tmpD, $restartD, qq(recall_genes_using_glimmer));
859 :     $meta->set_metadata('status.rp.recall_genes_using_glimmer',1);
860 :     }
861 :    
862 :    
863 :     # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
864 :     # #...Promote remaining ORFs to PEGs...
865 :     # #-----------------------------------------------------------------------
866 :     # if (! $meta->get_metadata('status.rp.promote_orfs_to_pegs')) {
867 :     # if ($keep) {
868 :     # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
869 :     # #...If just re-annotating an existing set of calls,
870 :     # # write assigned_functions for the PEG IDs listed in the found file.
871 :     # #-----------------------------------------------------------------------
872 :     # system "rm -fR $tmpD/Features/orf";
873 :     # my %func_of = map { m/^(\S+)\t([^\t]+)/o ? ($1 => $2) : () } &SeedUtils::file_read(qq($tmpD/found));
874 :     # foreach my $fid (sort { &SeedUtils::by_fig_id($a,$b) } (keys %func_of)) {
875 :     # $figV->assign_function($fid, q(master:rast), $func_of{$fid});
876 :     # }
877 :     # }
878 :     # else {
879 :     # #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
880 :     # #...Promote the remaining ORFs to PEGs (make assigned_functions entries)
881 :     # #-----------------------------------------------------------------------
882 :     # &my_run("$bin/promote_orfs_to_pegs $tmpD $tmpD/found > $errD/promote_orfs.stderr 2>&1");
883 :     # &log($tmpD, 'Finished promoting genes that could not be placed in any FIGfam');
884 :     # }
885 :     #
886 :     # &make_restart($chkD, $tmpD, $restartD, qq(promote_orfs_to_pegs));
887 :     # $meta->set_metadata('status.rp.promote_orfs_to_pegs',1);
888 :     # }
889 :    
890 :     return 1;
891 :     }
892 :    
893 :     sub add_pegs {
894 :     my ($keep, $meta,$figV,$sapO,$annoO, $taxonID,$code, $chkD,$tmpD,$restartD,$errD, $errfile_name,
895 :     $loc_of_protH, $seq_of_protH, $assignmentsH) = @_;
896 :    
897 :     open(FOUND, qq(>>$tmpD/found))
898 :     || die qq(could not append-open "found" file '$tmpD/found');
899 :    
900 :     foreach my $prot_id (sort { &SeedUtils::location_cmp($loc_of_protH->{$a}, $loc_of_protH->{$b}) } (keys %$loc_of_protH)) {
901 :     my $loc = $loc_of_protH->{$prot_id};
902 :     my $seq = $seq_of_protH->{$prot_id};
903 :    
904 :     my $fid = $figV->add_feature(q(master:rast), $taxonID, q(peg), $loc, q(), $seq,
905 :     qq($tmpD/called_by), q(protprocess_mga));
906 :     my $func = q();
907 :     if ($fid) {
908 :     if (defined($assignmentsH->{$prot_id}) && defined($func = $assignmentsH->{$prot_id}->[0])) {
909 :     print FOUND (join(qq(\t), ( $fid, @ { $assignmentsH->{$prot_id} })), qq(\n));
910 :    
911 :     if (not $figV->assign_function($fid, q(master:rast), $func, 2, qq($tmpD/proposed_functions))) {
912 :     die (qq(Could not add function for prot_id=$prot_id, fid=$fid, func=$func\n),
913 :     Dumper($loc_of_protH, $seq_of_protH, $assignmentsH)
914 :     );
915 :     }
916 :     }
917 :     }
918 :     else {
919 :     die (qq(Could not add feature=$prot_id\n), Dumper($loc_of_protH, $seq_of_protH, $assignmentsH));
920 :     }
921 :     }
922 :    
923 :     close(FOUND);
924 :     return 1;
925 :     }
926 :    
927 :    
928 :     sub kmer_approach {
929 :     die (q(Method 'kmer_approach' not yet implemented), qq(\n), Dumper(\@_));
930 :     }
931 :    
932 :    
933 :     sub get_kmer_based_functions {
934 :     my ($annoO, $kmer, $scoreThreshold, $seqHitThreshold, $protL, $kmerDataset) = @_;
935 :     my @kmerDataset = $kmerDataset ? (q(-kmerDataset) => $kmerDataset) : ();
936 :    
937 :     my $result_handle = $annoO->assign_function_to_prot(-input => $protL,
938 :     -kmer => $kmer,
939 :     -assignToAll => 0,
940 :     -scoreThreshold => $scoreThreshold,
941 :     -seqHitThreshold => $seqHitThreshold,
942 :     @kmerDataset
943 :     );
944 :     my $resultH = {};
945 :     while (my $result = $result_handle->get_next()) {
946 :     my ($prot_id, $function, $otu, $score, $nonoverlap_hits, $overlap_hits) = @$result;
947 :     if (defined($nonoverlap_hits) && ($nonoverlap_hits > 0)) {
948 :     $resultH->{$prot_id} = [$function, $otu, $score, $nonoverlap_hits, $overlap_hits];
949 :     }
950 :     }
951 :     return $resultH;
952 :     }
953 :    
954 :     sub find_nearest_neighbor {
955 :     my ($figV, $sapO, $orgdir, $num_neighbors) = @_;
956 :    
957 :     my $outfile = qq($orgdir/closest.genomes);
958 :     print STDERR qq(\nfind_nearest_neighbor: outfile=$outfile\n) if $ENV{VERBOSE};
959 :    
960 :     my @out = $figV->run_gathering_output(q(find_approx_neigh), $orgdir, $num_neighbors);
961 :     open(TMP, qq(>$outfile))
962 :     || die qq(could not write-open \'$outfile\');
963 :     print TMP @out;
964 :     close(TMP);
965 :     print STDERR (@out, qq(\n)) if $ENV{VERBOSE};
966 :    
967 :     my ($closest, $genus, $species) = ($out[0] =~ m/^(\S+)\t\d+\t(\S+)\s+(\S+)/o);
968 :     my $result = $sapO->genome_domain(-ids => [$closest]);
969 :     my $domain = $result->{$closest};
970 :    
971 :     return ($domain, $genus, $species);
972 :     }
973 :    
974 :    
975 :    
976 :     sub correct_frameshifts {
977 :     my ($keep,$fix_fs, $meta,$figV,$sapO,$annoO, $chkD,$tmpD,$restartD,$errD) = @_;
978 :    
979 :     if ($keep) {
980 :     #...Annotate FS only...
981 :     }
982 :     else {
983 :     if ($fix_fs || $meta->get_metadata("correction.frameshifts")) {
984 :     if (! $meta->get_metadata('status.rp.correct_frameshifts')) {
985 :     &my_run("correct_frameshifts -code=$code $tmpD > $errD/correct_frameshifts.stderr 2>&1");
986 :    
987 :     &make_restart($chkD, $tmpD, $restartD, qq(correct_frameshifts));
988 :     $meta->set_metadata('status.rp.correct_frameshifts',1);
989 :     }
990 :     }
991 :     }
992 :    
993 :     return 1;
994 :     }
995 :    
996 :    
997 :     sub find_and_backfill_missing_and_miscalled {
998 :     my ($keep,$backfill, $meta,$figV,$sapO,$annoO, $taxon_ID,$code, $chkD,$tmpD,$restartD,$errD) = @_;
999 :     #...Look for missing genes and miscalled genes
1000 :    
1001 :     if ($keep || (not $backfill)) {
1002 :     return 1;
1003 :     }
1004 :     else {
1005 :     if (not $meta->get_metadata("correction.backfill_gaps")) {
1006 :     return 1;
1007 :     }
1008 :     else {
1009 :     mkdir("$tmpD/CorrToReferenceGenomes")
1010 :     || die qq(Could not create $tmpD/CorrToReferenceGenomes);
1011 :    
1012 :     my $missing_dir = qq($tmpD/Missing_Genes);
1013 :     mkdir($missing_dir) || die qq(Could not create dir $missing_dir/);
1014 :     my @closest_genomes = map { m/^(\d+\.\d+)/o ? ($1) : () } &FIG::file_read(qq($tmpD/closest.genomes));
1015 :     foreach my $ref_genome (@closest_genomes) {
1016 :     my $corr_table = qq($tmpD/CorrToReferenceGenomes/$ref_genome);
1017 :     &my_run("svr_corresponding_genes -d $tmpD $taxonID $ref_genome > $corr_table");
1018 :    
1019 :     #
1020 :     # It is possible that corresponding genes generates zero-length output
1021 :     # in the event that there is asynchrony between the servers and the local
1022 :     # seed. Just skip them if that is the case.
1023 :     #
1024 :     if (-s $corr_table) {
1025 :     &my_run("find_missing_genes --orgdir $tmpD --corr $corr_table --ref $ref_genome > $missing_dir/$ref_genome 2> $missing_dir/$ref_genome\.err");
1026 :     }
1027 :     else {
1028 :     print STDERR "Zero-length correspondence table generated for $ref_genome\n";
1029 :     }
1030 :     }
1031 :     &my_run("merge_missing_gene_output $tmpD $missing_dir > $tmpD/missing_genes.out 2> $errD/missing_genes.err");
1032 :    
1033 :     my @missing_pegs = map { chomp; $_ =~ s/\%0A/\n/go; $_ } &FIG::file_read(qq($tmpD/missing_genes.out));
1034 :    
1035 :     foreach my $missing (@missing_pegs) {
1036 :     my ($loc, $frameshift, $length, $adjacent, $ref_pegs, $template_peg, $translation, $FS_evidence) = split(qq(\t), $missing);
1037 :     $loc =~ s/$taxonID\://g;
1038 :     my $new_fid = $figV->add_feature(q(rast), $taxonID, q(peg), $loc, undef, $translation);
1039 :     if (not $new_fid) {
1040 :     die qq(Could not add new PEG at loc=$loc, translation=$translation);
1041 :     }
1042 :     else {
1043 :     my $annot = (q(Gene prediction based on template PEG ) . $template_peg . qq(\n)
1044 :     . q(with support from PEGs) . $ref_pegs . qq(\n)
1045 :     );
1046 :    
1047 :     $figV->add_annotation($new_fid, q(rast), $annot);
1048 :    
1049 :     if ($frameshift) {
1050 :     $FS_evidence ||= qq(Sequence contains a probable frameshift);
1051 :     $figV->add_annotation($new_fid, q(rast), $FS_evidence);
1052 :     }
1053 :     }
1054 :     }
1055 :     }
1056 :    
1057 :     if (!-s "$tmpD/closest.genomes") {
1058 :     &log($tmpD, qq(Could not find any nearby genomes --- skipping \'backfill_gaps\'));
1059 :     }
1060 :     else {
1061 :     my $peg_tbl = "$tmpD/Features/peg/tbl";
1062 :     my $rna_tbl = (-s "$tmpD/Features/rna/tbl") ? "$tmpD/Features/rna/tbl" : "";
1063 :     my $extra_tbl = "$tmpD/Features/peg/tbl.extra";
1064 :    
1065 :     if (!-s "$tmpD/closest.genomes") {
1066 :     &log($tmpD, "Could not find any nearby genomes --- skipping \'backfill_gaps\'");
1067 :     }
1068 :     else {
1069 :     my $peg_tbl = "$tmpD/Features/peg/tbl";
1070 :     my $rna_tbl = (-s "$tmpD/Features/rna/tbl") ? "$tmpD/Features/rna/tbl" : "";
1071 :     my $extra_tbl = "$tmpD/Features/peg/tbl.extra";
1072 :    
1073 :     &my_run("backfill_gaps -orgdir=$tmpD -genetic_code=$code $tmpD/closest.genomes $taxonID $tmpD/contigs $rna_tbl $peg_tbl > $extra_tbl 2> $errD/backfill_gaps.stderr");
1074 :     if (-s $extra_tbl) {
1075 :     open( CALLED_BY, ">>$tmpD/called_by") || die "Could not append-open $tmpD/called_by";
1076 :     print CALLED_BY map { m/^(\S+)/ ? qq($1\tbackfill_gaps\n) : qq() } &FIG::file_read($extra_tbl);
1077 :     close(CALLED_BY);
1078 :     &my_run("cat $extra_tbl >> $peg_tbl");
1079 :     &my_run("get_fasta_for_tbl_entries -code=$code $tmpD/contigs < $extra_tbl >> $tmpD/Features/peg/fasta");
1080 :     system("rm -f $extra_tbl") && warn "WARNING: Could not remove $extra_tbl";
1081 :     }
1082 :     }
1083 :    
1084 :     &make_restart($chkD, $tmpD, $restartD, qq(backfill_gaps));
1085 :     $meta->set_metadata('status.rp.backfill_gaps',1);
1086 :     }
1087 :     }
1088 :    
1089 :     return 1;
1090 :     }
1091 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3