← 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: 2011-06-17 13:51:18 UTC
  • mfrom: (1.2.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20110617135118-xgpxhaanue57cwqs
Tags: 1.6.901-1
* New upstream release.
* Point debian/watch to search.cpan.org.
* Build using dh and overrides:
  - Use Debhelper 8 (debian/rules, debian/control).
  - Simplified debian/rules.
* Split into libbio-perl-perl, as discussed with the Debian Perl team.
  (debian/control, debian/bioperl.install, debian libbio-perl-perl.install)
* debian/control:
  - Incremented Standards-Version to reflect conformance with Policy 3.9.2.
    No other changes needed.
  - Vcs-Browser URL made redirectable to viewvc.
  - Removed useless ‘svn’ in the Vcs-Svn URL.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Assembly/Contig.pm
1
 
# $Id: Contig.pm 16123 2009-09-17 12:57:27Z cjfields $
2
1
#
3
2
# BioPerl module for Bio::Assembly::Contig
4
3
#   Mostly based on Bio::SimpleAlign by Ewan Birney
5
4
#
6
5
# Please direct questions and support issues to <bioperl-l@bioperl.org> 
7
6
#
8
 
# Cared for by Robson francisco de Souza <rfsouza@citri.iq.usp.br>
 
7
# Cared for by Robson Francisco de Souza <rfsouza@citri.iq.usp.br>
9
8
#
10
9
# Copyright Robson Francisco de Souza
11
10
#
25
24
 
26
25
    # Assembly loading methods
27
26
    $aio = Bio::Assembly::IO->new(-file=>"test.ace.1",
28
 
                               -format=>'phrap');
 
27
                                  -format=>'phrap');
29
28
 
30
29
    $assembly = $aio->next_assembly;
31
30
    foreach $contig ($assembly->all_contigs) {
202
201
the bugs and their resolution.  Bug reports can be submitted via the
203
202
web:
204
203
 
205
 
  http://bugzilla.open-bio.org/
 
204
  https://redmine.open-bio.org/projects/bioperl/
206
205
 
207
206
=head1 AUTHOR - Robson Francisco de Souza
208
207
 
220
219
 
221
220
use strict;
222
221
 
223
 
use Bio::SeqFeature::Collection;
224
 
use Bio::Seq::PrimaryQual;
 
222
use Bio::DB::SeqFeature::Store; # isa Bio::SeqFeature::CollectionI
 
223
use Bio::Seq::PrimaryQual;      # isa Bio::Seq::QualI
225
224
 
226
225
use Scalar::Util qw(weaken);
227
226
 
235
234
 Usage     : my $contig = Bio::Assembly::Contig->new();
236
235
 Function  : Creates a new contig object
237
236
 Returns   : Bio::Assembly::Contig
238
 
 Args      : -id => contig unique ID
239
 
             -source => string for the sequence assembly program used
240
 
             -collection => Bio::SeqFeature::Collection instance
 
237
 Args      : -id         => unique contig ID
 
238
             -source     => string for the sequence assembly program used
 
239
             -collection => Bio::SeqFeature::CollectionI instance
241
240
 
242
241
=cut
243
242
 
273
272
        $self->throw("Collection must implement Bio::SeqFeature::CollectionI") unless $collection->isa('Bio::SeqFeature::CollectionI');
274
273
        $self->{'_sfc'} = $collection;
275
274
    } else {
276
 
        $self->{'_sfc'} = Bio::SeqFeature::Collection->new()
 
275
        $self->{'_sfc'} = Bio::DB::SeqFeature::Store->new(
 
276
            -adaptor           => 'memory',
 
277
            -index_subfeatures => 1,
 
278
        );
277
279
    }
278
280
 
279
281
    # Assembly specifics
280
 
    $self->{'_assembly'} = undef; # Reference to a Bio::Assembly::Scaffold object, if contig belongs to one.
 
282
    $self->{'_assembly'} = undef; # Bio::Assembly::Scaffold the contig belongs to
281
283
    $self->{'_strand'} = 0; # Reverse (-1) or forward (1), if contig is in a scaffold. 0 otherwise
282
 
    $self->{'_neighbor_start'} = undef; # Will hold a reference to another contig
283
 
    $self->{'_neighbor_end'} = undef; # Will hold a reference to another contig
 
284
    $self->{'_neighbor_start'} = undef; # Neighbor Bio::Assembly::Contig
 
285
    $self->{'_neighbor_end'}   = undef; # Neighbor Bio::Assembly::Contig
284
286
 
285
287
    return $self; # success - we hope!
286
288
}
464
466
    }
465
467
 
