1 | #!/usr/bin/perl |
---|
2 | |
---|
3 | ############################################################################# |
---|
4 | # # |
---|
5 | # File: arb_sativa.pl # |
---|
6 | # Purpose: Adapter script for SATIVA taxonomy validation pipeline # |
---|
7 | # # |
---|
8 | # Author: Alexey Kozlov (alexey.kozlov@h-its.org) # |
---|
9 | # 2014 HITS gGmbH / Exelixis Lab # |
---|
10 | # http://h-its.org http://www.exelixis-lab.org/ # |
---|
11 | # # |
---|
12 | ############################################################################# |
---|
13 | |
---|
14 | use strict; |
---|
15 | use warnings; |
---|
16 | # use diagnostics; |
---|
17 | |
---|
18 | BEGIN { |
---|
19 | if (not exists $ENV{'ARBHOME'}) { die "Environment variable \$ARBHOME has to be defined"; } |
---|
20 | my $arbhome = $ENV{'ARBHOME'}; |
---|
21 | push @INC, "$arbhome/lib"; |
---|
22 | push @INC, "$arbhome/PERL_SCRIPTS/lib"; |
---|
23 | 1; |
---|
24 | } |
---|
25 | |
---|
26 | use ARB; |
---|
27 | use tools; |
---|
28 | |
---|
29 | my $sativa_home = $ENV{'ARBHOME'}.'/bin/sativa'; |
---|
30 | my ($start_time, $trainer_time, $mislabels_time); |
---|
31 | |
---|
32 | # ------------------------------------------------------------ |
---|
33 | |
---|
34 | sub expectEntry($$) { |
---|
35 | my ($gb_father,$key) = @_; |
---|
36 | my $gb = ARB::search($gb_father, $key, 'NONE'); |
---|
37 | defined $gb || die "Expected entry '$key' in '".ARB::get_db_path($gb_father)."'"; |
---|
38 | return $gb; |
---|
39 | } |
---|
40 | |
---|
41 | # ----------------------- |
---|
42 | |
---|
43 | sub exportTaxonomy($$$$$\@) { |
---|
44 | my ($seq_file,$tax_file,$tax_field,$sp_field,$marked_only,@marklist_r) = @_; |
---|
45 | |
---|
46 | open(TAXFILE,'>'.$tax_file) || die "can't open '$tax_file' (Reason: $!)"; |
---|
47 | # open(SEQFILE,'>'.$seq_file) || die "can't open '$seq_file' (Reason: $!)"; |
---|
48 | |
---|
49 | my $gb_main = ARB::open(":","r"); |
---|
50 | $gb_main || expectError('db connect (no running ARB?)'); |
---|
51 | |
---|
52 | dieOnError(ARB::begin_transaction($gb_main), 'begin_transaction'); |
---|
53 | |
---|
54 | for (my $gb_species = BIO::first_species($gb_main); |
---|
55 | $gb_species; |
---|
56 | $gb_species = BIO::next_species($gb_species)) |
---|
57 | { |
---|
58 | my $marked = ARB::read_flag($gb_species); |
---|
59 | if ($marked==1 || $marked_only==0) { |
---|
60 | my $acc_no = BIO::read_string($gb_species, "name"); |
---|
61 | $acc_no || expectError('read_string acc'); |
---|
62 | my $tax = BIO::read_string($gb_species, $tax_field); |
---|
63 | $tax || expectError('read_string '.$tax_field); |
---|
64 | my $species_name = BIO::read_string($gb_species, $sp_field); |
---|
65 | my $lineage = $tax; |
---|
66 | if (substr($lineage, -1) eq ';') { |
---|
67 | chop($lineage); |
---|
68 | } |
---|
69 | if (defined $species_name) { |
---|
70 | $lineage = $lineage.";".$species_name; |
---|
71 | } |
---|
72 | |
---|
73 | print TAXFILE $acc_no."\t".$lineage."\n"; |
---|
74 | |
---|
75 | # my $gb_ali = expectEntry($gb_species, 'ali_16s'); |
---|
76 | # my $gb_data = expectEntry($gb_ali, 'data'); |
---|
77 | # my $seq = ARB::read_string($gb_data); |
---|
78 | |
---|
79 | # print SEQFILE ">".$acc_no."\n"; |
---|
80 | # print SEQFILE $seq."\n"; |
---|
81 | |
---|
82 | push(@marklist_r, $acc_no); |
---|
83 | } |
---|
84 | } |
---|
85 | |
---|
86 | ARB::commit_transaction($gb_main); |
---|
87 | ARB::close($gb_main); |
---|
88 | |
---|
89 | close(TAXFILE); |
---|
90 | # close(SEQFILE); |
---|
91 | } |
---|
92 | |
---|
93 | # ------------------------------------------------------------ |
---|
94 | |
---|
95 | sub cpu_has_feature($) { |
---|
96 | my ($feature) = @_; |
---|
97 | my $cmd; |
---|
98 | my $uname = `uname`; |
---|
99 | chomp($uname); |
---|
100 | |
---|
101 | if ($uname eq "Darwin") { |
---|
102 | $cmd = "sysctl machdep.cpu.features"; |
---|
103 | } |
---|
104 | elsif ($uname eq "Linux") { |
---|
105 | $cmd = "grep flags /proc/cpuinfo"; |
---|
106 | } |
---|
107 | else { |
---|
108 | return 0; |
---|
109 | } |
---|
110 | |
---|
111 | `$cmd | grep -qi "$feature"`; |
---|
112 | return $? + 1; |
---|
113 | } |
---|
114 | |
---|
115 | sub cpu_get_cores() { |
---|
116 | my $cores = 2; |
---|
117 | my $uname = `uname`; |
---|
118 | chomp($uname); |
---|
119 | if ($uname eq "Darwin") { |
---|
120 | $cores = `sysctl -n hw.ncpu`; |
---|
121 | } |
---|
122 | elsif ($uname eq "Linux") { |
---|
123 | $cores = `grep -c "^processor" /proc/cpuinfo`; |
---|
124 | } |
---|
125 | |
---|
126 | return $cores; |
---|
127 | } |
---|
128 | |
---|
129 | sub runPythonPipeline($$$$$) { |
---|
130 | my ($seq_file,$tax_file,$dup_rank_names,$wrong_rank_count,$rank_test) = @_; |
---|
131 | |
---|
132 | my $refjson_file = 'arb_export.json'; |
---|
133 | my $cfg_file = $sativa_home.'/epac.cfg.sativa'; |
---|
134 | |
---|
135 | # auto-configure RAxML |
---|
136 | my $cores = cpu_get_cores(); |
---|
137 | |
---|
138 | my $bindir = $ENV{'ARBHOME'}."/bin"; |
---|
139 | my $tmpdir = `pwd`; |
---|
140 | chomp($tmpdir); |
---|
141 | |
---|
142 | my $stime = time; |
---|
143 | |
---|
144 | my $trainer_cmd = $sativa_home.'/epa_trainer.py'; |
---|
145 | system($trainer_cmd, '-t', $tax_file, '-s', $seq_file, '-r', $refjson_file, '-c', $cfg_file, '-T', $cores, '-no-hmmer', |
---|
146 | '-dup-rank-names', $dup_rank_names, '-wrong-rank-count', $wrong_rank_count, '-tmpdir', $tmpdir, '-C'); |
---|
147 | |
---|
148 | $trainer_time = time - $stime; |
---|
149 | |
---|
150 | if (-e $refjson_file) { |
---|
151 | $stime = time; |
---|
152 | my $checker_cmd = $sativa_home.'/find_mislabels.py'; |
---|
153 | my @args = ($checker_cmd, '-r', $refjson_file, '-c', $cfg_file, '-T', $cores, '-v', '-tmpdir', $tmpdir); |
---|
154 | if ($rank_test == 1) { push(@args, "-ranktest"); } |
---|
155 | system(@args); |
---|
156 | $mislabels_time = time - $stime; |
---|
157 | } |
---|
158 | else { |
---|
159 | die "There was a problem running SATIVA (see messages above). Exiting..." |
---|
160 | } |
---|
161 | } |
---|
162 | |
---|
163 | # ------------------------------------------------------------ |
---|
164 | |
---|
165 | sub deleteIfExists($$) { |
---|
166 | my ($father,$key) = @_; |
---|
167 | |
---|
168 | my $gb_field = ARB::entry($father,$key); |
---|
169 | if (defined $gb_field) { |
---|
170 | my $error = ARB::delete($gb_field); |
---|
171 | if ($error) { die $error; } |
---|
172 | } |
---|
173 | } |
---|
174 | |
---|
175 | sub importResults($$$) { |
---|
176 | my ($res_file,$marked_only,$mark_misplaced) = @_; |
---|
177 | |
---|
178 | open(FILE,'<'.$res_file) || die "can't open '$res_file' (Reason: $!)"; |
---|
179 | |
---|
180 | my %mis_map = (); |
---|
181 | |
---|
182 | # skip the header |
---|
183 | my $header = <FILE>; |
---|
184 | |
---|
185 | while (my $line = <FILE>) { |
---|
186 | chomp $line; |
---|
187 | my @fields = split "\t" , $line; |
---|
188 | my $seq_id = $fields[0]; |
---|
189 | my $lvl = $fields[1]; |
---|
190 | my $conf = $fields[4]; |
---|
191 | my $new_tax = $fields[6]; |
---|
192 | $mis_map{$seq_id} = {conf => $conf, new_tax => $new_tax, mis_lvl => $lvl}; |
---|
193 | if (scalar(@fields) == 9) { $mis_map{$seq_id}{'rank_conf'} = $fields[8]; } |
---|
194 | } |
---|
195 | close(FILE); |
---|
196 | |
---|
197 | my $gb_main = ARB::open(":","rw"); |
---|
198 | $gb_main || expectError('db connect (no running ARB?)'); |
---|
199 | |
---|
200 | dieOnError(ARB::begin_transaction($gb_main), 'begin_transaction'); |
---|
201 | |
---|
202 | my $count = 0; |
---|
203 | for (my $gb_species = BIO::first_species($gb_main); |
---|
204 | $gb_species; |
---|
205 | $gb_species = BIO::next_species($gb_species)) |
---|
206 | { |
---|
207 | my $error; |
---|
208 | my $name = BIO::read_string($gb_species, "name"); |
---|
209 | $name || expectError('read_string name'); |
---|
210 | my $marked = ARB::read_flag($gb_species); |
---|
211 | if (defined $mis_map{$name}) { |
---|
212 | $error = BIO::write_int($gb_species, "sativa_mislabel_flag", 1); |
---|
213 | if ($error) { die $error; } |
---|
214 | $error = BIO::write_string($gb_species, "sativa_tax", $mis_map{$name}{'new_tax'}); |
---|
215 | if ($error) { die $error; } |
---|
216 | $error = BIO::write_float($gb_species, "sativa_seqmis_conf", $mis_map{$name}{'conf'}); |
---|
217 | if ($error) { die $error; } |
---|
218 | $error = BIO::write_string($gb_species, "sativa_mislabel_level", $mis_map{$name}{'mis_lvl'}); |
---|
219 | if ($error) { die $error; } |
---|
220 | if (exists $mis_map{$name}{'rank_conf'}) { |
---|
221 | $error = BIO::write_string($gb_species, "sativa_rankmis_conf", $mis_map{$name}{'rank_conf'}); |
---|
222 | if ($error) { die $error; } |
---|
223 | } |
---|
224 | |
---|
225 | if ($mark_misplaced==1) { ARB::write_flag($gb_species, 1); } |
---|
226 | $count++; |
---|
227 | } |
---|
228 | elsif ($marked==1 || $marked_only==0) { |
---|
229 | $error = BIO::write_int($gb_species, "sativa_mislabel_flag", 0); |
---|
230 | if ($error) { die $error; } |
---|
231 | |
---|
232 | # deleteIfExists($gb_species, "sativa_tax"); |
---|
233 | # deleteIfExists($gb_species, "sativa_conf"); |
---|
234 | # deleteIfExists($gb_species, "sativa_seqmis_conf"); |
---|
235 | # deleteIfExists($gb_species, "sativa_mislabel_level"); |
---|
236 | # deleteIfExists($gb_species, "sativa_misrank_conf"); |
---|
237 | |
---|
238 | if ($mark_misplaced==1) { ARB::write_flag($gb_species, 0); } |
---|
239 | } |
---|
240 | } |
---|
241 | |
---|
242 | print "\nMislabels found: $count\n\n"; |
---|
243 | ARB::commit_transaction($gb_main); |
---|
244 | ARB::close($gb_main); |
---|
245 | } |
---|
246 | |
---|
247 | |
---|
248 | # ------------------------------------------------------------ |
---|
249 | |
---|
250 | sub die_usage($) { |
---|
251 | my ($err) = @_; |
---|
252 | print "Purpose: Run SATIVA taxonomy validation pipeline\n"; |
---|
253 | print "and import results back into ARB\n"; |
---|
254 | print "Usage: arb_sativa.pl [--marked-only] [--mark-misplaced] [--rank-test] tax_field [species_field]\n"; |
---|
255 | print " tax_field Field contatining full (original) taxonomic path (lineage)\n"; |
---|
256 | print " species_field Field containing species name\n"; |
---|
257 | print " --marked-only Process only species that are marked in ARB (default: process all)\n"; |
---|
258 | print " --mark-misplaced Mark putatively misplaced species in ARB upon completion\n"; |
---|
259 | print " --rank-test Test for misplaced higher ranks\n"; |
---|
260 | print "\n"; |
---|
261 | die "Error: $err\n"; |
---|
262 | } |
---|
263 | |
---|
264 | sub main() { |
---|
265 | my $args = scalar(@ARGV); |
---|
266 | if ($args<1) { die_usage('Missing arguments'); } |
---|
267 | |
---|
268 | my $tax_file = 'sativa_in.tax'; |
---|
269 | my $seq_file = 'sativa_in.phy'; |
---|
270 | my $res_file = 'result.mis'; |
---|
271 | |
---|
272 | my $marked_only = 0; |
---|
273 | my $mark_misplaced = 0; |
---|
274 | my $rank_test = 0; |
---|
275 | |
---|
276 | while (substr($ARGV[0],0,2) eq '--') { |
---|
277 | my $arg = shift @ARGV; |
---|
278 | if ($arg eq '--marked-only') { $marked_only = 1; } |
---|
279 | elsif ($arg eq '--mark-misplaced') { $mark_misplaced = 1; } |
---|
280 | elsif ($arg eq '--rank-test') { $rank_test = 1; } |
---|
281 | else { die_usage("Unknown switch '$arg'"); } |
---|
282 | $args--; |
---|
283 | } |
---|
284 | |
---|
285 | my $tax_field = shift @ARGV; |
---|
286 | my $species_field = shift @ARGV; |
---|
287 | my $dup_rank_names = shift @ARGV; |
---|
288 | my $wrong_rank_count = shift @ARGV; |
---|
289 | |
---|
290 | my @marklist = {}; |
---|
291 | |
---|
292 | $start_time = time; |
---|
293 | |
---|
294 | exportTaxonomy($seq_file,$tax_file,$tax_field,$species_field,$marked_only,@marklist); |
---|
295 | |
---|
296 | runPythonPipeline($seq_file,$tax_file,$dup_rank_names,$wrong_rank_count,$rank_test); |
---|
297 | |
---|
298 | importResults($res_file,$marked_only,$mark_misplaced); |
---|
299 | |
---|
300 | my $total_time = time - $start_time; |
---|
301 | print "Elapsed time: $total_time s (trainer: $trainer_time s, find_mislabels: $mislabels_time s)\n\n"; |
---|
302 | |
---|
303 | print "Done!\n\n"; |
---|
304 | } |
---|
305 | |
---|
306 | main(); |
---|
307 | |
---|