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

Diff of /FigKernelPackages/representative_sequences.pm

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

revision 1.6, Mon Mar 19 20:01:56 2007 UTC revision 1.7, Sat May 12 20:12:56 2007 UTC
# Line 627  Line 627 
627    
628      # ---------------------------------------# Default values for options      # ---------------------------------------# Default values for options
629    
630      my $max_sim   = 0.80;                  # Retain 80% identity of less      my $max_sim   = 0.80;                  # Retain 80% identity or less
631      my $logfile   = undef;                 # Log file of sequences processed      my $logfile   = undef;                 # Log file of sequences processed
632      my $max_e_val = 0.01;                  # Blast E-value to decrease output      my $max_e_val = 0.01;                  # Blast E-value to decrease output
633      my $sim_meas  = 'identity_fraction';   # Use sequence identity as measure      my $sim_meas  = 'identity_fraction';   # Use sequence identity as measure
# Line 1493  Line 1493 
1493      return ( $tmp, $save_tmp );      return ( $tmp, $save_tmp );
1494  }  }
1495    
1496    ###############################################################################
1497    
1498    sub n_rep_seqs {
1499        my(%args) = (ref($_[0]) eq 'HASH') ? %{$_[0]} : @_;
1500    
1501        my($seqs)          = $args{seqs}          || return undef;
1502        my($reps)          = $args{reps}          || undef;
1503        my($max_iden)      = $args{max_iden}      || 0.9;      # we don't keep seqs more than 90% identical
1504        my($max_rep)       = $args{max_rep}       || 50;       # maximum number of seqs in returned set
1505    
1506        my($lost) = {};
1507        my($repseqs,$representing) = &rep_seq_2($reps ? $reps : (), $seqs, { max_sim => $max_iden });
1508        if ($max_rep >= @$repseqs)
1509        {
1510            return ($repseqs,$representing);
1511        }
1512        my $n_rep = $reps ? @$reps : 0;
1513        my $incr = $max_iden / 2;
1514    #   print STDERR "max_iden=$max_iden, ", scalar @$repseqs,"\n";
1515        my $iterations_left = 7;
1516    
1517        my @seqs2;
1518        while ($iterations_left && ($max_rep != @$repseqs))
1519        {
1520            if ($max_rep > @$repseqs)
1521            {
1522                $max_iden += $incr;
1523            }
1524            else
1525            {
1526                @seqs2 = @$repseqs[$n_rep..(@$repseqs - 1)];
1527                &add_to_lost($lost,$representing);
1528                $max_iden -= $incr;
1529            }
1530            ($repseqs,$representing) = &rep_seq_2($reps ? $reps : (), \@seqs2, { max_sim => $max_iden });
1531    #       print STDERR "max_iden=$max_iden, ", scalar @$repseqs,"\n";
1532            $iterations_left--;
1533            $incr = $incr / 2;
1534        }
1535    
1536        foreach my $id (keys(%$lost))
1537        {
1538            my $rep_by = $lost->{$id};
1539            while ($lost->{$rep_by})
1540            {
1541                $rep_by = $lost->{$rep_by};
1542            }
1543            push(@{$representing->{$rep_by}},$id);
1544        }
1545        return ($repseqs,$representing);
1546    }
1547    
1548    sub add_to_lost {
1549        my($lost,$representing) = @_;
1550    
1551        foreach my $id (keys(%$representing))
1552        {
1553            my $x = $representing->{$id};
1554            foreach my $lost_id (@$x)
1555            {
1556                $lost->{$lost_id} = $id;
1557            }
1558        }
1559    }
1560    
1561    ###############################################################################
1562    
1563  1;  1;

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3