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

Annotation of /FigKernelPackages/Kmers.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 # -*- perl -*-
2 :     ########################################################################
3 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     ########################################################################
18 :    
19 :     package Kmers;
20 :     no warnings 'redefine';
21 :    
22 :     use strict;
23 :     use DB_File;
24 : olson 1.20 use File::Basename;
25 : overbeek 1.1 use FIG;
26 :    
27 :     use Tracer;
28 :    
29 : olson 1.15 use ProtSims;
30 :    
31 : overbeek 1.1 use Data::Dumper;
32 :     use Carp;
33 : overbeek 1.8 use FFs;
34 : overbeek 1.1
35 : olson 1.3 our $KmersC_available;
36 :     eval {
37 :     require KmersC;
38 :     $KmersC_available++;
39 :     };
40 :    
41 :    
42 : overbeek 1.1 sub new {
43 : olson 1.19 my $class = shift;
44 : overbeek 1.1
45 : olson 1.19 my $figfams = {};
46 :    
47 :     $KmersC_available or die "KmersC module not available in this perl build";
48 :    
49 :     #
50 :     # Support experiments with kmer builds by allowing
51 :     # arguments to specify fri.db and setI.db as well as
52 :     # a single kmer data file.
53 :     #
54 :    
55 :     my $dir;
56 :     my $FRIDB;
57 :     my $setIDB;
58 :    
59 :     if ($_[0] =~ /^-/)
60 :     {
61 :     my %args = @_;
62 : olson 1.3
63 : olson 1.19 $FRIDB = $args{-frIdb};
64 :     $setIDB = $args{-setIdb};
65 :     my $table = $args{-table};
66 : overbeek 1.11
67 : olson 1.19 if (!defined($FRIDB) || !defined($setIDB) || !defined($table))
68 :     {
69 :     warn "Kmer experimental interface must define -frIDb, -setIdb, and -table\n";
70 :     return undef;
71 :     }
72 : olson 1.20
73 :     #
74 :     # Using this interface, we expect that the
75 :     # actual kmer sdirectory is the one containing FRIDB.
76 :     #
77 :    
78 :     my $ffdir = dirname($FRIDB);
79 :     $figfams->{ffs} = FFs->new($ffdir);
80 :    
81 : olson 1.19 my $kmerc = new KmersC;
82 :     $kmerc->open_data($table);
83 :     my $sz = $kmerc->get_motif_len();
84 :     warn "Experimental Kmers $table of size $sz\n";
85 :     $figfams->{KmerC}->{$sz} = $kmerc;
86 :     $figfams->{default_kmer_size} = $sz;
87 :     }
88 :     else
89 :     {
90 :     $dir = shift;
91 :     -d $dir or return undef;
92 : olson 1.3
93 : olson 1.19 $FRIDB = "$dir/FRI.db";
94 :     $setIDB = "$dir/setI.db";
95 :    
96 :     #
97 :     # Open a KmersC for each Kmer size.
98 :     #
99 :    
100 :     for my $binary (<$dir/binary/table.binary.*>)
101 :     {
102 :     if ($binary =~ /(\d+)$/)
103 :     {
104 :     my $k = $1;
105 :     my $kc = new KmersC;
106 :     $kc->open_data($binary) or die "cannot load Kmer binary database $binary";
107 :    
108 :     $figfams->{KmerC}->{$k} = $kc;
109 :     }
110 :     }
111 :     $figfams->{default_kmer_size} = 8;
112 : olson 1.20
113 :     $figfams->{ffs} = FFs->new($dir);
114 : olson 1.19 }
115 : olson 1.3
116 :     my %fr_hash;
117 : olson 1.4 my $fr_hash_tie = tie %fr_hash, 'DB_File', $FRIDB, O_RDONLY, 0666, $DB_HASH;
118 : olson 1.12 $fr_hash_tie || die "tie failed for function index $FRIDB";
119 : olson 1.3
120 : olson 1.12 my %set_hash;
121 :     my $set_hash_tie = tie %set_hash, 'DB_File', $setIDB, O_RDONLY, 0666, $DB_HASH;
122 :     $set_hash_tie || warn "tie failed for function index $setIDB";
123 : olson 1.3
124 :     $figfams->{friH} = \%fr_hash;
125 : olson 1.12 $figfams->{setiH} = \%set_hash;
126 : overbeek 1.8 $figfams->{fig} = new FIG;
127 : olson 1.3
128 : olson 1.17
129 : olson 1.3 bless $figfams,$class;
130 :     return $figfams;
131 :     }
132 :    
133 : overbeek 1.1 sub DESTROY {
134 :     my ($self) = @_;
135 :     delete $self->{fig};
136 :     }
137 :    
138 : olson 1.12 sub match_seq {
139 : olson 1.17 my($self, $motif_sz, $seq) = @_;
140 : overbeek 1.6
141 : olson 1.17 if ($self->{KmerC}->{$motif_sz})
142 : olson 1.3 {
143 : olson 1.12 my $matches = [];
144 : parrello 1.16 Confess("No sequence specified.") if ! $seq;
145 : olson 1.17 $self->{KmerC}->{$motif_sz}->find_all_hits(uc $seq, $matches);
146 : olson 1.12 return $matches;
147 : olson 1.3 }
148 : olson 1.17 else
149 :     {
150 :     die "No KmerC found for size $motif_sz";
151 :     }
152 :     }
153 :    
154 :     =head3 assign_functions_to_prot_set
155 :    
156 :     my $result = $kmers->assign_functions_to_prot_set($args)
157 :    
158 :     =over 4
159 :    
160 :     =item args
161 :    
162 :     Reference to a hash of parameters.
163 :    
164 :     =item RETURN
165 :    
166 :     Returns a list of tuples. If -all => 1 was specified, there will be a tuple for
167 :     each input sequence, otherwise there will be a tuple only for each input sequence
168 :     that had a successful match.
169 :    
170 :     Each tuple is of the form
171 :    
172 :     [ sequence-id, assigned-function, OTU, score, non-overlapping hit-count, overlapping hit-count, detailed-hits]
173 :    
174 :     where detailed-hits is optional, depending on the value of the -detailed parameter.
175 :    
176 :     =back
177 :    
178 :     =head4 Parameter Hash Fields
179 :    
180 :     =over 4
181 : olson 1.3
182 : olson 1.17 =item -seqs list
183 :    
184 :     Reference to a list of triples [sequence-id, comment, sequence-data].
185 :    
186 :     =item -kmer N
187 :    
188 :     Specify the kmer size to use for analysis.
189 :    
190 :     =item -all 0|1
191 :    
192 :     If set, return an entry for each input sequence. If the scoring threshold was not
193 :     met, the assigned function will be undefined, but the scores and hit counts
194 :     will still be returned..
195 :    
196 :     =item -scoreThreshold N
197 :    
198 :     Require a Kmer score of at least N for a Kmer match to succeed.
199 :    
200 :     =item -hitThreshold N
201 :    
202 :     Require at least N (possibly overlapping) Kmer hits for a Kmer match to succeed.
203 :    
204 :     =item -seqHitThreshold N
205 :    
206 :     Require at least N sequential (non-overlapping) Kmer hits for a Kmer match to succeed.
207 :    
208 :     =item -normalizeScores 0|1
209 :    
210 :     Normalize the scores to the size of the protein.
211 :    
212 :     =item -detailed 0|1
213 :    
214 :     If true, return a detailed accounting of the kmers hit for each protein.
215 :    
216 :     =back
217 :    
218 :    
219 :     =cut
220 :    
221 :     sub assign_functions_to_prot_set {
222 :    
223 :     #
224 :     # This is ugly but needed for the transition to the new form of the code.
225 :     # If we were not invoked with the new args hash, reinvoke the previous
226 :     # version of the code which is named assign_functions_to_prot_set_compat.
227 :     # Otherwise fall through and use the new form.
228 :     #
229 :     my($self, @args) = @_;
230 :    
231 :     my $args;
232 :     if (@args == 1 && ref($args[0]) eq 'HASH')
233 :     {
234 :     $args = $args[0];
235 :     }
236 :     elsif (@args > 0 && $args[0] !~ /^-/)
237 :     {
238 :     return &assign_functions_to_prot_set_compat;
239 :     }
240 :     else
241 :     {
242 :     $args = { @args };
243 :     }
244 :    
245 :     my $seq_set = $args->{-seqs};
246 :     if (!defined($seq_set))
247 :     {
248 :     die "assign_functions_to_prot_set: No sequences provided via the -seqs argument";
249 :     }
250 :    
251 :     my $motif_sz = $args->{-kmer};
252 :     if (!ref($self->{KmerC}->{$motif_sz}))
253 : overbeek 1.1 {
254 : olson 1.17 die "No KmerC defined for kmer size $motif_sz";
255 : overbeek 1.1 }
256 :    
257 : olson 1.17 #
258 :     # Define the good-hit function.
259 :     #
260 :    
261 :     my $min_hits = $args->{-hitThreshold};
262 :     my $min_seqhits = $args->{-seqHitThreshold};
263 :     my $min_score = $args->{-scoreThreshold};
264 :    
265 :     my $good_hit = sub {
266 :     my($score, $hits, $seqhits) = @_;
267 : olson 1.3
268 : olson 1.17 return 0 if defined($min_hits) && $hits < $min_hits;
269 :     return 0 if defined($min_seqhits) && $seqhits < $min_seqhits;
270 :     return 0 if defined($min_score) && $score < $min_score;
271 :     return 1;
272 :     };
273 : olson 1.3
274 : olson 1.12 my $fr_hash = $self->{friH};
275 :     my $set_hash = $self->{setiH};
276 : olson 1.17 my $fig = $self->{fig};
277 :    
278 :     my @out;
279 :    
280 :     #
281 :     # Run the query.
282 :     #
283 :    
284 :     for my $seqent (@$seq_set)
285 : olson 1.12 {
286 : olson 1.17 my($id, $com, $seq) = @$seqent;
287 :     my $seq_len = length($seq);
288 :     my $matches = $self->match_seq($motif_sz, $seq);
289 :    
290 :     my(%hitsF,%hitsS);
291 :    
292 :     #
293 :     # Compute the non-overlapping hits.
294 :     # Do this by walking the list of hits in order of offset, and for each,
295 :     # compute the distance to the start of the last hit. If that distance
296 :     # is K or greater, increment the non-overlapping hit count.
297 :     #
298 :     my $last_hit_start;
299 :     my $last_frI;
300 :     my %non_overlapping_F;
301 :    
302 :     my @sorted_matches = sort { $a->[2] <=> $b->[2] or $a->[0] <=> $b->[0] } @$matches;
303 :     my @details;
304 :    
305 :     foreach my $match (@sorted_matches)
306 :     {
307 :     my($offset, $oligo, $frI, $setI) = @$match;
308 :    
309 :     if ($args->{-detailed})
310 :     {
311 :     push(@details, [$offset, $oligo, $fr_hash->{$frI}, $set_hash->{$setI}]);
312 :     }
313 :    
314 :     #
315 :     # $offset - offset in target of this hit
316 :     # $oligo - actual oligo that hit
317 :     # $frI - functional role index for this kmer hit
318 :     # $setI - OTU set index for this kmer hit
319 :     #
320 :    
321 :     #
322 :     # The sort above groups the hits by functional role and orders by
323 :     # offset within the target.
324 :     #
325 :    
326 :     if ($frI != $last_frI)
327 :     {
328 :     #
329 :     # We're searching in a new role now.
330 :     #
331 :    
332 :     undef $last_hit_start;
333 :     $last_frI = $frI;
334 :     }
335 :     # print "@$match lhs=$last_hit_start\n";
336 :    
337 :     #
338 :     # If this is the first hit, or if we are a kmer width away at least from the
339 :     # last non-overlapping hit, count another non-overlapping hit.
340 :     #
341 :    
342 :     if (!defined($last_hit_start) || ($offset - $last_hit_start) >= $motif_sz)
343 :     {
344 :     $non_overlapping_F{$frI}++;
345 :     $last_hit_start = $offset;
346 :     }
347 :    
348 :     #
349 :     # We count all functional role and OTU hits.
350 :     #
351 :     $hitsF{$frI}++;
352 :     if ($setI)
353 :     {
354 :     $hitsS{$setI}++ ;
355 :     }
356 :     }
357 :    
358 :     #
359 :     # Find the functional roles that had the best overlapping and non-overlapping hit counts.
360 :     #
361 :     my $nonoverlapFRI = &best_hit(\%non_overlapping_F, $min_hits);
362 :     my ($FRI, $bh2) = &best_2hits(\%hitsF, $min_seqhits);
363 :    
364 :     #
365 :     # Compute score and normalize if required. Score is based on the number of overlapping hits.
366 :     #
367 :    
368 :     my $nonoverlap_hits = $non_overlapping_F{$nonoverlapFRI};
369 :     my $overlap_hits = $hitsF{$FRI};
370 :     my $overlap_hits2 = $hitsF{$bh2};
371 :    
372 :     my $score = 0 + ($overlap_hits - $overlap_hits2);
373 :     if ($args->{-normalizeScores})
374 :     {
375 :     $score = $score / ($seq_len - $motif_sz + 1);
376 :     }
377 :    
378 :     my $setI = &best_hit(\%hitsS,$min_hits);
379 :    
380 :     my $fun;
381 :     my $set;
382 : olson 1.20 my $family;
383 :     my $family_sims;
384 : olson 1.17 if (($nonoverlapFRI == $FRI) && &$good_hit($score, $nonoverlap_hits, $overlap_hits))
385 :     {
386 :     $fun = $fr_hash->{$nonoverlapFRI};
387 :     $set = $set_hash->{$setI};
388 : olson 1.20
389 :     if (defined($args->{-determineFamily}) && $args->{-determineFamily} && ref($self->{ffs}))
390 :     {
391 :     my($fam, $sims) = $self->{ffs}->place_seq_and_function_in_family($seq, $fun);
392 :     $family = $fam;
393 :     $family_sims = $sims if $args->{-returnFamilySims};
394 :     }
395 : olson 1.12 }
396 : overbeek 1.14
397 : olson 1.17 if ($fun or $args->{-all})
398 : overbeek 1.14 {
399 : olson 1.17 push(@out, [$id, $fun, $set, $score, $nonoverlap_hits, $overlap_hits,
400 : olson 1.20 ($args->{-detailed} ? \@details : undef),
401 :     $family, $family_sims]);
402 : overbeek 1.14 }
403 :     }
404 : olson 1.17 return @out;
405 :     }
406 :    
407 :     sub assign_functions_using_similarity
408 :     {
409 :     my($self, $args) = _handle_args(@_);
410 :    
411 :     my $seqs = $args->{-seqs};
412 :     if (!defined($seqs))
413 : overbeek 1.14 {
414 : olson 1.17 die "assign_functions_using_similarity: No sequences provided via the -seqs argument";
415 :     }
416 : overbeek 1.14
417 : olson 1.17 my $db = $args->{-simDb};
418 :     if (!defined($db))
419 :     {
420 :     die "assign_functions_using_similarity: No similarity database provided by the -simDb argument";
421 :     }
422 : olson 1.15
423 : olson 1.18 if (@$seqs == 0)
424 :     {
425 :     return ();
426 :     }
427 :    
428 : olson 1.17 my @blastout = ProtSims::blastP($seqs, $db, 5);
429 : olson 1.15
430 : olson 1.17 my $cur;
431 :     my @set;
432 :     my @out;
433 :     for my $ent (@blastout)
434 :     {
435 :     if ($cur && $ent->id1 ne $cur)
436 : overbeek 1.14 {
437 : olson 1.17 my $val = $self->process_blast_set(@set);
438 :     push(@out, $val);
439 :     @set = ();
440 : overbeek 1.14 }
441 : olson 1.17 $cur = $ent->id1;
442 :     push(@set, $ent);
443 :     }
444 :     my $val = $self->process_blast_set(@set);
445 :     push(@out, $val);
446 :     return @out;
447 :     }
448 :    
449 :     sub process_blast_set
450 :     {
451 :     my($self, @blastout) = @_;
452 :    
453 :     my $id = $blastout[0]->id1;
454 :     my $fig = $self->{fig};
455 :    
456 :     if (@blastout > 5) { $#blastout = 4 }
457 :    
458 :     my %hit_pegs = map { $_->id2 => 1 } @blastout;
459 :     my @pegs = keys(%hit_pegs);
460 :     if (@pegs == 0)
461 :     {
462 :     return [$id,'hypothetical protein', undef, 0, 0, 0, undef];
463 :     }
464 :     else
465 :     {
466 :     my %funcs;
467 :     foreach my $peg (@pegs)
468 : overbeek 1.14 {
469 : olson 1.17 my $func = $fig->function_of($peg,1);
470 :     if (! &FIG::hypo($func))
471 : overbeek 1.14 {
472 : olson 1.17 $funcs{$func}++;
473 : overbeek 1.14 }
474 :     }
475 : olson 1.17 my @pos = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
476 :     my $proposed = (@pos > 0) ? $pos[0] : "hypothetical protein";
477 :     return [$id, $proposed, undef, 0, 0, 0, undef];
478 : olson 1.3 }
479 :     }
480 : overbeek 1.11
481 : olson 1.17 sub assign_functions_to_prot_set_compat {
482 : olson 1.15 my($self,$seq_set,$blast,$min_hits,$extra_blastdb) = @_;
483 :     $min_hits = 3 unless defined($min_hits);
484 : olson 1.17
485 :     my %match_set = map { my($id, $com, $seq) = @$_; $id => [$self->match_seq($self->{default_kmer_size},$seq), $seq] } @$seq_set;
486 :    
487 : olson 1.15 my $fr_hash = $self->{friH};
488 :     my $set_hash = $self->{setiH};
489 : olson 1.17
490 : olson 1.15 my $fig = $self->{fig};
491 : olson 1.17
492 : olson 1.15 my @missing;
493 :     while (my($id, $ent) = each %match_set)
494 :     {
495 :     my($matches, $seq) = @$ent;
496 :    
497 :     my(%hitsF,%hitsS);
498 : olson 1.17 foreach my $match (@$matches)
499 : olson 1.15 {
500 :     my($offset, $oligo, $frI, $setI) = @$match;
501 : olson 1.17 $hitsF{$frI}++;
502 : olson 1.15 if ($setI)
503 : olson 1.17 {
504 : olson 1.15 $hitsS{$setI}++ ;
505 :     }
506 :     }
507 :    
508 :     my $FRI = &best_hit(\%hitsF,$min_hits);
509 :     my $setI = &best_hit(\%hitsS,$min_hits);
510 :     push(@$ent, $FRI, $setI, \%hitsF);
511 : olson 1.17
512 : olson 1.15 if (!$fr_hash->{$FRI})
513 :     {
514 :     push(@missing, [$id, undef, $seq]);
515 :     }
516 :     }
517 : olson 1.17
518 : olson 1.15 #
519 :     # @missing now has the list of sequences that had no Kmer hits. If we have a
520 :     # blast db, blast 'em.
521 : olson 1.17
522 : olson 1.15 my @all_blastout;
523 :     if (@missing && -s $extra_blastdb)
524 :     {
525 :     #print Dumper(\@missing);
526 :     @all_blastout = ProtSims::blastP(\@missing, $extra_blastdb, 5);
527 :     #print Dumper(\@all_blastout);
528 :     }
529 : olson 1.17
530 : olson 1.15 #
531 :     # We now have Kmers output and blast output. Go through the original data and
532 :     # create the output.
533 :     #
534 : olson 1.17
535 : olson 1.15 my @out;
536 :    
537 :     for my $ent (@$seq_set)
538 :     {
539 :     my $id = $ent->[0];
540 :     my ($matches, $seq, $FRI, $setI, $hitsF) = @{$match_set{$id}};
541 : olson 1.17
542 : olson 1.15 my $blast_results = [];
543 :     if ($fr_hash->{$FRI})
544 :     {
545 :     if ($blast && ($fr_hash->{$FRI} || $set_hash->{$setI}))
546 :     {
547 :     $blast_results = &blast_data($self,$id,$seq,$fr_hash->{$FRI},$blast,'blastp');
548 :     }
549 : olson 1.17
550 : olson 1.15 push(@out, [$id, $fr_hash->{$FRI},$set_hash->{$setI}, $blast_results,$hitsF->{$FRI}]);
551 :     }
552 :     else
553 :     {
554 :     my @blastout = grep { $_->id1 eq $id } @all_blastout;
555 : olson 1.17
556 : olson 1.15 if (@blastout > 5) { $#blastout = 4 }
557 :    
558 :     my %hit_pegs = map { $_->id2 => 1 } @blastout;
559 :     my @pegs = keys(%hit_pegs);
560 :     if (@pegs == 0)
561 :     {
562 :     push(@out, [$id,'hypothetical protein','',[],0]);
563 :     }
564 :     else
565 :     {
566 :     my %funcs;
567 :     foreach my $peg (@pegs)
568 :     {
569 :     my $func = $fig->function_of($peg,1);
570 :     if (! &FIG::hypo($func))
571 :     {
572 :     $funcs{$func}++;
573 :     }
574 :     }
575 :     my @pos = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
576 :     my $proposed = (@pos > 0) ? $pos[0] : "hypothetical protein";
577 :     push(@out, [$id, $proposed,'',[],0]);
578 :     }
579 :     }
580 :     }
581 :     return @out;
582 :     }
583 :    
584 : olson 1.17
585 :    
586 :    
587 :    
588 : olson 1.12 sub best_hit {
589 : overbeek 1.14 my($hits,$min_hits) = @_;
590 : olson 1.12 my @poss = sort { $hits->{$b} <=> $hits->{$a} } keys(%$hits);
591 : overbeek 1.11
592 : olson 1.12 my $val;
593 :     if ((@poss > 0) && ($hits->{$poss[0]} >= $min_hits))
594 :     {
595 : overbeek 1.14 $val = $poss[0];
596 : olson 1.12 }
597 :     return $val;
598 :     }
599 : overbeek 1.11
600 : olson 1.17 sub best_2hits {
601 :     my($hits,$min_hits) = @_;
602 :     my @poss = sort { $hits->{$b} <=> $hits->{$a} } keys(%$hits);
603 :    
604 :     my $val;
605 :     if ((@poss > 0) && ($hits->{$poss[0]} >= $min_hits))
606 :     {
607 :     return @poss[0,1];
608 :     }
609 :     else
610 :     {
611 :     return undef;
612 :     }
613 :     }
614 :    
615 : olson 1.12 sub best_hit_in_group
616 :     {
617 :     my($group) = @_;
618 : overbeek 1.11
619 : olson 1.12 my %hash;
620 :     for my $tuple (@$group)
621 :     {
622 :     my($off,$oligo,$frI,$setI) = @$tuple;
623 :     if ($setI > 0)
624 :     {
625 :     $hash{$setI}++;
626 :     }
627 :     }
628 :     my @sorted = sort { $hash{$b} <=> $hash{$a} } keys %hash;
629 :     my $max = $sorted[0];
630 :     return $max;
631 : overbeek 1.11 }
632 :    
633 : olson 1.12 sub assign_functions_to_DNA_features {
634 : olson 1.17 my($self,$motif_sz,$seq,$min_hits,$max_gap,$blast) = @_;
635 : olson 1.7
636 :     $min_hits = 3 unless defined($min_hits);
637 : olson 1.12 $max_gap = 200 unless defined($max_gap);
638 : overbeek 1.1
639 : olson 1.12 my $fr_hash = $self->{friH};
640 :     my $set_hash = $self->{setiH};
641 : olson 1.3
642 : overbeek 1.1 my %hits;
643 :     my @ans;
644 : olson 1.12 my $matches = $self->process_dna_seq($seq);
645 : olson 1.3
646 : olson 1.17 push(@ans,&process_hits($self,$motif_sz, $matches,1,length($seq),$motif_sz, $min_hits, $max_gap,$blast,$seq));
647 : overbeek 1.1 undef %hits;
648 : overbeek 1.8
649 : olson 1.12 $matches = $self->process_dna_seq(&FIG::reverse_comp($seq));
650 : olson 1.17 push(@ans,&process_hits($self,$motif_sz, $matches,length($seq),1,$motif_sz, $min_hits, $max_gap,$blast,$seq));
651 : overbeek 1.8 return \@ans;
652 : olson 1.12 }
653 : overbeek 1.1
654 : olson 1.12 sub process_dna_seq {
655 : olson 1.3 my($self, $seq,$hits) = @_;
656 :    
657 : olson 1.12 my $matches = $self->match_seq($seq);
658 :     return $matches;
659 :     }
660 : overbeek 1.1
661 :    
662 : olson 1.12 sub process_hits {
663 : olson 1.17 my($self,$motif_sz, $matches,$beg,$end,$sz_of_match, $min_hits, $max_gap,$blast,$seq) = @_;
664 : olson 1.3
665 : olson 1.12 my $fr_hash = $self->{friH};
666 :     my $set_hash = $self->{setiH};
667 : overbeek 1.11
668 : olson 1.12 my $hits;
669 :     my %sets;
670 :     foreach my $tuple (@$matches)
671 : overbeek 1.11 {
672 : olson 1.12 my($off,$oligo,$frI,$setI) = @$tuple;
673 :     push(@{$hits->{$frI}},$tuple);
674 : overbeek 1.11 }
675 : overbeek 1.1
676 :     my @got = ();
677 : olson 1.12 my @poss = sort { (@{$hits->{$b}} <=> @{$hits->{$a}}) } keys(%$hits);
678 : overbeek 1.1 if (@poss != 0)
679 :     {
680 :     foreach my $frI (@poss)
681 :     {
682 : olson 1.12 my $hit_list = $hits->{$frI};
683 :     my @grouped = &group_hits($hit_list, $max_gap);
684 :     foreach my $group_ent (@grouped)
685 : overbeek 1.1 {
686 : olson 1.12 my($group, $group_hits) = @$group_ent;
687 :     my $N = @$group;
688 :     if ($N >= $min_hits) # consider only runs containing 3 or more hits
689 : overbeek 1.1 {
690 : overbeek 1.6 my $b1 = $group->[0];
691 : olson 1.12 my $e1 = $group->[-1] + ($sz_of_match-1);
692 : overbeek 1.6
693 :     my $loc;
694 :     if ($beg < $end)
695 :     {
696 : olson 1.12 $loc = [$beg+$b1,$beg+$e1];
697 : overbeek 1.6 }
698 :     else
699 :     {
700 : olson 1.12 $loc = [$beg-$b1,$beg-$e1];
701 : overbeek 1.6 }
702 :     my $func = $fr_hash->{$frI};
703 : olson 1.12
704 :     my $set = &best_hit_in_group($group_hits);
705 :     $set = $set_hash->{$set};
706 :    
707 :     my $blast_output = [];
708 :     if ($blast)
709 :     {
710 : overbeek 1.14 $blast_output = &blast_data($self,join("_",@$loc),$seq,$func,$blast,
711 : olson 1.12 ($motif_sz == $sz_of_match) ? 'blastn' : 'blastx');
712 :     }
713 :    
714 :     my $tuple = [$N,@$loc,$func,$set,$blast_output];
715 :    
716 : overbeek 1.8 push(@got,$tuple);
717 : overbeek 1.1 }
718 :     }
719 :     }
720 :     }
721 :     return @got;
722 :     }
723 :    
724 : overbeek 1.6 sub group_hits {
725 : olson 1.12 my($hits, $max_gap) = @_;
726 : overbeek 1.6
727 : olson 1.12 my @sorted = sort { $a->[0] <=> $b->[0] } @$hits;
728 : overbeek 1.6 my @groups = ();
729 :     my $position;
730 : olson 1.12 while (defined(my $hit = shift @sorted))
731 : overbeek 1.6 {
732 : olson 1.12 my($position,$oligo,$frI,$setI) = @$hit;
733 :    
734 : overbeek 1.6 my $group = [$position];
735 : olson 1.12 my $ghits = [$hit];
736 :     while ((@sorted > 0) && (($sorted[0]->[0] - $position) < $max_gap))
737 : overbeek 1.6 {
738 : olson 1.12 $hit = shift @sorted;
739 :     ($position,$oligo,$frI,$setI) = @$hit;
740 :     push(@$group,$position+1);
741 :     push(@$ghits, $hit);
742 : overbeek 1.6 }
743 :    
744 : olson 1.12 push(@groups,[$group, $ghits]);
745 : overbeek 1.6 }
746 :     return @groups;
747 :     }
748 :    
749 : olson 1.12 sub assign_functions_to_PEGs_in_DNA {
750 : olson 1.17 my($self,$motif_sz,$seq,$min_hits,$max_gap,$blast) = @_;
751 : olson 1.12
752 :     $blast = 0 unless defined($blast);
753 :     $min_hits = 3 unless defined($min_hits);
754 :     $max_gap = 200 unless defined($max_gap);
755 :    
756 :     my $fr_hash = $self->{friH};
757 :     my $set_hash = $self->{setiH};
758 :    
759 :     my %hits;
760 :     my @ans;
761 : olson 1.17 my $matches = $self->process_prot_seq($motif_sz, $seq);
762 :     push(@ans,&process_hits($self,$motif_sz, $matches,1,length($seq),3 * $motif_sz, $min_hits, $max_gap,$blast,$seq));
763 : olson 1.12 undef %hits;
764 : olson 1.17 $matches = $self->process_prot_seq($motif_sz, &FIG::reverse_comp($seq));
765 :     push(@ans,&process_hits($self,$motif_sz, $matches,length($seq),1,3 * $motif_sz, $min_hits, $max_gap,$blast,$seq));
766 : olson 1.12 return \@ans;
767 :     }
768 :    
769 :     sub process_prot_seq {
770 : olson 1.17 my($self, $motif_sz, $seq) = @_;
771 : olson 1.12
772 :     my $ans = [];
773 :     my $ln = length($seq);
774 :     my($i,$off);
775 :     for ($off=0; ($off < 3); $off++)
776 :     {
777 :     my $ln_tran = int(($ln - $off)/3) * 3;
778 :     my $tran = uc &FIG::translate(substr($seq,$off,$ln_tran));
779 :    
780 :    
781 : olson 1.17 my $matches = $self->match_seq($motif_sz, $tran);
782 : olson 1.12
783 :     push(@$ans, map { $_->[0] = ((3 * $_->[0]) + $off); $_ } @$matches);
784 :     }
785 :     return $ans;
786 :     }
787 :    
788 : overbeek 1.8 use Sim;
789 :    
790 :     sub blast_data {
791 : overbeek 1.14 my($self,$id,$seq,$func,$blast,$tool) = @_;
792 : overbeek 1.8
793 : overbeek 1.13 if ($tool eq "blastp")
794 :     {
795 : overbeek 1.14 return &blast_data1($self,$id,$seq,$func,$blast,$tool);
796 : overbeek 1.13 }
797 :    
798 :     if ($id =~ /^(\d+)_(\d+)$/)
799 :     {
800 :     my($b,$e) = ($1 < $2) ? ($1,$2) : ($2,$1);
801 :     my $b_adj = (($b - 5000) > 0) ? $b-5000 : 1;
802 :     my $e_adj = (($b + 5000) <= length($seq)) ? $b+5000 : length($seq);
803 :     my $seq1 = substr($seq,$b_adj-1, ($e_adj - $b_adj)+1);
804 : overbeek 1.14 my $blast_out = &blast_data1($self,$id,$seq1,$func,$blast,$tool);
805 : overbeek 1.13 foreach $_ (@$blast_out)
806 :     {
807 :     $_->[2] += $b_adj - 1;
808 :     $_->[3] += $b_adj - 1;
809 :     $_->[8] = length($seq);
810 :     }
811 :     return $blast_out;
812 :     }
813 :     else
814 :     {
815 : overbeek 1.14 return &blast_data1($self,$id,$seq,$func,$blast,$tool);
816 : overbeek 1.13 }
817 :     }
818 :    
819 :     sub blast_data1 {
820 : overbeek 1.14 my($self,$id,$seq,$func,$blast,$tool) = @_;
821 : overbeek 1.13
822 : overbeek 1.11
823 : overbeek 1.8 if (! $tool) { $tool = 'blastx' }
824 :     my $fig = $self->{fig};
825 : overbeek 1.11
826 : overbeek 1.8 my @blastout = ();
827 : overbeek 1.11 if ($tool ne 'blastn')
828 :     {
829 : olson 1.20 my $ffs = $self->{ffs};
830 : overbeek 1.11 my @fams = $ffs->families_implementing_role($func);
831 :     foreach my $fam (@fams)
832 :     {
833 :     my $subD = substr($fam,-3);
834 :     my $pegs_in_fam = "$FIG_Config::FigfamsData/FIGFAMS/$subD/$fam/PEGs.fasta";
835 :     push(@blastout,map { [$_->id2,$_->iden,$_->b1,$_->e1,$_->b2,$_->e2,$_->psc,$_->bsc,$_->ln1,$_->ln2,$fam] }
836 :     $fig->blast($id,$seq,$pegs_in_fam,0.1,"-FF -p $tool -b $blast"));
837 :     }
838 :     }
839 :     else
840 : overbeek 1.8 {
841 : overbeek 1.11 push(@blastout,map { [$_->id2,$_->iden,$_->b1,$_->e1,$_->b2,$_->e2,$_->psc,$_->bsc,$_->ln1,$_->ln2,$self->{what}] }
842 :     $fig->blast($id,$seq,$self->{blastdb},0.1,"-FF -p $tool -b $blast"));
843 : overbeek 1.8 }
844 :     @blastout = sort { $b->[7] <=> $a->[7] } @blastout;
845 :     if (@blastout > $blast) { $#blastout = $blast-1 }
846 :     return \@blastout;
847 :     }
848 :    
849 : olson 1.17 #
850 :     # Turn an argument list into a $self ref and an argument hash.
851 :     # Code lifted from ClientThing.
852 :     #
853 :     sub _handle_args
854 :     {
855 :     my $self = shift;
856 :     my $args = $_[0];
857 :     if (defined $args)
858 :     {
859 :     if (scalar @_ gt 1)
860 :     {
861 :     # Here we have multiple arguments. We check the first one for a
862 :     # leading hyphen.
863 :     if ($args =~ /^-/) {
864 :     # This means we have hash-form parameters.
865 :     my %args = @_;
866 :     $args = \%args;
867 :     } else {
868 :     # This means we have list-form parameters.
869 :     my @args = @_;
870 :     $args = \@args;
871 :     }
872 :     } else {
873 :     # Here we have a single argument. If it's a scalar, we convert it
874 :     # to a singleton list.
875 :     if (! ref $args) {
876 :     $args = [$args];
877 :     }
878 :     }
879 :     }
880 :     return($self, $args);
881 :     }
882 :    
883 :    
884 :    
885 :    
886 :    
887 : overbeek 1.1 1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3