[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.4, Wed Feb 28 21:55:44 2007 UTC revision 1.5, Mon Mar 19 18:58:10 2007 UTC
# Line 2  Line 2 
2    
3  use strict;  use strict;
4  use gjoparseblast;  use gjoparseblast;
 use Data::Dumper;  
 use Carp;  
5    
6  #  These four variables are used as globals below.  They would be changed  eval { use Data::Dumper };  # Not in all installations
7  #  for the SEED environment.  formatdb is run as a system command, and might  
8    #  These variables are used as globals below.  They have been made SEED-aware,
9    #  but not SEED-dependent.  formatdb is run as a system command, and might
10  #  be encapsulated in a run() construct for the SEED.  blastall is in an  #  be encapsulated in a run() construct for the SEED.  blastall is in an
11  #  input pipe.  #  input pipe.
12    
13  my $tmp_dir  = ".";  my $formatdb = $FIG_Config::ext_bin ? "$FIG_Config::ext_bin/formatdb"
14  my $formatdb = "formatdb";                                      : "formatdb";
15  my $blastall = "blastall";  my $blastall = $FIG_Config::ext_bin ? "$FIG_Config::ext_bin/blastall"
16                                        : "blastall";
17    
18  require Exporter;  require Exporter;
19    
20  our @ISA = qw(Exporter);  our @ISA = qw(Exporter);
21  our @EXPORT = qw( representative_sequences  our @EXPORT = qw( representative_sequences
22                    rep_seq_2                    rep_seq_2
23                      rep_seq
24                  );                  );
25    
26  #===============================================================================  #===============================================================================
27  #  Construct a representative set of related sequences:  #  Build or add to a set of representative sequences (if you to not want an
28    #  enrichment of sequences around a focus sequence (called the reference), this
29    #  is probably the subroutine that you want).
30  #  #
31  #    \@repseqs = representative_sequences( $ref, \@seqs, $max_sim, \%options );  #    \@reps = rep_seq( \@reps, \@new, \%options );
32    #    \@reps = rep_seq(         \@new, \%options );
33  #  #
34  #  or  #  or
35  #  #
36  #    ( \@repseqs, \%representing, \@low_sim ) = representative_sequences( $ref,  #    ( \@reps, \%representing ) = rep_seq( \@reps, \@new, \%options );
37  #                                                 \@seqs, $max_sim, \%options );  #    ( \@reps, \%representing ) = rep_seq(         \@new, \%options );
38  #  #
39  #  Build or add to a set of representative sequences (if you to not want an  #  or
 #  enrichment of sequences around a focus sequence (called the reference), this  
 #  is probably the subroutine that you want).  
40  #  #
41  #    \@reps = rep_seq_2( \@reps, \@new, \%options );  #    \@reps = rep_seq_2( \@reps, \@new, \%options );
42  #    \@reps = rep_seq_2(         \@new, \%options );  #    \@reps = rep_seq_2(         \@new, \%options );
# Line 43  Line 46 
46  #    ( \@reps, \%representing ) = rep_seq_2( \@reps, \@new, \%options );  #    ( \@reps, \%representing ) = rep_seq_2( \@reps, \@new, \%options );
47  #    ( \@reps, \%representing ) = rep_seq_2(         \@new, \%options );  #    ( \@reps, \%representing ) = rep_seq_2(         \@new, \%options );
48  #  #
49    #  Construct a representative set of related sequences:
50    #
51    #    \@repseqs = representative_sequences( $ref, \@seqs, $max_sim, \%options );
52    #
53    #  or
54    #
55    #    ( \@repseqs, \%representing, \@low_sim ) = representative_sequences( $ref,
56    #                                                 \@seqs, $max_sim, \%options );
57    #
58  #  Output:  #  Output:
59  #  #
60  #    \@repseqs  Reference to the list of retained (representative subset)  #    \@repseqs  Reference to the list of retained (representative subset)
# Line 121  Line 133 
133  #                   (a parameter for representative_sequences, but an option  #                   (a parameter for representative_sequences, but an option
134  #                   for rep_seq_2).  #                   for rep_seq_2).
135  #  #
136    #        save_tmp   Do not delete temporary files upon completion (for debug)
137    #
138  #        sim_meas   Measure similarity for inclusion or exclusion by  #        sim_meas   Measure similarity for inclusion or exclusion by
139  #                  'identity_fraction' (default), 'positive_fraction', or  #                  'identity_fraction' (default), 'positive_fraction', or
140  #                  'score_per_position'  #                  'score_per_position'
# Line 256  Line 270 
270  #  #
271  #        7c. If similarity of seq1 and seq2 > max_sim, veto seq2  #        7c. If similarity of seq1 and seq2 > max_sim, veto seq2
272  #  #
273  #        7d. Next seq2  #        7d. Next seq2
274    #
275    #    8. If there are more sequences to examine, go to 3.
276    #
277    #    9. Collect the sequences marked for keeping.
278    #
279    #===============================================================================
280    
281    #===============================================================================
282    #  Build or add to a set of representative sequences.  The difference of
283    #  rep_seq_2 and rep_seq is that rep_seq can have multiple representatives
284    #  in the blast database for a given group.  This helps prevent fragmentation
285    #  of clusters.
286    #
287    #    \@reps = rep_seq( \@reps, \@new, \%options );
288    #    \@reps = rep_seq(         \@new, \%options );
289    #
290    #  or
291    #
292    #    ( \@reps, \%representing ) = rep_seq( \@reps, \@new, \%options );
293    #    ( \@reps, \%representing ) = rep_seq(         \@new, \%options );
294    #
295    #===============================================================================
296    
297    sub rep_seq
298    {
299        ref $_[0] eq 'ARRAY'
300                or print STDERR "First parameter of rep_seq must be an ARRAY reference\n"
301                   and return undef;
302    
303        my ( $reps, $seqs, $options ) = ( [], undef, {} );
304    
305        if    ( @_ == 1 )
306        {
307            $seqs = shift;
308        }
309    
310        elsif ( @_ == 2 )
311        {
312            if    ( ref $_[1] eq 'ARRAY' ) { ( $reps, $seqs )    = @_ }
313            elsif ( ref $_[1] eq 'HASH'  ) { ( $seqs, $options ) = @_ }
314            else
315            {
316                print STDERR "Second parameter of rep_seq must be an ARRAY or HASH reference\n";
317                return undef;
318            }
319        }
320    
321        elsif ( @_ == 3 )
322        {
323            if ( ref $_[1] ne 'ARRAY' )
324            {
325                print STDERR "Second parameter of 3 parameter rep_seq must be an ARRAY reference\n";
326                return undef;
327            }
328            if ( ref $_[2] ne 'HASH' )
329            {
330                print STDERR "Third parameter of 3 parameter rep_seq must be a HASH reference\n";
331                return undef;
332            }
333    
334            ( $reps, $seqs, $options ) = @_;
335        }
336        else
337        {
338            print STDERR "rep_seq called with @{[scalar @_]} parameters\n";
339            return undef;
340        }
341    
342        # ---------------------------------------# Default values for options
343    
344        my $max_sim   = 0.80;                  # Retain 80% identity of less
345        my $logfile   = undef;                 # Log file of sequences processed
346        my $max_e_val = 0.01;                  # Blast E-value to decrease output
347        my $sim_meas  = 'identity_fraction';   # Use sequence identity as measure
348    
349        #  Two questionable decisions:
350        #     1. Be painfully flexible on option names.
351        #     2. Silently fix bad parameter values.
352    
353        foreach ( keys %$options )
354        {
355            my $value = $options->{ $_ };
356            if    ( m/^log/i )                  #  logfile
357            {
358                next if ! $value;
359                $logfile = ( ref $value eq "GLOB" ? $value : \*STDOUT );
360                select( ( select( $logfile ), $| = 1 )[0] );  #  autoflush on
361            }
362            elsif ( m/max/i && m/sim/i )        #  max(imum)_sim(ilarity)
363            {
364                $value += 0;
365                $value  = 0 if $value < 0;
366                $value  = 1 if $value > 1;
367                $max_sim = $value;
368            }
369            elsif ( m/max/i || m/[ep]_?val/i )  #  Other "max" tests must come first
370            {
371                $value += 0;
372                $value  = 0 if $value < 0;
373                $max_e_val = $value;
374            }
375            elsif ( m/sim/i || m/meas/i )       #  sim(ilarity)_meas(ure)
376            {
377                $sim_meas = standardize_similarity_measure( $value );
378            }
379            elsif ( m/save_?te?mp/i )           #  group temporary files
380            {
381                $options->{ savetmp } = 1;
382            }
383            else
384            {
385                print STDERR "WARNING: rep_seq bad option ignored: '$_' => '$value'\n";
386            }
387        }
388    
389        #  Check sequence ids for duplicates:
390    
391        my $reps2 = [];
392        my $seen  = {};
393    
394        my $id;
395        foreach ( @$reps )
396        {
397            $id = $_->[0];
398            if ( $seen->{ $id }++ )
399            {
400                print STDERR "Duplicate sequence id '$id' skipped by rep_seq\n";
401            }
402            else
403            {
404                push @$reps2, $_;
405            }
406        }
407    
408        my $seqs2 = [];
409        foreach ( @$seqs )
410        {
411            $id = $_->[0];
412            if ( $seen->{ $id }++ )
413            {
414                print STDERR "Duplicate sequence id '$_[0]' skipped by rep_seq\n";
415            }
416            else
417            {
418                push @$seqs2, $_;
419            }
420        }
421    
422        #  If no preexisting representatives, then take first sequence:
423    
424        @$reps2 or ( @$reps2 = shift @$seqs2 );
425    
426        if ( $logfile ) { foreach ( @$reps2 ) { print $logfile "$_->[0]\n" } }
427    
428        #  Search each rep sequence against itself to get max_bpp
429    
430        my ( $tmp_dir, $save_tmp ) = &tmp_dir( $options );
431        $tmp_dir or print STDERR "Unable to locate temporary file directory\n"
432                 and return;
433        my $db = "$tmp_dir/tmp_blast_db_$$";
434        my $protein = are_protein( $reps2 );
435        my %max_bpp;
436        my $entry;
437        foreach $entry ( @$reps2 )
438        {
439            $max_bpp{ $entry->[0] } = self_bpp( $db, $entry, $protein );
440        }
441    
442        my $naln = 10;  # Alignments to make
443        my $self =  1;  # Self match will never occur
444        my $prog = $protein ? 'blastp' : 'blastn';
445        my $blast_opt = "-e $max_e_val -v $naln -b $naln -F F -a 2";
446        $blast_opt   .= " -r 1 -q -1" if ! $protein;
447    
448        #  List of who is represented by a sequence:
449    
450        my %group = map { $_->[0] => [] } @$reps2;
451    
452        #  Groups can have more than one representative in the blast database:
453    
454        my $rep4blast = [ @$reps2 ];                        # initial reps
455        my %group_id  = map { $_->[0] => $_->[0] } @$reps2; # represent self
456    
457        #  Search each sequence against the database.
458    
459        my ( $seq, $bpp_max, $sid, $gid );
460        my $newdb = 1;
461    
462        foreach $entry ( @$seqs2 )
463        {
464            #  Is it time to rebuild a BLAST database?
465    
466            if ( $newdb )
467            {
468                make_blast_db( $db, $rep4blast, $protein );
469                $newdb = 0;
470            }
471    
472            #  Do the blast analysis.  Returned records are of the form:
473            #
474            #        0    1     2     3    4     5     6     7      8     9     10     11
475            #     [ qid, qdef, qlen, sid, sdef, slen, scr, e_val, n_mat, n_id, n_pos, n_gap ]
476            #
477            #  $tophit = [ $score, $blast_record, $surething ]
478    
479            $id  = $entry->[0];
480            $seq = $entry->[2];
481            my ( $tophit ) = sort { $b->[0] <=> $a->[0] }
482                             map  { in_group( $_, $max_sim, $sim_meas, $max_bpp{ $_->[3] } ) }
483                             top_blast_per_subject( $prog, $db, $id, $seq, $self, $blast_opt );
484    
485            # It matches an existing representative
486    
487            if ( $tophit )
488            {
489                $sid = $tophit->[1]->[3];      # id of the best matching sequence
490                $gid = $group_id{ $sid };      # look up representative for group
491                push @{ $group{ $gid } }, $id;  # add sequence to list in group
492                $group_id{ $id } = $gid;       # record group for id
493                print $logfile "$id\t$gid\n" if $logfile;
494    
495                # Add sequence to blast database if it is not a 'surething'
496    
497                if ( ! $tophit->[2] )
498                {
499                    push @$rep4blast, $entry;
500                    $max_bpp{ $id } = self_bpp( $db, $entry, $protein ) if $sim_meas =~ /^sc/;
501                    $newdb = 1;
502                }
503            }
504    
505            # It is a new representative
506    
507            else
508            {
509                push @$reps2, $entry;
510                push @$rep4blast, $entry;
511                $group{ $id } = [];
512                $group_id{ $id } = $id;   #  represent self
513                $max_bpp{ $id } = self_bpp( $db, $entry, $protein ) if $sim_meas =~ /^sc/;
514                $newdb = 1;
515                print $logfile "$id\n" if $logfile;
516            }
517        }
518    
519        if    ( $save_tmp ) {}  #  This is not used, might be in future
520        elsif ( $protein  ) { unlink $db, "$db.psq", "$db.pin", "$db.phr" }
521        else                { unlink $db, "$db.nsq", "$db.nin", "$db.nhr" }
522    
523        #  Return the surviving sequence entries, and optionally the hash of
524        #  ids represented by each survivor:
525    
526        wantarray ? ( $reps2, \%group ) : $reps2;
527    }
528    
529    
530    #===============================================================================
531    #  Caluculate sequence similarity according to the requested measure, and return
532    #  empty list if lower than max_sim.  Otherwise, return the hit and and
533    #  whether the hit is really strong:
534    #
535    #     [ $score, $hit, $surething ] = in_group( $hit, $max_sim, $measure, $bpp_max )
536    #     ()                           = in_group( $hit, $max_sim, $measure, $bpp_max )
537    #
538    #  $hit is a structure with blast information:
539    #
540    #  [ qid, qdef, qlen, sid, sdef, slen, scr, e_val, n_mat, n_id, n_pos, n_gap ]
541    #
542    #  The surething is the similarity for which $max_sim is 4 standard deviations
543    #  lower.
544    #===============================================================================
545    
546    sub in_group
547    {   my ( $hit, $max_sim, $measure, $bpp_max ) = @_;
548    
549        my $n = $hit->[8];           # aligned positions
550        return () if ( $n <= 0 );
551    
552        my $m;                       # matched positions
553    
554        if    ( $measure =~ /^sc/ ) { $m = $hit->[ 6] / ( $bpp_max || 2 ) } # score/pos
555        elsif ( $measure =~ /^po/ ) { $m = $hit->[10] }  # positives
556        else                        { $m = $hit->[ 9] }  # identities
557    
558        return () if $m < ( $max_sim * $n );
559    
560        my $u      = ( $n > $m ) ? ( $n - $m ) : 0;           #  differing positions
561        my $stddev = sqrt( $m * $u / $n );
562        my $conf   = 4;                      # standard deviations for "surething"
563        $max_sim   = 0.01 if $max_sim < 0.01;
564        my $surething = ( $u + $conf * $stddev ) <= ( ( 1 - $max_sim ) * $n ) ? 1 : 0;
565    
566        [ $m/$n, $hit, $surething ]
567    }
568    
569    
570    #===============================================================================
571    #  Build or add to a set of representative sequences.
572    #
573    #    \@reps = rep_seq_2( \@reps, \@new, \%options );
574    #    \@reps = rep_seq_2(         \@new, \%options );
575    #
576    #  or
577    #
578    #    ( \@reps, \%representing ) = rep_seq_2( \@reps, \@new, \%options );
579    #    ( \@reps, \%representing ) = rep_seq_2(         \@new, \%options );
580    #
581    #===============================================================================
582    
583    sub rep_seq_2
584    {
585        ref $_[0] eq 'ARRAY'
586                or print STDERR "First parameter of rep_seq_2 must be an ARRAY reference\n"
587                   and return undef;
588    
589        my ( $reps, $seqs, $options ) = ( [], undef, {} );
590    
591        if    ( @_ == 1 )
592        {
593            $seqs = shift;
594        }
595    
596        elsif ( @_ == 2 )
597        {
598            if    ( ref $_[1] eq 'ARRAY' ) { ( $reps, $seqs )    = @_ }
599            elsif ( ref $_[1] eq 'HASH'  ) { ( $seqs, $options ) = @_ }
600            else
601            {
602                print STDERR "Second parameter of rep_seq_2 must be an ARRAY or HASH reference\n";
603                return undef;
604            }
605        }
606    
607        elsif ( @_ == 3 )
608        {
609            if ( ref $_[1] ne 'ARRAY' )
610            {
611                print STDERR "Second parameter of 3 parameter rep_seq_2 must be an ARRAY reference\n";
612                return undef;
613            }
614            if ( ref $_[2] ne 'HASH' )
615            {
616                print STDERR "Third parameter of 3 parameter rep_seq_2 must be a HASH reference\n";
617                return undef;
618            }
619    
620            ( $reps, $seqs, $options ) = @_;
621        }
622        else
623        {
624            print STDERR "rep_seq_2 called with @{[scalar @_]} parameters\n";
625            return undef;
626        }
627    
628        # ---------------------------------------# Default values for options
629    
630        my $max_sim   = 0.80;                  # Retain 80% identity of less
631        my $logfile   = undef;                 # Log file of sequences processed
632        my $max_e_val = 0.01;                  # Blast E-value to decrease output
633        my $sim_meas  = 'identity_fraction';   # Use sequence identity as measure
634    
635        #  Two questionable decisions:
636        #     1. Be painfully flexible on option names.
637        #     2. Silently fix bad parameter values.
638    
639        foreach ( keys %$options )
640        {
641            my $value = $options->{ $_ };
642            if    ( m/^log/i )                  #  logfile
643            {
644                next if ! $value;
645                $logfile = ( ref $value eq "GLOB" ? $value : \*STDOUT );
646                select( ( select( $logfile ), $| = 1 )[0] );  #  autoflush on
647            }
648            elsif ( m/max/i && m/sim/i )        #  max(imum)_sim(ilarity)
649            {
650                $value += 0;
651                $value  = 0 if $value < 0;
652                $value  = 1 if $value > 1;
653                $max_sim = $value;
654            }
655            elsif ( m/max/i || m/[ep]_?val/i )  #  Other "max" tests must come first
656            {
657                $value += 0;
658                $value  = 0 if $value < 0;
659                $max_e_val = $value;
660            }
661            elsif ( m/sim/i || m/meas/i )       #  sim(ilarity)_meas(ure)
662            {
663                $sim_meas = standardize_similarity_measure( $value );
664            }
665            else
666            {
667                print STDERR "WARNING: rep_seq_2 bad option ignored: '$_' => '$value'\n";
668            }
669        }
670    
671        #  Check sequence ids for duplicates:
672    
673        my $reps2 = [];
674        my $seen  = {};
675    
676        my $id;
677        foreach ( @$reps )
678        {
679            $id = $_->[0];
680            if ( $seen->{ $id }++ )
681            {
682                print STDERR "Duplicate sequence id '$id' skipped by rep_seq_2\n";
683            }
684            else
685            {
686                push @$reps2, $_;
687            }
688        }
689    
690        my $seqs2 = [];
691        foreach ( @$seqs )
692        {
693            $id = $_->[0];
694            if ( $seen->{ $id }++ )
695            {
696                print STDERR "Duplicate sequence id '$_[0]' skipped by rep_seq_2\n";
697            }
698            else
699            {
700                push @$seqs2, $_;
701            }
702        }
703    
704        #  If no preexisting representatives, then take first sequence:
705    
706        @$reps2 or @$reps2 = ( shift @$seqs2 );
707    
708        if ( $logfile ) { foreach ( @$reps2 ) { print $logfile "$_->[0]\n" } }
709    
710        #  Search each rep sequence against itself for max_bpp
711    
712        my ( $tmp_dir, $save_tmp ) = &tmp_dir( $options );
713        $tmp_dir or print STDERR "Unable to locate temporary file directory\n"
714                 and return;
715        my $db = "$tmp_dir/tmp_blast_db_$$";
716        my $protein = are_protein( $reps2 );
717        my %max_bpp;
718        my $entry;
719        foreach $entry ( @$reps2 )
720        {
721            $max_bpp{ $entry->[0] } = self_bpp( $db, $entry, $protein );
722        }
723    
724        my $naln = 10;  # Alignments to make
725        my $self =  1;  # Self match will never occur
726        my $prog = $protein ? 'blastp' : 'blastn';
727        my $blast_opt = "-e $max_e_val -v $naln -b $naln -F F -a 2";
728        $blast_opt   .= " -r 1 -q -1" if ! $protein;
729    
730        #  List of who is represented by a sequence:
731    
732        my %group = map { $_->[0] => [] } @$reps2;
733    
734        #  Search each sequence against the database.
735    
736        my ( $seq, $bpp_max );
737        my $newdb = 1;
738    
739        foreach $entry ( @$seqs2 )
740        {
741            #  Is it time to rebuild a BLAST database?
742    
743            if ( $newdb )
744            {
745                make_blast_db( $db, $reps2, $protein );
746                $newdb = 0;
747            }
748    
749            #  Do the blast analysis.  Returned records are of the form:
750            #
751            #        0    1     2     3    4     5     6     7      8     9     10     11
752            #     [ qid, qdef, qlen, sid, sdef, slen, scr, e_val, n_mat, n_id, n_pos, n_gap ]
753    
754            $id  = $entry->[0];
755            $seq = $entry->[2];
756            my ( $tophit ) = sort { $b->[0] <=> $a->[0] }
757                             grep { $_->[0] >= $max_sim }
758                             map  { [ seq_similarity( $_, $sim_meas, $max_bpp{ $_->[3] } ), $_ ] }
759                             top_blast_per_subject( $prog, $db, $id, $seq, $self, $blast_opt );
760    
761    
762            # It matches an existing representative
763    
764            if ( $tophit )
765            {
766                # print STDERR join(", ", $tophit->[0], @{$tophit->[1]} ), "\n";
767                push @{ $group{ $tophit->[1]->[3] } }, $id;
768                print $logfile "$id\t$tophit->[1]->[3]\n" if $logfile;
769            }
770    
771            # It is a new representative
772    
773            else
774            {
775                push @$reps2, $entry;
776                $group{ $id } = [];
777                $max_bpp{ $id } = self_bpp( $db, $entry, $protein );
778                $newdb = 1;
779                print $logfile "$id\n" if $logfile;
780            }
781        }
782    
783        if    ( $save_tmp ) {}  #  This is not used, might be in future
784        elsif ( $protein  ) { unlink $db, "$db.psq", "$db.pin", "$db.phr" }
785        else                { unlink $db, "$db.nsq", "$db.nin", "$db.nhr" }
786    
787        #  Return the surviving sequence entries, and optionally the hash of
788        #  ids represented by each survivor:
789    
790        wantarray ? ( $reps2, \%group ) : $reps2;
791    }
792    
793    
794    #===============================================================================
795    #  Construct a representative set of related sequences:
796    #
797    #    \@repseqs = representative_sequences( $ref, \@seqs, $max_sim, \%options );
798  #  #
799  #    8. If there are more sequences to examine, go to 3.  #  or
800  #  #
801  #    9. Collect the sequences marked for keeping.  #    ( \@repseqs, \%representing, \@low_sim ) = representative_sequences( $ref,
802    #                                                 \@seqs, $max_sim, \%options );
803  #  #
804  #===============================================================================  #===============================================================================
   
805  sub representative_sequences {  sub representative_sequences {
806      my $seqs = ( shift @_ || shift @_ );  #  If $ref is undef, shift again      my $seqs = ( shift @_ || shift @_ );  #  If $ref is undef, shift again
807      ref( $seqs ) eq "ARRAY"      ref( $seqs ) eq "ARRAY"
# Line 390  Line 928 
928      #  Build blast database (it includes the reference):      #  Build blast database (it includes the reference):
929    
930      my $protein = are_protein( $seqs );      my $protein = are_protein( $seqs );
931      my $db = ( $tmp_dir ? "$tmp_dir/" : "" ) . "tmp_blast_db_$$";      my ( $tmp_dir, $save_tmp ) = &tmp_dir( $options );
932        $tmp_dir or print STDERR "Unable to locate temporary file directory\n"
933                 and return;
934        my $db = "$tmp_dir/tmp_blast_db_$$";
935      make_blast_db( $db, [ $ref, @$seqs ], $protein );      make_blast_db( $db, [ $ref, @$seqs ], $protein );
936    
937      #  Search query against new database      #  Search query against new database
# Line 419  Line 960 
960      my %hit = ();      my %hit = ();
961      @ref_hits = grep { my $sid = $_->[3]; $hit{ $sid } = 1; ( $sid ne $ref_id ) } @ref_hits;      @ref_hits = grep { my $sid = $_->[3]; $hit{ $sid } = 1; ( $sid ne $ref_id ) } @ref_hits;
962    
963      my %keep = ();      my %group = ();
964      $keep{ $ref_id } = [];      $group{ $ref_id } = [];
965      my %veto = ();      my %veto = ();
966      my $n_to_do = @ref_hits;      my $n_to_do = @ref_hits;
967      my $rebuild_d_n  = 40;      my $rebuild_d_n  = 40;
# Line 439  Line 980 
980          if ( $sim > ( $use_ref ? $max_ref_sim : $max_sim ) )          if ( $sim > ( $use_ref ? $max_ref_sim : $max_sim ) )
981          {          {
982              $veto{ $id } = 1;              $veto{ $id } = 1;
983              push @{ $keep{ $ref_id } }, $id;   #  Log the sequences represented              push @{ $group{ $ref_id } }, $id;   #  Log the sequences represented
984              $n_to_do--;              $n_to_do--;
985          }          }
986          else          else
# Line 452  Line 993 
993    
994      if ( $logfile )      if ( $logfile )
995      {      {
996          print $logfile join( "\t", $ref_id, @{ $keep{ $ref_id } } ), "\n";          print $logfile join( "\t", $ref_id, @{ $group{ $ref_id } } ), "\n";
997      }      }
998    
999      #  Search each sequence against the database.      #  Search each sequence against the database.
# Line 483  Line 1024 
1024          }          }
1025    
1026          $n_to_do--;          $n_to_do--;
1027          $keep{ $id1 } = [];          $group{ $id1 } = [];
1028    
1029          $max_sim_1 = $max_sim{ $id1 };          $max_sim_1 = $max_sim{ $id1 };
1030          $bpp_max = undef;          $bpp_max = undef;
# Line 492  Line 1033 
1033          {          {
1034              $bpp_max ||= $_->[6] / $_->[8];              $bpp_max ||= $_->[6] / $_->[8];
1035              $id2 = $_->[3];              $id2 = $_->[3];
1036              next if ( $veto{ $id2 } || $keep{ $id2 } );              next if ( $veto{ $id2 } || $group{ $id2 } );
1037              $max_sim_2 = $max_sim{ $id2 };              $max_sim_2 = $max_sim{ $id2 };
1038              $max_sim_2 = $max_sim_1 if ( $max_sim_1 > $max_sim_2 );              $max_sim_2 = $max_sim_1 if ( $max_sim_1 > $max_sim_2 );
1039              if ( seq_similarity( $_, $sim_meas, $bpp_max ) > $max_sim_2 )              if ( seq_similarity( $_, $sim_meas, $bpp_max ) > $max_sim_2 )
1040              {              {
1041                  $veto{ $id2 } = 1;                  $veto{ $id2 } = 1;
1042                  push @{ $keep{ $id1 } }, $id2;  #  Log the sequences represented                  push @{ $group{ $id1 } }, $id2;  #  Log the sequences represented
1043                  $n_to_do--;                  $n_to_do--;
1044              }              }
1045          }          }
1046    
1047          if ( $logfile )          if ( $logfile )
1048          {          {
1049              print $logfile join( "\t", $id1, @{ $keep{ $id1 } } ), "\n";              print $logfile join( "\t", $id1, @{ $group{ $id1 } } ), "\n";
1050          }          }
1051      }      }
1052    
1053      if ( $protein ) { unlink $db, "$db.psq", "$db.pin", "$db.phr" }      if    ( $save_tmp ) {}  #  This is not used, might be in future
1054        elsif ( $protein  ) { unlink $db, "$db.psq", "$db.pin", "$db.phr" }
1055      else            { unlink $db, "$db.nsq", "$db.nin", "$db.nhr" }      else            { unlink $db, "$db.nsq", "$db.nin", "$db.nhr" }
1056    
1057      #  Return the surviving sequence entries, and optionally the hash of      #  Return the surviving sequence entries, and optionally the hash of
1058      #  ids represented by each survivor:      #  ids represented by each survivor:
1059    
1060      my $kept = [ $ref, grep { $keep{ $_->[0] } } @$seqs ];      my $kept = [ $ref, grep { $group{ $_->[0] } } @$seqs ];
   
     wantarray ? ( $kept, \%keep, [ grep { ! $hit{ $_->[0] } } @$seqs ] ) : $kept;  
 }  
   
   
 #===============================================================================  
 #  Build or add to a set of representative sequences.  
 #  
 #    \@reps = rep_seq_2( \@reps, \@new, \%options );  
 #    \@reps = rep_seq_2(         \@new, \%options );  
 #  
 #  or  
 #  
 #    ( \@reps, \%representing ) = rep_seq_2( \@reps, \@new, \%options );  
 #    ( \@reps, \%representing ) = rep_seq_2(         \@new, \%options );  
 #  
 #===============================================================================  
   
 sub rep_seq_2  
 {  
     ref $_[0] eq 'ARRAY'  
             or print STDERR "First parameter of rep_seq_2 must be an ARRAY reference\n"  
                and return undef;  
   
     my ( $reps, $seqs, $options ) = ( [], undef, {} );  
   
     if    ( @_ == 1 )  
     {  
         $seqs = shift;  
     }  
   
     elsif ( @_ == 2 )  
     {  
         if    ( ref $_[1] eq 'ARRAY' ) { ( $reps, $seqs )    = @_ }  
         elsif ( ref $_[1] eq 'HASH'  ) { ( $seqs, $options ) = @_ }  
         else  
         {  
             print STDERR "Second parameter of rep_seq_2 must be an ARRAY or HASH reference\n";  
             return undef;  
         }  
     }  
   
     elsif ( @_ == 3 )  
     {  
         if ( ref $_[1] ne 'ARRAY' )  
         {  
             print STDERR "Second parameter of 3 parameter rep_seq_2 must be an ARRAY reference\n";  
             return undef;  
         }  
         if ( ref $_[2] ne 'HASH' )  
         {  
             print STDERR "Third parameter of 3 parameter rep_seq_2 must be a HASH reference\n";  
             return undef;  
         }  
   
         ( $reps, $seqs, $options ) = @_;  
     }  
     else  
     {  
         print STDERR "rep_seq_2 called with @{[scalar @_]} parameters\n";  
         return undef;  
     }  
   
     # ---------------------------------------# Default values for options  
   
     my $max_sim   = 0.80;                  # Retain 80% identity of less  
     my $logfile   = undef;                 # Log file of sequences processed  
     my $max_e_val = 0.01;                  # Blast E-value to decrease output  
     my $sim_meas  = 'identity_fraction';   # Use sequence identity as measure  
   
     #  Two questionable decisions:  
     #     1. Be painfully flexible on option names.  
     #     2. Silently fix bad parameter values.  
   
     foreach ( keys %$options )  
     {  
         my $value = $options->{ $_ };  
         if    ( m/^log/i )                  #  logfile  
         {  
             next if ! $value;  
             $logfile = ( ref $value eq "GLOB" ? $value : \*STDOUT );  
             select( ( select( $logfile ), $| = 1 )[0] );  #  autoflush on  
         }  
         elsif ( m/max/i && m/sim/i )        #  max(imum)_sim(ilarity)  
         {  
             $value += 0;  
             $value  = 0 if $value < 0;  
             $value  = 1 if $value > 1;  
             $max_sim = $value;  
         }  
         elsif ( m/max/i || m/[ep]_?val/i )  #  Other "max" tests must come first  
         {  
             $value += 0;  
             $value  = 0 if $value < 0;  
             $max_e_val = $value;  
         }  
         elsif ( m/sim/i || m/meas/i )       #  sim(ilarity)_meas(ure)  
         {  
             $sim_meas = standardize_similarity_measure( $value );  
         }  
         else  
         {  
             print STDERR "WARNING: rep_seq_2 bad option ignored: '$_' => '$value'\n";  
         }  
     }  
   
     #  Check sequence ids for duplicates:  
   
     my $reps2 = [];  
     my $seen  = {};  
   
     my $id;  
     foreach ( @$reps )  
     {  
         $id = $_->[0];  
         if ( $seen->{ $id }++ )  
         {  
             print STDERR "Duplicate sequence id '$id' skipped by rep_seq_2\n";  
         }  
         else  
         {  
             push @$reps2, $_;  
         }  
     }  
   
     my $seqs2 = [];  
     foreach ( @$seqs )  
     {  
         $id = $_->[0];  
         if ( $seen->{ $id }++ )  
         {  
             print STDERR "Duplicate sequence id '$_[0]' skipped by rep_seq_2\n";  
         }  
         else  
         {  
             push @$seqs2, $_;  
         }  
     }  
   
     #  If no preexisting representatives, then take first sequence:  
   
     @$reps2 or @$reps2 = ( shift @$seqs2 );  
   
     if ( $logfile ) { foreach ( @$reps2 ) { print $logfile "$_->[0]\n" } }  
   
     #  Search each rep sequence against itself for max_bpp  
   
     my $db = ( $tmp_dir ? "$tmp_dir/" : "" ) . "tmp_blast_db_$$";  
     my $protein = are_protein( $reps2 );  
   
     my %max_bpp;  
     my $entry;  
     foreach $entry ( @$reps2 )  
     {  
         $max_bpp{ $entry->[0] } = self_bpp( $db, $entry, $protein );  
     }  
   
     my $naln = 10;  # Alignments to make  
     my $self =  1;  # Self match will never occur  
     my $prog = $protein ? 'blastp' : 'blastn';  
     my $blast_opt = "-e $max_e_val -v $naln -b $naln -F F -a 2";  
     $blast_opt   .= " -r 1 -q -1" if ! $protein;  
   
     #  List of who is represented by a sequence:  
   
     my %keep = map { $_->[0] => [] } @$reps2;  
   
     #  Search each sequence against the database.  
   
     my ( $seq, $bpp_max );  
     my $newdb = 1;  
   
     foreach $entry ( @$seqs2 )  
     {  
         #  Is it time to rebuild a BLAST database?  
   
         if ( $newdb )  
         {  
             make_blast_db( $db, $reps2, $protein );  
             $newdb = 0;  
         }  
         #  Do the blast analysis.  Returned records are of the form:  
         #  
         #        0    1     2     3    4     5     6     7      8     9     10     11  
         #     [ qid, qdef, qlen, sid, sdef, slen, scr, e_val, n_mat, n_id, n_pos, n_gap ]  
   
         $id  = $entry->[0];  
         $seq = $entry->[2];  
         my ( $tophit ) = sort { $b->[0] <=> $a->[0] }  
                          grep { $_->[0] >= $max_sim }  
                          map  { [ seq_similarity( $_, $sim_meas, $max_bpp{ $_->[3] } ), $_ ] }  
                          top_blast_per_subject( $prog, $db, $id, $seq, $self, $blast_opt );  
   
   
         # It matches an existing representative  
   
         if ( $tophit )  
         {  
             # print STDERR join(", ", $tophit->[0], @{$tophit->[1]} ), "\n";  
             push @{ $keep{ $tophit->[1]->[3] } }, $id;  
             print $logfile "$id\t$tophit->[1]->[3]\n" if $logfile;  
         }  
   
         # It is a new representative  
   
         else  
         {  
             push @$reps2, $entry;  
             $keep{ $id } = [];  
             $max_bpp{ $id } = self_bpp( $db, $entry, $protein );  
             $newdb = 1;  
             print $logfile "$id\n" if $logfile;  
         }  
     }  
   
     if ( $protein ) { unlink $db, "$db.psq", "$db.pin", "$db.phr" }  
     else            { unlink $db, "$db.nsq", "$db.nin", "$db.nhr" }  
1061    
1062      #  Return the surviving sequence entries, and optionally the hash of      wantarray ? ( $kept, \%group, [ grep { ! $hit{ $_->[0] } } @$seqs ] ) : $kept;
     #  ids represented by each survivor:  
   
     wantarray ? ( $reps2, \%keep ) : $reps2;  
1063  }  }
1064    
1065    
# Line 817  Line 1139 
1139    
1140      #  First hit is always a perfect match, so we get bits per position:      #  First hit is always a perfect match, so we get bits per position:
1141      #  This is only used if the measure is bits per position      #  This is only used if the measure is bits per position
1142      if (! $hit->[8]) { return 0 }  
1143      $hit->[6] / $hit->[8];      $hit->[6] / $hit->[8];
1144  }  }
1145    
# Line 884  Line 1206 
1206  {   my ( $prog, $db, $qid, $qseq, $self, $blast_opt, $sort, $no_merge ) = @_;  {   my ( $prog, $db, $qid, $qseq, $self, $blast_opt, $sort, $no_merge ) = @_;
1207    
1208      $cnt++;      $cnt++;
1209      my $query_file = ( $tmp_dir ? "$tmp_dir/" : "" ) . "tmp_blast_seq_${$}_$cnt";      my ( $tmp_dir, $save_tmp ) = &tmp_dir( );
1210        $tmp_dir or print STDERR "Unable to locate temporary file directory\n"
1211                 and return;
1212        my $query_file = "$tmp_dir/tmp_blast_seq_${$}_$cnt";
1213      my    $QFILE;      my    $QFILE;
1214      open  $QFILE, ">$query_file" or die "Could not create sequence file \"$query_file\"";      open  $QFILE, ">$query_file" or die "Could not create sequence file \"$query_file\"";
1215      print $QFILE  ">$qid\n$qseq\n";      print $QFILE  ">$qid\n$qseq\n";
1216      close $QFILE;      close $QFILE;
1217    
1218      my $blast_cmd = "$blastall -p $prog -d '$db' -i '$query_file' $blast_opt";      my $blast_cmd = "$blastall -p $prog -d '$db' -i '$query_file' $blast_opt";
1219    
1220      my   $BPIPE;      my   $BPIPE;
1221      open $BPIPE, "$blast_cmd |" or die "Could open blast pipe\n";      open $BPIPE, "$blast_cmd |" or die "Could open blast pipe\n";
1222      my $sims = integrate_blast_segments( $BPIPE, $sort, $no_merge, $self );      my $sims = integrate_blast_segments( $BPIPE, $sort, $no_merge, $self );
1223      close $BPIPE;      close $BPIPE;
1224        unlink $query_file if ! $save_tmp;
     unlink $query_file;  
1225    
1226      my $pq = "";  #  Previous query id      my $pq = "";  #  Previous query id
1227      my $ps = "";  #  Previous subject id      my $ps = "";  #  Previous subject id
# Line 1146  Line 1471 
1471  }  }
1472    
1473    
1474    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1475    #  Locate the directory for temporary files in a SEED-aware, but not SEED-
1476    #  dependent manner:
1477    #
1478    #     ( $tmp_dir, $save_tmp ) = tmp_dir( \%options )
1479    #
1480    #  $save_tmp is true if the option 'savetmp' is true.
1481    #- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1482    sub tmp_dir
1483    {
1484        my $options = ref( $_[0] ) eq 'HASH' ? shift : {};
1485    
1486        my $tmp = $options->{ tmp } && -d  $options->{ tmp } ?  $options->{ tmp }
1487                : $FIG_Config::temp && -d  $FIG_Config::temp ?  $FIG_Config::temp
1488                :                      -d '/tmp'             ? '/tmp'
1489                :                                              '.';
1490    
1491        my $save_tmp = $options->{ savetmp } || 0;
1492    
1493        return ( $tmp, $save_tmp );
1494    }
1495    
1496    
1497    
1498  1;  1;

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3