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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3