466
468
    # Add feature to feature collection
467
 
    my $nof_added = $self->{'_sfc'}->add_features($args);
 
469
    my $nof_added = $self->get_features_collection->add_features($args);
468
470
 
469
471
    return $nof_added;
470
472
}
474
476
 Title     : remove_features
475
477
 Usage     : $contig->remove_features(@feat)
476
478
 Function  : Remove an array of contig features
477
 
 Returns   : number of features removed.
478
 
 Argument  : An array of Bio::SeqFeatureI
 
479
 Returns   : true if successful
 
480
 Argument  : An array of Bio::SeqFeature::Generic (Bio::SeqFeatureI)
479
481
 
480
482
=cut
481
483
 
483
485
    my ($self, @args) = @_;
484
486
 
485
487
    # Removing shortcuts for aligned sequence features
486
 
    foreach my $feat (@args) {
 
488
    for my $feat (@args) {
487
489
        if (my $seq = $feat->entire_seq()) {
488
 
            my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
 
490
            my $seqID = $seq->id || $seq->display_id || $seq->primary_id;
489
491
            my $tag = $feat->primary_tag;
490
492
            $tag =~ s/:$seqID$/$1/g;
491
493
            delete( $self->{'_elem'}{$seqID}{'_feat'}{$tag} )
493
495
                $self->{'_elem'}{$seqID}{'_feat'}{$tag} eq $feat);
494
496
        }
495
497
    }
496
 
    
497
 
    # Removing Bio::SeqFeature::Collection features
498
 
    return $self->{'_sfc'}->remove_features(\@args);
 
498
   
 
499
    # Removing Bio::SeqFeature objects
 
500
    return $self->get_features_collection->delete(@args);
499
501
}
500
502
 
501
503
=head2 get_features_collection
502
504
 
503
505
 Title     : get_features_collection
504
506
 Usage     : $contig->get_features_collection()
505
 
 Function  : Get the collection of all contig features
506
 
 Returns   : Bio::SeqFeature::Collection
 
507
 Function  : Get the collection of all contig features and seqfeatures
 
508
 Returns   : Bio::DB::SeqFeature::Store (Bio::SeqFeature::CollectionI)
507
509
 Argument  : none
508
510
 
509
511
=cut
723
725
 Function  : Get "gapped consensus" location for aligned sequence
724
726
 Returns   : Bio::SeqFeature::Generic for coordinates or undef.
725
727
             A warning is printed if sequence coordinates were not set.
726
 
 Argument  : Bio::LocatabaleSeq object
 
728
 Argument  : Bio::LocatableSeq object
727
729
 
728
730
=cut
729
731
 
769
771
             Note: the original feature primary tag will
770
772
                   be lost.
771
773
 
772
 
             $seq   : a Bio::LocatabaleSeq object
 
774
             $seq   : a Bio::LocatableSeq object
773
775
 
774
776
=cut
775
777
 
789
791
        unless (defined $feat->start);
790
792
 
791
793
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
792
 
    if (exists( $self->{'_elem'}{$seqID} ) &&
793
 
    exists( $self->{'_elem'}{$seqID}{'_seq'} ) &&
794
 
    defined( $self->{'_elem'}{$seqID}{'_seq'} ) &&
795
 
    ($seq ne $self->{'_elem'}{$seqID}{'_seq'}) ) {
 
794
    if ( exists( $self->{'_elem'}{$seqID} ) &&
 
795
         exists( $self->{'_elem'}{$seqID}{'_seq'} ) &&
 
796
         defined( $self->{'_elem'}{$seqID}{'_seq'} ) &&
 
797
         ($seq ne $self->{'_elem'}{$seqID}{'_seq'}) ) {
796
798
        $self->warn("Replacing sequence $seqID\n");
797
799
        $self->remove_seq($self->{'_elem'}{$seqID}{'_seq'});
 
800
        $self->remove_features($feat);
798
801
    }
 
802
 
 
803
    # Add new sequence and Bio::Generic::SeqFeature
799
804
    $self->add_seq($seq);
800
805
 
801
 
    # Remove previous coordinates, if any
802
 
    $self->remove_features($feat);
803
 
 
804
 
    # Add new Bio::Generic::SeqFeature
805
 
    $feat->add_tag_value('contig',$self->id)
806
 
        unless ( $feat->has_tag('contig') );
807
 
    $feat->primary_tag("_aligned_coord:$seqID");
 
806
    $feat->add_tag_value('contig',$self->id) unless ( $feat->has_tag('contig') );
 
807
    $feat->primary_tag("_aligned_coord");
 
808
    $feat->source_tag($seqID);
808
809
    $feat->attach_seq($seq);
 
810
 
809
811
    $self->{'_elem'}{$seqID}{'_feat'}{"_aligned_coord:$seqID"} = $feat;
810
812
    $self->add_features([ $feat ]);
811
813
}
852
854
=cut
853
855
 
