← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty-proposed

« back to all changes in this revision

Viewing changes to Bio/Assembly/Contig.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2009-03-10 07:19:11 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20090310071911-fukqzw54pyb1f0bd
Tags: 1.6.0-2
https://launchpad.net/bugs/324001
* Removed patch system (not used):
  - removed instuctions in debian/rules;
  - removed quilt from Build-Depends in debian/control.
* Re-enabled tests:
  - uncommented test command in debian/rules;
  - uncommented previously missing build-dependencies in debian/control.
  - Re-enabled tests and uncommented build-dependencies accordingly.
* Removed libmodule-build-perl and libtest-harness-perl from
  Build-Depends-Indep (provided by perl-modules).
* Better cleaning of empty directories using find -type d -empty -delete
  instead of rmdir in debian/rules (LP: #324001).

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Assembly/Contig.pm
1
 
# $Id: Contig.pm,v 1.11.4.3 2006/10/02 23:10:12 sendu Exp $
 
1
# $Id: Contig.pm 14942 2008-10-21 04:01:57Z fangly $
2
2
#
3
3
# BioPerl module for Bio::Assembly::Contig
4
4
#   Mostly based on Bio::SimpleAlign by Ewan Birney
22
22
    use Bio::Assembly::IO;
23
23
 
24
24
    # Assembly loading methods
25
 
    $aio = new Bio::Assembly::IO(-file=>"test.ace.1",
 
25
    $aio = Bio::Assembly::IO->new(-file=>"test.ace.1",
26
26
                               -format=>'phrap');
27
27
 
28
28
    $assembly = $aio->next_assembly;
210
210
use Bio::SeqFeature::Collection;
211
211
use Bio::Seq::PrimaryQual;
212
212
 
 
213
use Scalar::Util qw(weaken);
 
214
 
213
215
use base qw(Bio::Root::Root Bio::Align::AlignI);
214
216
 
215
217
=head1 Object creator
217
219
=head2 new
218
220
 
219
221
 Title     : new
220
 
 Usage     : my $contig = new Bio::Assembly::Contig();
 
222
 Usage     : my $contig = Bio::Assembly::Contig->new();
221
223
 Function  : Creates a new contig object
222
224
 Returns   : Bio::Assembly::Contig
223
 
 Args      : -source => string representing the source
224
 
                        program where this contig came
225
 
                        from
226
 
             -id => contig unique ID
 
225
 Args      : -id => contig unique ID
 
226
             -source => string for the sequence assembly program used
 
227
             -collection => Bio::SeqFeature::Collection instance
227
228
 
228
229
=cut
229
230
 
230
231
#-----------
231
232
sub new {
232
233
#-----------
233
 
    my ($class,@args) = @_;
 
234
    my ($class, @args) = @_;
234
235
 
235
236
    my $self = $class->SUPER::new(@args);
236
237
 
237
 
    my ($src, $id) = $self->_rearrange([qw(SOURCE ID)], @args);
 
238
    my ($src, $id, $collection) = $self->_rearrange([qw(SOURCE ID COLLECTION)], @args);
238
239
    $src && $self->source($src);
239
240
    ($id && $self->id($id)) || ($self->{'_id'} = 'NoName'); # Alignment (contig) name
240
241
    ($id && $self->id($id)) || ($self->{'_source'} = 'Unknown'); # Program used to build the contig
243
244
    # Bio::SimpleAlign derived fields (check which ones are needed for AlignI compatibility)
244
245
    $self->{'_elem'} = {}; # contig elements: aligned sequence objects (keyed by ID)
245
246
    $self->{'_order'} = {}; # store sequence order
246
 
#    $self->{'start_end_lists'} = {}; # References to entries in {'_seq'}. Keyed by seq ids.
247
 
#    $self->{'_dis_name'} = {}; # Display names for each sequence
 
247
    # $self->{'start_end_lists'} = {}; # References to entries in {'_seq'}. Keyed by seq ids.
 
248
    # $self->{'_dis_name'} = {}; # Display names for each sequence
248
249
    $self->{'_symbols'} = {}; # List of symbols
249
250
 
250
251
    #Contig specific slots
252
253
    $self->{'_consensus_quality'} = undef;
253
254
    $self->{'_nof_residues'} = 0;
254
255
    $self->{'_nof_seqs'} = 0;
255
 
#    $self->{'_nof_segments'} = 0; # Let's not make it heavier than needed by now...
256
 
    $self->{'_sfc'} = Bio::SeqFeature::Collection->new();
 
256
    # $self->{'_nof_segments'} = 0; # Let's not make it heavier than needed by now...
 
257
    
 
258
    # for cases where SF::Collection is shared between Bio::Assembly::Contig 
 
259
    if ($collection) {
 
260
        $self->throw("Collection must implement Bio::SeqFeature::CollectionI") unless $collection->isa('Bio::SeqFeature::CollectionI');
 
261
        $self->{'_sfc'} = $collection;
 
262
    } else {
 
263
        $self->{'_sfc'} = Bio::SeqFeature::Collection->new()
 
264
    }
257
265
 
258
 
    # Assembly specifcs
 
266
    # Assembly specifics
259
267
    $self->{'_assembly'} = undef; # Reference to a Bio::Assembly::Scaffold object, if contig belongs to one.
260
268
    $self->{'_strand'} = 0; # Reverse (-1) or forward (1), if contig is in a scaffold. 0 otherwise
261
269
    $self->{'_neighbor_start'} = undef; # Will hold a reference to another contig
305
313
    my $assembly = shift;
306
314
 
307
315
    $self->throw("Using non Bio::Assembly::Scaffold object when assign contig to assembly")
308
 
        if (defined $assembly && ! $assembly->isa("Bio::Assembly::Scaffold"));
 
316
    if (defined $assembly && ! $assembly->isa("Bio::Assembly::Scaffold"));
 
317
    # We create a circular reference to a Scaffold object. It is made weak
 
318
    # to prevent memory leaks.
 
319
    $self->{'_assembly'} = $assembly if (defined $assembly); 
 
320
    weaken($self->{'_assembly'});
309
321
 
310
 
    $self->{'_assembly'} = $assembly if (defined $assembly);
311
322
    return $self->{'_assembly'};
312
323
}
313
324
 
330
341
    my $self = shift;
331
342
    my $ori = shift;
332
343
 
333
 
        if (defined $ori) {
334
 
    $self->throw("Contig strand must be either 1, -1 or 0")
 
344
    if (defined $ori) {
 
345
        $self->throw("Contig strand must be either 1, -1 or 0")
335
346
            unless $ori == 1 || $ori == 0 || $ori == -1;
336
 
 
337
 
    $self->{'_strand'} = $ori;
 
347
        $self->{'_strand'} = $ori;
338
348
    }
339
349
 
340
350
    return $self->{'_strand'};
357
367
    my $ref = shift;
358
368
 
359
369
    $self->throw("Trying to assign a non Bio::Assembly::Contig object to upstream contig")
360
 
        if (defined $ref && ! $ref->isa("Bio::Assembly::Contig"));
 
370
        if (defined $ref && ! $ref->isa("Bio::Assembly::Contig"));
361
371
 
362
372
    $self->{'_neighbor_start'} = $ref if (defined $ref);
363
373
    return $self->{'_neighbor_start'};
380
390
    my $ref = shift;
381
391
 
382
392
    $self->throw("Trying to assign a non Bio::Assembly::Contig object to downstream contig")
383
 
        if (defined $ref && ! $ref->isa("Bio::Assembly::Contig"));
 
393
        if (defined $ref && ! $ref->isa("Bio::Assembly::Contig"));
384
394
    $self->{'_neighbor_end'} = $ref if (defined $ref);
385
395
    return $self->{'_neighbor_end'};
386
396
}
424
434
    # Adding shortcuts for aligned sequence features
425
435
    $flag = 0 unless (defined $flag);
426
436
    if ($flag && defined $self->{'_consensus_sequence'}) {
427
 
        foreach my $feat (@$args) {
428
 
            next if (defined $feat->seq);
429
 
            $feat->attach_seq($self->{'_consensus_sequence'});
430
 
        }
 
437
        foreach my $feat (@$args) {
 
438
            next if (defined $feat->seq);
 
439
            $feat->attach_seq($self->{'_consensus_sequence'});
 
440
        }
431
441
    } elsif (!$flag) { # Register aligned sequence features
432
 
        foreach my $feat (@$args) {
433
 
            if (my $seq = $feat->entire_seq()) {
434
 
                my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
435
 
                $self->warn("Adding contig feature attached to unknown sequence $seqID!")
436
 
                    unless (exists $self->{'_elem'}{$seqID});
437
 
                my $tag = $feat->primary_tag;
438
 
                $self->{'_elem'}{$seqID}{'_feat'}{$tag} = $feat;
439
 
            }
440
 
        }
 
442
        foreach my $feat (@$args) {
 
443
            if (my $seq = $feat->entire_seq()) {
 
444
                my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
 
445
                $self->warn("Adding contig feature attached to unknown sequence $seqID!")
 
446
                    unless (exists $self->{'_elem'}{$seqID});
 
447
                my $tag = $feat->primary_tag;
 
448
                $self->{'_elem'}{$seqID}{'_feat'}{$tag} = $feat;
 
449
            }
 
450
        }
441
451
    }
442
452
 
443
453
    # Add feature to feature collection
461
471
 
462
472
    # Removing shortcuts for aligned sequence features
463
473
    foreach my $feat (@args) {
464
 
        if (my $seq = $feat->entire_seq()) {
465
 
            my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
466
 
            my $tag = $feat->primary_tag;
467
 
            $tag =~ s/:$seqID$/$1/g;
468
 
            delete( $self->{'_elem'}{$seqID}{'_feat'}{$tag} )
469
 
                if (exists $self->{'_elem'}{$seqID}{'_feat'}{$tag} &&
470
 
                    $self->{'_elem'}{$seqID}{'_feat'}{$tag} eq $feat);
471
 
        }
 
474
        if (my $seq = $feat->entire_seq()) {
 
475
            my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
 
476
            my $tag = $feat->primary_tag;
 
477
            $tag =~ s/:$seqID$/$1/g;
 
478
            delete( $self->{'_elem'}{$seqID}{'_feat'}{$tag} )
 
479
                if (exists $self->{'_elem'}{$seqID}{'_feat'}{$tag} &&
 
480
                $self->{'_elem'}{$seqID}{'_feat'}{$tag} eq $feat);
 
481
        }
472
482
    }
473
 
 
 
483
    
 
484
    # Removing Bio::SeqFeature::Collection features
474
485
    return $self->{'_sfc'}->remove_features(\@args);
475
486
}
476
487
 
486
497
 
487
498
sub get_features_collection {
488
499
    my $self = shift;
489
 
 
490
500
    return $self->{'_sfc'};
491
501
}
492
502
 
 
503
=head2 remove_features_collection
 
504
 
 
505
 Title     : remove_features_collection
 
506
 Usage     : $contig->remove_features_collection()
 
507
 Function  : Remove the collection of all contig features. It is useful
 
508
             to save some memory (when contig features are not needed).
 
509
 Returns   : none
 
510
 Argument  : none
 
511
 
 
512
=cut
 
513
 
 
514
sub remove_features_collection {
 
515
    my $self = shift;
 
516
    # Removing shortcuts for aligned sequence features
 
517
    for my $seqID (keys %{$self->{'_elem'}}) {
 
518
        delete $self->{'_elem'}{$seqID};
 
519
    }
 
520
    # Removing Bio::SeqFeature::Collection features
 
521
    $self->{'_sfc'} = {};
 
522
    return;
 
523
}
 
524
 
493
525
=head1 Coordinate system's related methods
494
526
 
495
527
See L<Coordinate_Systems> above.
533
565
    my $out_ID = ( split(' ',$type_out) )[1];
534
566
 
535
567
    if ($in_ID  ne 'consensus') {
536
 
        $read_in  = $self->get_seq_coord( $self->get_seq_by_name($in_ID)  );
537
 
        $self->throw("Can't change coordinates without sequence location for $in_ID")
538
 
            unless (defined $read_in);
 
568
        $read_in  = $self->get_seq_coord( $self->get_seq_by_name($in_ID)  );
 
569
        $self->throw("Can't change coordinates without sequence location for $in_ID")
 
570
            unless (defined $read_in);
539
571
    }
540
572
    if ($out_ID ne 'consensus') {
541
 
        $read_out = $self->get_seq_coord( $self->get_seq_by_name($out_ID) );
542
 
        $self->throw("Can't change coordinates without sequence location for $out_ID")
543
 
            unless (defined $read_out);
 
573
        $read_out = $self->get_seq_coord( $self->get_seq_by_name($out_ID) );
 
574
        $self->throw("Can't change coordinates without sequence location for $out_ID")
 
575
            unless (defined $read_out);
544
576
    }
545
577
 
546
578
    # Performing transformation between coordinates
547
 
  SWITCH1: {
548
 
 
549
 
      # Transformations between contig padded and contig unpadded
550
 
      (($type_in eq 'gapped consensus') && ($type_out eq 'ungapped consensus')) && do {
551
 
          $self->throw("Can't use ungapped consensus coordinates without a consensus sequence")
552
 
              unless (defined $self->{'_consensus_sequence'});
553
 
          $query = &_padded_unpadded($self->{'_consensus_gaps'}, $query);
554
 
          last SWITCH1;
555
 
      };
556
 
      (($type_in eq 'ungapped consensus') && ($type_out eq 'gapped consensus')) && do {
557
 
          $self->throw("Can't use ungapped consensus coordinates without a consensus sequence")
558
 
              unless (defined $self->{'_consensus_sequence'});
559
 
          $query = &_unpadded_padded($self->{'_consensus_gaps'},$query);
560
 
          last SWITCH1;
561
 
      };
562
 
 
563
 
      # Transformations between contig (padded) and read (padded)
564
 
      (($type_in  eq 'gapped consensus') &&
565
 
       ($type_out =~ /^aligned /) && defined($read_out)) && do {
566
 
           $query = $query - $read_out->start() + 1;
567
 
           last SWITCH1;
568
 
       };
569
 
      (($type_in =~ /^aligned /) && defined($read_in) &&
570
 
       ($type_out  eq 'gapped consensus')) && do {
571
 
           $query = $query + $read_in->start() - 1;
572
 
           last SWITCH1;
573
 
       };
574
 
 
575
 
      # Transformations between contig (unpadded) and read (padded)
576
 
      (($type_in eq 'ungapped consensus') &&
577
 
       ($type_out =~ /^aligned /) && defined($read_out)) && do {
578
 
           $query = $self->change_coord('ungapped consensus','gapped consensus',$query);
579
 
           $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
580
 
           last SWITCH1;
581
 
       };
582
 
      (($type_in =~ /^aligned /) && defined($read_in) &&
583
 
       ($type_out eq 'ungapped consensus')) && do {
584
 
           $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
585
 
           $query = $self->change_coord('gapped consensus','ungapped consensus',$query);
586
 
           last SWITCH1;
587
 
       };
588
 
 
589
 
      # Transformations between seq $read_in padded and seq $read_out padded
590
 
      (defined($read_in)  && ($type_in  =~ /^aligned /)  &&
591
 
       defined($read_out) && ($type_out =~ /^aligned /)) && do {
592
 
           $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
593
 
           $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
594
 
           last SWITCH1;
595
 
       };
596
 
 
597
 
      # Transformations between seq $read_in padded and seq $read_out unpadded
598
 
      (defined($read_in)  && ($type_in  =~ /^aligned /)    &&
599
 
       defined($read_out) && ($type_out =~ /^unaligned /)) && do {
600
 
           if ($read_in ne $read_out) {
601
 
               $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
602
 
               $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
603
 
           }
604
 
           my $list_out = $self->{'_elem'}{$out_ID}{'_gaps'};
605
 
           $query = &_padded_unpadded($list_out,$query);
606
 
           # Changing read orientation if read was reverse complemented when aligned
607
 
           if ($read_out->strand == -1) {
608
 
               my ($length) = $read_out->length();
609
 
               $length = $length - &_nof_gaps($list_out,$length);
610
 
               $query  = $length - $query + 1;
611
 
           }
612
 
           last SWITCH1;
613
 
       };
614
 
      (defined($read_in)  && ($type_in  =~ /^unaligned /) &&
615
 
       defined($read_out) && ($type_out =~ /^aligned /))  && do {
616
 
           my $list_in = $self->{'_elem'}{$in_ID}{'_gaps'};
617
 
           # Changing read orientation if read was reverse complemented when aligned
618
 
           if ($read_in->strand == -1) {
619
 
               my ($length) = $read_in->length();
620
 
               $length = $length - &_nof_gaps($list_in,$length);
621
 
               $query  = $length - $query + 1;
622
 
           }
623
 
           $query = &_unpadded_padded($list_in,$query);
624
 
           if ($read_in ne $read_out) {
625
 
               $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
626
 
               $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
627
 
           }
628
 
           last SWITCH1;
629
 
       };
630
 
 
631
 
      # Transformations between seq $read_in unpadded and seq $read_out unpadded
632
 
      (defined($read_in)  && ($type_in  =~ /^unaligned /)    &&
633
 
       defined($read_out) && ($type_out =~ /^unaligned /)) && do {
634
 
           $query = $self->change_coord("unaligned $in_ID","aligned $out_ID",$query);
635
 
           $query = $self->change_coord("aligned $out_ID","unaligned $out_ID",$query);
636
 
           last SWITCH1;
637
 
       };
638
 
 
639
 
      # Transformations between contig (padded) and read (unpadded)
640
 
      (($type_in eq 'gapped consensus') &&
641
 
       ($type_out =~ /^unaligned /) && defined($read_out)) && do {
642
 
           $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
643
 
           $query = $self->change_coord("aligned $out_ID","unaligned $out_ID",$query);
644
 
           last SWITCH1;
645
 
       };
646
 
      (($type_in =~ /^unaligned /) && defined($read_in) &&
647
 
       ($type_out eq 'gapped consensus')) && do {
648
 
           $query = $self->change_coord("unaligned $in_ID","aligned $in_ID",$query);
649
 
           $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
650
 
           last SWITCH1;
651
 
       };
652
 
 
653
 
      # Transformations between contig (unpadded) and read (unpadded)
654
 
      (($type_in eq 'ungapped consensus') &&
655
 
       ($type_out =~ /^unaligned /) && defined($read_out)) && do {
656
 
           $query = $self->change_coord('ungapped consensus','gapped consensus',$query);
657
 
           $query = $self->change_coord('gapped consensus',"unaligned $out_ID",$query);
658
 
           last SWITCH1;
659
 
       };
660
 
      (($type_in =~ /^unaligned /) && defined($read_in) &&
661
 
       ($type_out eq 'ungapped consensus')) && do {
662
 
           $query = $self->change_coord("unaligned $in_ID",'gapped consensus',$query);
663
 
           $query = $self->change_coord('gapped consensus','ungapped consensus',$query);
664
 
           last SWITCH1;
665
 
       };
666
 
 
667
 
      $self->throw("Unknow coordinate system. Args: $type_in, $type_out.");
668
 
      $query = undef; # If a coordinate systems just requested is unknown
669
 
  }
 
579
    SWITCH1: {
 
580
 
 
581
        # Transformations between contig padded and contig unpadded
 
582
        (($type_in eq 'gapped consensus') && ($type_out eq 'ungapped consensus')) && do {
 
583
            $self->throw("Can't use ungapped consensus coordinates without a consensus sequence")
 
584
                unless (defined $self->{'_consensus_sequence'});
 
585
            $query = &_padded_unpadded($self->{'_consensus_gaps'}, $query);
 
586
            last SWITCH1;
 
587
        };
 
588
        (($type_in eq 'ungapped consensus') && ($type_out eq 'gapped consensus')) && do {
 
589
            $self->throw("Can't use ungapped consensus coordinates without a consensus sequence")
 
590
                unless (defined $self->{'_consensus_sequence'});
 
591
            $query = &_unpadded_padded($self->{'_consensus_gaps'},$query);
 
592
            last SWITCH1;
 
593
        };
 
594
 
 
595
        # Transformations between contig (padded) and read (padded)
 
596
        (($type_in  eq 'gapped consensus') &&
 
597
        ($type_out =~ /^aligned /) && defined($read_out)) && do {
 
598
            $query = $query - $read_out->start() + 1;
 
599
            last SWITCH1;
 
600
        };
 
601
        (($type_in =~ /^aligned /) && defined($read_in) &&
 
602
        ($type_out  eq 'gapped consensus')) && do {
 
603
            $query = $query + $read_in->start() - 1;
 
604
            last SWITCH1;
 
605
        };
 
606
 
 
607
        # Transformations between contig (unpadded) and read (padded)
 
608
        (($type_in eq 'ungapped consensus') &&
 
609
        ($type_out =~ /^aligned /) && defined($read_out)) && do {
 
610
            $query = $self->change_coord('ungapped consensus','gapped consensus',$query);
 
611
            $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
 
612
            last SWITCH1;
 
613
        };
 
614
        (($type_in =~ /^aligned /) && defined($read_in) &&
 
615
        ($type_out eq 'ungapped consensus')) && do {
 
616
            $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
 
617
            $query = $self->change_coord('gapped consensus','ungapped consensus',$query);
 
618
            last SWITCH1;
 
619
        };
 
620
 
 
621
        # Transformations between seq $read_in padded and seq $read_out padded
 
622
        (defined($read_in)  && ($type_in  =~ /^aligned /)  &&
 
623
        defined($read_out) && ($type_out =~ /^aligned /)) && do {
 
624
            $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
 
625
            $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
 
626
            last SWITCH1;
 
627
        };
 
628
 
 
629
        # Transformations between seq $read_in padded and seq $read_out unpadded
 
630
        (defined($read_in)  && ($type_in  =~ /^aligned /)    &&
 
631
        defined($read_out) && ($type_out =~ /^unaligned /)) && do {
 
632
            if ($read_in ne $read_out) {
 
633
                $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
 
634
                $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
 
635
            }
 
636
            my $list_out = $self->{'_elem'}{$out_ID}{'_gaps'};
 
637
            $query = &_padded_unpadded($list_out,$query);
 
638
            # Changing read orientation if read was reverse complemented when aligned
 
639
            if ($read_out->strand == -1) {
 
640
                my ($length) = $read_out->length();
 
641
                $length = $length - &_nof_gaps($list_out,$length);
 
642
                $query  = $length - $query + 1;
 
643
            }
 
644
            last SWITCH1;
 
645
        };
 
646
        (defined($read_in)  && ($type_in  =~ /^unaligned /) &&
 
647
        defined($read_out) && ($type_out =~ /^aligned /))  && do {
 
648
            my $list_in = $self->{'_elem'}{$in_ID}{'_gaps'};
 
649
            # Changing read orientation if read was reverse complemented when aligned
 
650
            if ($read_in->strand == -1) {
 
651
                my ($length) = $read_in->length();
 
652
                $length = $length - &_nof_gaps($list_in,$length);
 
653
                $query  = $length - $query + 1;
 
654
            }
 
655
            $query = &_unpadded_padded($list_in,$query);
 
656
            if ($read_in ne $read_out) {
 
657
                $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
 
658
                $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
 
659
            }
 
660
            last SWITCH1;
 
661
        };
 
662
 
 
663
        # Transformations between seq $read_in unpadded and seq $read_out unpadded
 
664
        (defined($read_in)  && ($type_in  =~ /^unaligned /)    &&
 
665
        defined($read_out) && ($type_out =~ /^unaligned /)) && do {
 
666
            $query = $self->change_coord("unaligned $in_ID","aligned $out_ID",$query);
 
667
            $query = $self->change_coord("aligned $out_ID","unaligned $out_ID",$query);
 
668
            last SWITCH1;
 
669
        };
 
670
 
 
671
        # Transformations between contig (padded) and read (unpadded)
 
672
        (($type_in eq 'gapped consensus') &&
 
673
        ($type_out =~ /^unaligned /) && defined($read_out)) && do {
 
674
            $query = $self->change_coord('gapped consensus',"aligned $out_ID",$query);
 
675
            $query = $self->change_coord("aligned $out_ID","unaligned $out_ID",$query);
 
676
            last SWITCH1;
 
677
        };
 
678
        (($type_in =~ /^unaligned /) && defined($read_in) &&
 
679
        ($type_out eq 'gapped consensus')) && do {
 
680
            $query = $self->change_coord("unaligned $in_ID","aligned $in_ID",$query);
 
681
            $query = $self->change_coord("aligned $in_ID",'gapped consensus',$query);
 
682
            last SWITCH1;
 
683
        };
 
684
 
 
685
        # Transformations between contig (unpadded) and read (unpadded)
 
686
        (($type_in eq 'ungapped consensus') &&
 
687
        ($type_out =~ /^unaligned /) && defined($read_out)) && do {
 
688
            $query = $self->change_coord('ungapped consensus','gapped consensus',$query);
 
689
            $query = $self->change_coord('gapped consensus',"unaligned $out_ID",$query);
 
690
            last SWITCH1;
 
691
        };
 
692
        (($type_in =~ /^unaligned /) && defined($read_in) &&
 
693
        ($type_out eq 'ungapped consensus')) && do {
 
694
            $query = $self->change_coord("unaligned $in_ID",'gapped consensus',$query);
 
695
            $query = $self->change_coord('gapped consensus','ungapped consensus',$query);
 
696
            last SWITCH1;
 
697
        };
 
698
 
 
699
        $self->throw("Unknow coordinate system. Args: $type_in, $type_out.");
 
700
        $query = undef; # If a coordinate systems just requested is unknown
 
701
    }
670
702
 
671
703
    return $query;
672
704
}
686
718
    my ($self,$seq) = @_;
