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

Annotation of /FigKernelPackages/KmersOld.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 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 :     use FIG;
25 :    
26 :     use Tracer;
27 :    
28 :     use ProtSims;
29 :    
30 :     use Data::Dumper;
31 :     use Carp;
32 :     use FFs;
33 :    
34 :     our $KmersC_available;
35 :     eval {
36 :     require KmersC;
37 :     $KmersC_available++;
38 :     };
39 :    
40 :    
41 :     # This is the constructor. Presumably, $class is 'Kmers'.
42 :     #
43 :    
44 :     sub new {
45 :     my($class,$KmerDB,$FRIDB,$setIDB) = @_;
46 :    
47 :     my $figfams = {};
48 :     $figfams->{what} = '';
49 :     $figfams->{blastdb} = "./blastdb";
50 :    
51 :     my $dir;
52 :     if (defined($KmerDB) && (! defined($FRIDB)))
53 :     {
54 :     $dir = $KmerDB;
55 :     ($KmerDB,$FRIDB,$setIDB) = ("$dir/kmer.db","$dir/FRI.db","$dir/setI.db");
56 :    
57 :     if (open(WHAT,"<$dir/what") && defined($_ = <WHAT>))
58 :     {
59 :     chomp;
60 :     $figfams->{what} = $_;
61 :     close(WHAT);
62 :     }
63 :     $figfams->{blastdb} = "$dir/blastdb";
64 :     }
65 :     else
66 :     {
67 :     $dir = "."; # look for 'blastdb' and 'what' in current directory
68 :     }
69 :     $figfams->{dir} = $dir;
70 :    
71 :     if ((! defined($KmerDB)) || (! defined($FRIDB)) || (! defined($setIDB))) { return undef }
72 :    
73 :     my %fr_hash;
74 :     my %set_hash;
75 :     my %kmer_hash;
76 :    
77 :     my $fr_hash_tie = tie %fr_hash, 'DB_File', $FRIDB, O_RDONLY, 0666, $DB_HASH;
78 :     my $set_hash_tie = tie %set_hash, 'DB_File', $setIDB, O_RDONLY, 0666, $DB_HASH;
79 :     my $kmer_hash_tie = tie %kmer_hash, 'DB_File', $KmerDB, O_RDONLY, 0666, $DB_HASH;
80 :    
81 :     $fr_hash_tie || die "tie failed for function index $FRIDB";
82 :     $set_hash_tie || die "tie failed for function index $FRIDB";
83 :     $kmer_hash_tie || die "tie failed for kmer hash $KmerDB";
84 :    
85 :     my($motif,undef) = each %kmer_hash;
86 :     $figfams->{size} = length($motif);
87 :    
88 :     $figfams->{KmerH} = \%kmer_hash;
89 :     $figfams->{friH} = \%fr_hash;
90 :     $figfams->{setiH} = \%set_hash;
91 :    
92 :     $figfams->{fig} = new FIG;
93 :     bless $figfams,$class;
94 :     return $figfams;
95 :     }
96 :    
97 :     sub new_using_C {
98 :     my($class,$KmerBinaryDB,$FRIDB,$setIDB) = @_;
99 :    
100 :     if (defined($KmerBinaryDB) && (! defined($FRIDB)))
101 :     {
102 :     my $dir = $KmerBinaryDB;
103 :     ($KmerBinaryDB,$FRIDB,$setIDB) = ("$dir/table.binary","$dir/FRI.db","$dir/setI.db");
104 :     }
105 :    
106 :     my $figfams = {};
107 :     if ((! defined($KmerBinaryDB)) || (! defined($FRIDB)) || (! defined("$setIDB"))) { return undef }
108 :    
109 :    
110 :     $KmersC_available or die "KmersC module not available in this perl build";
111 :    
112 :     my %fr_hash;
113 :     my $fr_hash_tie = tie %fr_hash, 'DB_File', $FRIDB, O_RDONLY, 0666, $DB_HASH;
114 :     $fr_hash_tie || die "tie failed for function index $FRIDB";
115 :    
116 :     my %set_hash;
117 :     my $set_hash_tie = tie %set_hash, 'DB_File', $setIDB, O_RDONLY, 0666, $DB_HASH;
118 :     $set_hash_tie || warn "tie failed for function index $setIDB";
119 :    
120 :     my $kc = new KmersC;
121 :     $kc->open_data($KmerBinaryDB) or die "cannot load Kmer binary database $KmerBinaryDB";
122 :    
123 :     $figfams->{size} = $kc->get_motif_len();
124 :     $figfams->{KmerC} = $kc;
125 :     $figfams->{friH} = \%fr_hash;
126 :     $figfams->{setiH} = \%set_hash;
127 :     $figfams->{fig} = new FIG;
128 :    
129 :     bless $figfams,$class;
130 :     return $figfams;
131 :     }
132 :    
133 :     sub DESTROY {
134 :     my ($self) = @_;
135 :     delete $self->{fig};
136 :     }
137 :    
138 :     sub match_seq {
139 :     my($self,$seq) = @_;
140 :    
141 :     if ($self->{KmerC})
142 :     {
143 :     my $matches = [];
144 :     Confess("No sequence specified.") if ! $seq;
145 :     $self->{KmerC}->find_all_hits(uc $seq, $matches);
146 :     return $matches;
147 :     }
148 :    
149 :     my $kmer_hash = $self->{KmerH};
150 :     my $motif_sz = $self->{size};
151 :     my $matches = [];
152 :     my $ln = length($seq);
153 :     my $i;
154 :     for ($i=0; ($i < ($ln - $motif_sz)); $i++)
155 :     {
156 :     my $oligo = uc substr($seq,$i,$motif_sz);
157 :     my $x = $kmer_hash->{$oligo};
158 :     if (defined($x))
159 :     {
160 :     push(@$matches,[$i,$oligo,split(/\t/,$x)]);
161 :     }
162 :     }
163 :     return $matches;
164 :     }
165 :    
166 :     sub assign_function_to_prot {
167 :     my($self,$seq,$blast,$min_hits,$extra_blastdb) = @_;
168 :     $min_hits = 3 unless defined($min_hits);
169 :    
170 :     my $matches = $self->match_seq($seq);
171 :    
172 :     my $fr_hash = $self->{friH};
173 :     my $set_hash = $self->{setiH};
174 :    
175 :     my(%hitsF,%hitsS);
176 :     foreach my $match (@$matches)
177 :     {
178 :     my($offset, $oligo, $frI, $setI) = @$match;
179 :     $hitsF{$frI}++;
180 :     if ($setI)
181 :     {
182 :     $hitsS{$setI}++ ;
183 :     }
184 :     }
185 :    
186 :     my $FRI = &best_hit(\%hitsF,$min_hits);
187 :     my $setI = &best_hit(\%hitsS,$min_hits);
188 :     my $blast_results = [];
189 :     if ($fr_hash->{$FRI})
190 :     {
191 :     if ($blast && ($fr_hash->{$FRI} || $set_hash->{$setI}))
192 :     {
193 :     $blast_results = &blast_data($self,'query',$seq,$fr_hash->{$FRI},$blast,'blastp');
194 :     }
195 :     return [$fr_hash->{$FRI},$set_hash->{$setI}, $blast_results,$hitsF{$FRI}];
196 :     }
197 :     elsif ((-s $extra_blastdb) &&
198 :     (-s "$extra_blastdb.psq") &&
199 :     (-M $extra_blastdb >= -M "$extra_blastdb.pdq"))
200 :     {
201 :     my $fig = $self->{fig};
202 :    
203 :     # my $tmpF = "$FIG_Config::temp/tmpseq.$$.fasta";
204 :     # open(TMP,">$tmpF") || die "could not open $tmpF";
205 :     # print TMP ">query\n$seq\n";
206 :     # close(TMP);
207 :    
208 :     my $seq_inp = [['query', '', $seq]];
209 :     my @blastout = ProtSims::blastP($seq_inp, $extra_blastdb, 5);
210 :    
211 :     #my @blastout = `$FIG_Config::ext_bin/blastall -p blastp -FF -m 8 -e 1.0-20 -d $extra_blastdb -i $tmpF`;
212 :     #unlink $tmpF;
213 :    
214 :     # if (@blastout > 5) { $#blastout = 4 }
215 :     # my %hit_pegs = map { $_ =~ /^\S+\t(\S+)/; $1 => 1 } @blastout;
216 :     my %hit_pegs = map { $_->id2 => 1 } @blastout;
217 :     my @pegs = keys(%hit_pegs);
218 :     if (@pegs == 0)
219 :     {
220 :     return ['hypothetical protein','',[],0];
221 :     }
222 :     else
223 :     {
224 :     my %funcs;
225 :     foreach my $peg (@pegs)
226 :     {
227 :     my $func = $fig->function_of($peg,1);
228 :     if (! &FIG::hypo($func))
229 :     {
230 :     $funcs{$func}++;
231 :     }
232 :     }
233 :     my @pos = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
234 :     my $proposed = (@pos > 0) ? $pos[0] : "hypothetical protein";
235 :     return [$proposed,'',[],0];
236 :     }
237 :     }
238 :     else
239 :     {
240 :     return ['','',[],0];
241 :     }
242 :     }
243 :    
244 :     sub assign_functions_to_prot_set {
245 :     my($self,$seq_set,$blast,$min_hits,$extra_blastdb) = @_;
246 :     $min_hits = 3 unless defined($min_hits);
247 :    
248 :     my %match_set = map { my($id, $com, $seq) = @$_; $id => [$self->match_seq($seq), $seq] } @$seq_set;
249 :    
250 :     my $fr_hash = $self->{friH};
251 :     my $set_hash = $self->{setiH};
252 :    
253 :     my $fig = $self->{fig};
254 :    
255 :     my @missing;
256 :     while (my($id, $ent) = each %match_set)
257 :     {
258 :     my($matches, $seq) = @$ent;
259 :    
260 :     my(%hitsF,%hitsS);
261 :     foreach my $match (@$matches)
262 :     {
263 :     my($offset, $oligo, $frI, $setI) = @$match;
264 :     $hitsF{$frI}++;
265 :     if ($setI)
266 :     {
267 :     $hitsS{$setI}++ ;
268 :     }
269 :     }
270 :    
271 :     my $FRI = &best_hit(\%hitsF,$min_hits);
272 :     my $setI = &best_hit(\%hitsS,$min_hits);
273 :     push(@$ent, $FRI, $setI, \%hitsF);
274 :    
275 :     if (!$fr_hash->{$FRI})
276 :     {
277 :     push(@missing, [$id, undef, $seq]);
278 :     }
279 :     }
280 :    
281 :     #
282 :     # @missing now has the list of sequences that had no Kmer hits. If we have a
283 :     # blast db, blast 'em.
284 :    
285 :     my @all_blastout;
286 :     if (@missing && -s $extra_blastdb)
287 :     {
288 :     #print Dumper(\@missing);
289 :     @all_blastout = ProtSims::blastP(\@missing, $extra_blastdb, 5);
290 :     #print Dumper(\@all_blastout);
291 :     }
292 :    
293 :     #
294 :     # We now have Kmers output and blast output. Go through the original data and
295 :     # create the output.
296 :     #
297 :    
298 :     my @out;
299 :    
300 :     for my $ent (@$seq_set)
301 :     {
302 :     my $id = $ent->[0];
303 :     my ($matches, $seq, $FRI, $setI, $hitsF) = @{$match_set{$id}};
304 :    
305 :     my $blast_results = [];
306 :     if ($fr_hash->{$FRI})
307 :     {
308 :     if ($blast && ($fr_hash->{$FRI} || $set_hash->{$setI}))
309 :     {
310 :     $blast_results = &blast_data($self,$id,$seq,$fr_hash->{$FRI},$blast,'blastp');
311 :     }
312 :    
313 :     push(@out, [$id, $fr_hash->{$FRI},$set_hash->{$setI}, $blast_results,$hitsF->{$FRI}]);
314 :     }
315 :     else
316 :     {
317 :     my @blastout = grep { $_->id1 eq $id } @all_blastout;
318 :    
319 :     if (@blastout > 5) { $#blastout = 4 }
320 :    
321 :     my %hit_pegs = map { $_->id2 => 1 } @blastout;
322 :     my @pegs = keys(%hit_pegs);
323 :     if (@pegs == 0)
324 :     {
325 :     push(@out, [$id,'hypothetical protein','',[],0]);
326 :     }
327 :     else
328 :     {
329 :     my %funcs;
330 :     foreach my $peg (@pegs)
331 :     {
332 :     my $func = $fig->function_of($peg,1);
333 :     if (! &FIG::hypo($func))
334 :     {
335 :     $funcs{$func}++;
336 :     }
337 :     }
338 :     my @pos = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
339 :     my $proposed = (@pos > 0) ? $pos[0] : "hypothetical protein";
340 :     push(@out, [$id, $proposed,'',[],0]);
341 :     }
342 :     }
343 :     }
344 :     return @out;
345 :     }
346 :    
347 :     sub best_hit {
348 :     my($hits,$min_hits) = @_;
349 :     my @poss = sort { $hits->{$b} <=> $hits->{$a} } keys(%$hits);
350 :    
351 :     my $val;
352 :     if ((@poss > 0) && ($hits->{$poss[0]} >= $min_hits))
353 :     {
354 :     $val = $poss[0];
355 :     }
356 :     return $val;
357 :     }
358 :    
359 :     sub best_hit_in_group
360 :     {
361 :     my($group) = @_;
362 :    
363 :     my %hash;
364 :     for my $tuple (@$group)
365 :     {
366 :     my($off,$oligo,$frI,$setI) = @$tuple;
367 :     if ($setI > 0)
368 :     {
369 :     $hash{$setI}++;
370 :     }
371 :     }
372 :     my @sorted = sort { $hash{$b} <=> $hash{$a} } keys %hash;
373 :     my $max = $sorted[0];
374 :     return $max;
375 :     }
376 :    
377 :     sub assign_functions_to_DNA_features {
378 :     my($self,$seq,$min_hits,$max_gap,$blast) = @_;
379 :    
380 :     $min_hits = 3 unless defined($min_hits);
381 :     $max_gap = 200 unless defined($max_gap);
382 :    
383 :     my $fr_hash = $self->{friH};
384 :     my $set_hash = $self->{setiH};
385 :     my $motif_sz = $self->{size};
386 :    
387 :     my %hits;
388 :     my @ans;
389 :     my $matches = $self->process_dna_seq($seq);
390 :    
391 :     push(@ans,&process_hits($self,$matches,1,length($seq),$motif_sz, $min_hits, $max_gap,$blast,$seq));
392 :     undef %hits;
393 :    
394 :     $matches = $self->process_dna_seq(&FIG::reverse_comp($seq));
395 :     push(@ans,&process_hits($self,$matches,length($seq),1,$motif_sz, $min_hits, $max_gap,$blast,$seq));
396 :     return \@ans;
397 :     }
398 :    
399 :     sub process_dna_seq {
400 :     my($self, $seq,$hits) = @_;
401 :    
402 :     my $matches = $self->match_seq($seq);
403 :     return $matches;
404 :     }
405 :    
406 :    
407 :     sub process_hits {
408 :     my($self,$matches,$beg,$end,$sz_of_match, $min_hits, $max_gap,$blast,$seq) = @_;
409 :    
410 :     my $fr_hash = $self->{friH};
411 :     my $set_hash = $self->{setiH};
412 :     my $motif_sz = $self->{size};
413 :    
414 :     my $hits;
415 :     my %sets;
416 :     foreach my $tuple (@$matches)
417 :     {
418 :     my($off,$oligo,$frI,$setI) = @$tuple;
419 :     push(@{$hits->{$frI}},$tuple);
420 :     }
421 :    
422 :     my @got = ();
423 :     my @poss = sort { (@{$hits->{$b}} <=> @{$hits->{$a}}) } keys(%$hits);
424 :     if (@poss != 0)
425 :     {
426 :     foreach my $frI (@poss)
427 :     {
428 :     my $hit_list = $hits->{$frI};
429 :     my @grouped = &group_hits($hit_list, $max_gap);
430 :     foreach my $group_ent (@grouped)
431 :     {
432 :     my($group, $group_hits) = @$group_ent;
433 :     my $N = @$group;
434 :     if ($N >= $min_hits) # consider only runs containing 3 or more hits
435 :     {
436 :     my $b1 = $group->[0];
437 :     my $e1 = $group->[-1] + ($sz_of_match-1);
438 :    
439 :     my $loc;
440 :     if ($beg < $end)
441 :     {
442 :     $loc = [$beg+$b1,$beg+$e1];
443 :     }
444 :     else
445 :     {
446 :     $loc = [$beg-$b1,$beg-$e1];
447 :     }
448 :     my $func = $fr_hash->{$frI};
449 :    
450 :     my $set = &best_hit_in_group($group_hits);
451 :     $set = $set_hash->{$set};
452 :    
453 :     my $blast_output = [];
454 :     if ($blast)
455 :     {
456 :     $blast_output = &blast_data($self,join("_",@$loc),$seq,$func,$blast,
457 :     ($motif_sz == $sz_of_match) ? 'blastn' : 'blastx');
458 :     }
459 :    
460 :     my $tuple = [$N,@$loc,$func,$set,$blast_output];
461 :    
462 :     push(@got,$tuple);
463 :     }
464 :     }
465 :     }
466 :     }
467 :     return @got;
468 :     }
469 :    
470 :     sub group_hits {
471 :     my($hits, $max_gap) = @_;
472 :    
473 :     my @sorted = sort { $a->[0] <=> $b->[0] } @$hits;
474 :     my @groups = ();
475 :     my $position;
476 :     while (defined(my $hit = shift @sorted))
477 :     {
478 :     my($position,$oligo,$frI,$setI) = @$hit;
479 :    
480 :     my $group = [$position];
481 :     my $ghits = [$hit];
482 :     while ((@sorted > 0) && (($sorted[0]->[0] - $position) < $max_gap))
483 :     {
484 :     $hit = shift @sorted;
485 :     ($position,$oligo,$frI,$setI) = @$hit;
486 :     push(@$group,$position+1);
487 :     push(@$ghits, $hit);
488 :     }
489 :    
490 :     push(@groups,[$group, $ghits]);
491 :     }
492 :     return @groups;
493 :     }
494 :    
495 :     sub assign_functions_to_PEGs_in_DNA {
496 :     my($self,$seq,$min_hits,$max_gap,$blast) = @_;
497 :    
498 :     $blast = 0 unless defined($blast);
499 :     $min_hits = 3 unless defined($min_hits);
500 :     $max_gap = 200 unless defined($max_gap);
501 :    
502 :     my $fr_hash = $self->{friH};
503 :     my $set_hash = $self->{setiH};
504 :     my $motif_sz = $self->{size};
505 :    
506 :     my %hits;
507 :     my @ans;
508 :     my $matches = $self->process_prot_seq($seq);
509 :     push(@ans,&process_hits($self,$matches,1,length($seq),3 * $motif_sz, $min_hits, $max_gap,$blast,$seq));
510 :     undef %hits;
511 :     $matches = $self->process_prot_seq(&FIG::reverse_comp($seq));
512 :     push(@ans,&process_hits($self,$matches,length($seq),1,3 * $motif_sz, $min_hits, $max_gap,$blast,$seq));
513 :     return \@ans;
514 :     }
515 :    
516 :     sub process_prot_seq {
517 :     my($self, $seq) = @_;
518 :    
519 :     my $ans = [];
520 :     my $ln = length($seq);
521 :     my($i,$off);
522 :     for ($off=0; ($off < 3); $off++)
523 :     {
524 :     my $ln_tran = int(($ln - $off)/3) * 3;
525 :     my $tran = uc &FIG::translate(substr($seq,$off,$ln_tran));
526 :    
527 :    
528 :     my $matches = $self->match_seq($tran);
529 :    
530 :     push(@$ans, map { $_->[0] = ((3 * $_->[0]) + $off); $_ } @$matches);
531 :     }
532 :     return $ans;
533 :     }
534 :    
535 :     use Sim;
536 :    
537 :     sub blast_data {
538 :     my($self,$id,$seq,$func,$blast,$tool) = @_;
539 :    
540 :     if ($tool eq "blastp")
541 :     {
542 :     return &blast_data1($self,$id,$seq,$func,$blast,$tool);
543 :     }
544 :    
545 :     if ($id =~ /^(\d+)_(\d+)$/)
546 :     {
547 :     my($b,$e) = ($1 < $2) ? ($1,$2) : ($2,$1);
548 :     my $b_adj = (($b - 5000) > 0) ? $b-5000 : 1;
549 :     my $e_adj = (($b + 5000) <= length($seq)) ? $b+5000 : length($seq);
550 :     my $seq1 = substr($seq,$b_adj-1, ($e_adj - $b_adj)+1);
551 :     my $blast_out = &blast_data1($self,$id,$seq1,$func,$blast,$tool);
552 :     foreach $_ (@$blast_out)
553 :     {
554 :     $_->[2] += $b_adj - 1;
555 :     $_->[3] += $b_adj - 1;
556 :     $_->[8] = length($seq);
557 :     }
558 :     return $blast_out;
559 :     }
560 :     else
561 :     {
562 :     return &blast_data1($self,$id,$seq,$func,$blast,$tool);
563 :     }
564 :     }
565 :    
566 :     sub blast_data1 {
567 :     my($self,$id,$seq,$func,$blast,$tool) = @_;
568 :    
569 :    
570 :     if (! $tool) { $tool = 'blastx' }
571 :     my $fig = $self->{fig};
572 :    
573 :     my @blastout = ();
574 :     if ($tool ne 'blastn')
575 :     {
576 :     my $ffs = new FFs($FIG_Config::FigfamsData);
577 :     my @fams = $ffs->families_implementing_role($func);
578 :     foreach my $fam (@fams)
579 :     {
580 :     my $subD = substr($fam,-3);
581 :     my $pegs_in_fam = "$FIG_Config::FigfamsData/FIGFAMS/$subD/$fam/PEGs.fasta";
582 :     push(@blastout,map { [$_->id2,$_->iden,$_->b1,$_->e1,$_->b2,$_->e2,$_->psc,$_->bsc,$_->ln1,$_->ln2,$fam] }
583 :     $fig->blast($id,$seq,$pegs_in_fam,0.1,"-FF -p $tool -b $blast"));
584 :     }
585 :     }
586 :     else
587 :     {
588 :     push(@blastout,map { [$_->id2,$_->iden,$_->b1,$_->e1,$_->b2,$_->e2,$_->psc,$_->bsc,$_->ln1,$_->ln2,$self->{what}] }
589 :     $fig->blast($id,$seq,$self->{blastdb},0.1,"-FF -p $tool -b $blast"));
590 :     }
591 :     @blastout = sort { $b->[7] <=> $a->[7] } @blastout;
592 :     if (@blastout > $blast) { $#blastout = $blast-1 }
593 :     return \@blastout;
594 :     }
595 :    
596 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3