Bio LocatableSeq
SummaryIncluded librariesPackage variablesSynopsisDescriptionGeneral documentationMethods
Toolbar
WebCvs
Summary
Bio::LocatableSeq - A Bio::PrimarySeq object with start/end points on it
that can be projected into a MSA or have coordinates relative to
another seq.
Package variables
Globals (from "use vars" definitions)
$OTHER_SYMBOLS
$GAP_SYMBOLS
$MATCHPATTERN
$RESIDUE_SYMBOLS
$FRAMESHIFT_SYMBOLS
Included modules
Bio::Location::Fuzzy
Bio::Location::Simple
Inherit
Bio::PrimarySeq Bio::RangeI
Synopsis
    use Bio::LocatableSeq;
my $seq = Bio::LocatableSeq->new(-seq => "CAGT-GGT",
-id => "seq1",
-start => 1,
-end => 7);
# a normal sequence object $locseq->seq(); $locseq->id(); # has start,end points $locseq->start(); $locseq->end(); # inherits off RangeI, so range operations possible
Description
The LocatableSeq sequence object was developed mainly because the SimpleAlign
object requires this functionality, and in the rewrite of the Sequence object we
had to decide what to do with this.
It is, to be honest, not well integrated with the rest of bioperl. For example,
the trunc() function does not return a LocatableSeq object, as some might have
thought. Also, the sequence is not a Bio::SeqI, so the location is simply
inherited from Bio::RangeI and is not stored in a Bio::Location.
There are all sorts of nasty gotcha's about interactions between coordinate
systems when these sort of objects are used. Some mapping now occurs to deal
with HSP data, however it can probably be integrated in better and most methods
do not implement it correctly yet. Also, several PrimarySeqI methods (subseq(),
trunc(), etc.) do not behave as expected and must be used with care. Due to this,
LocatableSeq functionality is to be refactored in a future BioPerl release.
However, for alignment functionality it works adequately for the time being.
If you do not need alignment functionality, Bio::SeqfeatureI-implementing
modules may be a suitable alternative to Bio::LocatableSeq. For example,
Bio::SeqFeature::Generic and Bio::SeqFeature::Lite provide methods to
attach a sequence to a specific region of a parent sequence and to set other
useful attributes.
Methods
new
No description
Code
startDescriptionCode
endDescriptionCode
_ungapped_len
No description
Code
strandDescriptionCode
mappingDescriptionCode
frameshiftsDescriptionCode
get_nseDescriptionCode
force_nseDescriptionCode
num_gapsDescriptionCode
column_from_residue_numberDescriptionCode
location_from_columnDescriptionCode
revcomDescriptionCode
truncDescriptionCode
validate_seqDescriptionCode
no_gaps
No description
Code
no_sequencesDescriptionCode
Methods description
startcode    nextTop
 Title   : start
Usage : $obj->start($newval)
Function: Get/set the 1-based start position of this sequence in the original
sequence. '0' means before the original sequence starts.
Returns : value of start
Args : newvalue (optional)
endcodeprevnextTop
 Title   : end
Usage : $obj->end($newval)
Function: Get/set the 1-based end position of this sequence in the original
sequence. '0' means before the original sequence starts.
Returns : value of end
Args : newvalue (optional)
Note : although this is a get/set, it checks passed values against the
calculated end point ( derived from the sequence and based on
$GAP_SYMBOLS and possible frameshifts() ). If there is no match,
it will warn and set the proper value. Probably best used for
debugging proper sequence calculations.
strandcodeprevnextTop
 Title   : strand
Usage : $obj->strand($newval)
Function: return or set the strandedness
Returns : the value of the strandedness (-1, 0 or 1)
Args : the value of the strandedness (-1, 0 or 1)
mappingcodeprevnextTop
 Title   : mapping