687
719
 
688
720
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
689
 
        $self->throw("$seq is not a Bio::LocatableSeq");
 
721
        $self->throw("$seq is not a Bio::LocatableSeq");
690
722
    }
691
723
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
692
724
 
693
725
    unless (exists( $self->{'_elem'}{$seqID} )) {
694
 
        $self->warn("No such sequence ($seqID) in contig ".$self->id);
695
 
        return;
 
726
        $self->warn("No such sequence ($seqID) in contig ".$self->id);
 
727
        return;
696
728
    }
697
729
    unless (exists( $self->{'_elem'}{$seqID}{'_feat'}{"_aligned_coord:$seqID"} )) {
698
 
        # $self->warn("Chad. Location not set for sequence ($seqID) in contig ".$self->id);
699
 
        return;
 
730
        # $self->warn("Chad. Location not set for sequence ($seqID) in contig ".$self->id);
 
731
        return;
700
732
    }
701
733
 
702
734
    return $self->{'_elem'}{$seqID}{'_feat'}{"_aligned_coord:$seqID"};
732
764
    my ($self,$feat,$seq) = @_;
733
765
 
734
766
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
735
 
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
 
767
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
736
768
    }
737
769
 
738
770
    # Complaining about inadequate feature object
739
 
     $self->throw("Coordinates must be a Bio::SeqFeature::Generic object!")
