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

Diff of /FigKernelPackages/FIG.pm

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

revision 1.237, Tue Mar 8 21:14:53 2005 UTC revision 1.238, Fri Mar 11 01:04:32 2005 UTC
# Line 3987  Line 3987 
3987  =cut  =cut
3988    
3989  sub sims {  sub sims {
3990      my ($self,$id,$maxN,$maxP,$select,$max_expand) = @_;      my ( $self, $id, $maxN, $maxP, $select, $max_expand, $filters ) = @_;
3991      my($sim);      my($sim);
3992    
3993      $max_expand = defined($max_expand) ? $max_expand : $maxN;      $max_expand = defined($max_expand) ? $max_expand : $maxN;
3994    
3995      if ($self->is_deleted_fid($id)) { return () }      return () if $self->is_deleted_fid( $id );
3996    
     my @sims = ();  
3997      my @maps_to = $self->mapped_prot_ids($id);      my @maps_to = $self->mapped_prot_ids($id);
3998      if (@maps_to > 0)      ( @maps_to > 0 ) or return ();
3999      {  
4000          my $rep_id   = $maps_to[0]->[0];          my $rep_id   = $maps_to[0]->[0];
4001          my @entry = grep { $_->[0] eq $id } @maps_to;      if ( ! defined( $maps_to[0]->[1] ) )
         if ((@entry == 1) && defined($entry[0]->[1]))  
         {  
             if ((! defined($maps_to[0]->[1])) ||  
                 (! defined($entry[0]->[1])))  
4002              {              {
4003                  print STDERR &Dumper(\@maps_to,\@entry);          print STDERR &Dumper( \@maps_to );
4004                  confess "bad";                  confess "bad";
4005              }              }
4006    
4007        my @entry = grep { $_->[0] eq $id } @maps_to;
4008        ( @entry == 1 ) and defined( $entry[0]->[1] ) or return ();
4009    
4010        #  Get the similarities
4011    
4012        my @raw_sims = get_raw_sims( $self, $rep_id, $maxP, $filters );
4013    
4014        #  If the query is not the representative, make sims look like it is
4015        #  by replacing id1 and fixing match coordinates if lengths differ.
4016    
4017              my $delta = $maps_to[0]->[1] - $entry[0]->[1];              my $delta = $maps_to[0]->[1] - $entry[0]->[1];
             my @raw_sims = &get_raw_sims($self,$rep_id,$maxN,$maxP);  
4018              if ($id ne $rep_id)              if ($id ne $rep_id)
4019              {              {
4020                  foreach $sim (@raw_sims)                  foreach $sim (@raw_sims)
4021                  {                  {
   
4022                      $sim->[0] = $id;                      $sim->[0] = $id;
4023                      $sim->[6] -= $delta;                      $sim->[6] -= $delta;
4024                      $sim->[7] -= $delta;                      $sim->[7] -= $delta;
4025                  }                  }
4026              }              }
4027    
4028        #  The query must be present for expanding matches to identical sequences.
4029    
4030              if (($max_expand > 0) && ($select ne "raw"))              if (($max_expand > 0) && ($select ne "raw"))
4031              {              {
4032                  unshift(@raw_sims,bless([$id,                  unshift(@raw_sims,bless([$id,
4033                                           $rep_id,                                           $rep_id,
4034                                           100.00,                                       "100.00",
4035                                           $entry[0]->[1],                                           $entry[0]->[1],
4036                                           0,                                           0,
4037                                           0,                                           0,
# Line 4033  Line 4041 
4041                                           2 * $entry[0]->[1],                                           2 * $entry[0]->[1],
4042                                           $entry[0]->[1],                                           $entry[0]->[1],
4043                                           $maps_to[0]->[1],                                           $maps_to[0]->[1],
4044                                           "blastp",                                       "blastp"
4045                                           undef,                                     ], 'Sim'
4046                                           undef                                   )
4047                                          ],'Sim'));                 );
4048                  $max_expand++;                  $max_expand++;
4049              }              }
4050              @sims = grep { $_->id1 ne $_->id2 } &expand_raw_sims($self,\@raw_sims,$maxP,$select,0,$max_expand);  
4051          }      # print STDERR "\n\n"; for ( @raw_sims ) { print STDERR join( ", ", @{ $_ } ), "\n" }
4052      }  
4053      return grep { ! $self->is_deleted_fid($_->id2) } @sims;      #  expand_raw_sims now handles sanity checks on id1 eq id2 and id2
4054        #  is not deleted.  This lets it keep count of the actual number of
4055        #  sims reported!
4056    
4057        return expand_raw_sims( $self, \@raw_sims, $maxP, $select, 0, $max_expand, $filters );
4058  }  }
4059    
4060    
4061    #  maxP is not used.  It is checked by both functions that call here.
4062    
4063  sub expand_raw_sims {  sub expand_raw_sims {
4064      my($self,$raw_sims,$maxP,$select,$dups,$max_expand) = @_;      my( $self, $raw_sims, $maxP, $select, $dups, $max_expand, $filters ) = @_;
4065      my($sim,$id2,%others,$x);      my( $sim, $id1, $id2, %others, $x );
4066    
4067        #  Set up behavior defaults (pretty wide open):
4068    
4069        my ( $maxN, $show_env );
4070        if ( $filters && ref( $filters ) eq "HASH" )
4071        {
4072            defined( $filters->{ maxN }       ) and $maxN       = $filters->{ maxN };
4073            defined( $filters->{ select }     ) and $select     = $filters->{ select };
4074            defined( $filters->{ max_expand } ) and $max_expand = $filters->{ max_expand };
4075            defined( $filters->{ show_env }   ) and $show_env   = $filters->{ show_env };
4076            defined( $filters->{ dups }       ) and $dups       = $filters->{ dups };
4077        }
4078        defined( $maxN )       or $maxN       = 1000000;   # Unlimited sims
4079        defined( $select )     or $select     =    'all';  # Show all expansions
4080        defined( $max_expand ) or $max_expand =       0;   # But none by default
4081        defined( $show_env )   or $show_env   =       1;   # Show environmental by default
4082    
4083      my @sims = ();      my @sims = ();
4084      foreach $sim (@$raw_sims)      foreach $sim (@$raw_sims)
4085      {      {
         next if ($sim->psc > $maxP);  
4086          $id2 = $sim->id2;          $id2 = $sim->id2;
4087          next if ($others{$id2} && (! $dups));          if ( ! $dups )
4088            {
4089                next if $others{ $id2 };
4090          $others{$id2} = 1;          $others{$id2} = 1;
4091          if (($select && ($select eq "raw")) || ($max_expand <= 0))          }
4092    
4093            $id1 = $sim->id1;
4094            if ( ( $select eq "raw" ) || ( $max_expand <= 0 ) )
4095          {          {
4096                next if ( ! $show_env && ( $id2 =~ /^fig\|9999999/ ) );
4097                next if ( $id1 eq $id2 ) || $self->is_deleted_fid( $id2 );
4098              push(@sims,$sim);              push(@sims,$sim);
4099                last if ( @sims >= $maxN );
4100          }          }
4101          else          else
4102          {          {
# Line 4067  Line 4104 
4104              $max_expand--;              $max_expand--;
4105    
4106              my @maps_to = $self->mapped_prot_ids($id2);              my @maps_to = $self->mapped_prot_ids($id2);
4107              if ((! $select) || ($select eq "fig"))              if ( $select eq "fig" )            # Only fig
4108              {              {
4109                  @relevant = grep { $_->[0] =~ /^fig/ }  @maps_to;                  @relevant = grep { $_->[0] =~ /^fig/ }  @maps_to;
4110              }              }
4111              elsif ($select && ($select =~ /^ext/i))              elsif ( $select =~ /^fig_?pref/i ) # FIG preferred
4112                {
4113                    @relevant = grep { $_->[0] =~ /^fig/ } @maps_to;
4114                    if ( ! @relevant )
4115                    {
4116                        push @sims, $sim;
4117                        last if ( @sims >= $maxN );
4118                        next;
4119                    }
4120                }
4121                elsif ( $select =~ /^ext/i )       # Not fig
4122              {              {
4123                  @relevant = grep { $_->[0] !~ /^fig/ } @maps_to;                  @relevant = grep { $_->[0] !~ /^fig/ } @maps_to;
4124              }              }
4125              else              else                               # All
4126              {              {
4127                  @relevant = @maps_to;                  @relevant = @maps_to;
4128              }              }
4129    
4130              foreach $x (@relevant)              foreach $x (@relevant)
4131              {              {
                 my $sim1 = [@$sim];  
4132                  my($x_id,$x_ln) = @$x;                  my($x_id,$x_ln) = @$x;
4133                  defined($x_ln) || confess "x_ln id2=$id2 x_id=$x_id";                  defined( $x_ln ) || confess "x_ln id2='$id2' x_id='$x_id'";
4134                  defined($maps_to[0]->[1]) || confess "maps_to";                  defined($maps_to[0]->[1]) || confess "maps_to";
4135                  my $delta2 = $maps_to[0]->[1] - $x_ln;                  next if ( ! $show_env && ( $x_id =~ /^fig\|9999999/ ) );
4136                    next if ( $id1 eq $x_id ) || $self->is_deleted_fid( $x_id );
4137    
4138                    my $delta2  = $maps_to[0]->[1] - $x_ln;   # Coordinate shift
4139                    my $sim1    = [ @$sim ];                  # Make a copy
4140                  $sim1->[1] = $x_id;                  $sim1->[1] = $x_id;
4141                  $sim1->[8] -= $delta2;                  $sim1->[8] -= $delta2;
4142                  $sim1->[9] -= $delta2;                  $sim1->[9] -= $delta2;
4143                  bless($sim1,"Sim");                  bless($sim1,"Sim");
4144                  push(@sims,$sim1);                  push(@sims,$sim1);
4145                    last if ( @sims >= $maxN );
4146              }              }
4147          }          }
4148      }      }
4149    
4150      return @sims;      return @sims;
4151  }  }
4152    
4153    
4154  sub get_raw_sims {  sub get_raw_sims {
4155      my($self,$rep_id,$maxN,$maxP) = @_;      my ( $self, $rep_id, $maxP, $filters ) = @_;
4156      my(@sims,$seek,$fileN,$ln,$fh,$file,$readC,@lines,$i,$sim);      my ( $sim_chunk, $seek, $fileN, $ln, $fh, $file, @lines, $sim );
     my($sim_chunk,$psc,$id2);  
4157    
4158      $maxN = $maxN ? $maxN : 500;      #  Set up behavior defaults (pretty wide open):
4159    
4160        my ( $show_env, $min_sim, $sim_meas, $min_q_cov, $min_s_cov, $sort_by );
4161        if ( $filters && ref( $filters ) eq "HASH" )
4162        {
4163            defined( $filters->{ maxP }      ) and $maxP      = $filters->{ maxP };
4164            defined( $filters->{ show_env }  ) and $show_env  = $filters->{ show_env };
4165            defined( $filters->{ min_sim }   ) and $min_sim   = $filters->{ min_sim };
4166            defined( $filters->{ sim_meas }  ) and $sim_meas  = $filters->{ sim_meas };
4167            defined( $filters->{ min_q_cov } ) and $min_q_cov = $filters->{ min_q_cov };
4168            defined( $filters->{ min_s_cov } ) and $min_s_cov = $filters->{ min_s_cov };
4169            defined( $filters->{ sort_by }   ) and $sort_by   = $filters->{ sort_by };
4170        }
4171        defined( $maxP )      or $maxP       =    10;
4172        defined( $show_env )  or $show_env   =     1;
4173        defined( $min_sim )   or $min_sim    =     0;
4174        defined( $sim_meas )  or $sim_meas   =   'id';
4175        defined( $min_q_cov ) or $min_q_cov  =     0;
4176        defined( $min_s_cov ) or $min_s_cov  =     0;
4177        defined( $sort_by )   or $sort_by = 'bits';
4178    
     @sims = ();  
4179      my $rdbH = $self->db_handle;      my $rdbH = $self->db_handle;
4180      my $relational_db_response = $rdbH->SQL("SELECT seek, fileN, len FROM sim_seeks WHERE id = \'$rep_id\' ");      my $relational_db_response = $rdbH->SQL("SELECT seek, fileN, len FROM sim_seeks WHERE id = \'$rep_id\' ");
4181    
4182        #  Gather all of the acceptable "lines" from the sim chunks
4183    
4184      foreach $sim_chunk (@$relational_db_response)      foreach $sim_chunk (@$relational_db_response)
4185      {      {
4186          ($seek,$fileN,$ln) = @$sim_chunk;          ($seek,$fileN,$ln) = @$sim_chunk;
4187          $file = $self->N2file($fileN);          $file = $self->N2file($fileN);
4188          $fh   = $self->openF($file);          $fh   = $self->openF($file);
4189          if (! $fh)          $fh or confess "could not open sims for $file";
4190    
4191            #  Read file, parse lines, sanity check values, and filter E-value
4192            #   0.  The query peg
4193            #   1.  The similar peg
4194            #   2.  The percent id
4195            #   3.  Alignment length
4196            #   4.  Mismatches
4197            #   5.  Gap openings
4198            #   6.  The start of the match in the query peg
4199            #   7.  The end of the match in the query peg
4200            #   8.  The start of the match in the similar peg
4201            #   9.  The end of the match in the similar peg
4202            #  10.  E-value
4203            #  11.  Bit score
4204            #  12.  Length of query peg
4205            #  13.  Length of similar peg
4206            #  14.  Method
4207    
4208            push @lines, grep { ( @$_ >= 15 ) &&
4209                                ( $_->[10] =~ /^[0-9.e-]+$/ ) &&  # E-value
4210                                ( $_->[10] <= $maxP )   &&        # E-value test
4211                                ( $_->[11] =~ /^[0-9.]+$/ ) &&    # bit score
4212                                ( $_->[12] =~ /^\d+$/ ) &&        # query len
4213                                ( $_->[13] =~ /^\d+$/ ) &&        # subj len
4214                                ( $_->[6]  =~ /^\d+$/ ) &&        # q-match start
4215                                ( $_->[7]  =~ /^\d+$/ ) &&        # q-match end
4216                                ( $_->[8]  =~ /^\d+$/ ) &&        # s-match start
4217                                ( $_->[9]  =~ /^\d+$/ ) &&        # s-match end
4218                                ( $_->[2]  =~ /^[0-9.]+$/ )       # percent id
4219                              }
4220                         map  { [ split( /\t/, $_ ), "blastp" ] }
4221                         @{ read_block( $fh, $seek, $ln-1 ) };
4222        }
4223    
4224        #  Similarity filter
4225    
4226        if ( $min_sim > 0 )
4227        {
4228            if    ( $sim_meas eq 'id'  )
4229            {
4230                @lines = grep { $_->[2] >= $min_sim } @lines;
4231            }
4232            elsif ( $sim_meas eq 'bpp' )
4233          {          {
4234              confess "could not open sims for $file";              @lines = grep { $_->[11] >= $min_sim * ( $_->[7] - $_->[6] + 1 ) } @lines;
4235          }          }
         $readC = &read_block($fh,$seek,$ln-1);  
         @lines = grep    {  
                              (@$_ == 15) &&  
                              ($_->[12] =~ /^\d+$/) &&  
                              ($_->[13] =~ /^\d+$/) &&  
                              ($_->[6] =~ /^\d+$/) &&  
                              ($_->[7] =~ /^\d+$/) &&  
                              ($_->[8] =~ /^\d+$/) &&  
                              ($_->[9] =~ /^\d+$/) &&  
                              ($_->[2] =~ /^[0-9.]+$/) &&  
                              ($_->[10] =~ /^[0-9.e-]+$/)  
4236                           }                           }
                          map  { [split(/\t/,$_),"blastp"] }  
                          @$readC;  
4237    
4238          @lines = sort { $b->[11] <=> $a->[11] } @lines;      #  Query coverage filter
4239    
4240          for ($i=0; ($i < @lines); $i++)      if ( $min_q_cov > 0 )
4241          {          {
4242              $psc    = $lines[$i]->[10];          my $thresh = 0.01 * $min_q_cov;
4243              $id2    = $lines[$i]->[1];          @lines = grep { ( abs( $_->[7] - $_->[6] ) + 1 ) >= ( $thresh * $_->[12] ) } @lines;
4244              if ($maxP >= $psc)      }
4245    
4246        #  Subject coverage filter
4247    
4248        if ( $min_s_cov > 0 )
4249              {              {
4250                  $sim = $lines[$i];          my $thresh = 0.01 * $min_s_cov;
4251                  bless($sim,"Sim");          @lines = grep { ( abs( $_->[9] - $_->[8] ) + 1 ) >= ( $thresh * $_->[13] ) } @lines;
4252                  push(@sims,$sim);      }
4253                  if (@sims == $maxN) { return @sims }  
4254        #  Order the surviving raw sims by requested score:
4255    
4256        if    ( $sort_by eq 'id' )                        # Percent identity
4257        {
4258            @lines = sort { $b->[2] <=> $a->[2] } @lines;
4259              }              }
4260        elsif ( $sort_by eq 'bpp' )                       # Bits per position
4261        {
4262            @lines = map  { $_->[0] }
4263                     sort { $b->[1] <=> $a->[1] }
4264                     map  { [ $_, $_->[11] / abs( $_->[7] - $_->[6] ) ] }
4265                     @lines;
4266          }          }
4267        else                                                 # Bit score (bits)
4268        {
4269            @lines = sort { $b->[11] <=> $a->[11] } @lines;
4270      }      }
4271      return @sims;  
4272        #  Bless the raw sims:
4273    
4274        return map { bless( $_, 'Sim' ); $_ } @lines;
4275  }  }
4276    
4277    
4278  sub read_block {  sub read_block {
4279      shift if UNIVERSAL::isa($_[0],__PACKAGE__);      shift if UNIVERSAL::isa($_[0],__PACKAGE__);
4280      my($fh,$seek,$ln) = @_;      my($fh,$seek,$ln) = @_;
# Line 4308  Line 4433 
4433    
4434      @hits = sort { ($a->psc <=> $b->psc) or ($a->iden cmp $b->iden) } grep { $_->id2 ne $id } @hits;      @hits = sort { ($a->psc <=> $b->psc) or ($a->iden cmp $b->iden) } grep { $_->id2 ne $id } @hits;
4435      if ($maxN && ($maxN < @hits)) { $#hits = $maxN - 1 }      if ($maxN && ($maxN < @hits)) { $#hits = $maxN - 1 }
4436      return &expand_raw_sims($self,\@hits,$maxP,$select);      return expand_raw_sims( $self, \@hits, $maxP, $select );
4437  }  }
4438    
4439  sub blastit {  sub blastit {

Legend:
Removed from v.1.237  
changed lines
  Added in v.1.238

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3