Usage : $obj->mapping($newval)
Function: return or set the mapping indices (indicates # symbols/positions in
the source string mapping to # of coordinate positions)
Returns : two-element array (# symbols => # coordinate pos)
Args : two elements (# symbols => # coordinate pos); this can also be
passed in as an array reference of the two elements (as might be
passed upon Bio::LocatableSeq instantiation, for instance).
frameshiftscodeprevnextTop
 Title   : frameshifts
Usage : $obj->frameshifts($newval)
Function: get/set the frameshift hash, which contains sequence positions as
keys and the shift (-2, -1, 1, 2) as the value
Returns : hash
Args : hash or hash reference
get_nsecodeprevnextTop
 Title   : get_nse
Usage :
Function: read-only name of form id/start-end
Example :
Returns :
Args :
force_nsecodeprevnextTop
 Title   : force_nse
Usage : $ls->force_nse()
Function: Boolean which forces get_nse() to build an NSE, regardless
of whether id(), start(), or end() is set
Returns : Boolean value
Args : (optional) Boolean (1 or 0)
Note : This will convert any passed value evaluating as TRUE/FALSE to 1/0
respectively
num_gapscodeprevnextTop
 Title   : num_gaps
Usage :$self->num_gaps('.')
Function:Gets number of gaps in the sequence. The count excludes
leading or trailing gap characters.
Valid bioperl sequence characters are [A-Za-z\-\.\*]. Of these, '.' and '-' are counted as gap characters unless an optional argument specifies one of them. Returns : number of internal gaps in the sequence. Args : a gap character (optional) Status : Stable Note : replaces no_gaps
column_from_residue_numbercodeprevnextTop
 Title   : column_from_residue_number
Usage : $col = $seq->column_from_residue_number($resnumber)
Function:
This function gives the position in the alignment (i.e. column number) of the given residue number in the sequence. For example, for the sequence Seq1/91-97 AC..DEF.GH column_from_residue_number(94) returns 6. An exception is thrown if the residue number would lie outside the length of the aligment (e.g. column_from_residue_number( "Seq2", 22 ) Returns : A column number for the position of the given residue in the given sequence (1 = first column) Args : A residue number in the whole sequence (not just that segment of it in the alignment)
location_from_columncodeprevnextTop
 Title   : location_from_column
Usage : $loc = $ali->location_from_column($column_number)
Function:
This function gives the residue number for a given position in the alignment (i.e. column number) of the given. Gaps complicate this process and force the output to be a Bio::Location::Simple where values can be undefined.
For example, for the sequence:
Seq/91-96 .AC..DEF.G. location_from_column( 3 ) position 92 location_from_column( 4 ) position 92^93 location_from_column( 9 ) position 95^96 location_from_column( 1 ) position undef An exact position returns a Bio::Location::Simple object where where location_type() returns 'EXACT', if a position is between bases location_type() returns 'IN-BETWEEN'. Column before the first residue returns undef. Note that if the position is after the last residue in the alignment, that there is no guarantee that the original sequence has residues after that position. An exception is thrown if the column number is not within the sequence. Returns : Bio::Location::Simple or undef Args : A column number Throws : If column is not within the sequence
See Bio::Location::Simple for more.
revcomcodeprevnextTop
 Title   : revcom
Usage : $rev = $seq->revcom()
Function: Produces a new Bio::LocatableSeq object which
has the reversed complement of the sequence. For protein
sequences this throws an exception of "Sequence is a
protein. Cannot revcom"
Returns : A new Bio::LocatableSeq object Args : none
trunccodeprevnextTop
 Title   : trunc
Usage : $subseq = $myseq->trunc(10,100);
Function: Provides a truncation of a sequence,
Example : Returns : a fresh Bio::PrimarySeqI implementing object Args : Two integers denoting first and last columns of the sequence to be included into sub-sequence.
validate_seqcodeprevnextTop
 Title   : validate_seq
Usage : if(! $seq->validate_seq($seq_str) ) {
print "sequence $seq_str is not valid for an object of
alphabet ",$seq->alphabet, "\n";
}
Function: Validates a given sequence string. A validating sequence string
must be accepted by seq(). A string that does not validate will
lead to an exception if passed to seq().
The implementation provided here does not take alphabet() into account. Allowed are all letters (A-Z), numbers [0-9] and common symbols used for gaps, stop codons, unknown residues, and frameshifts, including '-','.','*','?','=',and '~'. Example : Returns : 1 if the supplied sequence string is valid for the object, and 0 otherwise. Args : The sequence string to be validated.
no_sequencescodeprevnextTop
 Title     : no_sequences
Usage : $gaps = $seq->no_sequences
Function : number of sequence in the sequence alignment
Returns : integer
Argument :
Status : Deprecated (in favor of num_sequences())
Methods code
newdescriptionprevnextTop
sub new {
    my ($class, @args) = @_;
    my $self = $class->SUPER::new(@args);

    my ($start,$end,$strand, $mapping, $fs, $nse) =
    $self->_rearrange( [qw(START
                        END
                        STRAND
                        MAPPING
                        FRAMESHIFTS
                        FORCE_NSE
                        )],
               @args);
    
    $mapping ||= [1,1];
    $self->mapping($mapping);
    $nse || 0;
    $self->force_nse($nse);
    defined $fs    && $self->frameshifts($fs);
    defined $start && $self->start($start);
    defined $end   && $self->end($end);
    defined $strand && $self->strand($strand);

    return $self; # success - we hope!
}
startdescriptionprevnextTop
sub start {
    my $self = shift;
    if( @_ ) {
        my $value = shift;
        $self->{'start'} = $value;
    }
    return $self->{'start'} if defined $self->{'start'};
    return 1                if $self->seq;
    return;
}
enddescriptionprevnextTop
sub end {
    my $self = shift;
    if( @_ ) {
        my $value = shift;
        my $st = $self->start;
        # start of 0 usually means the sequence is all gaps but maps to
# other sequences in an alignment
if ($self->seq && $st != 0 ) { my $len = $self->_ungapped_len; my $calend = $st + $len - 1; my $id = $self->id || 'unknown'; if ($calend != $value) { $self->warn("In sequence $id residue count gives end value ". "$calend.\n Overriding value [$value] with value $calend for ". "Bio::LocatableSeq::end().\n".$self->seq); $value = $calend; } } $self->{'end'} = $value; } if (defined $self->{'end'}) { return $self->{'end'} } elsif ( my $len = $self->_ungapped_len) { return $len + $self->start - 1; } else { return; } } # changed 08.10.26 to return ungapped length, not the calculated end
# of the sequence
}
_ungapped_lendescriptionprevnextTop
sub _ungapped_len {
    my $self = shift;
    return unless my $string = $self->seq;
    my ($map_res, $map_coord) = $self->mapping;
    my $offset = 0;
    if (my %data = $self->frameshifts) {
        map {$offset += $_} values %data;
    }
    $string =~ s{[$GAP_SYMBOLS$FRAMESHIFT_SYMBOLS]+}{}g;
    return CORE::length($string)/($map_res/$map_coord) + $offset/($map_coord/$map_res);
}

#sub length {
# my $self = shift;
# return unless my $string = $self->seq;
# $string =~ s{[$GAP_SYMBOLS$FRAMESHIFT_SYMBOLS]+}{}g;
# return CORE::length($string);
#
}
stranddescriptionprevnextTop
sub strand {
   my $self = shift;
   if( @_ ) {
        my $value = shift;
        $self->{'strand'} = $value;
    }
    return $self->{'strand'};
}
mappingdescriptionprevnextTop
sub mapping {
    my $self = shift;
    if( @_ ) {
        my @mapping = (ref($_[0]) eq 'ARRAY') ? @{$_[0]} : @_;
        $self->throw("Must pass two values (# residues mapped to # positions)")
            if @mapping != 2;
        if ((grep {$_ != 1 && $_ != 3} @mapping) || ($mapping[0] == 3 && $mapping[1] == 3)) {
            $self->throw("Mapping values other than 1 or 3 are not currently supported")
        }
        $self->{'_mapping'} =\@ mapping;
    }
    $self->throw('Mapping for LocatableSeq not set') if !exists $self->{'_mapping'};
    return @{ $self->{'_mapping'} };
}
frameshiftsdescriptionprevnextTop
sub frameshifts {
    my $self = shift;
    if( @_ ) {
        if (ref $_[0] eq 'HASH') {
            $self->{_frameshifts} = $_[0];
        } else {
            # assume this is a full list to be converted to a hash
$self->{_frameshifts} =\% {@_} # coerce into hash ref
} } (defined $self->{_frameshifts} && ref $self->{_frameshifts} eq 'HASH') ? return %{$self->{_frameshifts}} : return ();
}
get_nsedescriptionprevnextTop
sub get_nse {
   my ($self,$char1,$char2) = @_;

   $char1 ||= "/";
   $char2 ||= "-";
   
   my ($id, $st, $end, $strand)  = ($self->id(), $self->start(),
                                    $self->end(), $self->strand || 0);
   
   if ($self->force_nse) {
        $id  ||= '';
        $st  ||= 0;
        $end ||= 0;
   }
   
   $self->throw("Attribute id not set") unless defined($id);
   $self->throw("Attribute start not set") unless defined($st);
   $self->throw("Attribute end not set") unless defined($end);
   
   if ($strand && $strand == -1) {
      ($st, $end) = ($end, $st);
   }
   
   #Stockholm Rfam includes version if present so it is optional
my $v = $self->version ? '.'.$self->version : ''; return join('',$id, $v, $char1, $st, $char2, $end);
}
force_nsedescriptionprevnextTop
sub force_nse {
    my ($self, $flag) = @_;
    if (defined $flag) {
        $flag ? (return $self->{'_force_nse'} = 1) : (return $self->{'_force_nse'} = 0);
    }
    return $self->{'_force_nse'};
}
num_gapsdescriptionprevnextTop
sub num_gaps {
    my ($self,$char) = @_;
    my ($seq, $count) = (undef, 0);

    # default gap characters
$char ||= $GAP_SYMBOLS; $self->warn("I hope you know what you are doing setting gap to [$char]") unless $char =~ /[$GAP_SYMBOLS]/; $seq = $self->seq; return 0 unless $seq; # empty sequence does not have gaps
$seq =~ s/^([$char]+)//; $seq =~ s/([$char]+)$//; while ( $seq =~ /[$char]+/g ) { $count++; } return $count;
}
column_from_residue_numberdescriptionprevnextTop
sub column_from_residue_number {
    my ($self, $resnumber) = @_;

    $self->throw("Residue number has to be a positive integer, not [$resnumber]")
    unless $resnumber =~ /^\d+$/ and $resnumber > 0;

    if ($resnumber >= $self->start() and $resnumber <= $self->end()) {
	my @chunks;
	my $column_incr;
	my $current_column;
	my $current_residue = $self->start - 1;
	my $seq = $self->seq;
	my $strand = $self->strand || 0;

	if ($strand == -1) {
#	    @chunks = reverse $seq =~ m/[^\.\-]+|[\.\-]+/go;
@chunks = reverse $seq =~ m/[$RESIDUE_SYMBOLS]+|[$GAP_SYMBOLS]+/go;
$column_incr = -1; $current_column = (CORE::length $seq) + 1; } else { # @chunks = $seq =~ m/[^\.\-]+|[\.\-]+/go;
@chunks = $seq =~ m/[$RESIDUE_SYMBOLS]+|[$GAP_SYMBOLS]+/go;
$column_incr = 1; $current_column = 0; } while (my $chunk = shift @chunks) { # if ($chunk =~ m|^[\.\-]|o) {
if ($chunk =~ m|^[$GAP_SYMBOLS]|o) {
$current_column += $column_incr * CORE::length($chunk);
} else { if ($current_residue + CORE::length($chunk) < $resnumber) { $current_column += $column_incr * CORE::length($chunk); $current_residue += CORE::length($chunk); } else { if ($strand == -1) { $current_column -= $resnumber - $current_residue; } else { $current_column += $resnumber - $current_residue; } return $current_column; } } } } $self->throw("Could not find residue number $resnumber");
}
location_from_columndescriptionprevnextTop
sub location_from_column {
    my ($self, $column) = @_;

    $self->throw("Column number has to be a positive integer, not [$column]")
    unless $column =~ /^\d+$/ and $column > 0;
    $self->throw("Column number [$column] is larger than".
         " sequence length [". $self->length. "]")
    unless $column <= $self->length;

    my ($loc);
    my $s = $self->subseq(1,$column);
    $s =~ s/[^a-zA-Z\*]//g;

    my $pos = CORE::length $s;

    my $start = $self->start || 0 ;
    my $strand = $self->strand() || 1;
    my $relative_pos = ($strand == -1)
        ? ($self->end - $pos + 1)
    : ($pos + $start - 1);
    if ($self->subseq($column, $column) =~ /[a-zA-Z\*]/ ) {
    $loc = Bio::Location::Simple->new
        (-start  => $relative_pos,
         -end    => $relative_pos,
         -strand => 1,
         );
    } elsif ($pos == 0 and $self->start == 1) {
    } else {
      my ($start,$end) = ($relative_pos, $relative_pos + $strand);
      if ($strand == -1) {
    ($start,$end) = ($end,$start);
      }
    $loc = Bio::Location::Simple->new
        (-start         => $start,
         -end           => $end,
         -strand        => 1,
         -location_type => 'IN-BETWEEN'
         );
    }
    return $loc;
}
revcomdescriptionprevnextTop
sub revcom {
    my ($self) = @_;
    # since we don't know whether sequences without 1 => 1 correlation can be
# revcom'd, kick back
if (grep {$_ != 1} $self->mapping) { $self->warn('revcom() not supported for sequences with mapped values of > 1'); return; } my $new = $self->SUPER::revcom; $new->strand($self->strand * -1) if $self->strand; $new->start($self->start) if $self->start; $new->end($self->end) if $self->end; return $new;
}
truncdescriptionprevnextTop
sub trunc {
    my ($self, $start, $end) = @_;
    my $new = $self->SUPER::trunc($start, $end);
    $new->strand($self->strand);

    # end will be automatically calculated
$start = $end if $self->strand && $self->strand == -1; $start = $self->location_from_column($start); $start ? ($start = $start->end) : ($start = 1); $new->start($start) if $start; return $new;
}
validate_seqdescriptionprevnextTop
sub validate_seq {
    my ($self,$seqstr) = @_;
    if( ! defined $seqstr ){ $seqstr = $self->seq(); }
    return 0 unless( defined $seqstr);
    
    if((CORE::length($seqstr) > 0) &&
       ($seqstr !~ /^([$MATCHPATTERN]+)$/)) {
        $self->warn("seq doesn't validate with [$MATCHPATTERN], mismatch is " .
            join(",",($seqstr =~ /([^$MATCHPATTERN]+)/g)));
        return 0;
    }
    return 1;
}

################## DEPRECATED METHODS ##################
}
no_gapsdescriptionprevnextTop
sub no_gaps {
	my $self = shift;
	$self->deprecated(-warn_version => 1.0069,
					  -throw_version => 1.0075,
                      -message => 'Use of method no_gaps() is deprecated, use num_gaps() instead');
    $self->num_gaps(@_);
}
no_sequencesdescriptionprevnextTop
sub no_sequences {
	my $self = shift;
	$self->deprecated(-warn_version => 1.0069,
					  -throw_version => 1.0075,
                      -message => 'Use of method no_sequences() is deprecated, use num_sequences() instead');
    $self->num_sequences(@_);
}

1;
}
General documentation
FEEDBACKTop
Mailing ListsTop
User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to one
of the Bioperl mailing lists. Your participation is much appreciated.
  bioperl-l@bioperl.org                  - General discussion
http://bioperl.org/wiki/Mailing_lists - About the mailing lists
Support Top
Please direct usage questions or support issues to the mailing list:
bioperl-l@bioperl.org
rather than to the module maintainer directly. Many experienced and
reponsive experts will be able look at the problem and quickly
address it. Please include a thorough description of the problem
with code and data examples if at all possible.
Reporting BugsTop
Report bugs to the Bioperl bug tracking system to help us keep track
the bugs and their resolution. Bug reports can be submitted via the
web:
  https://redmine.open-bio.org/projects/bioperl/
APPENDIXTop
The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a _
no_gapTop
 Title     : no_gaps
Usage : $self->no_gaps('.')
Function : Gets number of gaps in the sequence. The count excludes
leading or trailing gap characters.
Valid bioperl sequence characters are [A-Za-z\-\.\*]. Of these, '.' and '-' are counted as gap characters unless an optional argument specifies one of them. Returns : number of internal gaps in the sequence. Args : a gap character (optional) Status : Deprecated (in favor of num_gaps())