source: branches/lib/GDE/SATIVA/arb_sativa.pl

Last change on this file was 14663, checked in by akozlov, 10 years ago
  • Property svn:executable set to *
File size: 10.2 KB
Line 
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#         2015 HITS gGmbH / Exelixis Lab                                    #
10#         http://h-its.org  http://www.exelixis-lab.org/                    #
11#                                                                           #
12#############################################################################
13
14use strict;
15use warnings;
16# use diagnostics;
17
18BEGIN {
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
26use ARB;
27use tools;
28
29my $sativa_home = $ENV{'ARBHOME'}.'/lib/sativa';
30my ($start_time, $trainer_time, $mislabels_time);
31
32# ------------------------------------------------------------
33
34sub 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
43sub 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 $gb_field = ARB::entry($gb_species, $tax_field);
63           if (not defined $gb_field) {     
64             print "WARNING: field ".$tax_field." is missing for sequence ".$acc_no."\n";
65             next;
66           }     
67           my $tax = BIO::read_string($gb_species, $tax_field);
68           $tax || expectError('read_string '.$tax_field);
69           my $lineage = $tax;
70           if (substr($lineage, -1) eq ';') {
71             chop($lineage);
72           }
73           if ($sp_field) {
74             my $species_name = BIO::read_string($gb_species, $sp_field);
75             if (defined $species_name) {
76               $lineage = $lineage.";".$species_name;
77             }
78           }
79           
80           print TAXFILE $acc_no."\t".$lineage."\n";
81
82#           my $gb_ali = expectEntry($gb_species, 'ali_16s');
83#           my $gb_data = expectEntry($gb_ali, 'data');
84#           my $seq = ARB::read_string($gb_data);
85
86#           print SEQFILE ">".$acc_no."\n";
87#           print SEQFILE $seq."\n";
88
89           push(@marklist_r, $acc_no);
90         }
91       }
92
93  ARB::commit_transaction($gb_main);
94  ARB::close($gb_main);
95
96  close(TAXFILE);
97#  close(SEQFILE);
98}
99
100# ------------------------------------------------------------
101
102sub cpu_has_feature($) {
103    my ($feature) = @_;
104    my $cmd;
105    my $uname = `uname`;
106    chomp($uname);
107
108    if ($uname eq "Darwin") {
109       $cmd = "sysctl machdep.cpu.features";
110    }
111    elsif ($uname eq "Linux") {
112       $cmd = "grep flags /proc/cpuinfo";
113    }
114    else {
115       return 0;
116    }
117
118    `$cmd | grep -qi "$feature"`;
119    return $? + 1;
120}
121
122sub cpu_get_cores() {
123    my $cores = 2;
124    my $uname = `uname`;
125    chomp($uname);
126    if ($uname eq "Darwin") {
127       $cores = `sysctl -n hw.ncpu`;
128    }
129    elsif ($uname eq "Linux") {
130       $cores = `grep -c "^processor" /proc/cpuinfo`;
131    }
132
133    return $cores;
134}
135
136sub runPythonPipeline($$$$$$$$$$) {
137  my ($seq_file,$tax_file,$tax_code,$raxml_method,$num_reps,$conf_cutoff,$rand_seed,$rank_test,$verbose,$cores) = @_;
138
139  my $refjson_file = 'arb_export.json';
140  my $cfg_file = $sativa_home.'/epac.cfg.sativa';
141
142  # auto-configure RAxML
143  if ($cores == 0) { $cores = cpu_get_cores(); }
144
145  my $bindir = $ENV{'ARBHOME'}."/bin";
146  my $tmpdir = `pwd`;
147  chomp($tmpdir);
148
149  my $stime = time;
150
151  my $sativa_cmd = $sativa_home.'/sativa.py';
152  my @args = ($sativa_cmd, '-t', $tax_file, '-s', $seq_file, '-c', $cfg_file, '-T', $cores, '-tmpdir', $tmpdir, 
153               '-N', $num_reps, '-x', $tax_code, '-m', $raxml_method, '-p', $rand_seed, '-C', $conf_cutoff);
154  if ($rank_test == 1) { push(@args, "-ranktest"); }
155  if ($verbose == 1) { push(@args, "-v"); }
156  system(@args);
157
158#  $trainer_time = time - $stime;
159}
160
161# ------------------------------------------------------------
162
163sub deleteIfExists($$) {
164  my ($father,$key) = @_;
165
166  my $gb_field = ARB::entry($father,$key);
167  if (defined $gb_field) { 
168    my $error = ARB::delete($gb_field);
169    if ($error) { die $error; }
170  }
171}
172
173sub importResults($$$$) {
174  my ($res_file,$marked_only,$mark_misplaced,$field_suffix) = @_;
175
176  open(FILE,'<'.$res_file) || die "can't open '$res_file' (Reason: $!)";
177
178  my %mis_map = ();
179
180  # skip the header
181  my $header = <FILE>;
182
183  while (my $line = <FILE>) {
184    chomp $line;
185    my @fields = split "\t" , $line;
186    my $seq_id = $fields[0];
187    my $lvl = $fields[1];
188    my $conf = $fields[4];
189    my $new_tax = $fields[6];
190    $mis_map{$seq_id} = {conf => $conf, new_tax => $new_tax, mis_lvl => $lvl};
191    if (scalar(@fields) == 9) { $mis_map{$seq_id}{'rank_conf'} = $fields[8]; }
192  }
193  close(FILE);
194
195  my $gb_main = ARB::open(":","rw");
196  $gb_main || expectError('db connect (no running ARB?)');
197
198  dieOnError(ARB::begin_transaction($gb_main), 'begin_transaction');
199
200  my $count = 0;
201  for (my $gb_species = BIO::first_species($gb_main);
202       $gb_species;
203       $gb_species = BIO::next_species($gb_species))
204       {
205         my $error;
206         my $name = BIO::read_string($gb_species, "name");         
207         $name || expectError('read_string name');
208         my $marked = ARB::read_flag($gb_species);
209         if (defined $mis_map{$name}) {
210           $error = BIO::write_int($gb_species, "sativa_mislabel_flag".$field_suffix, 1);
211           if ($error) { die $error; }
212           $error = BIO::write_string($gb_species, "sativa_tax".$field_suffix, $mis_map{$name}{'new_tax'});
213           if ($error) { die $error; }
214           $error = BIO::write_float($gb_species, "sativa_seqmis_conf".$field_suffix, $mis_map{$name}{'conf'});
215           if ($error) { die $error; }
216           $error = BIO::write_string($gb_species, "sativa_mislabel_level".$field_suffix, $mis_map{$name}{'mis_lvl'});
217           if ($error) { die $error; }
218           if (exists $mis_map{$name}{'rank_conf'}) {
219               $error = BIO::write_string($gb_species, "sativa_rankmis_conf".$field_suffix, $mis_map{$name}{'rank_conf'});
220               if ($error) { die $error; }
221           }
222           
223           if ($mark_misplaced==1) { ARB::write_flag($gb_species, 1); }
224           $count++; 
225         }
226         elsif ($marked==1 || $marked_only==0) {
227           $error = BIO::write_int($gb_species, "sativa_mislabel_flag".$field_suffix, 0);
228           if ($error) { die $error; }
229 
230#           deleteIfExists($gb_species, "sativa_tax");
231#           deleteIfExists($gb_species, "sativa_conf");
232#           deleteIfExists($gb_species, "sativa_seqmis_conf");
233#           deleteIfExists($gb_species, "sativa_mislabel_level");
234#           deleteIfExists($gb_species, "sativa_misrank_conf");
235
236           if ($mark_misplaced==1) { ARB::write_flag($gb_species, 0); }
237         }
238       }
239
240  print "\nMislabels found: $count\n\n";
241  ARB::commit_transaction($gb_main);
242  ARB::close($gb_main);
243}
244
245
246# ------------------------------------------------------------
247
248sub die_usage($) {
249  my ($err) = @_;
250  print "Purpose: Run SATIVA taxonomy validation pipeline\n";
251  print "and import results back into ARB\n";
252  print "Usage: arb_sativa.pl [--marked-only] [--mark-misplaced] [--rank-test] [--verbose] tax_field tax_code num_reps raxml_method conf_cutoff rand_seed cores species_field\n";
253  print "       tax_field         Field contatining full (original) taxonomic path (lineage)\n";
254  print "       species_field     Field containing species name\n";
255  print "       --marked-only     Process only species that are marked in ARB (default: process all)\n";
256  print "       --mark-misplaced  Mark putatively misplaced species in ARB upon completion\n";
257  print "       --rank-test       Test for misplaced higher ranks\n";
258  print "\n";
259  die "Error: $err\n";
260}
261
262sub main() {
263  my $args = scalar(@ARGV);
264  if ($args<1) { die_usage('Missing arguments'); }
265
266  my $tax_file      = 'sativa_in.tax';
267  my $seq_file      = 'sativa_in.phy';
268  my $res_file      = 'sativa_in.mis';
269
270  my $marked_only = 0;
271  my $mark_misplaced = 0;
272  my $rank_test = 0;
273  my $verbose = 0;
274
275  while (substr($ARGV[0],0,2) eq '--') {
276    my $arg = shift @ARGV;
277    if ($arg eq '--marked-only') { $marked_only = 1; }
278    elsif ($arg eq '--mark-misplaced') { $mark_misplaced = 1; }
279    elsif ($arg eq '--rank-test') { $rank_test = 1; }
280    elsif ($arg eq '--verbose') { $verbose = 1; }
281    else { die_usage("Unknown switch '$arg'"); }
282    $args--;
283  }
284
285  my $tax_field  = shift @ARGV;
286  my $tax_code = shift @ARGV;
287  my $num_reps = shift @ARGV;
288  my $raxml_method = shift @ARGV;
289  my $conf_cutoff = shift @ARGV;
290  my $rand_seed = shift @ARGV;
291  my $cores = shift @ARGV;
292  my $species_field = shift @ARGV;
293 
294  my $field_suffix = $tax_field;
295  $field_suffix =~ s/^tax\_//g;
296  $field_suffix = "_".$field_suffix;
297 
298  my @marklist = {};
299 
300  $start_time = time;
301
302  exportTaxonomy($seq_file,$tax_file,$tax_field,$species_field,$marked_only,@marklist);
303
304  runPythonPipeline($seq_file,$tax_file,$tax_code,$num_reps,$raxml_method,$conf_cutoff,$rand_seed,$rank_test,$verbose,$cores);
305
306  importResults($res_file,$marked_only,$mark_misplaced,$field_suffix);
307 
308  my $total_time = time - $start_time;
309#  print "Elapsed time: $total_time s (trainer: $trainer_time s, find_mislabels: $mislabels_time s)\n\n";
310 
311  print "Done!\n\n";
312}
313
314main();
315
Note: See TracBrowser for help on using the repository browser.