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

Annotation of /FigKernelPackages/ProtFamsLite.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : arodri7 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 ProtFamsLite;
20 :    
21 :     use strict;
22 :     use DB_File;
23 :    
24 :     #use ProtFamLite;
25 :     use Tracer;
26 :    
27 :     use Data::Dumper;
28 :     use Carp;
29 :     use Digest::MD5;
30 :    
31 :     # This is the constructor. Presumably, $class is 'ProtFamsLite'.
32 :     #
33 :    
34 :     sub new {
35 :     my($class,$fam_data) = @_;
36 :    
37 :     my $protfams = {};
38 :    
39 :     defined($fam_data) || return undef;
40 :     $protfams->{dir} = $fam_data;
41 :    
42 :     bless $protfams,$class;
43 :     return $protfams;
44 :     }
45 :    
46 :     sub create_tie
47 :     {
48 :     my($self, $hash, $file, $type) = @_;
49 :    
50 :     my $mode = -w $file ? O_RDWR : O_RDONLY;
51 :    
52 :     return tie %$hash, 'DB_File', $file, $mode, 0666, $type;
53 :     }
54 :    
55 :     sub check_db_role {
56 :     my($self) = @_;
57 :     my $fam_data = $self->{dir};
58 :    
59 :     my $db_role = "$fam_data/role.db";
60 :     if (! -s $db_role) { return undef }
61 :    
62 :     my %role_hash;
63 :     my $role_hash_tie = $self->create_tie(\%role_hash, $db_role, $DB_HASH);
64 :     $role_hash_tie || die "tie $db_role failed: $!";
65 :     $self->{role_db} = \%role_hash;
66 :     }
67 :    
68 :     sub check_db_function {
69 :     my($self) = @_;
70 :     my $fam_data = $self->{dir};
71 :    
72 :     my $db_function = "$fam_data/function.db";
73 :     if (! -s $db_function) { return undef }
74 :     my %function_hash;
75 :     my $function_hash_tie = $self->create_tie(\%function_hash, $db_function, $DB_HASH);
76 :     $function_hash_tie || die "tie $db_function failed: $!";
77 :     $self->{function_db} = \%function_hash;
78 :     }
79 :    
80 :     sub check_db_prot_to_fams {
81 :     my($self) = @_;
82 :     my $fam_data = $self->{dir};
83 :    
84 :     my $db_prot_to_fams = "$fam_data/protein.db";
85 :     if (! -s $db_prot_to_fams) { return undef }
86 :     my %prot_to_fams_hash;
87 :     my $prot_to_fams_hash_tie = $self->create_tie(\%prot_to_fams_hash, $db_prot_to_fams, $DB_HASH);
88 :     $prot_to_fams_hash_tie || die "tie $db_prot_to_fams failed: $!";
89 :     $self->{prot_to_fams_db} = \%prot_to_fams_hash;
90 :     }
91 :    
92 :     sub check_db_genome_to_fams {
93 :     my($self) = @_;
94 :     my $fam_data = $self->{dir};
95 :    
96 :     my $db_genome_to_fams = "$fam_data/genome.db";
97 :     if (! -s $db_genome_to_fams) { return undef }
98 :     my %genome_to_fams_hash;
99 :     my $genome_to_fams_hash_tie = $self->create_tie(\%genome_to_fams_hash, $db_genome_to_fams, $DB_HASH);
100 :     $genome_to_fams_hash_tie || die "tie $db_genome_to_fams failed: $!";
101 :     $self->{genome_to_fams_db} = \%genome_to_fams_hash;
102 :     }
103 :    
104 :     sub check_db_relevant_prot_data {
105 :     my($self) = @_;
106 :     my $fam_data = $self->{dir};
107 :    
108 :     my $db_relevant_prot_data = "$fam_data/relevant.prot.data.db";
109 :     if (! -s $db_relevant_prot_data) { return undef }
110 :     my %relevant_prot_data_hash;
111 :     my $relevant_prot_data_hash_tie = $self->create_tie(\%relevant_prot_data_hash, $db_relevant_prot_data, $DB_HASH);
112 :     $relevant_prot_data_hash_tie || die "tie $db_relevant_prot_data failed: $!";
113 :     $self->{relevant_prot_data_db} = \%relevant_prot_data_hash;
114 :     }
115 :    
116 :     sub check_db_PDB_connections {
117 :     my($self) = @_;
118 :    
119 :     my $fam_data = $self->{dir};
120 :     my $db_PDB_connections = "$fam_data/PDB.connections.db";
121 :     if (! -s $db_PDB_connections) { return undef }
122 :     my %PDB_hash;
123 :     my $PDB_hash_tie = $self->create_tie(\%PDB_hash, $db_PDB_connections, $DB_HASH);
124 :     $PDB_hash_tie || die "tie $db_PDB_connections failed: $!";
125 :     $self->{PDB_connections_db} = \%PDB_hash;
126 :     }
127 :    
128 :     sub check_db_md5_to_fams {
129 :     my($self) = @_;
130 :    
131 :     my $fam_data = $self->{dir};
132 :     my $db_md5_fams = "$fam_data/md5.protfams.db";
133 :     if (! -s $db_md5_fams) { return undef }
134 :     my %md5_hash;
135 :     my $md5_hash_tie = $self->create_tie(\%md5_hash, $db_md5_fams, $DB_HASH);
136 :     $md5_hash_tie || die "tie $db_md5_fams failed: $!";
137 :     $self->{md5_fams_db} = \%md5_hash;
138 :     }
139 :    
140 :     sub check_db_family_id_map {
141 :     my($self) = @_;
142 :     my $fam_data = $self->{dir};
143 :    
144 :     my $db_map = "$fam_data/family_id_map.db";
145 :     if (! -s $db_map) { return undef }
146 :    
147 :     my %map_hash;
148 :     my $map_hash_tie = $self->create_tie(\%map_hash, $db_map, $DB_HASH);
149 :     $map_hash_tie || die "tie $db_map failed: $!";
150 :     $self->{map_db} = \%map_hash;
151 :     }
152 :    
153 :     sub PDB_connections {
154 :     my($self,$fam,$raw) = @_;
155 :    
156 :     $self->check_db_PDB_connections;
157 :     my $sims = $self->{PDB_connections_db}->{$fam};
158 :     my @sims = map { $_ =~ /pdb\|([0-9a-zA-Z]+)/; [$1,[split(/\t/,$_)]] } split(/\n/,$sims);
159 :     if (! $raw) { @sims = map { $_->[0] } grep { ($_->[1]->[11] > 0.5) && ((($_->[1]->[4] - $_->[1]->[3]) / $_->[1]->[5]) > 0.8) } @sims}
160 :     return \@sims;
161 :     }
162 :    
163 :     sub function_of {
164 :     my($self,$prot,$ignore_comments) = @_;
165 :    
166 :     $self->check_db_relevant_prot_data; # lazily tie DB
167 :    
168 :     my $prot_data;
169 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
170 :     {
171 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
172 :     if ($ignore_comments)
173 :     {
174 :     $func =~ s/\s*\#.*$//;
175 :     }
176 :     return $func;
177 :     }
178 :     return "";
179 :     }
180 :    
181 :     sub get_translation {
182 :     my($self,$prot,$ignore_comments) = @_;
183 :    
184 :     $self->check_db_relevant_prot_data; # lazily tie DB
185 :    
186 :     my $prot_data;
187 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
188 :     {
189 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
190 :     return $seq;
191 :     }
192 :     return "";
193 :     }
194 :    
195 :    
196 :     sub md5_function_of {
197 :     my ($self, $md5, $preferred_db) = @_;
198 :    
199 :     # find if md5 is in a family
200 :    
201 :     # if in_family get the function
202 :    
203 :     # else go to seed to get function
204 :     my (@tuples, $function);
205 :     if (@tuples = $self->{fig}->mapped_prot_ids(qq(gnl|md5|$md5))) {
206 :     my @tmp;
207 :     if (@tmp = grep { $_->[0] =~ m/^fig\|\d+\.\d+\.peg\.\d+$/o } @tuples) {
208 :     #...Pick the function of one of the FIDs
209 :     }
210 :     elsif (@tmp = grep { $_->[0] =~ m/^sp\|\S+$/o } @tuples) {
211 :     #...Pick the function of one of the SP's
212 :     }
213 :     elsif (@tmp = grep { $_->[0] =~ m/^pir\|\S+$/o } @tuples) {
214 :     #...DWIM
215 :     }
216 :     elsif (@tmp = grep { $_->[0] =~ m/^kegg\|\S+$/o } @tuples) {
217 :     }
218 :     elsif (@tmp = grep { $_->[0] =~ m/^gi\|\S+$/o } @tuples) {
219 :     }
220 :     else {
221 :     $function = qq(Hypothetical protein);
222 :     }
223 :     return $function;
224 :     }
225 :     }
226 :    
227 :     sub org_of {
228 :     my($self,$prot) = @_;
229 :    
230 :     $self->check_db_relevant_prot_data; # lazily tie DB
231 :    
232 :     my $prot_data;
233 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
234 :     {
235 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
236 :     return $org;
237 :     }
238 :     return "";
239 :     }
240 :    
241 :     sub path_to {
242 :     my ($self,$family_name) = @_;
243 :    
244 :     $self->check_db_family_id_map; # lazily tie DB
245 :    
246 :     my $internal_id;
247 :     if ($internal_id = $self->{family_id_map_db}->{$family_name})
248 :     {
249 :     #my($internal_id) = split(/\t/,$family_name);
250 :    
251 :     my $group = substr($internal_id, -3);
252 :     $group = (qq(0) x (3-length($group))) . $group;
253 :    
254 :     my $path_to_family = qq($self->{dir}/FAMS/$group/$internal_id);
255 :    
256 :     return $path_to_family;
257 :     }
258 :     return "";
259 :     }
260 :    
261 :     sub seq_of {
262 :     my($self,$prot) = @_;
263 :    
264 :     $self->check_db_relevant_prot_data; # lazily tie DB
265 :    
266 :     my $prot_data;
267 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
268 :     {
269 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
270 :     return $seq;
271 :     }
272 :     return "";
273 :     }
274 :    
275 :     sub aliases_of {
276 :     my($self,$prot) = @_;
277 :    
278 :     $self->check_db_relevant_prot_data; # lazily tie DB
279 :    
280 :     my $prot_data;
281 :     if ($prot_data = $self->{relevant_prot_data_db}->{$prot})
282 :     {
283 :     my($org,$func,$seq,$aliases) = split(/\t/,$prot_data);
284 :     return wantarray() ? split(/,/,$aliases) : $aliases;
285 :     }
286 :     return "";
287 :     }
288 :    
289 :     sub families_implementing_role {
290 :     my($self,$role) = @_;
291 :    
292 :     $self->check_db_role; # lazily tie DB
293 :    
294 :     my $fams = $self->{role_db}->{$role};
295 :     return $fams ? split("\t",$fams) : ();
296 :     }
297 :    
298 :     sub families_with_function {
299 :     my($self,$function) = @_;
300 :    
301 :     $self->check_db_function; # lazily tie DB
302 :    
303 :     my $fams = $self->{function_db}->{$function};
304 :     return $fams ? split("\t",$fams) : ();
305 :     }
306 :    
307 :     sub family_containing_prot {
308 :     my($self,$prot) = @_;
309 :     my @fams = $self->families_containing_prot($prot);
310 :     return (@fams > 0) ? $fams[0] : undef;
311 :     }
312 :    
313 :     sub families_containing_prot {
314 :     my($self,$prot) = @_;
315 :    
316 :     $self->check_db_prot_to_fams; # lazily tie DB
317 :     my $fams = $self->{prot_to_fams_db}->{$prot};
318 :     return $fams ? split("\t",$fams) : ();
319 :     }
320 :    
321 :     sub families_in_genome {
322 :     my($self,$genome) = @_;
323 :    
324 :     $self->check_db_genome_to_fams; # lazily tie DB
325 :    
326 :     my $fams = $self->{genome_to_fams_db}->{$genome};
327 :     return $fams ? split("\t",$fams) : ();
328 :     }
329 :    
330 :     sub families_containing_md5 {
331 :     my ($self,$md5) = @_;
332 :    
333 :     $self->check_db_md5_to_fams; # lazily tie DB
334 :    
335 :     my $fams = $self->{md5_fams_db}->{$md5};
336 :     if ($fams){
337 :     return $fams;
338 :     }
339 :     else{
340 :     return undef;
341 :     }
342 :     }
343 :    
344 :     sub all_families {
345 :     my($self) = @_;
346 :    
347 :     return sort map { chomp; $_ } `cut -f1 $self->{dir}/family.functions`;
348 :     }
349 :    
350 :     sub is_prots_in_family {
351 :     my($self,$seqs) = @_;
352 :     my $id_2_protfam = {};
353 :    
354 :     my $protD = $self->{dir};
355 :     my $protfam_file = "$protD/families.3c";
356 :    
357 :     my $contents;
358 :     sysopen(FAM_PROTS, $protfam_file, 0) or die "could not open file '$protfam_file': $!";
359 :     sysread(FAM_PROTS, $contents, 1000000000);
360 :     close(FAM_PROTS) or die "could not close file '$protfam_file': $!";
361 :     my %protfam_prots = map { $_ =~ /^(\S+)\t(\S+)\t/; $2 => $1 } split("\n", $contents);
362 :    
363 :     foreach (@$seqs){
364 :     if ($protfam_prots{$_}){
365 :     $id_2_protfam->{$_} = $protfam_prots{$_};
366 :     }
367 :     }
368 :     return $id_2_protfam;
369 :     }
370 :    
371 :     sub place_in_family {
372 :     my($self,$seq,$debug,$loose,$debug_prot,$nuc) = @_;
373 :     my($protfam,$should_be,$sims);
374 :    
375 :     my $old_sep = $/;
376 :     $/ = "\n";
377 :    
378 :     my $dir = $self->{dir};
379 :     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";
380 :     open(TMP,">$tmpF")
381 :     || die "could not open $tmpF";
382 :     print TMP ">query\n$seq\n";
383 :     close(TMP);
384 :    
385 :     # check if sequence is in any protfam by md5 checksum
386 :     my $md5 = Digest::MD5::md5_hex( uc $seq );
387 :     my $fams = $self->families_containing_md5($md5);
388 :     if (($fams) && ($fams !~ /\,/)){ # do regular check if the checksum returns more than one protfam
389 :     my $got = new ProtFamLite($fams,$self);
390 :     ($should_be, $sims) = $got->should_be_member($seq,$debug,$loose,$debug_prot,$nuc);
391 :     return ($got,$sims);
392 :     }
393 :    
394 :     if ($nuc){
395 :     # open(BLAST,"blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g F | head -n 20 |")
396 :     Open(\*BLAST,"$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T |");
397 :     print STDERR "blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T\n" if ($debug);
398 :     }
399 :     else{
400 :     # open(BLAST,"blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastp | head -n 20 |")
401 :     Open(\*BLAST, "$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastp |");
402 :    
403 :     }
404 :    
405 :     my(%seen);
406 :     my $got = undef;
407 :     my $min_sc = 1;
408 :     my $checked = 0;
409 :     while ((! $got) && ($checked < 30) && (defined($_ = <BLAST>)))
410 :     {
411 :     if ($debug) { print STDERR $_ }
412 :     chop;
413 :     my $sim = [split(/\t/,$_)];
414 :     if ($min_sc > $sim->[10]) { $min_sc = $sim->[10] }
415 :     my $fam_id = ($sim->[1] =~ /(.*)-/) ? $1 : "";
416 :     if (! $seen{$fam_id})
417 :     {
418 :     $seen{$fam_id} = 1;
419 :    
420 :     # skip if the e-value is 1 or higher
421 :     if ($sim->[10] >= 1){
422 :     print STDERR "SIM evalue: " . $sim->[10] . "\n" if ($debug);
423 :     $checked++;
424 :     next;
425 :     }
426 :    
427 :     $protfam = new PrtoFamLite($fam_id,$dir);
428 :    
429 :     if (not defined($protfam))
430 :     {
431 :     next;
432 :     }
433 :     else
434 :     {
435 :     $checked++;
436 :     if ($debug) { print STDERR "checking family $fam_id ",$protfam->family_function,"\n" }
437 :     ($should_be, $sims) = $protfam->should_be_member($seq,$debug,$loose,$debug_prot,$nuc);
438 :     if ($debug) { print STDERR " should_be=$should_be\n" }
439 :     my $psc;
440 :     if ($should_be && defined($psc = &psc($sims->[0])) && ($psc <= $min_sc))
441 :     {
442 :     if ($debug) { print STDERR "min_sc=$min_sc best=$psc\n" }
443 :     $got = $protfam;
444 :     }
445 :     }
446 :     }
447 :     }
448 :    
449 :     if ($debug)
450 :     {
451 :     while ($_ && (defined($_ = <BLAST>))) ### JUST PRINT OUT REMAINING BLASTS AGAINST REPDB
452 :     {
453 :     print STDERR $_;
454 :     }
455 :     }
456 :     $/ = $old_sep;
457 :     close(BLAST);
458 :    
459 :     unlink($tmpF);
460 :     return $got ? ($got,$sims) : (undef,undef);
461 :     }
462 :    
463 :     =head3
464 :     usage: my $out = $protfams->place_in_family_bulk($seqs);
465 :    
466 :     The genereal overview of this function is that it does an initial blast of all
467 :     input sequences against the repdb with one call, and then does individual blast
468 :     calls for each input sequence to determine to which family it should be placed.
469 :    
470 :     input
471 :     the only difference to place_in_family is that the $seqs parameter is a
472 :     list of sequences to be placed in a family.
473 :     The list can be either a number of sequences in fasta format or an array
474 :     of just the sequences to be placed in a family
475 :    
476 :     i.e. ->>>> [">seq1\nADTGGHHH\n", ">seq2\nRARGHTGK\n", ....]
477 :     or
478 :     ["ADTGGHHH", "RARGHTGK", ....]
479 :    
480 :     output
481 :     return reference to a hash where the keys are the input sequence
482 :     and the values are tuples same as what place in family returns: [$family, $sims]
483 :    
484 :     $family is a hash reference to the Protfam it was placed in. If the sequence
485 :     was not placed in a family, then it will be an empty tuple.
486 :    
487 :     $sims is a reference to the similarities of the Protfam selected
488 :    
489 :    
490 :     =cut
491 :    
492 :     sub place_in_family_bulk {
493 :     my($self,$seqs,$debug,$loose,$debug_prot,$nuc) = @_;
494 :     my($protfam,$should_be,$sims);
495 :     my $out = {};
496 :     my @run_blast = ();
497 :    
498 :     my $dir = $self->{dir};
499 :     my $old_sep = $/;
500 :     $/ = "\n";
501 :     my $seq_list = {};
502 :     my $count = 0;
503 :    
504 :     foreach my $fasta (@$seqs){
505 :     my ($id, $header, $seq);
506 :     if ($fasta =~ m/^>(.*?)\n(.*?)\n/){
507 :     #($header,$seq) = $fasta =~ m/^>(.*?)\n(.*?)\n/;
508 :     $header = $1;
509 :     $seq = $2;
510 :     ($id) = split (/\s+/, $header);
511 :     }
512 :     else{
513 :     $id = "fasta".$count;
514 :     $seq = $fasta;
515 :     $count++;
516 :     $fasta = ">$id\n$seq\n";
517 :     }
518 :    
519 :    
520 :     # check if sequence is in any protein fam by md5 checksum
521 :     my $fams;
522 :     if (!$nuc){
523 :     my $md5 = Digest::MD5::md5_hex( uc $seq );
524 :     $fams = $self->families_containing_md5($md5);
525 :     }
526 :     if (($fams) && ($fams !~ /\,/)){ # do regular check if the checksum returns more than one prtofam
527 :     my $got = new ProtFamLite($fams,$self);
528 :     ($should_be, $sims) = $got->should_be_member($seq,$debug,$loose,$debug_prot,$nuc);
529 :     $out->{$seq} = [$got, $sims];
530 :     #push (@$out, [$seq, $got, $sims]);
531 :     }
532 :     else{
533 :     $seq_list->{$id} = $seq;
534 :     push (@run_blast, $fasta);
535 :     }
536 :     }
537 :    
538 :     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";
539 :     open(TMP,">$tmpF")
540 :     || die "could not open $tmpF";
541 :     print TMP @run_blast;
542 :     close(TMP);
543 :     my @tmp;
544 :    
545 :     if ($nuc){
546 :     @tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T` || die "could not open blast";
547 :     print STDERR "$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T\n" if ($debug);
548 :     }
549 :     else{
550 :     @tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastp` || die "could not open blast";
551 :     print STDERR "$FIG_Config::ext_bin/blastall -i $tmpF -d $dir/repdb -m 8 -FF -p blastx -g T\n" if ($debug);
552 :     }
553 :    
554 :     my $prev_id;
555 :     my $seen = {};
556 :     my $got = {};
557 :     my $min_sc = {};
558 :     my $checked = {};
559 :    
560 :     while ( (defined( $_ = shift @tmp)) )
561 :     {
562 :     my $sim = [split(/\t/,$_)];
563 :    
564 :     next if ($got->{$sim->[0]});
565 :     next if ($checked->{$sim->[0]} > 30);
566 :     if (!$min_sc->{$sim->[0]}){
567 :     $min_sc->{$sim->[0]} = 1;
568 :     }
569 :    
570 :     if ($debug) { print STDERR $_ }
571 :     chop;
572 :     if ($min_sc->{$sim->[0]} > $sim->[10]) { $min_sc->{$sim->[0]} = $sim->[10] }
573 :     my $fam_id = ($sim->[1] =~ /(.*)-/) ? $1 : "";
574 :     if (! $seen->{$sim->[0]}->{$fam_id})
575 :     {
576 :     $seen->{$sim->[0]}->{$fam_id} = 1;
577 :    
578 :     # skip if the e-value is 1 or higher
579 :     if ($sim->[10] >= 1){
580 :     #print STDERR "SIM evalue: " . $sim->[10] . "\n" if ($debug);
581 :     $checked->{$sim->[0]}++;
582 :     next;
583 :     }
584 :    
585 :     $protfam = new ProtFamLite($fam_id,$self);
586 :    
587 :     if (not defined($protfam))
588 :     {
589 :     next;
590 :     }
591 :     else
592 :     {
593 :     $checked->{$sim->[0]}++;
594 :     if ($debug) { print STDERR "checking family $fam_id ",$protfam->family_function,"\n" }
595 :     my ($should_be, $sims) = $protfam->should_be_member($seq_list->{$sim->[0]},$debug,$loose,$debug_prot,$nuc);
596 :     if ($debug) { print STDERR " should_be=$should_be\n" }
597 :     my $psc;
598 :     if ($should_be && defined($psc = &psc($sims->[0])) && ($psc <= $min_sc->{$sim->[0]}))
599 :     {
600 :     if ($debug) { print STDERR "min_sc=$min_sc best=$psc\n" }
601 :     $got->{$sim->[0]} = $protfam;
602 :     #push (@$out, [$seq_list->{$sim->[0]}, $got->{$sim->[0]}, $sims]);
603 :     $out->{$seq_list->{$sim->[0]}} = [$got->{$sim->[0]}, $sims];
604 :     }
605 :     }
606 :     }
607 :    
608 :     # if ($debug)
609 :     # {
610 :     # my $blast_line;
611 :     # while ($blast_line && (defined($blast_line = shift @blast_tmp))) ### JUST PRINT OUT REMAINING BLASTS AGAINST REPDB
612 :     # {
613 :     # print STDERR $blast_line;
614 :     # }
615 :     # }
616 :     #
617 :     # if (!$got){
618 :     # push (@$out, [$prev_id, undef, undef]);
619 :     # }
620 :    
621 :     }
622 :    
623 :     $/ = $old_sep;
624 :     # close(BLAST);
625 :    
626 :     foreach my $id (keys %$seq_list)
627 :     {
628 :     if (!$got->{$id})
629 :     {
630 :     #push (@$out, [$seq_list->{$id}, undef, undef]);
631 :     $out->{$seq_list->{$id}} = [undef, undef];
632 :     }
633 :     }
634 :    
635 :     unlink($tmpF);
636 :     return $out;
637 :     }
638 :    
639 :     sub psc {
640 :     my ($sim) = @_;
641 :     return ($sim->[10] =~ /^e-/) ? "1.0" . $sim->[10] : $sim->[10];
642 :     }
643 :    
644 :    
645 :     =head3
646 :     usage: $protfams->family_functions();
647 :    
648 :     returns a hash of all the functions for all protein fams from the family.functions file
649 :    
650 :     =cut
651 :    
652 :     sub family_functions {
653 :     my($self) = @_;
654 :     my $ffD = $self->{dir};
655 :     my $ff_file = "$ffD/family.functions";
656 :    
657 :     my $contents;
658 :     sysopen(FAM_FUNC, $ff_file, 0) or die "could not open file '$ff_file': $!";
659 :     sysread(FAM_FUNC, $contents, 1000000000);
660 :     close(FAM_FUNC) or die "could not close file '$ff_file': $!";
661 :     my %ff_name = map {split(/\t/)} split("\n", $contents);
662 :    
663 :     return \%ff_name;
664 :     }
665 :    
666 :    
667 :    
668 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3