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

Annotation of /FigKernelPackages/ProtFamLite.pm

Parent Directory Parent Directory | Revision Log Revision Log


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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3