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

Annotation of /FigKernelPackages/ProtFamLite.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 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 ProtFamLite;
20 :    
21 :     use strict;
22 :     use Carp;
23 :     use Data::Dumper;
24 :    
25 :     =head1 Module to access a Protein Family
26 :    
27 :     =head3 new
28 :    
29 :     usage:
30 :     my $ProtFam_Obj = ProtFamLite->new($ProtFams_Obj, $fam_id);
31 :    
32 :     C<$fam_id> is the ID of the family, of the form C<n+> where C<n+> is one or more digits;
33 :     it is required.
34 :    
35 :     C<$ProtFams_Obj> is required, since it is the Families-Object
36 :     that provides access to the desired collection of protein-family data.
37 :    
38 :     =cut
39 :    
40 :     sub new {
41 :     my ($class, $ProtFams_Obj, $fam_id) = @_;
42 :    
43 :     defined($ProtFams_Obj) || confess "ProtFams_Obj is undefined";
44 :    
45 :     my $fig = $ProtFams_Obj->{fig} ||
46 :     confess "ProtFams_Obj does not contain a FIG object";
47 :    
48 : arodri7 1.2 # ($fam_id =~ /\d+$/) || confess "invalid family id: $fam_id";
49 : gdpusch 1.1
50 :     my $fam = {};
51 :     $fam->{id} = $fam_id;
52 :     $fam->{root} = $ProtFams_Obj->{root};
53 : arodri7 1.2 $fam->{fig} = $fig;
54 : gdpusch 1.1
55 :     my $fams_dir = qq($fam->{root}/FAMS);
56 : arodri7 1.2 # my $dir = &fam_dir($ProtFams_Obj, $fam_id);
57 :     my $dir = $ProtFams_Obj->path_to($fam_id);
58 : gdpusch 1.1 $fam->{dir} = $dir;
59 :     (-d $dir) || return undef;
60 :    
61 :     $fam->{function} = $fig->file_read( qq($dir/function), 1) || qq();
62 : arodri7 1.2 chomp $fam->{function};
63 :     $fam->{protfams_obj} = $ProtFams_Obj;
64 : gdpusch 1.1
65 :     my ($prot, $prots);
66 :     my $protsL = [ map { $_ =~ /^(\S+)/ ? ($1) : () }
67 :     $fig->file_read( qq($dir/PROTS), qq(*))
68 :     ];
69 :    
70 :     if (@$protsL < 2) {
71 :     return undef;
72 :     }
73 :    
74 :     $fam->{protsL} = $protsL;
75 :     my $protsH = {};
76 :     foreach $prot (@$protsL) {
77 :     $protsH->{$prot} = 1;
78 :     }
79 :     $fam->{protsH} = $protsH;
80 :    
81 :     if (-s "$dir/PROTS.fasta") {
82 :     open(FASTA,"<$dir/PROTS.fasta") || die "could not read-open $dir/PROTS.fasta";
83 :     while (my ($prot, $seqP) = &read_fasta_record(\*FASTA)) {
84 :     $fam->{prot_lengths}->{$prot} = length($$seqP);
85 :     }
86 :     close(FASTA);
87 :     }
88 :     else {
89 :     confess "$fam_id is missing PROTS.fasta";
90 :     }
91 :    
92 : arodri7 1.2 if ((-f $fam->{root}."/FIG") && (&use_ross_bounds($class,$dir)) && (-s "$dir/bounds")) {
93 : gdpusch 1.1 $fam->{bounds} = &load_bounds("$dir/bounds");
94 :     }
95 :    
96 :     bless $fam,$class;
97 :     return $fam;
98 :     }
99 :    
100 :     sub use_blast {
101 :     my($dir) = @_;
102 :    
103 :     return ((-s "$dir/decision.procedure") && &which_dec("$dir/dec",'blast'));
104 :     }
105 :    
106 :     sub use_ross_bounds {
107 :     my ($self, $dir) = @_;
108 :    
109 : arodri7 1.2 return ( (! (-s "$dir/decision.procedure"))
110 : gdpusch 1.1 ||
111 :     &which_dec("$dir/dec",'ross')
112 :     );
113 :     }
114 :    
115 :     sub which_dec {
116 :     my($dir,$pat) = @_;
117 :    
118 :     if (open(DEC,"<$dir/decision.procedure"))
119 :     {
120 :     my $x = <DEC>;
121 :     close(DEC);
122 :     if ($x)
123 :     {
124 :    
125 :     return ($x =~ /^$pat/);
126 :     }
127 :     }
128 :     return 0;
129 :     }
130 :    
131 :     sub load_bounds {
132 :     my($file) = @_;
133 :    
134 :     my $bounds;
135 :     if (open(BOUNDS,"<$file"))
136 :     {
137 :     $bounds = {};
138 :     my $x;
139 :     while (defined($x = <BOUNDS>))
140 :     {
141 :     chomp $x;
142 :     my @flds = split(/\t/,$x);
143 :     my $prot = shift @flds;
144 :     $bounds->{$prot} = [@flds];
145 :     }
146 :     close(BOUNDS);
147 :     }
148 :     return $bounds;
149 :     }
150 :    
151 :     =head3 prots_of
152 :    
153 :     usage:
154 :     print $protfam_obj->prots_of();
155 :    
156 :     Returns a list of just prots.
157 :    
158 :     =cut
159 :    
160 :     sub prots_of {
161 :     my($self) = @_;
162 :     return [$self->list_members];
163 :     }
164 :    
165 :     =head3 list_members
166 :    
167 :     usage:
168 :     @ids = $protfam_obj->list_members();
169 :    
170 :     Returns a list of the PROT FIDs in a family.
171 :    
172 :     =cut
173 :    
174 :     sub list_members {
175 :     my ($self) = @_;
176 :    
177 :     my $fam_dir = $self->{dir};
178 :     my @prots = map { ($_ =~ /^(\S+)/) ? ($1) : ()
179 :     } $self->{fig}->file_read(qq($fam_dir/PROTS), qq(*));
180 :    
181 :     return @prots;
182 :     }
183 :    
184 :     =head3 representatives
185 :    
186 :     usage:
187 :     C<< @rep_seqs = $protfam_obj->representatives(); >>
188 :    
189 :     Returns a list of the "representative sequences" characterizing a FIGfam.
190 :    
191 :     =cut
192 :    
193 :     sub representatives {
194 :     my($self) = @_;
195 :    
196 :     my $reps = $self->{reps};
197 :     if (! $reps)
198 :     {
199 :     my $rep_file = "$self->{dir}/reps";
200 :     $reps = $self->{reps} = (-s $rep_file) ? [map { $_ =~ /(\S+)/; $1 } `fgrep -v ">" $rep_file`] : [];
201 :     }
202 :     return @$reps;
203 :     }
204 :    
205 :    
206 :     =head3 should_be_member
207 :    
208 :     usage:
209 :     if ( $protfam_obj->should_be_member( $seq ) ) { #...do something... }
210 :    
211 :     Returns ($placed,$sims). $placed will be
212 :     C<TRUE> if the protein sequence in C<$seq> is judged to be
213 :     "similar enough" to the members of a family to potentially be included.
214 :    
215 :     I have added the "loose" argument as an optional last argument. This means that
216 :    
217 :     if ( $protfam_obj->should_be_member( $seq,0,1 ) ) { #...do something... }
218 :    
219 :     will return true, even if the input sequence is truncated (i.e., we do not force the
220 :     similarity to go across 80% of matched sequences in the family).
221 :    
222 :     =cut
223 :    
224 :     sub should_be_member {
225 :     my($self,$seq,$debug,$loose,$debug_prot,$nuc) = @_;
226 :    
227 :     my $old_eol = $/;
228 :     $/ = "\n";
229 :    
230 :     my $dir = $self->{dir};
231 :     my($in,@rc);
232 :     if ( ((open(DEC,"<$dir/decision.procedure") && ($in = <DEC>) && ($in =~ /^(\S+)(\s+(\S.*\S))?/) )) ||
233 :     (($nuc) && open(DEC,"<$dir/decision.procedure.blast") && ($in = <DEC>) && ($in =~ /^(\S+)(\s+(\S.*\S))?/) ) )
234 :     # if ( ((open(DEC,"<$dir/decision.procedure") && ($in = <DEC>) && ($in =~ /^(\S+)(\s+(\S.*\S))?/) )))
235 :     {
236 :     no strict 'refs';
237 :    
238 :     my $procedure = $1;
239 :     my @args = $3 ? split(/\s+/,$3) : ();
240 :     close(DEC);
241 :     @rc = &{$procedure}($self,$debug,$loose,$seq,$dir,$debug_prot,$nuc,@args);
242 :     }
243 :     else
244 :     {
245 : arodri7 1.2 # @rc = &ross_hack($self,$debug,$loose,$seq,$dir,$debug_prot,$nuc);
246 :    
247 :     open(DEC,"<$dir/decision.procedure.blast");
248 :     $in = <DEC>;
249 :     $in =~ /^(\S+)(\s+(\S.*\S))?/;
250 :     my $procedure = $1;
251 :     my @args = $3 ? split(/\s+/,$3) : ();
252 :     close(DEC);
253 :    
254 :     @rc = &blast_vote($self,$debug,$loose,$seq,$dir,$debug_prot,$nuc,@args);
255 : gdpusch 1.1 }
256 :     $/ = $old_eol;
257 :     return @rc;
258 :     }
259 :    
260 :     sub min {
261 :     my($x,$y) = @_;
262 :    
263 :     return ($x <= $y) ? $x : $y;
264 :     }
265 :    
266 :     use ProtFamsLite;
267 :     sub blast_vote {
268 :     my($self,$debug, $loose, $seq,$dir,$debug_prot,$nuc,$partition,$min_bsc) = @_;
269 :    
270 :     if ($ENV{DEBUG}) { $debug = 1 }
271 :     (-s "$dir/PROTS") || return undef;
272 :    
273 :     my $PFsO = ProtFamsLite->new($self->{root});
274 : arodri7 1.2 my $fig = $self->{fig};
275 :    
276 : gdpusch 1.1 if ($debug) { print STDERR "checking: ",$self->family_id," min_bsc=$min_bsc\n" }
277 : arodri7 1.2 my %prots_in = map { $_ =~ /(\S+)/; $1 => 1 } $self->{fig}->file_read( qq($dir/PROTS), qq(*) );
278 : gdpusch 1.1
279 :     my $N = &min(10,scalar keys(%prots_in) - ($debug_prot ? 1 : 0));
280 :     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";
281 :     open(TMP,">$tmpF")
282 :     || die "could not open tmp$$.fasta";
283 :     print TMP ">query\n$seq\n";
284 :     close(TMP);
285 :    
286 :     my $query_ln = length($seq);
287 :     my $partitionF = "$self->{root}/Partitions/" . ($partition % 1000) . "/$partition/fasta";
288 :     my @tmp;
289 :     if ($nuc){
290 :     #@tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $partitionF -m 8 -FF -p blastx -g F`;
291 :     @tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $partitionF -m 8 -FF -p blastx -g T`;
292 :     print STDERR "blastall -i $tmpF -d $partitionF -m 8 -FF -p blastx -g T\n" if ($debug);
293 :     }
294 :     else{
295 :     @tmp = `$FIG_Config::ext_bin/blastall -i $tmpF -d $partitionF -m 8 -FF -p blastp`;
296 :     }
297 :     unlink($tmpF);
298 :    
299 :     my $sims = [];
300 :     my $in = 0;
301 :     my $out = 0;
302 :     my(%seen);
303 :     for (my $simI=0; ($simI < @tmp); $simI++)
304 :     {
305 :     $_ = $tmp[$simI];
306 :     if ($debug) { print STDERR $_ }
307 :     chop;
308 :    
309 :     my $sim = [split(/\t/,$_)];
310 :     my $prot = $sim->[1];
311 : arodri7 1.2 #next if ((! -f "$self->{root}/FIG") && ($prot =~ /fig\|/));
312 : gdpusch 1.1 next if ($debug_prot && ($debug_prot eq $prot));
313 :     my $bit_score = $sim->[11];
314 :     my $matched1 = abs($sim->[7] - $sim->[6]) + 1;
315 :     my $matched2 = abs($sim->[9] - $sim->[8]) + 1;
316 :     if ($debug) { print STDERR "normalized bit-score=",sprintf("%3.2f",$bit_score / $matched2),"\n" }
317 :     if (! $seen{$prot})
318 :     {
319 :     $seen{$prot} = 1;
320 : arodri7 1.2 my $count_in = &count_in($self,$PFsO,$fig,\%prots_in,$loose,$prot);
321 :     my $ln2 = &get_len2($self,$PFsO,$fig,$prot,\%prots_in);
322 : gdpusch 1.1
323 :     if (sprintf("%3.2f",($bit_score / $matched2)) >= $min_bsc)
324 :     {
325 :     if ($nuc){
326 :     if ($count_in && # (print "in-ok\n") &&
327 :     $ln2 &&
328 :     ($matched1 >= (0.7 * $query_ln)) # (print "mat1-ok\n") &&
329 :     )
330 :     {
331 :     push @$sim, $query_ln, $self->{prot_lengths}->{$prot};
332 :     bless $sim, 'Sim';
333 :     push @$sims, $sim;
334 :     if ($N > 0)
335 :     {
336 :     $in++;
337 :     }
338 :     }
339 :     else
340 :     {
341 :     if ($N > 0)
342 :     {
343 :     $out++;
344 :     }
345 :     }
346 :     if ($debug) {print STDERR &Dumper(["in=$in",
347 :     "out=$out",
348 :     $count_in,
349 :     "ln1=$query_ln",
350 :     "ln2=$ln2",
351 :     "matched1=$matched1",
352 :     "matched2=$matched2",
353 :     "bsc=$bit_score",
354 :     $ln2 ? sprintf("%3.2f",$bit_score/$matched2) : ""]); }
355 :     if ($N > 0)
356 :     {
357 :     $N--;
358 :     last if ($N <= 0);
359 :     }
360 :     else
361 :     {
362 :     last;
363 :     }
364 :     }
365 :     else{
366 :     if ($count_in && # (print "in-ok\n") &&
367 :     $ln2 &&
368 :     ($matched1 >= (0.7 * $query_ln)) && # (print "mat1-ok\n") &&
369 :     ($matched2 >= (0.7 * $ln2)) # (print "mat2-ok\n") &&
370 :     )
371 :     {
372 :     push @$sim, $query_ln, $self->{prot_lengths}->{$prot};
373 :     bless $sim, 'Sim';
374 :     push @$sims, $sim;
375 :     if ($N > 0)
376 :     {
377 :     $in++;
378 :     }
379 :     }
380 :     else
381 :     {
382 :     if ($N > 0)
383 :     {
384 :     $out++;
385 :     }
386 :     }
387 :     if ($debug) {print STDERR &Dumper(["in=$in",
388 :     "out=$out",
389 :     $count_in,
390 :     "ln1=$query_ln",
391 :     "ln2=$ln2",
392 :     "matched1=$matched1",
393 :     "matched2=$matched2",
394 :     "bsc=$bit_score",
395 :     $ln2 ? sprintf("%3.2f",$bit_score/$matched2) : ""]); }
396 :     if ($N > 0)
397 :     {
398 :     $N--;
399 :     last if ($N <= 0);
400 :     }
401 :     else
402 :     {
403 :     last;
404 :     }
405 :     }
406 :     }
407 :     else{
408 :     $out++;
409 :     $N--;
410 :     last if ($N <= 0);
411 :     }
412 :     }
413 :     }
414 :     if ($debug) { print STDERR "in=$in out=$out, FINAL VOTE\n" }
415 :     return (($in > $out),$sims);
416 :     }
417 :    
418 :     sub get_len2 {
419 : arodri7 1.2 my($self,$PFsO,$fig,$prot,$prots_in) = @_;
420 : gdpusch 1.1
421 :     if ($prots_in->{$prot})
422 :     {
423 :     return $self->{prot_lengths}->{$prot};
424 :     }
425 : arodri7 1.2 elsif ($prot =~ /fig\|\d+/)
426 :     {
427 :     return $fig->translation_length($prot);
428 :     }
429 : gdpusch 1.1 else
430 :     {
431 :     return length($PFsO->seq_of($prot));
432 :     }
433 :     }
434 :    
435 : arodri7 1.2 use Digest::MD5;
436 : gdpusch 1.1 sub count_in {
437 : arodri7 1.2 my($self,$PFsO,$fig,$prots_in,$loose,$prot) = @_;
438 : gdpusch 1.1
439 :     if ($prots_in->{$prot}) { return 1 }
440 : arodri7 1.2
441 :     # figure out if it should be count_in by md5 sequence
442 :     my $seq;
443 :     if ($prot =~ /^fig\|/)
444 :     {
445 :     $seq = $fig->get_translation($prot);
446 :     }
447 :     else
448 :     {
449 :     $seq = $PFsO->get_translation($prot);
450 :     }
451 :    
452 :     my $md5 = Digest::MD5::md5_hex( uc $seq );
453 :     my @proteins_with_md5;
454 :     push (@proteins_with_md5, split (/\n/, $PFsO->proteins_containing_md5($md5)));
455 :    
456 :     foreach my $md5_prot (@proteins_with_md5)
457 :     {
458 :     if ($prots_in->{$md5_prot}) { return 1 }
459 :     }
460 :    
461 : gdpusch 1.1 if (! $loose) { return 0 }
462 :     my $fam_func = $self->family_function;
463 : arodri7 1.2 $fam_func =~ s/.*: //;
464 :     chomp $fam_func;
465 :     my $prot_func;
466 :     #my $fig = $self->{fig};
467 :    
468 :     if ($prot =~ /^fig\|/)
469 :     {
470 :     $prot_func = $fig->function_of($prot);
471 :     }
472 :     else
473 :     {
474 :     $prot_func = $PFsO->function_of($prot);
475 :     }
476 :    
477 : gdpusch 1.1 if ($prot_func)
478 :     {
479 :     $prot_func =~ s/\s*\#.*$//;
480 :     return ($prot_func eq $fam_func);
481 :     }
482 :     return 0;
483 :     }
484 :    
485 :     sub ross_hack {
486 :     my($self,$debug,$loose,$seq,$dir,$debug_prot,$nuc,$boundsFile) = @_;
487 :    
488 :     my $all_bounds;
489 :     if ($boundsFile)
490 :     {
491 :     if ($boundsFile eq "bounds") { $all_bounds = $self->{bounds} }
492 :     elsif (-s "$boundsFile") { $all_bounds = &load_bounds($boundsFile) }
493 :     elsif (-s "$dir/$boundsFile") { $all_bounds = &load_bounds("$dir/$boundsFile") }
494 :     else { $all_bounds = $self->{bounds} }
495 :     }
496 :     else
497 :     {
498 :     $all_bounds = $self->{bounds};
499 :     }
500 :    
501 :     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";
502 :     open(TMP,">$tmpF")
503 :     || die "could not open tmp$$.fasta";
504 :     print TMP ">query\n$seq\n";
505 :     close(TMP);
506 :    
507 :     my $query_ln = length($seq);
508 :     my @tmp;
509 :     if ($nuc){
510 :     @tmp = `blastall -i $tmpF -d $dir/PROTS.fasta -m 8 -FF -p blastx -g T`;
511 :     }
512 :     else{
513 :     @tmp = `blastall -i $tmpF -d $dir/PROTS.fasta -m 8 -FF -p blastp`;
514 :     }
515 :    
516 :     unlink($tmpF);
517 :    
518 :     my %seen;
519 :     my $should_be = 0;
520 :     my $yes = 0;
521 :     my $no = 0;
522 :    
523 :     my $ln1 = length($seq);
524 :     my $good = 0;
525 :     my $bad = 0;
526 :    
527 :     my $sims = [];
528 :     my $required_better = ((@{$self->{protsL}} - ($debug_prot ? 1 : 0)) > 1) ? 1 : 0;
529 :     # print STDERR "required_better=$required_better\n";
530 :    
531 :     for (my $simI=0; ($simI < @tmp) && (! $good) && (! $bad); $simI++)
532 :     {
533 :     $_ = $tmp[$simI];
534 :     if ($debug) { print STDERR $_ }
535 :     chop;
536 :    
537 :     my $sim = [split(/\t/,$_)];
538 :     my $prot = $sim->[1];
539 :     next if ($debug_prot && ($debug_prot eq $prot));
540 :    
541 :     my $sc = $sim->[10];
542 :     my $bit_score = $sim->[11];
543 :    
544 :     my $matched1 = abs($sim->[7] - $sim->[6]) + 1;
545 :     my $matched2 = abs($sim->[9] - $sim->[8]) + 1;
546 :     my $ln2 = $self->{prot_lengths}->{$prot};
547 :    
548 :     my $bounds;
549 :     if ($nuc) {
550 :     if ((! $seen{$prot}) &&
551 :     ($loose ||
552 :     (
553 :     ($matched1 >= (0.7 * $ln1)) )
554 :     )
555 :     )
556 :     {
557 :     $seen{$prot} = 1;
558 :     push @$sim, $query_ln, $self->{prot_lengths}->{$prot};
559 :     bless $sim, 'Sim';
560 :     push @$sims, $sim;
561 :    
562 :     $bounds = $all_bounds->{$prot};
563 :     if ($bounds && (@$sims <= 10))
564 :     {
565 :     if ((($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score)) ||
566 :     ($loose &&
567 :     ((($bit_score/$ln1) >= ($bounds->[1] / $ln2)) &&
568 :     ((! $bounds->[2]) || (($bounds->[3]/$ln2) < ($bit_score/$ln1))))))
569 :     {
570 :     if ($debug) { print STDERR " yes - $prot matched1=$matched1 ln1=$ln1 matched2=$matched2 ln2=$ln2\n" }
571 :     ++$yes;
572 :     if ($yes > ($no+$required_better)) { $good = 1 }
573 :     }
574 :     else
575 :     {
576 :     if ($debug) { print STDERR " no - $prot ", join(",",@$bounds),"\n" }
577 :     ++$no;
578 :     if ($no > ($yes+$required_better)) { $bad = 1 }
579 :     }
580 :     }
581 :     else {
582 :     if ($debug) {
583 :     print STDERR " bounds=", ($bounds ? qq(non-nil) : qq(nil))
584 :     , ", num_sims=", (scalar @$sims), "\n";
585 :     }
586 :     }
587 :     }
588 :     else {
589 :     if ($debug) {
590 :     print STDERR
591 :     " seen=", ($seen{$prot} ? $seen{$prot} : 0), " score=$sc,"
592 :     , " matched1=$matched1, ln1=$ln1,"
593 :     , " matched2=$matched2, ln2=$ln2,"
594 :     , "\n";
595 :     }
596 :     }
597 :     }
598 :     else{
599 :     if ((! $seen{$prot}) &&
600 :     ($loose ||
601 :     (($sc <= 1.0e-10) &&
602 :     ($matched1 >= (0.7 * $ln1)) &&
603 :     ($matched2 >= (0.7 * $ln2)))
604 :     )
605 :     )
606 :     {
607 :     $seen{$prot} = 1;
608 :     push @$sim, $query_ln, $self->{prot_lengths}->{$prot};
609 :     bless $sim, 'Sim';
610 :     push @$sims, $sim;
611 :    
612 :     $bounds = $all_bounds->{$prot};
613 :     if ($bounds && (@$sims <= 10))
614 :     {
615 :     if ((($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score)) ||
616 :     ($loose &&
617 :     ((($bit_score/$ln1) >= ($bounds->[1] / $ln2)) &&
618 :     ((! $bounds->[2]) || (($bounds->[3]/$ln2) < ($bit_score/$ln1))))))
619 :     {
620 :     if ($debug) { print STDERR " yes - $prot matched1=$matched1 ln1=$ln1 matched2=$matched2 ln2=$ln2\n" }
621 :     ++$yes;
622 :     if ($yes > ($no+$required_better)) { $good = 1 }
623 :     }
624 :     else
625 :     {
626 :     if ($debug) { print STDERR " no - $prot ", join(",",@$bounds),"\n" }
627 :     ++$no;
628 :     if ($no > ($yes+$required_better)) { $bad = 1 }
629 :     }
630 :     }
631 :     else {
632 :     if ($debug) {
633 :     print STDERR " bounds=", ($bounds ? qq(non-nil) : qq(nil))
634 :     , ", num_sims=", (scalar @$sims), "\n";
635 :     }
636 :     }
637 :     }
638 :     else {
639 :     if ($debug) {
640 :     print STDERR
641 :     " seen=", ($seen{$prot} ? $seen{$prot} : 0), " score=$sc,"
642 :     , " matched1=$matched1, ln1=$ln1,"
643 :     , " matched2=$matched2, ln2=$ln2,"
644 :     , "\n";
645 :     }
646 :     }
647 :     }
648 :     }
649 :     if ($yes > $no) { $good = 1 }
650 :     if ($debug) { print STDERR " yes=$yes no=$no good=$good\n"; }
651 :     return ($good,$sims);
652 :     }
653 :    
654 :     =head3 family_function
655 :    
656 :     usage:
657 :     $func = $protfam_obj->family_function();
658 :    
659 :     Returns the "consensus function" assigned to a FIGfam object.
660 :    
661 :     =cut
662 :    
663 :     sub family_function {
664 :     my($self,$full) = @_;
665 :    
666 :     my $fam_id = $self->{id};
667 :     my $func = $self->{function};
668 :     if (! $full) { $func =~ s/^$fam_id \(not subsystem-based\): // }
669 :     return $func;
670 :     }
671 :    
672 :    
673 :    
674 :     =head3 family_id
675 :    
676 :     usage:
677 :     $fam_id = $figfam_obj->family_id();
678 :    
679 :     Returns the FIGfam ID of a FIGfam object.
680 :    
681 :     =cut
682 :    
683 :     sub family_id {
684 :     my($self) = @_;
685 :    
686 :     return $self->{id};
687 :     }
688 :    
689 :    
690 :    
691 :     =head3 family_dir
692 :    
693 :     usage:
694 :    
695 :     $dir = &ProtFamLite::family_dir( $ProtFams_Obj, $fam_id );
696 :    
697 :     Returns the path to the subdirectory of C<$ProtFams_Obj>
698 :     that the protein Family data for a Family with ID C<$fam_id>
699 :     would be stored in.
700 :    
701 :     =cut
702 :    
703 :     sub fam_dir {
704 :     my ($ProtFams_Obj, $fam_id) = @_;
705 :     my $protfam_root = $ProtFams_Obj->{root};
706 :    
707 :     my $group = substr($fam_id, -3) ||
708 :     confess "Could not extract group-number from Family-ID $fam_id";
709 :     $group = (qq(0) x (3-length($group))) . $group;
710 :    
711 :     return qq($protfam_root/FAMS/$group/$fam_id);
712 :     }
713 :    
714 :     sub by_fig_id {
715 :     my($a,$b) = @_;
716 :     my($g1,$g2,$t1,$t2,$n1,$n2);
717 :     if (($a =~ /^fig\|(\d+\.\d+).([^\.]+)\.(\d+)$/) && (($g1,$t1,$n1) = ($1,$2,$3)) &&
718 :     ($b =~ /^fig\|(\d+\.\d+).([^\.]+)\.(\d+)$/) && (($g2,$t2,$n2) = ($1,$2,$3))) {
719 :     ($g1 <=> $g2) or ($t1 cmp $t2) or ($n1 <=> $n2);
720 :     } else {
721 :     $a cmp $b;
722 :     }
723 :     }
724 :    
725 :     sub read_fasta_record {
726 :    
727 :     shift if UNIVERSAL::isa($_[0],__PACKAGE__);
728 :     my ($file_handle) = @_;
729 :     my ($old_end_of_record, $fasta_record, @lines, $head, $sequence, $seq_id, $comment, @parsed_fasta_record);
730 :    
731 :     if (not defined($file_handle)) { $file_handle = \*STDIN; }
732 :    
733 :     $old_end_of_record = $/;
734 :     $/ = "\n>";
735 :    
736 :     if (defined($fasta_record = <$file_handle>)) {
737 :     chomp $fasta_record;
738 :     @lines = split( /\n/, $fasta_record );
739 :     $head = shift @lines;
740 :     $head =~ s/^>?//;
741 :     $head =~ m/^(\S+)/;
742 :     $seq_id = $1;
743 :     if ($head =~ m/^\S+\s+(.*)$/) { $comment = $1; } else { $comment = ""; }
744 :     $sequence = join( "", @lines );
745 :     @parsed_fasta_record = ( $seq_id, \$sequence, $comment );
746 :     } else {
747 :     @parsed_fasta_record = ();
748 :     }
749 :    
750 :     $/ = $old_end_of_record;
751 :    
752 :     return @parsed_fasta_record;
753 :     }
754 :    
755 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3