740
 
        unless ( $feat->isa("Bio::SeqFeature::Generic") );
 
771
    $self->throw("Coordinates must be a Bio::SeqFeature::Generic object!")
 
772
        unless ( $feat->isa("Bio::SeqFeature::Generic") );
741
773
    $self->throw("Sequence coordinates must have an end!")
742
 
        unless (defined $feat->end);
 
774
        unless (defined $feat->end);
743
775
    $self->throw("Sequence coordinates must have a start!")
744
 
        unless (defined $feat->start);
 
776
        unless (defined $feat->start);
745
777
 
746
778
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
747
779
    if (exists( $self->{'_elem'}{$seqID} ) &&
748
 
        exists( $self->{'_elem'}{$seqID}{'_seq'} ) &&
749
 
        defined( $self->{'_elem'}{$seqID}{'_seq'} ) &&
750
 
        ($seq ne $self->{'_elem'}{$seqID}{'_seq'}) ) {
751
 
        $self->warn("Replacing sequence $seqID\n");
752
 
        $self->remove_seq($self->{'_elem'}{$seqID}{'_seq'});
 
780
    exists( $self->{'_elem'}{$seqID}{'_seq'} ) &&
 
781
    defined( $self->{'_elem'}{$seqID}{'_seq'} ) &&
 
782
    ($seq ne $self->{'_elem'}{$seqID}{'_seq'}) ) {
 
783
        $self->warn("Replacing sequence $seqID\n");
 
784
        $self->remove_seq($self->{'_elem'}{$seqID}{'_seq'});
753
785
    }
