[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.112, Tue Jun 8 17:19:43 2004 UTC revision 1.113, Tue Jun 8 20:26:23 2004 UTC
# Line 1508  Line 1508 
1508    
1509  =cut  =cut
1510    
1511  sub genome_version {  sub genome_version :scalar {
1512      my($self,$genome) = @_;      my($self,$genome) = @_;
1513    
1514      my(@tmp);      my(@tmp);
# Line 1648  Line 1648 
1648  =cut  =cut
1649    
1650    
1651  sub genome_of {  sub genome_of :scalar {
1652      shift if UNIVERSAL::isa($_[0],__PACKAGE__);      shift if UNIVERSAL::isa($_[0],__PACKAGE__);
1653      my $prot_id = (@_ == 1) ? $_[0] : $_[1];      my $prot_id = (@_ == 1) ? $_[0] : $_[1];
1654    
# Line 4964  Line 4964 
4964    
4965  =cut  =cut
4966    
4967  sub crude_estimate_of_distance {  sub crude_estimate_of_distance :scalar {
4968      my($self,$genome1,$genome2) = @_;      my($self,$genome1,$genome2) = @_;
4969      my($i,$v,$d,$dist);      my($i,$v,$d,$dist);
4970    
# Line 4981  Line 4981 
4981      return $self->crude_estimate_of_distance1($genome1,$genome2);      return $self->crude_estimate_of_distance1($genome1,$genome2);
4982  }  }
4983    
4984  sub crude_estimate_of_distance1 {  sub crude_estimate_of_distance1 :scalar {
4985      my($self,$genome1,$genome2) = @_;      my($self,$genome1,$genome2) = @_;
4986      my($i,$v,$d,$dist);      my($i,$v,$d,$dist);
4987    
# Line 5841  Line 5841 
5841      }      }
5842  }  }
5843    
5844    
5845    ###########
5846    #
5847    # Some routines for dealing with peg search and similarities.
5848    #
5849    # This is code lifted from pom.cgi and reformatted for more general use.
5850    #
5851    #
5852    # Find the given role in the given (via CGI params) organism.
5853    #
5854    # We do this by finding a list of pegs that are annotated to have
5855    # this role in other organisms that are "close enough" to our organism
5856    #
5857    # We then find pegs in this organism that are similar to
5858    # these pegs.
5859    #
5860    
5861    sub find_role_in_org
5862    {
5863        my($self,$role, $org, $user, $sims_cutoff) = @_;
5864    
5865        my($id2,$psc,$col_hdrs,$tab,$peg,$curr_func,$id2_func);
5866        my($seen,$peg);
5867    
5868        if (!$org)
5869        {
5870            return undef;
5871        }
5872    
5873        #
5874        # Create a list of candidates.
5875        #
5876        # These are the list of sequences that contain the given role,
5877        # sorted by the crude_estimate_of_distance from the given peg.
5878        #
5879    
5880        my @cand = map { $_->[0] }
5881        sort { $a->[1] <=> $b->[1] }
5882        map {
5883            $peg = $_;
5884            [$peg,$self->crude_estimate_of_distance($org,&FIG::genome_of($peg))]
5885            }
5886        $self->seqs_with_role($role,$user);
5887    
5888        my $hits = {};
5889        $seen = {};
5890    
5891        #
5892        # Pick the top 10 hits if there are more than 10.
5893        #
5894        my $how_many = (@cand > 10) ? 10 : scalar @cand;
5895    
5896        $self->try_to_locate($org,$hits,[@cand[1..$how_many]],$seen, $sims_cutoff);
5897    
5898        if (keys(%$hits) == 0)
5899        {
5900            splice(@cand,0,$how_many+1);
5901            &try_to_locate($self,$org,$hits,\@cand,$seen, $sims_cutoff);
5902        }
5903    
5904        #
5905        # At this point %$hits contains the pegs in our organism that
5906        # may have the given role. The key is the peg, the value
5907        # is a pair [score, similar-peg]
5908        #
5909        #
5910        # We reformat this into a list of entries
5911        # [ $psc, $peg-in-this-org, $length, $current-fn, $matched-protein, $matched-len, $matched-fun]
5912        #
5913    
5914    
5915        $col_hdrs = ["P-Sc","PEG","Ln1","Current Function", "Protein Hit","Ln2","Function"];
5916    
5917        my @ret;
5918    
5919        foreach $peg ( sort {$hits->{$a}->[0] <=> $hits->{$b}->[0]} keys(%$hits))
5920        {
5921            ($psc,$id2) = @{$hits->{$peg}};
5922            $curr_func = $self->function_of($peg,$user);
5923            $id2_func  = $self->function_of($id2,$user);
5924    
5925            push(@ret, [$psc, $peg, $self->translation_length($peg),
5926                        $curr_func, $id2, $self->translation_length($id2),$id2_func]);
5927        }
5928        return @ret;
5929    }
5930    
5931    #
5932    # Mark in $hits the pegs in $org that are similar to
5933    # pegs in other organisms that have the given role.
5934    #
5935    sub try_to_locate {
5936        my($self,$org,$hits,$to_try,$seen, $cutoff) = @_;
5937        my($prot,$id2,$psc,$id2a,$x,$sim);
5938    
5939        if (! $cutoff) { $cutoff = 1.0e-5 }
5940    
5941        #
5942        # @$to_try is a list of pegs
5943        #
5944        foreach $prot (@$to_try)
5945        {
5946            #
5947            # If we've not looked at it before ...
5948            #
5949            if (! $seen->{$prot})
5950            {
5951                #
5952                # Retrieve the top 1000 sims for this peg.  raw
5953                # means don't expand.
5954                #
5955                foreach $sim ($self->sims($prot,1000,$cutoff,"raw",0))
5956                {
5957                    $id2 = $sim->id2;
5958                    $psc = $sim->psc;
5959    
5960                    #
5961                    # Retrieve the proteins that the sims map to.
5962                    #
5963    
5964                    foreach $id2a (map { $_->[0] } $self->mapped_prot_ids($id2))
5965                    {
5966                        #
5967                        # If it's a protein in the organism we're looking in,
5968                        # and if it's a better hit than the hit we had before,
5969                        # mark it in $hits->{id} with the score and the
5970                        # protein id.
5971                        #
5972                        if (($id2a =~ /^fig\|(\d+\.\d+)/) && ($1 eq $org))
5973                        {
5974                            $x = $hits->{$id2a};
5975                            if ((! $x) || ($x->[0] > $psc))
5976                            {
5977                                $hits->{$id2a} = [$sim->psc,$prot];
5978                            }
5979                        }
5980                        #
5981                        # Otherwise, mark it as having been seen if the score is good enough.
5982                        #
5983                        elsif ($psc < 1.0e-20)
5984                        {
5985                            {
5986                                $seen->{$id2a} = 1;
5987                            }
5988                        }
5989    
5990                    }
5991                }
5992            }
5993        }
5994    }
5995    
5996    
5997    
5998    
5999  1;  1;
6000    
6001    

Legend:
Removed from v.1.112  
changed lines
  Added in v.1.113

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3