[Bio] / FigKernelPackages / ProtSims.pm Repository:
ViewVC logotype

Annotation of /FigKernelPackages/ProtSims.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1
2 : olson 1.10 # This is a SAS component
3 :    
4 : overbeek 1.1 #
5 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
6 :     # for Interpretations of Genomes. All Rights Reserved.
7 :     #
8 :     # This file is part of the SEED Toolkit.
9 :     #
10 :     # The SEED Toolkit is free software. You can redistribute
11 :     # it and/or modify it under the terms of the SEED Toolkit
12 :     # Public License.
13 :     #
14 :     # You should have received a copy of the SEED Toolkit Public License
15 :     # along with this program; if not write to the University of Chicago
16 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
17 :     # Genomes at veronika@thefig.info or download a copy from
18 :     # http://www.theseed.org/LICENSE.TXT.
19 :     #
20 :     package ProtSims;
21 :    
22 :     use strict;
23 : olson 1.4 use Sim;
24 : overbeek 1.1 use Data::Dumper;
25 :     use Carp;
26 :     use gjoseqlib;
27 : olson 1.15 use DB_File;
28 : overbeek 1.1
29 : olson 1.13 use SeedAware;
30 :    
31 :     my $blat_cmd = SeedAware::executable_for("blat");
32 :     my $blastall_cmd = SeedAware::executable_for("blastall");
33 :     my $formatdb_cmd = SeedAware::executable_for("formatdb");
34 : olson 1.10
35 : olson 1.17 #
36 :     # Need to call temps different names on each invocation -
37 :     # if a subprocess hung on windows, we get problems with
38 :     # "file in use" errors.
39 :     #
40 :     my $tmp_serial = 0;
41 :    
42 : overbeek 1.1 sub blastP {
43 : overbeek 1.12 my($q,$db,$min_hits,$use_blast) = @_;
44 : gdpusch 1.21
45 : overbeek 1.2 my(@q,@db);
46 : gdpusch 1.21
47 : olson 1.13 my $tmp_dir = SeedAware::location_of_tmp();
48 : olson 1.18 my $tmp_suffix = $$ . "." . $tmp_serial++ . "." . time;
49 : gdpusch 1.21
50 : overbeek 1.1 my $qF;
51 :     if (ref $q)
52 :     {
53 : gdpusch 1.21 $qF = "$tmp_dir/tmp.query.$tmp_suffix";
54 : overbeek 1.1 &gjoseqlib::print_alignment_as_fasta($qF,$q);
55 :     }
56 :     else
57 :     {
58 :     $qF = $q;
59 :     if (-e $qF)
60 :     {
61 : overbeek 1.2 @q = &gjoseqlib::read_fasta($qF);
62 : overbeek 1.1 $q = \@q;
63 :     }
64 :     else
65 :     {
66 :     die "$qF is missing";
67 :     }
68 :     }
69 : gdpusch 1.21
70 : olson 1.15 my $db_lengths = {};
71 : gdpusch 1.21
72 : overbeek 1.1 my $dbF;
73 :     if (ref $db)
74 :     {
75 : gdpusch 1.21 $dbF = "$tmp_dir/tmp.db.$tmp_suffix";
76 :     my $tmp_formatdb_logfile = "$tmp_dir/tmp.db.$tmp_suffix.formatdb.log";
77 :    
78 : overbeek 1.1 &gjoseqlib::print_alignment_as_fasta($dbF,$db);
79 : gdpusch 1.21 system($formatdb_cmd, '-i', $dbF, '-l', $tmp_formatdb_logfile);
80 :     unlink($tmp_formatdb_logfile) unless $ENV{DEBUG};
81 :    
82 : olson 1.15 $db_lengths->{$_->[0]} = length($_->[2]) for @$db;
83 : overbeek 1.1 }
84 :     else
85 :     {
86 :     $dbF = $db;
87 :     if (-e $dbF)
88 :     {
89 : gdpusch 1.21 my $formatdb_logfile = "$dbF.formatdb.log";
90 :    
91 : overbeek 1.1 if (! ((-e "$dbF.psq") && ((-M "$dbF.psq") < (-M $dbF))))
92 :     {
93 : gdpusch 1.21 system($formatdb_cmd, '-i', $dbF, '-l', $formatdb_logfile);
94 : overbeek 1.1 }
95 : gdpusch 1.21
96 : olson 1.15 my $db_tie;
97 :     if (-f "$dbF.lengths")
98 :     {
99 :     $db_tie = tie %$db_lengths, 'DB_File', "$dbF.lengths", O_RDONLY, 0, $DB_BTREE;
100 :     $db_tie or warn "Cannot tie $dbF.lengths: $!";
101 :     }
102 :     if (!defined($db_tie))
103 :     {
104 :     if (open(my $dbfh, "<", $dbF))
105 :     {
106 :     while (my($id, $def, $seq) = gjoseqlib::read_next_fasta_seq($dbfh))
107 :     {
108 :     $db_lengths->{$id} = length($seq);
109 :     }
110 :     close($dbfh);
111 :     }
112 :     else
113 :     {
114 :     warn "Cannot open $dbF: $!";
115 :     }
116 :     }
117 : overbeek 1.1 }
118 :     else
119 :     {
120 :     die "$dbF is missing";
121 :     }
122 :     }
123 : gdpusch 1.21
124 :    
125 : olson 1.17 my $tmpF = "$tmp_dir/sim.out.$tmp_suffix";
126 : gdpusch 1.21
127 : olson 1.13 if ($use_blast)
128 :     {
129 :     open(my $fh, ">", $tmpF);
130 :     close($fh);
131 :     }
132 :     else
133 : olson 1.4 {
134 : olson 1.13 my @cmd = ($blat_cmd, $dbF, $qF, "-prot", "-out=blast8", $tmpF);
135 : gdpusch 1.21
136 : olson 1.16 #
137 :     # When running under FCGI, the system_with_redirect fails due to
138 :     # the FCGI library redefining open.
139 :     #
140 : gdpusch 1.21
141 : olson 1.16 my $rc;
142 :     if (defined($ENV{FCGI_ROLE}))
143 :     {
144 :     $rc = system(@cmd);
145 :     }
146 :     else
147 :     {
148 :     $rc = system_with_redirect(\@cmd, { stdout => '/dev/null' } );
149 :     }
150 : gdpusch 1.20
151 :     if ($ENV{DEBUG} || $ENV{VERBOSE} || $ENV{FIG_VERBOSE})
152 :     {
153 : gdpusch 1.22 print STDERR "Blat run succeeds: @cmd\n";
154 : gdpusch 1.20 }
155 :    
156 : olson 1.13 if ($rc != 0)
157 :     {
158 : gdpusch 1.22 warn "Blat run failed with rc=$rc: @cmd\n";
159 : olson 1.13 open(my $fh, ">", $tmpF);
160 :     close($fh);
161 :     }
162 : olson 1.4 }
163 : gdpusch 1.21
164 : olson 1.13 # my $cmd = $use_blast ? "touch $tmpF" : "$blat_cmd $dbF $qF -prot -out=blast8 $tmpF > /dev/null";
165 :     # #print STDERR "$cmd\n";
166 :     # my $rc = system $cmd;
167 :     # if ($rc != 0)
168 :     # {
169 :     # die "ProtSims::blastP: blat run failed with rc=$rc: $cmd\n";
170 :     # }
171 : gdpusch 1.21
172 : overbeek 1.5 my @sims1 = ();
173 :     open(BLAT,"<$tmpF") || die "could not open $tmpF";
174 :     my $sim = <BLAT>;
175 :     while ($sim && ($sim=~ /^(\S+\t\S+)/))
176 :     {
177 :     my $pegs = $1;
178 :     my @pieces = ();
179 :     while ($sim && ($sim=~ /^((\S+\t\S+).*\S)/) && ($2 eq $pegs))
180 :     {
181 :     push(@pieces,[split(/\t/,$1)]);
182 :     $sim = <BLAT>;
183 :     }
184 :     push(@sims1,&condense(\@pieces));
185 :     }
186 :     close(BLAT);
187 : overbeek 1.1 unlink $tmpF;
188 : olson 1.4 #print STDERR &Dumper(sims1 => \@sims1);
189 : gdpusch 1.21
190 : overbeek 1.1 my @rerun = ();
191 :     my @sims = ();
192 :     my $qI = 0;
193 :     my $simI = 0;
194 :     while ($qI < @$q)
195 :     {
196 :     my $peg = $q->[$qI]->[0];
197 : olson 1.4 #print STDERR "processing $peg\n";
198 : overbeek 1.1 my $hits = 0;
199 :     my $simI1 = $simI;
200 :     while (($simI1 < @sims1) && ($sims1[$simI1]->[0] eq $peg))
201 :     {
202 :     if (($simI1 == $#sims1) || ($sims1[$simI1]->[1] ne $sims1[$simI1+1]->[1]))
203 :     {
204 :     $hits++;
205 :     }
206 :     $simI1++;
207 :     }
208 : olson 1.4 #print STDERR "hits=$hits min_hits=$min_hits\n";
209 : overbeek 1.1 if ($hits >= $min_hits)
210 :     {
211 : olson 1.4 push(@sims,@sims1[$simI..$simI1-1]);
212 : overbeek 1.1 }
213 :     else
214 :     {
215 :     push(@rerun,$q->[$qI]);
216 :     }
217 :     $simI = $simI1;
218 :     $qI++;
219 :     }
220 : olson 1.4 #print STDERR &Dumper(rerun => \@rerun);
221 : gdpusch 1.21
222 : overbeek 1.1 if (@rerun > 0)
223 :     {
224 : olson 1.17 my $tmpQ = "$tmp_dir/tmpQ.$tmp_suffix";
225 : overbeek 1.1 &gjoseqlib::print_alignment_as_fasta($tmpQ,\@rerun);
226 : gdpusch 1.21
227 : olson 1.16 #
228 :     # If we're under FCGI (and thus under the servers), and the loadavg is low, do a small parallel run.
229 :     #
230 : gdpusch 1.21
231 : olson 1.16 my @par = ();
232 :     if (defined($ENV{FCGI_ROLE}) && open(my $la, "/proc/loadavg"))
233 :     {
234 :     my $l = <$la>;
235 :     chomp $l;
236 :     my @vals = split(/\s+/, $l);
237 :     my @procs = split(/\//, $vals[3]);
238 :     if ($vals[0] < 4 && $procs[0] < 8)
239 :     {
240 :     @par = ("-a", 4);
241 :     }
242 :     }
243 : gdpusch 1.21
244 : olson 1.16
245 : olson 1.13 # my $cmd = "$blastall_cmd -m 8 -i $tmpQ -d $dbF -FF -p blastp -e 1e-5";
246 : olson 1.16 my @cmd = ($blastall_cmd, '-m', 8, '-i', $tmpQ, '-d', $dbF, '-FF', '-p', 'blastp', '-e', '1e-5', @par);
247 : olson 1.4 #print STDERR "$cmd\n";
248 : olson 1.13 #open(BL, "$cmd|") or die "ProtSims::blastP: pipe to blast failed with $!: $cmd\n";
249 : gdpusch 1.21
250 : olson 1.14 #
251 :     # It'd be nice to do this but windows doesn't support it.
252 :     #
253 :     #open(BL, "-|", @cmd) or die "ProtSims::blastP: pipe to blast failed with $!: @cmd\n";
254 : olson 1.17 my $out_tmp = "$tmp_dir/blast_out.$tmp_suffix";
255 : olson 1.14 push(@cmd, "-o", $out_tmp);
256 : overbeek 1.19 #print STDERR "@cmd\n";
257 : olson 1.14 my $rc = system(@cmd);
258 :     if ($rc != 0)
259 :     {
260 :     warn "Error rc=$rc running blast: @cmd\n";
261 :     }
262 :     else
263 :     {
264 :     if (open(BL, "<", $out_tmp))
265 :     {
266 :     while (<BL>)
267 :     {
268 :     chomp;
269 :     push @sims, [ split(/\s+/, $_) ];
270 :     }
271 :     #my @blastout = map { chomp; [split(/\s+/,$_)] } <BL>;
272 :     # push(@sims,@blastout);
273 :     close(BL);
274 :     }
275 :     }
276 :     unlink $out_tmp;
277 : overbeek 1.1 unlink $tmpQ;
278 :     }
279 : gdpusch 1.21
280 : overbeek 1.3 my %lnQ = map { $_->[0] => length($_->[2]) } @$q;
281 : gdpusch 1.21
282 : olson 1.15 @sims = map { push(@$_, $lnQ{$_->[0]}, $db_lengths->{$_->[1]}); bless($_,'Sim') } @sims;
283 : gdpusch 1.21
284 : gdpusch 1.22 if ($qF eq "$tmp_dir/tmp.query.$tmp_suffix") { unlink $qF }
285 :     if ($dbF eq "$tmp_dir/tmp.db.$tmp_suffix") { unlink($dbF,"$dbF.psq","$dbF.pin","$dbF.phr") }
286 : overbeek 1.9 return sort { ($a->id1 cmp $b->id1) or ($a->psc <=> $b->psc) or ($a->id2 cmp $b->id2) } @sims;
287 : overbeek 1.5 }
288 :    
289 :     sub condense {
290 :     my($pieces) = @_;
291 :    
292 : overbeek 1.6 my $best_sc = $pieces->[0]->[10];
293 : overbeek 1.5 my @sims = sort { ($a->[6] <=> $b->[6]) } @$pieces;
294 : overbeek 1.7 while (@sims > 1)
295 : overbeek 1.5 {
296 : overbeek 1.7 my $gap1 = $sims[1]->[6] - $sims[0]->[7];
297 :     my $gap2 = $sims[1]->[8] - $sims[0]->[9];
298 :     my $diff = abs($gap1 - $gap2);
299 :     if (($gap1 <= 300) && ($gap2 <= 300) && ($diff <= (0.1 * &min($gap1,$gap2))))
300 :     {
301 :     $sims[0]->[7] = $sims[1]->[7];
302 :     $sims[0]->[9] = $sims[1]->[9];
303 :     }
304 : overbeek 1.5 splice(@sims,1,1);
305 :     }
306 : overbeek 1.6 $sims[0]->[10] = $best_sc;
307 : overbeek 1.5 return $sims[0];
308 : overbeek 1.1 }
309 :    
310 : overbeek 1.7 sub min {
311 :     my($x,$y) = @_;
312 :     return ($x <= $y) ? $x : $y;
313 :     }
314 :    
315 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3