754
786
    $self->add_seq($seq);
755
787
 
758
790
 
759
791
    # Add new Bio::Generic::SeqFeature
760
792
    $feat->add_tag_value('contig',$self->id)
761
 
        unless ( $feat->has_tag('contig') );
 
793
        unless ( $feat->has_tag('contig') );
762
794
    $feat->primary_tag("_aligned_coord:$seqID");
763
795
    $feat->attach_seq($seq);
764
796
    $self->{'_elem'}{$seqID}{'_feat'}{"_aligned_coord:$seqID"} = $feat;
782
814
    my $seq  = shift;
783
815
 
784
816
    $self->throw("Consensus sequence must be a Bio::LocatableSeq!")
785
 
        unless ($seq->isa("Bio::LocatableSeq"));
786
 
 
787
 
    my $con_len = $seq->length;
788
 
    $seq->start(1); $seq->end($con_len);
 
817
        unless ($seq->isa("Bio::LocatableSeq"));
789
818
 
790
819
    $self->{'_consensus_gaps'} = []; # Consensus Gap registry
791
 
    $self->_register_gaps($seq->seq,
792
 
                          $self->{'_consensus_gaps'});
 
820
    $self->_register_gaps( $seq->seq, $self->{'_consensus_gaps'} );
793
821
    $self->{'_consensus_sequence'} = $seq;
