[Bio] / Babel / bin / sims2lca.pl Repository:
ViewVC logotype

Annotation of /Babel/bin/sims2lca.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : tharriso 1.1 #!/usr/bin/env perl
2 :    
3 :     use strict;
4 :     use warnings;
5 :    
6 :     use Data::Dumper;
7 :     use Getopt::Long;
8 :     use Cache::Memcached;
9 :    
10 :     my $verbose = 0;
11 :     my $simf = '';
12 :     my $outf = '';
13 :     my $eval_rng = 0.2;
14 :     my $fcache = '';
15 :     my $mcache = ''; # _ach
16 :     my $memhost = "140.221.76.21:11211";
17 :     my $usage = "$0 [--verbose] [--file_cache <md52lca file> || --mem_cache <key extension>] [--eval <evalue range: default '$eval_rng'>] --sims <sim file> --out <out file>\n";
18 :    
19 :     if ( (@ARGV > 0) && ($ARGV[0] =~ /-h/) ) { print STDERR $usage; exit 1; }
20 :     if ( ! GetOptions('verbose!' => \$verbose,
21 :     'sims=s' => \$simf,
22 :     'out=s' => \$outf,
23 :     'eval:f' => \$eval_rng,
24 :     'file_cache:s' => \$fcache,
25 :     'mem_cache:s' => \$mcache
26 :     ) ) {
27 :     print STDERR $usage; exit 1;
28 :     }
29 :     unless ($simf && $outf && (-s $simf)) { print STDERR "Missing input and/or output files:\n$usage"; exit 1; }
30 :    
31 :     my ($md5_mem, $md5_file);
32 :     print STDERR time . "\n";
33 :    
34 :     # load md5 lcas
35 :     if ($fcache && (-s $fcache)) {
36 :     $md5_file = {};
37 :     print STDERR "Loading mapping file ... " if ($verbose);
38 :     open(FILE, "<$fcache") || die "Can't open md5 to lca mapping file $fcache: $!\n";
39 :     while (my $line = <FILE>) {
40 :     chomp $line;
41 :     my ($md5, @taxa) = split(/\t/, $line);
42 :     my $rank = pop @taxa;
43 :     $md5_file->{$md5} = \@taxa;
44 :     }
45 :     close FILE;
46 :     print STDERR "Done\n" if ($verbose);
47 :     }
48 :     elsif ($mcache) {
49 :     $md5_mem = new Cache::Memcached {'servers' => [$memhost], 'debug' => 0, 'compress_threshold' => 10_000};
50 :     }
51 :     else {
52 :     print STDERR "Provide one of file_cache or mem_cache:\n$usage"; exit 1;
53 :     }
54 :    
55 :     my $curr = '';
56 :     my ($top_exp, @md5s, @idents, @lens, @evals);
57 :    
58 :     print STDERR time . "\n";
59 :     print STDERR "Reading sims file:\n" if ($verbose);
60 :     open(SIMF, "<$simf") || die "Can't open file $simf: $!\n";
61 :     open(OUTF, ">$outf") || die "Can't open file $outf: $!\n";
62 :     while (my $line = <SIMF>) {
63 :     chomp $line;
64 :     my ($fid,$md5,$ident,$len,undef,undef,undef,undef,undef,undef,$eval,$score) = split(/\t/, $line);
65 :    
66 :     if ($eval =~ /^(\d\.\d)e([-+])(\d+)$/) {
67 :     my ($int, $pos, $exp) = ($1, $2, $3);
68 :     my $is_zero = (($int eq '0.0') && ($exp eq '00')) ? 1 : 0;
69 :     $exp = $is_zero ? 0 : ($pos eq '-') ? -1 * $exp : $exp;
70 :     next if (($pos eq '+') && (! $is_zero));
71 :     next if ((! $fid) || (! $line));
72 :    
73 :     if (! defined $top_exp) { $top_exp = $exp; }
74 :     if ($curr eq '') { $curr = $fid; }
75 :     if ($curr ne $fid) {
76 :     # get lca for md5s
77 :     if (@md5s > 0) {
78 :     my @taxa = &taxa_for_md5s(\@md5s);
79 :     my @lca = &get_lca($curr, \@taxa);
80 :     if (@lca == 9) {
81 :     print OUTF join("\t", ($curr, join(";",@md5s), join(";",@idents), join(";",@lens), join(";",@evals), join(";",@lca))) . "\n";
82 :     }
83 :     }
84 :     # reset variables
85 :     $curr = $fid;
86 :     $top_exp = $exp;
87 :     @md5s = ();
88 :     @idents = ();
89 :     @lens = ();
90 :     @evals = ();
91 :     }
92 :     # add data
93 :     if (($top_exp == 0) || ($exp < $top_exp)) {
94 :     $top_exp = $exp;
95 :     }
96 :     elsif ($exp <= ($top_exp - ($top_exp * $eval_rng))) {
97 :     next;
98 :     }
99 :     push @md5s , $md5;
100 :     push @idents, $ident;
101 :     push @lens , $len;
102 :     push @evals , $exp;
103 :     }
104 :     }
105 :    
106 :     if (@md5s > 0) {
107 :     my @taxa = &taxa_for_md5s(\@md5s);
108 :     my @lca = &get_lca($curr, \@taxa);
109 :     if (@lca == 9) {
110 :     print OUTF join("\t", ($curr, join(",",@md5s), join(",",@idents), join(",",@lens), join(",",@evals), @lca)) . "\n";
111 :     }
112 :     }
113 :     close OUTF;
114 :     print STDERR "Done\n" if ($verbose);
115 :     print STDERR time . "\n";
116 :    
117 :     exit 0;
118 :    
119 :     sub taxa_for_md5s {
120 :     my ($md5s) = @_;
121 :    
122 :     my @taxa = ();
123 :     if ($md5_file && ref($md5_file)) {
124 :     @taxa = map { $md5_file->{$_} } grep { exists $md5_file->{$_} } @$md5s;
125 :     }
126 :     elsif ($md5_mem && ref($md5_mem)) {
127 :     my @keys = map { $_ . $mcache } @$md5s;
128 :     my $data = $md5_mem->get_multi(@keys);
129 :     if ($data && ref($data)) {
130 :     @taxa = map { $data->{$_}{lca} } grep { exists $data->{$_}{lca} } @keys;
131 :     }
132 :     }
133 :     return @taxa;
134 :     }
135 :    
136 :     sub get_lca {
137 :     my ($frag, $taxa) = @_;
138 :    
139 :     my $coverage = {};
140 :     foreach my $t (@$taxa) {
141 :     for (my $i = 0; $i < scalar(@$t); $i++) {
142 :     $coverage->{$i+1}->{ $t->[$i] }++ if ($t->[$i]);
143 :     }
144 :     }
145 :    
146 :     if ( scalar(keys %$coverage) < 8 ) {
147 :     print STDERR "Incomplete Taxonomy:\t$frag\t" . join("\t", map { join(";", @$_) } @$taxa) . "\n";
148 :     return ();
149 :     }
150 :     if ( scalar(keys %{$coverage->{1}}) > 1 ) {
151 :     print STDERR "No LCA possible:\t$frag\t" . join("\t", keys %{$coverage->{1}}) . "\n";
152 :     return ();
153 :     }
154 :    
155 :     my @lca = ();
156 :     my $pos = 0;
157 :     my $max = 0;
158 :    
159 :     foreach my $key (sort { $a <=> $b } keys %$coverage) {
160 :     my $num = scalar keys %{$coverage->{$key}};
161 :     if (($num <= $max) || ($max == 0)) {
162 :     $max = $num;
163 :     $pos = $key;
164 :     }
165 :     }
166 :    
167 :     if ( scalar(keys %{$coverage->{$pos}}) == 1 ) {
168 :     @lca = map { keys %{$coverage->{$_}} } (1 .. $pos);
169 :     push @lca, ( map {'-'} ($pos + 1 .. 8) ) if ($pos < 8);
170 :     push @lca, $pos;
171 :     }
172 :     else {
173 :     print STDERR "LCA error ($pos)\t$frag\n";
174 :     }
175 :    
176 :     return @lca;
177 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3