[Bio] / FortyEightMeta / mg_postproc_taxa_sims.pl Repository:
ViewVC logotype

Annotation of /FortyEightMeta/mg_postproc_taxa_sims.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1
2 :     #
3 :     # Postprocess the metagenome LSU/SSU/Greengene sims runs.
4 :     #
5 :    
6 :     use DB_File;
7 :     use Data::Dumper;
8 :    
9 : olson 1.13 use JobStage;
10 : olson 1.1 use strict;
11 :     use FIG;
12 :     use FIGV;
13 :     use FIG_Config;
14 :     use File::Basename;
15 :     use File::Copy;
16 :     use GenomeMeta;
17 :     use Carp 'croak';
18 :     use FortyEightMeta::ExpandSims;
19 :    
20 : olson 1.13 my $STAGE = "sims_postprocess";
21 :    
22 : olson 1.1 @ARGV == 1 or die "Usage: $0 job-dir\n";
23 :    
24 : olson 1.13 my $jobdir = shift;
25 : olson 1.1
26 : olson 1.13 my $stage = new JobStage('Job48', $STAGE, $jobdir);
27 :     $stage or die "Cannot create job for $jobdir\n";
28 : olson 1.1
29 : olson 1.13 my $job_id = basename($jobdir);
30 :     my $job = $stage->job();
31 : olson 1.1
32 : olson 1.13 my $meta = $job->meta;
33 : olson 1.1
34 : olson 1.13 print "Running job! $jobdir\n";
35 : olson 1.1
36 : olson 1.13 my $blast2gff_exe = "$FIG_Config::bin/blast2gff";
37 :     -x $blast2gff_exe or $stage->fatal("Executable $blast2gff_exe missing");
38 : olson 1.1
39 : olson 1.13 my $correct_contigs_exe = "$FIG_Config::bin/correct_gff_missing_contigs";
40 :     -x $correct_contigs_exe or $stage->fatal("Executable $correct_contigs_exe missing");
41 : olson 1.1
42 : olson 1.13 my $gff2seed_exe = "$FIG_Config::bin/gff2seed";
43 :     -x $gff2seed_exe or $stage->fatal("Executable $gff2seed_exe missing");
44 : olson 1.1
45 : olson 1.13 my $get_ss_conns_exe = "$FIG_Config::bin/get_ss_connections";
46 :     -x $get_ss_conns_exe or $stage->fatal("Executable $get_ss_conns_exe missing");
47 : olson 1.1
48 : olson 1.13 $stage->set_status("in_progress");
49 : olson 1.1
50 :     &process_rdp();
51 :     &create_seed_gff();
52 :     &install_seed_gff();
53 :    
54 :     $meta->set_metadata("status.sims_postprocess", "complete");
55 :     $meta->set_metadata("sims_postprocess.running", "no");
56 :    
57 :     exit(0);
58 :    
59 :     #
60 :     # Process the seed blastx results into an organism directory.
61 :     #
62 :     sub create_seed_gff
63 :     {
64 :    
65 :     #
66 :     # consolidate all the blastx generated sims into a single file
67 :     #
68 :    
69 :     my $sims = "$jobdir/proc/sims.seed.raw";
70 : mkubal 1.12 #if (! -e $sims){
71 :     # consolidate_sims("$jobdir/proc/sims.seed", $sims);
72 :     #}
73 : mkubal 1.11
74 : olson 1.1 my $gff = "$jobdir/proc/$genome.gff";
75 :     my $out_sims = "$jobdir/proc/$genome.sims";
76 :    
77 :     #
78 : olson 1.5 # Use the btree-indexed NR from the data directory if possible.
79 :     #
80 :    
81 : olson 1.6 my $nr = "$FIG_Config::metagenome_data/nr.2007-0208.btree";
82 :     my $nrlen = "$FIG_Config::metagenome_data/nr.2007-0208.len.btree";
83 : olson 1.5 my @nrarg;
84 :     if (-f $nr)
85 :     {
86 :     @nrarg = ('-nr', $nr);
87 :     }
88 : olson 1.6 if (-f $nrlen)
89 :     {
90 :     push(@nrarg, "-nrlen", $nrlen);
91 :     }
92 : olson 1.5
93 :    
94 :     #
95 : olson 1.1 # Verify we have the inputs we need for generating the orgdir.
96 :     #
97 :    
98 :     my $fasta = $meta->get_metadata("preprocess.fasta_file");
99 : olson 1.13 ($fasta and -f $fasta) or $stage->fatal("fasta not found: '$fasta'");
100 : olson 1.1
101 :     my $taxonomy = &FIG::file_head("$jobdir/TAXONOMY", 1);
102 :     my $project = &FIG::file_head("$jobdir/PROJECT", 1);
103 :     chomp $taxonomy;
104 :     chomp $project;
105 :    
106 :     #
107 :     # Start the process.
108 :     #
109 :    
110 :     my @args = (-f => $fasta,
111 :     -s => $out_sims,
112 :     -o => $gff,
113 :     -gs => $job->genome_name(),
114 :     -taxonomy => $taxonomy,
115 :     -proj => $project,
116 :     -g => $genome,
117 : olson 1.5 @nrarg,
118 : olson 1.1 $sims);
119 :    
120 :     my $pid = open(BOUT, "-|");
121 :     if ($pid == 0)
122 :     {
123 :     exec { $blast2gff_exe } $blast2gff_exe, @args;
124 :     die "Failed: $!: $blast2gff_exe @args";
125 :     }
126 :    
127 :     print "Started pid=$pid $blast2gff_exe @args\n";
128 :    
129 :     while (<BOUT>)
130 :     {
131 :     print "Blast2gff: $_";
132 :     }
133 :    
134 :     if (!close(BOUT))
135 :     {
136 : olson 1.13 $stage->fatal("Failed with \$!=$! \$?=$?: $blast2gff_exe @args");
137 : olson 1.1 }
138 :     }
139 :    
140 :     #
141 :     # Unpack the gff file into the organism directory, do the rest of the post-install cleanup.
142 :    
143 :     #
144 :     sub install_seed_gff
145 :     {
146 :     my $genome = $job->genome_id();
147 :     my $gff = "$jobdir/proc/$genome.gff";
148 :     my $out_sims = "$jobdir/proc/$genome.sims";
149 :    
150 : olson 1.13 -f $gff or $stage->fatal("gff file $gff does not exist");
151 : olson 1.1
152 :     makedir("$jobdir/rp");
153 :     my $orgdir = "$jobdir/rp/$genome";
154 :    
155 :     my $fasta = $meta->get_metadata("preprocess.fasta_file");
156 : olson 1.13 ($fasta and -f $fasta) or $stage->fatal("fasta not found: '$fasta'");
157 : olson 1.1
158 :     #
159 :     # seed2gff won't run if the orgdir exists
160 :     #
161 :    
162 :     if (-d $orgdir)
163 :     {
164 :     my $new = "$orgdir." . time;
165 :     print "Moving $orgdir to $new\n";
166 : olson 1.13 rename $orgdir, $new or $stage->fatal("Cannot rename $orgdir to $new: $!");
167 : olson 1.1 }
168 :    
169 :     #
170 :     # Create seed dir
171 :     #
172 :     my @args = ("-debug", $gff, "$jobdir/rp");
173 :     &run_pipe(\*STDOUT, $gff2seed_exe, @args);
174 :    
175 :     #
176 :     # Copy in missing contigs
177 :     #
178 :    
179 :     &run_pipe(\*STDOUT, $correct_contigs_exe, "-orgdir", $orgdir, $fasta, $genome);
180 :    
181 :     #
182 : olson 1.9 # Index contigs.
183 :     #
184 :    
185 :     &run_pipe(\*STDOUT, "$FIG_Config::bin/make_fasta_btree",
186 :     "$orgdir/contigs",
187 :     "$orgdir/contigs.btree",
188 :     "$orgdir/contig_len.btree");
189 :    
190 :     #
191 : olson 1.1 # Copy sims.
192 :     #
193 :     copy($out_sims, "$orgdir/similarities");
194 :    
195 :     #
196 :     # Expand sims.
197 :     #
198 :    
199 : olson 1.7 my $peg_syns = "$FIG_Config::fortyeight_data/peg.synonyms";
200 : olson 1.8 my $exp_sims = "$orgdir/expanded_similarities";
201 : olson 1.7 &run_pipe(\*STDOUT, "$FIG_Config::bin/expand_sims",
202 :     "-orgdir", $orgdir,
203 :     $peg_syns,
204 : olson 1.8 "$orgdir/similarities", $exp_sims);
205 : olson 1.7
206 :     # my $fig = new FIGV($orgdir);
207 :    
208 : olson 1.1
209 :    
210 : olson 1.7 # my $syns_fh = new FileHandle($peg_syns);
211 : olson 1.13 # $syns_fh or $stage->fatal("cannot open $peg_syns: $!");
212 : olson 1.1
213 : olson 1.7 # my $syns = load_peg_syns($syns_fh);
214 : olson 1.1
215 : olson 1.7 # eval {
216 :     # expand_sims($fig, $out_sims, $exp_sims, $syns);
217 :     # };
218 :     # if ($@)
219 :     # {
220 : olson 1.13 # $stage->fatal("error in expand_sims: $@");
221 : olson 1.7 # }
222 :    
223 : olson 1.1
224 :     #
225 :     # Set up to compute subsys mappings.
226 :     #
227 :    
228 :     makedir("$orgdir/Subsystems");
229 :     my $bindings = "$orgdir/Subsystems/bindings";
230 :     my $subsystems = "$orgdir/Subsystems/subsystems";
231 :    
232 :     &touch($bindings);
233 :     &touch($subsystems);
234 :    
235 :     my $fh = new FileHandle(">$orgdir/ss.attributes");
236 : olson 1.13 $fh or $stage->fatal("Cannot write $orgdir/ss.attributes: $!");
237 : olson 1.1
238 :     &run_pipe($fh,
239 :     $get_ss_conns_exe,
240 :     "-figv", $bindings, $subsystems,
241 :     "-orgdir", $orgdir,
242 :     "-g", $genome);
243 :     close($fh);
244 :    
245 :     my $fh = new FileHandle(">$orgdir/taxa_summary_by_blast");
246 : olson 1.13 $fh or $stage->fatal("Cannot write $orgdir/taxa_summary_by_blast: $!");
247 : olson 1.1
248 :     &run_pipe($fh,
249 :     $summarize_exe,
250 :     "-orgdir", $orgdir,
251 : olson 1.9 "-b", "$orgdir/taxa_summary_by_blast.btree",
252 : olson 1.1 $genome);
253 :     close($fh);
254 : olson 1.4
255 :     #
256 :     # Copy in the taxa hits that we computed earlier.
257 :     #
258 :     system("cp $jobdir/proc/taxa*hits $orgdir/.");
259 : olson 1.1 }
260 :    
261 :     #
262 :     # If we are running here, we should have generated sims in the proc/sims.{16s,gg,lsu,ssu} directories.
263 :     # Scan the task.list files in each of them, and concatenate the output into files proc/sims.WHAT.raw.
264 :     #
265 :     sub process_rdp
266 :     {
267 : olson 1.3 my $orgdir = "$jobdir/rp/$genome";
268 :    
269 : olson 1.1 for my $name (qw(16s gg lsu ssu))
270 :     {
271 : mkubal 1.11 my $dir = "$jobdir/proc/sims.$name";
272 : olson 1.1 my $sims = "$jobdir/proc/sims.$name.raw";
273 :    
274 : mkubal 1.12 #if (! -e $sims){
275 :     # consolidate_sims($dir, $sims);
276 :     #}
277 : olson 1.1 # Now we can invoke the appropriate analysis.
278 : mkubal 1.11
279 : olson 1.3
280 : olson 1.4 my $besthits = "$jobdir/proc/taxa.$name.besthits";
281 :     my $allhits = "$jobdir/proc/taxa.$name.allhits";
282 : olson 1.1
283 :     if ($name eq 'gg')
284 :     {
285 : olson 1.3 my @cmd = ("$FIG_Config::bin/blast2taxa", $sims, $besthits, $allhits);
286 : olson 1.2 my $rc = system(@cmd);
287 : olson 1.1 if ($rc != 0)
288 :     {
289 : olson 1.13 $stage->fatal("error running $FIG_Config::bin/blast2taxa $sims");
290 : olson 1.1 }
291 :     }
292 :     elsif ($name eq 'lsu' or $name eq 'ssu' or $name eq '16s')
293 :     {
294 :     my $db = $name;
295 :     $db = "rdp" if $name eq '16s';
296 :    
297 :     my @cmd = ("$FIG_Config::bin/16s_taxa_mysql",
298 :     '-d', $db,
299 :     '-b', $sims,
300 :     '-c', '1e-5',
301 :     '-l', '50',
302 : olson 1.3 '-ah', $allhits,
303 :     '-bh', $besthits);
304 : olson 1.1 print "Run @cmd\n";
305 :     my $rc = system(@cmd);
306 :     if ($rc != 0)
307 :     {
308 : olson 1.13 $stage->fatal("error rc=$rc running @cmd");
309 : olson 1.1 }
310 :     }
311 :     }
312 :     }
313 :    
314 :     sub consolidate_sims
315 :     {
316 :     my($dir, $sims) = @_;
317 :     if (! -d $dir)
318 :     {
319 : olson 1.13 $stage->fatal("Sim directory $dir not found");
320 : olson 1.1 }
321 :    
322 : olson 1.13 open(TL, "<$dir/task.list") or $stage->fatal("Cannot open $dir/task.list: $!");
323 : olson 1.1
324 :     #
325 :     # Open output file.
326 :     #
327 : olson 1.13 open(O, ">$sims") or $stage->fatal("Cannot create $sims: $!");
328 : olson 1.1
329 :     while (<TL>)
330 :     {
331 :     chomp;
332 :     my ($id, $in, $nr, $flags, $out, $err) = split(/\t/);
333 :    
334 :     if (! -f $out)
335 :     {
336 : olson 1.13 $stage->fatal("Sims output file $out missing for sims dir $dir");
337 : olson 1.1 }
338 :    
339 :     copy($out, \*O);
340 :     }
341 :     close(O);
342 :     }
343 :    
344 :     sub run_pipe
345 :     {
346 :     my($outfh, $cmd, @args) = @_;
347 :    
348 :     my $pid = open(BOUT, "-|");
349 :     if ($pid == 0)
350 :     {
351 :     exec { $cmd } $cmd, @args;
352 :     die "Failed: $!: $cmd @args";
353 :     }
354 :    
355 :     print "Started pid=$pid $cmd @args\n";
356 :    
357 :     while (<BOUT>)
358 :     {
359 :     print $outfh $_;
360 :     }
361 :    
362 :     if (!close(BOUT))
363 :     {
364 : olson 1.13 $stage->fatal("Failed with \$!=$! \$?=$?: $cmd @args");
365 : olson 1.1 }
366 :     }
367 :    
368 :     sub makedir
369 :     {
370 :     my($dir) = @_;
371 :    
372 :     if (! -d $dir)
373 :     {
374 : olson 1.13 mkdir($dir, 0777) or $stage->fatal("mkdir $dir failed: $!");
375 : olson 1.1 }
376 :     chmod 0777, $dir;
377 :     }
378 :    
379 :     sub touch
380 :     {
381 :     my($file) = @_;
382 : olson 1.13 my $fh = new FileHandle(">$file") or $stage->fatal("touch: cannot open $file for writing: $!");
383 : olson 1.1 close($fh);
384 :     }
385 :    
386 :    

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3