794
822
 
 
823
    $seq->start(1);
 
824
    $seq->end($seq->_ungapped_len);
 
825
 
 
826
    my $con_len = $seq->length;
 
827
 
795
828
    return $con_len;
796
829
}
797
830
 
810
843
    my $qual  = shift;
811
844
 
812
845
    $self->throw("Consensus quality must be a Bio::Seq::Quality object!")
813
 
        unless ( $qual->isa("Bio::Seq::Quality") );
 
846
        unless ( $qual->isa("Bio::Seq::Quality") );
814
847
 
815
848
    $self->throw("Consensus quality can't be added before you set the consensus sequence!")
816
 
        unless (defined $self->{'_consensus_sequence'});
 
849
        unless (defined $self->{'_consensus_sequence'});
817
850
 
818
851
    $self->{'_consensus_quality'} = $qual;
819
852
}
887
920
    my ($self,$seq,$qual) = @_;
888
921
 
889
922
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
890
 
        $self->throw("Unable to process non locatable sequences [", ref($seq), "]");
 
923
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
891
924
    }
892
925
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
893
926
 
894
927
    $self->throw("Consensus quality must be a Bio::Seq::QualI object!")
895
 
        unless ( $qual->isa("Bio::Seq::QualI") );
 
928
        unless ( $qual->isa("Bio::Seq::QualI") );
896
929
    $self->throw("Use add_seq first: aligned sequence qualities can't be added before you load the sequence!")
897
 
        unless (exists $self->{'_elem'}{$seqID}{'_seq'});
 
930
        unless (exists $self->{'_elem'}{$seqID}{'_seq'});
898
931
    $self->throw("Use set_seq_coord first: aligned sequence qualities can't be added before you add coordinates for the sequence!") unless (defined( $self->get_seq_coord($seq) ));
899
932
 
900
933
    # Adding gaps to quality object
906
939
    my $next     = 0;