854
856
sub set_consensus_quality {
855
 
    my $self = shift;
856
 
    my $qual  = shift;
 
857
    my ($self, $qual) = @_;
857
858
 
858
 
    $self->throw("Consensus quality must be a Bio::Seq::Quality object!")
859
 
        unless ( $qual->isa("Bio::Seq::Quality") );
 
859
    $self->throw("Consensus quality must be a Bio::Seq::QualI object!")
 
860
        unless ( $qual->isa("Bio::Seq::QualI") );
860
861
 
861
862
    $self->throw("Consensus quality can't be added before you set the consensus sequence!")
862
863
        unless (defined $self->{'_consensus_sequence'});
903
904
 Usage     : $contig->get_consensus_quality()
904
905
 Function  : Get a reference to the consensus quality object
905
906
             for this contig.
906
 
 Returns   : A Bio::QualI object
 
907
 Returns   : A Bio::Seq::QualI object
907
908
 Argument  : none
908
909
 
909
910
=cut
951
952
    my $previous = 0;
952
953
    my $next     = 0;
953
954
    my $i = 0; my $j = 0;
954
 
    while ($i<=$#{$tmp}) {
 
955
    while ($i <= $#{$tmp}) {
955
956
        # IF base is a gap, quality is the average for neighbouring sites
956
 
        if (substr($sequence,$j,1) eq '-') {
 
957
        if ($j > $i && substr($sequence,$j,1) eq '-') {
957
958
            $previous = $tmp->[$i-1] unless ($i == 0);
958
959
            if ($i < $#{$tmp}) {
959
960
                $next = $tmp->[$i+1];
975
976
=head2 get_seq_ids
976
977
 
977
978
 Title     : get_seq_ids
978
 
 Usage     : $contig->get_seq_ids(-start=>$start,
979
 
                  -end=>$end,
980
 
                  -type=>"gapped A0QR67B08.b");
 
979
 Usage     : $contig->get_seq_ids( -start => $start,
 
980
                                   -end   => $end,
 
981
                                   -type  => "gapped A0QR67B08.b" );
981
982
 Function  : Get list of sequence IDs overlapping interval [$start, $end]
982
983
             The default interval is [1,$contig->length]
983
984
             Default coordinate system is "gapped contig"
986
987
             -start : consensus subsequence start
987
988
             -end   : consensus subsequence end
988
989
             -type  : the coordinate system type for $start and $end arguments
989
 
                      Coordinate system avaliable are:
 
990
                      Coordinate system available are:
990
991
                      "gapped consensus"   : consensus coordinates with gaps
991
992
                      "ungapped consensus" : consensus coordinates without gaps
992
993
                      "aligned $ReadID"    : read $ReadID coordinates with gaps
998
999
sub get_seq_ids {
999
1000
    my ($self, @args) = @_;
1000
1001
 
1001
 
    my ($type,$start,$end) =
1002
 
    $self->_rearrange([qw(TYPE START END)], @args);
 
1002
    my ($type, $start, $end) = $self->_rearrange([qw(TYPE START END)], @args);
1003
1003
 
 
1004
    my @list;
1004
1005
    if (defined($start) && defined($end)) {
1005
1006
        if (defined($type) && ($type ne 'gapped consensus')) {
1006
1007
            $start = $self->change_coord($type,'gapped consensus',$start);
1007
1008
            $end   = $self->change_coord($type,'gapped consensus',$end);
1008
1009
        }
1009
 
 
1010
 
        my @list = grep { $_->isa("Bio::SeqFeature::Generic") &&
1011
 
        ($_->primary_tag =~ /^_aligned_coord:/) }
1012
 
        $self->{'_sfc'}->features_in_range( -start=>$start,
1013
 
                                            -end=>$end,
1014
 
                                            -contain=>0,
1015
 
                                            -strandmatch=>'ignore' );
 
1010
        @list = $self->get_features_collection->features(
 
1011
           -type         => '_aligned_coord', # primary tag
 
1012
           -start        => $start,
 
1013
           -end          => $end,
 
1014
           #-contain     => 0,
 
1015
           #-strandmatch => 'ignore',
 
1016
        );
1016
1017
        @list = map { $_->entire_seq->id } @list;
1017
 
        return @list;
 
1018
    } else {
 
1019
        # Entire aligned sequences list
 
1020
        @list = map { $self->{'_order'}{$_} } sort { $a cmp $b } keys %{ $self->{'_order'} };
1018
1021
    }
1019
1022
 
1020
 
    # Entire aligned sequences list
1021
 
    return map { $self->{'_order'}{$_} } sort { $a cmp $b } keys %{ $self->{'_order'} };
 
1023
    return @list;
1022
1024
}
1023
1025
 
1024
1026
=head2 get_seq_feat_by_tag
1025
1027
 
1026
1028
 Title     : get_seq_feat_by_tag
1027
1029
 Usage     : $seq = $contig->get_seq_feat_by_tag($seq,"_aligned_coord:$seqID")
1028
 
 Function  :
1029
 
 
1030
 
             Get a sequence feature based on its primary_tag.
1031
 
             When you add
1032
 
 
 
1030
 Function  : Get a sequence feature based on its primary_tag.
1033
1031
 Returns   : a Bio::SeqFeature object
1034
1032
 Argument  : a Bio::LocatableSeq and a string (feature primary tag)
1035
1033
 
1041
1039
    if( !ref $seq || ! $seq->isa('Bio::LocatableSeq') ) {
1042
1040
        $self->throw("Unable to process non locatable sequences [".ref($seq)."]");
1043
1041
    }
1044
 
    my $seqID = $seq->id() || $seq->display_id || $seq->primary_id;
 
1042
    my $seqID = $seq->id || $seq->display_id || $seq->primary_id;
1045
1043
 
1046
1044
    return $self->{'_elem'}{$seqID}{'_feat'}{$tag};
1047
1045
}
1139
1137
    # Our locatable sequences are always considered to be complete sequences
1140
1138
    $seq->start(1);
1141
1139
    $seq->end($seq->_ungapped_len);
1142
 
 
 
1140
    
 
1141
    my $alphabet = $seq->alphabet;
 
1142
    
 
1143
    $alphabet = lc($alphabet) if defined $alphabet;
 
1144
    
1143
1145
    $self->warn("Adding non-nucleotidic sequence ".$seqID)
1144
 
        if (lc($seq->alphabet) ne 'dna' && lc($seq->alphabet) ne 'rna');
 
1146
        if (!$alphabet || ($alphabet ne 'dna' && $alphabet ne 'rna'));
1145
1147
 
1146
1148
    # build the symbol list for this sequence,
1147
1149
    # will prune out the gap and missing/match chars
1407
1409
 
1408
1410
             Creates a new alignment from a subset of
1409
1411
             sequences.  Numbering starts from 1.  Sequence positions
1410
 
             larger than num_sequences() will thow an error.
 
1412
             larger than num_sequences() will throw an error.
1411
1413
 
1412
1414
 Returns   : a Bio::Assembly::Contig object
1413
1415
 Args      : array of integers for the sequences
1649
1651
 Usage     : $str = $contig->consensus_string($threshold_percent)
1650
1652
 Function  : Makes a strict consensus
1651
1653
 Returns   :
1652
 
 Argument  : Optional treshold ranging from 0 to 100.
 
1654
 Argument  : Optional threshold ranging from 0 to 100.
1653
1655
             The consensus residue has to appear at least threshold %
1654
1656
             of the sequences at a given location, otherwise a '?'
1655
1657
             character will be placed at that location.
1726
1728
    $self->throw_not_implemented();
1727
1729
}
1728
1730
 
1729
 
=head2 maxdname_length
 
1731
=head2 maxname_length
1730
1732
 
1731
1733
 Title     : maxname_length
1732
1734
 Usage     : $contig->maxname_length()
2172
2174
=cut
2173
2175
 
2174
2176
sub no_residues {
2175
 
        my $self = shift;
2176
 
        $self->deprecated(-warn_version => 1.0069,
2177
 
                                          -throw_version => 1.0075);
 
2177
    my $self = shift;
 
2178
    $self->deprecated(-warn_version  => 1.0069,
 
2179
                      -throw_version => 1.0075);
2178
2180
    $self->num_residues(@_);
2179
2181
}
2180
2182
 
2190
2192
=cut
2191
2193
 
2192
2194
sub no_sequences {
2193
 
        my $self = shift;
2194
 
        $self->deprecated(-warn_version => 1.0069,
2195
 
                                          -throw_version => 1.0075);
 
2195
    my $self = shift;
 
2196
    $self->deprecated(-warn_version => 1.0069,
 
2197
                      -throw_version => 1.0075);
2196
2198
    $self->num_sequences(@_);
2197
2199
}
2198
2200