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

Diff of /FigKernelPackages/FigFam.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.64, Wed Oct 24 14:25:28 2007 UTC revision 1.65, Thu Nov 15 14:30:50 2007 UTC
# Line 24  Line 24 
24  use gjoseqlib;  use gjoseqlib;
25  use gjoalignment;  use gjoalignment;
26  use PHOB;  use PHOB;
27    use FF;
28    
29  use Tracer;  use Tracer;
30    
31  use Carp;  use Carp;
# Line 70  Line 72 
72      my($class,$fig,$fam_id,$fam_data) = @_;      my($class,$fig,$fam_id,$fam_data) = @_;
73    
74      ($fam_id =~ /^FIG\d{6}$/) || confess "invalid family id: $fam_id";      ($fam_id =~ /^FIG\d{6}$/) || confess "invalid family id: $fam_id";
75      if ($ENV{'VERBOSE'})  { print STDERR "building FigFam for $fam_id\n"; }      Trace "building FigFam for $fam_id\n" if T(3);
76    
77      my $fam = {};      my $fam = {};
78      $fam->{id}  = $fam_id;      $fam->{id}  = $fam_id;
79      $fam->{fig} = $fig;      $fam->{fig} = $fig;
# Line 107  Line 110 
110              return undef;              return undef;
111          }          }
112      }      }
113      if ($ENV{'VERBOSE'})  { print STDERR "PEGs built\n"; }      Trace "PEGs built\n" if T(3);
114    
115      if ((! -s "$dir/function") && (! -s "$dir/built"))      if ((! -s "$dir/function") && (! -s "$dir/built"))
116      {      {
# Line 120  Line 123 
123      chomp $func;      chomp $func;
124      close(FUNC);      close(FUNC);
125      $fam->{function} = $func;      $fam->{function} = $func;
126      if ($ENV{'VERBOSE'})  { print STDERR "function=$func\n"; }      Trace "function=$func\n" if T(3);
127    
128      my($peg,$pegs);      my($peg,$pegs);
129      my $pegsL = [      my $pegsL = [
130                   sort { &FIG::by_fig_id($a,$b) }                   sort { &FIG::by_fig_id($a,$b) }
131  #                grep { scalar $fig->function_of($_) eq $func }  #                grep { scalar $fig->function_of($_) eq $func }
132                   grep { $fig->is_real_feature($_) }                   grep { $fig->is_real_feature($_) }
133                   map { $_ =~ /^\S+\t(\S+)/; $1 }                   map { $_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/; $1 }
134                   `cat \"$dir/PEGs\"`                   `cat \"$dir/PEGs\"`
135                  ];                  ];
136      if ($ENV{'VERBOSE'}) { print STDERR &Dumper($pegsL) }  
137        Trace &Dumper($pegsL) if T(3);
138    
139      if (@$pegsL < 2)      if (@$pegsL < 2)
140      {      {
141          if (-w $dir) {          if (-w $dir) {
# Line 237  Line 242 
242              my @sims = sort { $b->[2] <=> $a->[2] } map { [split(/\t/,$_)] } split(/\n/,$sims);              my @sims = sort { $b->[2] <=> $a->[2] } map { [split(/\t/,$_)] } split(/\n/,$sims);
243              if ($peg = $sims[0]->[0])              if ($peg = $sims[0]->[0])
244              {              {
245  #               print STDERR "processing $peg\n";                  print STDERR "processing $peg\n";
246                  my($i,$hits,$bad,$bad_sc,$last_good,$func2);                  my($i,$hits,$bad,$bad_sc,$last_good,$func2);
247                  for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)                  for ($i=0, $hits=0, $bad="", $bad_sc=0, $last_good=1; ($i < @sims) && (! $bad); $i++)
248                  {                  {
# Line 265  Line 270 
270                  }                  }
271                  else                  else
272                  {                  {
273  #                   print STDERR "    failed\n";                      print STDERR "    failed\n";
274                  }                  }
275              }              }
276          }          }
# Line 317  Line 322 
322      my $rep_file   = "$dir/reps";      my $rep_file   = "$dir/reps";
323      if ((! -f $blast_file) && (! -s "$dir/built"))      if ((! -f $blast_file) && (! -s "$dir/built"))
324      {      {
325          &FIG::run("split_and_trim_sequences $dir/Split blastout=$blast_file < $dir/PEGs.fasta");          &FIG::run("split_and_trim_sequences $dir/Split 1.0e-10 0.7 blastout=$blast_file < $dir/PEGs.fasta");
326            if (! -s $blast_file)
327            {
328                my $fullF = "$dir/PEGs.fasta";
329                open(BLASTOUT,"blastall -i $fullF -d $fullF -FF -m 8 -p blastp -b 20000 -v 20000 |")
330                    || die "could not run blastall -i $fullF -d $fullF -FF -m 8 -p blastp -b 20000 -v 20000";
331                my @all = ();
332                while (defined($_ = <BLASTOUT>))
333                {
334                    if (($_ =~ /^(\S+)\t(\S+)\t(\S+\t){4}(\S+)\t(\S+)\t(\S+)\t(\S+)\s+(\S+)/) &&
335                        ($1 ne $2) &&
336                        ($8 <= 1.0e-10) &&
337                        ((@all == 0) || (($1 ne $all[$#all]->[0]) || ($2 ne $all[$#all]->[1]))))
338                    {
339                        push(@all,[$1,$2,$4,$5,$6,$7,$8,$fam->{peg_lengths}->{$1},$fam->{peg_lengths}->{$2}]);
340                    }
341                }
342                close(BLASTOUT);
343    
344                open(BLASTOUT,">$blast_file") || die "could not open $blast_file";
345                foreach $_ (@all)
346                {
347                    print BLASTOUT join("\t",@$_),"\n";
348                }
349                close(BLASTOUT);
350            }
351    
352          if (-s $blast_file)          if (-s $blast_file)
353          {          {
354              my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;              my @out = map { chop; [split(/\t/,$_) ] } `cat $blast_file`;
# Line 564  Line 595 
595    
596  sub pegs_of {  sub pegs_of {
597      my($self) = @_;      my($self) = @_;
598      return [$self->list_members];  
599        return &FF::pegs_of($self);
600  }  }
601    
602    
# Line 767  Line 799 
799  sub representatives {  sub representatives {
800      my($self) = @_;      my($self) = @_;
801    
802      my $reps = $self->{reps};      return &FF::representatives($self);
     return @$reps;  
803  }  }
804    
805    
# Line 794  Line 825 
825  sub should_be_member {  sub should_be_member {
826      my($self,$seq,$debug,$loose,$debug_peg) = @_;      my($self,$seq,$debug,$loose,$debug_peg) = @_;
827    
828      if ($debug) {      return &FF::should_be_member($self,$seq,$debug,$loose,$debug_peg);
         $ENV{DEBUG}   = 1;  
         $ENV{VERBOSE} = $ENV{VERBOSE} ? $ENV{VERBOSE} : 1;  
         print STDERR "\nTesting for membership:\n";  
     }  
   
     my $fig = $self->{fig};  
     my $dir = $self->{dir};  
   
     if ($ENV{DEBUG}) { print STDERR &Dumper([$self->family_id,$seq]) }  
   
   
     my($in);  
     if (open(DEC,"<$dir/decision.procedure") && ($in = <DEC>) && ($in =~ /^(\S+)(\s+(\S.*\S))?/))  
     {  
         no strict 'refs';  
   
         my $procedure = $1;  
         my @args      = $3 ? split(/\s+/,$3) : ();  
         return &{$procedure}($self,$fig,$loose,$seq,$dir,$debug_peg,@args);  
     }  
     else  
     {  
         return &ross_hack($self,$fig,$loose,$seq,$dir,$debug_peg);  
     }  
 }  
   
 sub blast_vote {  
     my($self,$fig,$loose,$seq,$dir,$debug_peg,$partition,$min_bsc) = @_;  
   
     (-s "$dir/PEGs") || return undef;  
     if ($ENV{DEBUG}) { print STDERR "checking: min_bsc=$min_bsc\n" }  
     my %pegs_in = map { $_ =~ /(fig\|\d+\.\d+\.peg\.\d+)/; $1 => 1 } `cat $dir/PEGs`;  
   
     my $N = &FIG::min(10,scalar keys(%pegs_in));  
     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";  
     open(TMP,">$tmpF")  
         || die "could not open $FIG_Config::temp/tmp$$.fasta";  
     print TMP ">query\n$seq\n";  
     close(TMP);  
     my $query_ln = length($seq);  
     my $partitionF = "$FIG_Config::data/FigfamsData/Partitions/" . ($partition % 1000) . "/$partition/fasta";  
     my @tmp = `blastall -i $tmpF -d $partitionF -m 8 -FF -p blastp`;  
   
     my $sims = [];  
     my $in = 0;  
     my $out = 0;  
     my(%seen);  
     for (my $simI=0; $N && ($simI < @tmp); $simI++)  
     {  
         $_ = $tmp[$simI];  
         if ($ENV{DEBUG}) { print STDERR $_ }  
         chop;  
   
         my $sim = [split(/\t/,$_)];  
         my $peg = $sim->[1];  
         next if ($debug_peg && ($debug_peg eq $peg));  
         my $bit_score = $sim->[11];  
         my $matched1 = abs($sim->[7] - $sim->[6]) + 1;  
         my $matched2 = abs($sim->[9] - $sim->[8]) + 1;  
         if ($ENV{DEBUG}) { print STDERR sprintf("%3.2f",$bit_score / $matched2),"\n" }  
         if (! $seen{$peg})  
         {  
             $seen{$peg} = 1;  
             my $ln2;  
             if ($pegs_in{$peg} &&                        # (print "in-ok\n") &&  
                 ($ln2 = $self->{peg_lengths}->{$peg}) && # (print "ln2-ok\n") &&  
                 ($matched1 >= (0.7 * $query_ln)) &&      # (print "mat1-ok\n") &&  
                 ($matched2 >= (0.7 * $ln2)) &&           # (print "mat2-ok\n") &&  
                 (($bit_score / $matched2) >= $min_bsc)        # && (print "bsc-ok\n")  
                 )  
             {  
                 push @$sim, $query_ln, $self->{peg_lengths}->{$peg};  
                 bless $sim, 'Sim';  
                 push @$sims, $sim;  
                 $in++;  
             }  
             else  
             {  
                 $out++;  
             }  
             if ($ENV{DEBUG}) {print STDERR &Dumper(["in=$in",  
                                                     "out=$out",  
                                                     $pegs_in{$peg},  
                                                     "ln1=$query_ln",  
                                                     "ln2=$ln2",  
                                                     "matched1=$matched1",  
                                                     "matched2=$matched2",  
                                                     "bsc=$bit_score",  
                                                     $ln2 ? sprintf("%3.2f",$bit_score/$matched2) : ""]); }  
             $N--;  
         }  
     }  
     return (($in > $out),$sims);  
 }  
   
 sub ross_hack {  
     my($self,$fig,$loose,$seq,$dir,$debug_peg) = @_;  
   
     my $tmpF = "$FIG_Config::temp/tmp$$.fasta";  
     open(TMP,">$tmpF")  
         || die "could not open $FIG_Config::temp/tmp$$.fasta";  
     print TMP ">query\n$seq\n";  
     close(TMP);  
   
     my $query_ln = length($seq);  
     my @tmp = `blastall -i $tmpF -d $dir/PEGs.fasta -m 8 -FF -p blastp`;  
   
     my %seen;  
     my $should_be = 0;  
     my $yes = 0;  
     my $no  = 0;  
   
     my $ln1 = length($seq);  
     my $good = 0;  
     my $bad = 0;  
   
     my $sims = [];  
     for (my $simI=0; ($simI < @tmp) && (! $good) && (! $bad); $simI++)  
     {  
         $_ = $tmp[$simI];  
         if ($ENV{DEBUG}) { print STDERR $_ }  
         chop;  
   
         my $sim = [split(/\t/,$_)];  
         my $peg = $sim->[1];  
         next if ($debug_peg && ($debug_peg eq $peg));  
   
         my $sc = $sim->[10];  
         my $bit_score = $sim->[11];  
   
         my $matched1 = abs($sim->[7] - $sim->[6]) + 1;  
         my $matched2 = abs($sim->[9] - $sim->[8]) + 1;  
         my $ln2 = $self->{peg_lengths}->{$peg};  
   
         if ((! $seen{$peg}) &&  
             ($loose ||  
              (($sc <= 1.0e-10) &&  
               ($matched1 >= (0.8 * $ln1)) &&  
               ($matched2 >= (0.8 * $ln2)))  
             )  
            )  
         {  
             $seen{$peg} = 1;  
             push @$sim, $query_ln, $self->{peg_lengths}->{$peg};  
             bless $sim, 'Sim';  
             push @$sims, $sim;  
   
             my $bounds = $self->{bounds}->{$peg};  
             if ($bounds && (@$sims <= 10))  
             {  
                 if ((($bit_score >= $bounds->[1]) && ((! $bounds->[2]) || $bounds->[3] < $bit_score)) ||  
                     ($loose &&  
                      ((($bit_score/$ln1) >= ($bounds->[1] / $ln2)) &&  
                       ((! $bounds->[2]) || (($bounds->[3]/$ln2) < ($bit_score/$ln1))))))  
                 {  
                     if ($ENV{DEBUG}) { print STDERR "    yes - $peg matched1=$matched1 ln1=$ln1 matched2=$matched2 ln2=$ln2\n" }  
                     ++$yes;  
                     if ($yes > ($no+1)) { $good = 1 }  
                 }  
                 else  
                 {  
                     my $bad_func = $bounds->[2] ? $fig->function_of($bounds->[2]) : "";  
                     if ($ENV{DEBUG}) { print STDERR "    no - $peg ", join(",",@$bounds)," $bad_func\n" }  
                     ++$no;  
                     if ($no > $yes) { $bad = 1 }  
                 }  
             }  
             else {  
                 if ($ENV{DEBUG}) {  
                     print STDERR "    bounds=", ($bounds ? qq(non-nil) : qq(nil))  
                         , ", num_sims=", (scalar @$sims), "\n";  
                 }  
             }  
         }  
         else {  
             if ($ENV{DEBUG}) {  
                 print STDERR  
                     "    seen=", ($seen{$peg} ? $seen{$peg} : 0), " score=$sc,"  
                     , " matched1=$matched1, ln1=$ln1,"  
                     , " matched2=$matched2, ln2=$ln2,"  
                     , "\n";  
             }  
         }  
     }  
     if ($ENV{DEBUG}) { print STDERR "        yes=$yes no=$no good=$good\n"; }  
     return ($good,$sims);  
829  }  }
830    
831  =head3 family_function  =head3 family_function
# Line 995  Line 840 
840  sub family_function {  sub family_function {
841      my($self) = @_;      my($self) = @_;
842    
843      return $self->{function}      return &FF::family_function($self);
844  }  }
845    
846    
# Line 1012  Line 857 
857  sub family_id {  sub family_id {
858      my($self) = @_;      my($self) = @_;
859    
860      return $self->{id};      return &FF::family_function($self);
861  }  }
862    
   
   
863  =head3 family_dir  =head3 family_dir
864    
865  usage:  usage:

Legend:
Removed from v.1.64  
changed lines
  Added in v.1.65

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3