907
940
    my $i = 0; my $j = 0;
908
941
    while ($i<=$#{$tmp}) {
909
 
        # IF base is a gap, quality is the average for neighbouring sites
910
 
        if (substr($sequence,$j,1) eq '-') {
911
 
            $previous = $tmp->[$i-1] unless ($i == 0);
912
 
            if ($i < $#{$tmp}) {
913
 
                $next = $tmp->[$i+1];
914
 
            } else {
915
 
                $next = 0;
916
 
            }
917
 
            push(@quality,int( ($previous+$next)/2 ));
918
 
        } else {
919
 
            push(@quality,$tmp->[$i]);
920
 
            $i++;
921
 
        }
922
 
        $j++;
 
942
        # IF base is a gap, quality is the average for neighbouring sites
 
943
        if (substr($sequence,$j,1) eq '-') {
 
944
            $previous = $tmp->[$i-1] unless ($i == 0);
 
945
            if ($i < $#{$tmp}) {
 
946
                $next = $tmp->[$i+1];
 
947
            } else {
 
948
                $next = 0;
 
949
            }
 
950
            push(@quality,int( ($previous+$next)/2 ));
 
951
        } else {
 
952
            push(@quality,$tmp->[$i]);
 
953
            $i++;
 
954
        }
 
955
        $j++;
923
956
    }
924
957
 
925
 
    $self->{'_elem'}{$seqID}{'_qual'} = Bio::Seq::PrimaryQual->new(-qual=>join(" ",@quality),
926
 
                                                                   -id=>$seqID);
 
958
    $self->{'_elem'}{$seqID}{'_qual'} = Bio::Seq::PrimaryQual->new(
 
959
        -qual=>join(" ",@quality), -id=>$seqID );
927
960
}
928
961
 
929
962
=head2 get_seq_ids
930
963
 
931
964
 Title     : get_seq_ids
932
965
 Usage     : $contig->get_seq_ids(-start=>$start,
933
 
                                  -end=>$end,
934
 
                                  -type=>"gapped A0QR67B08.b");
935
 
 Function  : Get list of sequence IDs overlapping inteval [$start, $end]
 
966
                  -end=>$end,
 
967
                  -type=>"gapped A0QR67B08.b");
 
968
 Function  : Get list of sequence IDs overlapping interval [$start, $end]
936
969
             The default interval is [1,$contig->length]
937
970
             Default coordinate system is "gapped contig"
938
971
 Returns   : An array
953
986
    my ($self, @args) = @_;
954
987
 
955
988
    my ($type,$start,$end) =
956
 
        $self->_rearrange([qw(TYPE START END)], @args);
 
989
    $self->_rearrange([qw(TYPE START END)], @args);
957
990
 
958
991
    if (defined($start) && defined($end)) {
959
 
        if (defined($type) && ($type ne 'gapped consensus')) {
960
 
            $start = $self->change_coord($type,'gapped consensus',$start);
961
 
            $end   = $self->change_coord($type,'gapped consensus',$end);
962
 
        }
 
992
        if (defined($type) && ($type ne 'gapped consensus')) {
 
993
            $start = $self->change_coord($type,'gapped consensus',$start);
 
994
            $end   = $self->change_coord($type,'gapped consensus',$end);
 
995
        }
963
996
 
964
 
        my @list = grep { $_->isa("Bio::SeqFeature::Generic") &&
965
 
                              ($_->primary_tag =~ /^_aligned_coord:/) }
966
 
        $self->{'_sfc'}->features_in_range(-start=>$start,
967
 
                                           -end=>$end,
968
 
                                           -contain=>0,
969
 
                                           -strandmatch=>'ignore');
970
 
        @list = map { $_->entire_seq->id } @list;
971
 
        return @list;
 
997
        my @list = grep { $_->isa("Bio::SeqFeature::Generic") &&
 
998
        ($_->primary_tag =~ /^_aligned_coord:/) }
 
999
        $self->{'_sfc'}->features_in_range( -start=>$start,
 
1000
                                            -end=>$end,
 
1001
                                            -contain=>0,
 
1002
                                            -strandmatch=>'ignore' );
 
1003
        @list = map { $_->entire_seq->id } @list;
 
1004
        return @list;
972
1005
    }
973
1006
 
974
1007
    # Entire aligned sequences list
993
1026
    my ($self,$seq,$tag) = @_;
994
1027
 
995
1028
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
996
 
        $self->throw("Unable to process non locatable sequences [", ref($seq), "]");
 
1029
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
997
1030
    }
998
1031
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
999
1032
 
1016
1049
    my ($seqID) = @_;
1017
1050
 
1018
1051
    unless (exists $self->{'_elem'}{$seqID}{'_seq'}) {
1019
 
        $self->throw("Could not find sequence $seqID in contig ".$self->id);
1020
 
        return;
 
1052
        $self->throw("Could not find sequence $seqID in contig ".$self->id);
 
1053
    return;
1021
1054
    }
1022
1055
 
1023
1056
    return $self->{'_elem'}{$seqID}{'_seq'};
1043
1076
    my ($seqID) = @_;
1044
1077
 
1045
1078
    unless (exists $self->{'_elem'}{$seqID}{'_qual'}) {
1046
 
        $self->warn("Could not find quality for $seqID in contig!");
1047
 
        return;
 
1079
        $self->warn("Could not find quality for $seqID in contig!");
 
1080
        return;
1048
1081
    }
1049
1082
 
1050
1083
    return $self->{'_elem'}{$seqID}{'_qual'};
1079
1112
    my $seq = shift;
1080
1113
 
1081
1114
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
1082
 
        $self->throw("Unable to process non locatable sequences [", ref($seq), "]");
 
1115
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
1083
1116
    }
1084
1117
 
1085
1118
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
1086
 
    $self->{'_elem'}{$seqID} = {} unless (exists $self->{'elem'}{$seqID});
 
1119
    $self->{'_elem'}{$seqID} = {} unless (exists $self->{'_elem'}{$seqID});
1087
1120
 
1088
1121
    if (exists( $self->{'_elem'}{$seqID}{'_seq'} ) &&
1089
 
        ($seq eq $self->{'_elem'}{$seqID}{'_seq'}) ) {
1090
 
        $self->warn("Adding sequence $seqID, which has already been added");
 
1122
    ($seq eq $self->{'_elem'}{$seqID}{'_seq'}) ) {
 
1123
        $self->warn("Adding sequence $seqID, which has already been added");
1091
1124
    }
1092
1125
 
