view release on metacpan or  search on metacpan

Bio/Map/GeneMap.pm  view on Meta::CPAN

    if (defined($value)) {
        $value >= 0 || $self->throw("Supplied value must be a positive integer");
        $pos->start($value + 1);
    }
    
    return $pos->start - 1;
}

=head2 downstream

 Title   : downstream
 Usage   : my $distance = $map->downstream;
           $map->downstream($distance);
 Function: Get/set the nominal end of the map relative to the end of the
           Bio::Map::Gene object on this map.
 Returns : int
 Args    : none to get, OR int to set (the number of bases the map extends
           beyond the end of the gene)

=cut

sub downstream {
    my $self = shift;
    if (@_) { $self->{_downstream} = shift }
    return $self->{_downstream} || 0;
}

=head2 length

 Title   : length
 Usage   : my $length = $map->length();
 Function: Retrieves the length of the map. This is normally the length of the
           upstream region + length of the gene + length of the downstream
           region, but may be longer if positions have been placed on the map
           beyond the end of the nominal downstream region.
 Returns : int
 Args    : none

=cut

sub length {
	my $self = shift;
	my $expected_length = $self->gene->position($self)->length + $self->upstream + $self->downstream;
    my $actual_length = $self->SUPER::length;
    return $actual_length > $expected_length ? $actual_length : $expected_length;
}

=head2 seq

 Title   : seq
 Usage   : $string = $obj->seq()
 Function: Get/set the sequence as a string of letters. When getting, If the
           GeneMap object didn't have sequence attached directly to it for the
           region requested, the map's gene's database will be asked for the
           sequence, and failing that, the map's gene's positions will be asked
           for their sequences. Areas for which no sequence could be found will
           be filled with Ns, unless no sequence was found anywhere, in which
           case undef is returned.
 Returns : string
 Args    : Optionally on set the new value (a string). An optional second
           argument presets the alphabet (otherwise it will be guessed).

=cut

sub seq {
    my ($self, @args) = @_;
    my $seq = $self->SUPER::seq(@args);
    my $expected_length = $self->length;
    if (! $seq || CORE::length($seq) < $expected_length) {
        my @have = split('', $seq || '');
        my @result;
        for (0..($expected_length - 1)) {
            $result[$_] = shift(@have) || 'N';
        }
        
        # build map sequence by asking gene or positions
        my @slice_stuff = $self->gene->_get_slice($self);
        if (@slice_stuff) {
            my ($slice_adaptor, $slice, $strand) = @slice_stuff;
            my ($start, $end, $gene_start) = (CORE::length($seq || '') + 1, $expected_length, $self->upstream + 1);
            
            # convert map coords to genomic coords
            my $adjust = $strand == -1 ? $slice->end : $slice->start;
            my $adjustment = sub { return $strand == -1 ? $adjust - shift() : shift() + $adjust; };
            my $converted_start = &$adjustment($start - $gene_start);
            my $converted_end = &$adjustment($end - $gene_start);
            ($converted_start, $converted_end) = ($converted_end, $converted_start) if $converted_start > $converted_end;
            
            # get sequence from a new slice of desired region
            #*** what happens if desired region starts or ends off end of chromo?...
            my $new_slice = $slice_adaptor->fetch_by_region($slice->coord_system_name, $slice->seq_region_name, $converted_start, $converted_end);
            if ($new_slice && (my $seq_str = $new_slice->seq)) {
                if ($strand == -1) {
                    $seq_str = $self->_revcom($seq_str);
                }
                splice(@result, CORE::length($seq || ''), CORE::length($seq_str), split('', $seq_str));
            }
        }
        else {
            foreach my $pos ($self->get_positions) {
                next unless $pos->can('seq');
                my @pos_seq = split('', $pos->seq(undef, undef, 1) || next);
                for my $i ($pos->start($pos->absolute_relative)..$pos->end($pos->absolute_relative)) {
                    $i--;
                    my $base = shift(@pos_seq);
                    if ($result[$i] eq 'N') {
                        $result[$i] = $base;
                    }
                }
            }
        }
        
        $seq = join('', @result);
    }
    return $seq;
}

=head2 subseq

 Title   : subseq
 Usage   : $substring = $obj->subseq(10, 40);