1093
1126
    # Our locatable sequences are always considered to be complete sequences
1094
 
    $seq->start(1); $seq->end($seq->length());
 
1127
    $seq->start(1);
 
1128
    $seq->end($seq->_ungapped_len);
1095
1129
 
1096
1130
    $self->warn("Adding non-nucleotidic sequence ".$seqID)
1097
 
        if (lc($seq->alphabet) ne 'dna' && lc($seq->alphabet) ne 'rna');
 
1131
        if (lc($seq->alphabet) ne 'dna' && lc($seq->alphabet) ne 'rna');
1098
1132
 
1099
1133
    # build the symbol list for this sequence,
1100
1134
    # will prune out the gap and missing/match chars
1101
1135
    # when actually asked for the symbol list in the
1102
1136
    # symbol_chars
1103
1137
    if (defined $seq->seq) {
1104
 
        map { $self->{'_symbols'}->{$_} = 1; } split(//,$seq->seq);
 
1138
        map { $self->{'_symbols'}->{$_} = 1; } split(//,$seq->seq);
1105
1139
    } else {
1106
 
        $self->{'_symbols'} = {};
 
1140
        $self->{'_symbols'} = {};
1107
1141
    }
1108
1142
 
1109
1143
    my $seq_no = ++$self->{'_nof_seqs'};
1110
1144
 
1111
1145
    if (ref( $self->{'_elem'}{$seqID}{'_seq'} )) {
1112
 
        $self->warn("Replacing one sequence [$seqID]\n");
 
1146
        $self->warn("Replacing one sequence [$seqID]\n");
1113
1147
    } else {
1114
 
        #print STDERR "Assigning $seqID to $order\n";
1115
 
        $self->{'_order'}->{$seq_no} = $seqID;
1116
 
#       $self->{'_start_end_lists'}->{$id} = []
1117
 
#           unless(exists $self->{'_start_end_lists'}->{$id});
1118
 
#       push @{$self->{'_start_end_lists'}->{$id}}, $seq;
 
1148
        #print STDERR "Assigning $seqID to $order\n";
 
1149
        $self->{'_order'}->{$seq_no} = $seqID;
 
1150
        # $self->{'_start_end_lists'}->{$id} = []
 
1151
        # unless(exists $self->{'_start_end_lists'}->{$id});
 
1152
        # push @{$self->{'_start_end_lists'}->{$id}}, $seq;
1119
1153
    }
1120
1154
 
1121
1155
    $self->{'_elem'}{$seqID}{'_seq'}  = $seq;
1134
1168
 
1135
1169
 Title     : remove_seq
1136
1170
 Usage     : $contig->remove_seq($seq);
1137
 
 Function  : Removes a single sequence from an alignment
 
1171
 Function  : Removes a single sequence from a contig
1138
1172
 Returns   : 1 on success, 0 otherwise
1139
1173
 Argument  : a Bio::LocatableSeq object
1140
1174
 
1144
1178
    my ($self,$seq) = @_;
1145
1179
 
1146
1180
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
1147
 
        $self->throw("Unable to process non locatable sequences [", ref($seq), "]");
 
1181
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
1148
1182
    }
1149
1183
 
1150
1184
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
1151
1185
    unless (exists $self->{'_elem'}{$seqID} ) {
1152
 
        $self->warn("No sequence named $seqID  [$seq]");
1153
 
        return 0;
 
1186
        $self->warn("No sequence named $seqID  [$seq]");
 
1187
        return 0;
1154
1188
    }
1155
1189
 
1156
1190
    # Updating residue count
1157
1191
    $self->{'_nof_residues'} -= $seq->length() +
1158
 
        &_nof_gaps( $self->{'_elem'}{$seqID}{'_gaps'}, $seq->length );
 
1192
    &_nof_gaps( $self->{'_elem'}{$seqID}{'_gaps'}, $seq->length );
 
1193
    
 
1194
    # Update number of sequences
 
1195
    $self->{'_nof_seqs'}--; 
 
1196
    
 
1197
    # Update order of sequences (order starts at 1)
 
1198
    my $max_order = $self->{'_nof_seqs'} + 1;
 
1199
    my $target_order = $max_order + 1;
 
1200
    for (my $order = 1 ; $order <= $max_order ; $order++) {
 
1201
      if ($self->{'_order'}->{$order} eq $seqID) {
 
1202
        # Found the wanted sequence order
 
1203
        $target_order = $order;
 
1204
      }
 
1205
      if ($order > $target_order) {
 
1206
        # Decrement this sequence order by one order
 
1207
        $self->{'_order'}->{$order-1} = $self->{'_order'}->{$order};
 
1208
      }
 
1209
      if ($order == $max_order) {
 
1210
        # Remove last order
 
1211
        delete $self->{'_order'}->{$order};
 
1212
      }
 
1213
    }
1159
1214
 
1160
1215
    # Remove all references to features of this sequence
1161
1216
    my @feats = ();
1162
 
    foreach my $tag (keys %{ $self->{'_elem'}{$seqID}{'_feat'} }) {
1163
 
        push(@feats, $self->{'_elem'}{$seqID}{'_feat'}{$tag});
 
1217
    for my $tag (keys %{ $self->{'_elem'}{$seqID}{'_feat'} }) {
 
1218
        push(@feats, $self->{'_elem'}{$seqID}{'_feat'}{$tag});
1164
1219
    }
1165
1220
    $self->{'_sfc'}->remove_features(\@feats);
1166
1221
    delete $self->{'_elem'}{$seqID};
1230
1285
    my (@arr,$seqID);
1231
1286
 
1232
1287
    foreach $seqID ( map { $self->{'_order'}{$_} } sort { $a <=> $b } keys %{$self->{'_order'}} ) {
1233
 
        push(@arr,$self->{'_elem'}{$seqID}{'_seq'});
 
1288
        push(@arr,$self->{'_elem'}{$seqID}{'_seq'});
1234
1289
    }
1235
1290
 
1236
1291
    return @arr;
1296
1351
    my ($pos) = @_;
1297
1352
 
1298
1353
    $self->throw("Sequence position has to be a positive integer, not [$pos]")
1299
 
        unless $pos =~ /^\d+$/ and $pos > 0;
 
1354
        unless $pos =~ /^\d+$/ and $pos > 0;
1300
1355
    $self->throw("No sequence at position [$pos]")
1301
 
        unless $pos <= $self->no_sequences ;
 
1356
        unless $pos <= $self->no_sequences ;
1302
1357
 
1303
1358
    my $seqID = $self->{'_order'}->{--$pos};
1304
1359
    return $self->{'_elem'}{$seqID}{'_seq'};
1502
1557
    my ($self, $contig_name) = @_;
1503
1558
 
1504
1559
    if (defined( $contig_name )) {
1505
 
        $self->{'_id'} = $contig_name;
 
1560
        $self->{'_id'} = $contig_name;
1506
1561
    }
1507
1562
 
1508
1563
    return $self->{'_id'};
1658
1713
    $self->throw_not_implemented();
1659
1714
}
1660
1715
 
1661
 
=head2 maxdisplayname_length
 
1716
=head2 maxdname_length
1662
1717
 
1663
 
 Title     : maxdisplayname_length
1664
 
 Usage     : $contig->maxdisplayname_length()
 
1718
 Title     : maxname_length
 
1719
 Usage     : $contig->maxname_length()
1665
1720
 Function  :
1666
1721
 
1667
1722
             Gets the maximum length of the displayname in the
1715
1770
 Usage   : $id = $contig->percentage_identity
1716
1771
 Function: The function calculates the percentage identity of the alignment
1717
1772
 Returns : The percentage identity of the alignment (as defined by the
1718
 
                                                     implementation)
 
1773
                             implementation)
1719
1774
 Argument: None
1720
1775
 
1721
1776
=cut
1782
1837
           sequence with the given name. For example, for the
1783
1838
           alignment
1784
1839
 
1785
 
             Seq1/91-97 AC..DEF.GH
1786
 
             Seq2/24-30 ACGG.RTY..
1787
 
             Seq3/43-51 AC.DDEFGHI
 
1840
           Seq1/91-97 AC..DEF.GH
 
1841
           Seq2/24-30 ACGG.RTY..
 
1842
           Seq3/43-51 AC.DDEFGHI
1788
1843
 
1789
1844
           column_from_residue_number( "Seq1", 94 ) returns 5.
1790
1845
           column_from_residue_number( "Seq2", 25 ) returns 2.
1794
1849
           outside the length of the aligment
1795
1850
           (e.g. column_from_residue_number( "Seq2", 22 )
1796
1851
 
1797
 
          Note: If the the parent sequence is represented by more than
1798
 
          one alignment sequence and the residue number is present in
1799
 
          them, this method finds only the first one.
 
1852
      Note: If the the parent sequence is represented by more than
 
1853
      one alignment sequence and the residue number is present in
 
1854
      them, this method finds only the first one.
1800
1855
 
1801
1856
 Returns : A column number for the position in the alignment of the
1802
1857
           given residue in the given sequence (1 = first column)
1920
1975
    my $middle = 0;
1921
1976
    while ($end - $start > 1) {
1922
1977
        $middle = int(($end+$middle)/2);
1923
 
        (&_compare($query,$list->[$middle]) == 0) && return ($middle);
1924
 
        (&_compare($query,$list->[$middle]) <  0) && do { $end   = $middle ; $middle = 0; next };
1925
 
        $start = $middle; # If &_compare() > 0, move region beggining
 
1978
        (&_compare($query,$list->[$middle]) == 0) && return ($middle);
 
1979
        (&_compare($query,$list->[$middle]) <  0) && do { $end   = $middle ; $middle = 0; next };
 
1980
        $start = $middle; # If &_compare() > 0, move region beggining
1926
1981
    }
1927
1982
    return ($start,$end);
1928
1983
}
1970
2025
    # If before any alignment, return 0
1971
2026
    elsif ($index[0] == -1) { $query = 0 }
1972
2027
    elsif ($index[0] >= 0) {
1973
 
        # If query is between alignments, translate coordinates
1974
 
        if ($#index > 0) { $query = $index[0] + 1 }
1975
 
        # If query sits upon an alignment, do another correction
1976
 
        elsif ($#index == 0) { $query = $index[0] }
 
2028
    # If query is between alignments, translate coordinates
 
2029
    if ($#index > 0) { $query = $index[0] + 1 }
 
2030
    # If query sits upon an alignment, do another correction
 
2031
    elsif ($#index == 0) { $query = $index[0] }
1977
2032
    }
1978
2033
    #
1979
2034
    return $query;
2032
2087
    $query = $query + $align;
2033
2088
    my $new_align = &_nof_gaps($list,$query);
2034
2089
    while ($new_align - $align > 0) {
2035
 
        $query = $query + $new_align - $align;
2036
 
        $align  = $new_align;
2037
 
        $new_align = &_nof_gaps($list,$query);
 
2090
        $query = $query + $new_align - $align;
 
2091
        $align  = $new_align;
 
2092
        $new_align = &_nof_gaps($list,$query);
2038
2093
    }
2039
2094
    # If current position is also a align, look for the first upstream base
2040
2095
    while (defined($list->[$align]) && ($list->[$align] == $query)) {
2041
 
        $query++; $align++;
 
2096
        $query++; $align++;
2042
2097
    }
2043
2098
    #
2044
2099
    return $query;
2064
2119
    my $dbref    = shift;
2065
2120
 
2066
2121
    $self->throw("Not an aligned sequence string to register gaps")
2067
 
        if (ref($sequence));
 
2122
        if (ref($sequence));
2068
2123
 
2069
2124
    $self->throw("Not an array reference for gap registry")
2070
 
        unless (ref($dbref) eq 'ARRAY');
 
2125
        unless (ref($dbref) eq 'ARRAY');
2071
2126
 
2072
2127
    # Registering alignments
2073
2128
    @{$dbref} = (); # Cleaning registry
2074
2129
    if (defined $sequence) {
2075
 
        my $i = -1;
2076
 
        while(1) {
2077
 
            $i = index($sequence,"-",$i+1);
2078
 
            last if ($i == -1);
2079
 
            push(@{$dbref},$i+1);
2080
 
        }
 
2130
        my $i = -1;
 
2131
        while(1) {
 
2132
            $i = index($sequence,"-",$i+1);
 
2133
            last if ($i == -1);
 
2134
            push(@{$dbref},$i+1);
 
2135
        }
2081
2136
    } else {
2082
 
#       $self->warn("Found undefined sequence while registering gaps");
2083
 
        return 0;
 
2137
        # $self->warn("Found undefined sequence while registering gaps");
 
2138
        return 0;
2084
2139
    }
2085
2140
 
2086
2141
    return scalar(@{$dbref});
2087
2142
}
2088
2143
 
 
2144
 
2089
2145
1;