File Coverage

blib/lib/Ace/Sequence.pm
Criterion Covered Total %
statement 278 393 70.7
branch 107 204 52.5
condition 38 81 46.9
subroutine 40 58 69.0
pod 19 34 55.9
total 482 770 62.6


line stmt bran cond sub pod time code
1             package Ace::Sequence;
2 1     1   14 use strict;
  1         14  
  1         13  
3              
4 1     1   14 use Carp;
  1         8  
  1         22  
5 1     1   15 use strict;
  1         9  
  1         11  
6 1     1   30 use Ace 1.50 qw(:DEFAULT rearrange);
  1         24  
  1         34  
7 1     1   41 use Ace::Sequence::FeatureList;
  1         12  
  1         22  
8 1     1   37 use Ace::Sequence::Feature;
  1         10  
  1         27  
9 1     1   172 use AutoLoader 'AUTOLOAD';
  1         22  
  1         17  
10 1     1   18 use vars '$VERSION';
  1         10  
  1         14  
11             my %CACHE;
12              
13             $VERSION = '1.51';
14              
15 1     1   16 use constant CACHE => 1;
  1         10  
  1         15  
16              
17             use overload
18 1         16   '""' => 'asString',
19               cmp => 'cmp',
20 1     1   16 ;
  1         9  
21              
22             # synonym: stop = end
23             *stop = \&end;
24             *abs = \&absolute;
25             *source_seq = \&source;
26             *source_tag = \&subtype;
27             *primary_tag = \&type;
28              
29             my %plusminus = ( '+' => '-',
30             '-' => '+',
31             '.' => '.');
32              
33             # internal keys
34             # parent => reference Sequence in "+" strand
35             # p_offset => our start in the parent
36             # length => our length
37             # strand => our strand (+ or -)
38             # refseq => reference Sequence for coordinate system
39              
40             # object constructor
41             # usually called like this:
42             # $seq = Ace::Sequence->new($object);
43             # but can be called like this:
44             # $seq = Ace::Sequence->new(-db=>$db,-name=>$name);
45             # or
46             # $seq = Ace::Sequence->new(-seq => $object,
47             # -offset => $offset,
48             # -length => $length,
49             # -ref => $refseq
50             # );
51             # $refseq, if provided, will be used to establish the coordinate
52             # system. Otherwise the first base pair will be set to 1.
53             sub new {
54 25     25 0 245   my $pack = shift;
55 25         889   my ($seq,$start,$end,$offset,$length,$refseq,$db) =
56                 rearrange([
57             ['SEQ','SEQUENCE','SOURCE'],
58             'START',
59             ['END','STOP'],
60             ['OFFSET','OFF'],
61             ['LENGTH','LEN'],
62             'REFSEQ',
63             ['DATABASE','DB'],
64             ],@_);
65              
66             # Object must have a parent sequence and/or a reference
67             # sequence. In some cases, the parent sequence will be the
68             # object itself. The reference sequence is used to set up
69             # the frame of reference for the coordinate system.
70              
71             # fetch the sequence object if we don't have it already
72 25 50 33     470   croak "Please provide either a Sequence object or a database and name"
      66        
73                 unless ref($seq) || ($seq && $db);
74              
75             # convert start into offset
76 25 50 33     253   $offset = $start - 1 if defined($start) and !defined($offset);
77              
78             # convert stop/end into length
79 25 0 33     244   $length = ($end > $start) ? $end - $offset : $end - $offset - 2
    50          
80                 if defined($end) && !defined($length);
81              
82             # if just a string is passed, try to fetch a Sequence object
83 25 100       281   my $obj = ref($seq) ? $seq : $db->fetch('Sequence'=>$seq);
84 25 50       291   unless ($obj) {
85 0         0     Ace->error("No Sequence named $obj found in database");
86 0         0     return;
87               }
88              
89             # get parent coordinates and length of this sequence
90             # the parent is an Ace Sequence object in the "+" strand
91 25         433   my ($parent,$p_offset,$p_length,$strand) = find_parent($obj);
92 25 50       1625   return unless $parent;
93              
94             # handle negative strands
95 25         271   my $r_strand = $strand;
96 25         234   my $r_offset = $p_offset;
97 25   100     286   $offset ||= 0;
98 25 100       241   $offset *= -1 if $strand < 0;
99              
100             # handle feature objects
101 25 100       632   $offset += $obj->offset if $obj->can('smapped');
102              
103             # get source
104 25 100       340   my $source = $obj->can('smapped') ? $obj->source : $obj;
105              
106             # store the object into our instance variables
107 25   66     1025   my $self = bless {
108             obj        => $source,
109             offset     => $offset,
110             length     => $length || $p_length,
111             parent     => $parent,
112             p_offset   => $p_offset,
113             refseq     => [$source,$r_offset,$r_strand],
114             strand     => $strand,
115             absolute   => 0,
116             automerge  => 1,
117             },$pack;
118              
119             # set the reference sequence
120 25 100 50     269   eval { $self->refseq($refseq) } or return if defined $refseq;
  17         175  
121              
122             # wheww!
123 25         391   return $self;
124             }
125              
126             # return the "source" object that the user offset from
127             sub source {
128 141     141 0 1920   $_[0]->{obj};
129             }
130              
131             # return the parent object
132 79     79 0 890 sub parent { $_[0]->{parent} }
133              
134             # return the length
135             #sub length { $_[0]->{length} }
136             sub length {
137 6     6 1 56   my $self = shift;
138 6         57   my ($start,$end) = ($self->start,$self->end);
139 6 100       132   return $end - $start + ($end > $start ? 1 : -1); # for stupid 1-based adjustments
140             }
141              
142 0     0 1 0 sub reversed { return shift->strand < 0; }
143              
144             sub automerge {
145 10     10 1 88   my $self = shift;
146 10         196   my $d = $self->{automerge};
147 10 50       97   $self->{automerge} = shift if @_;
148 10         106   $d;
149             }
150              
151             # return reference sequence
152             sub refseq {
153 1018     1018 1 8496   my $self = shift;
154 1018         10917   my $prev = $self->{refseq};
155 1018 100       14215   if (@_) {
156 18         155     my $refseq = shift;
157 18         148     my $arrayref;
158              
159 18 50       172   BLOCK: {
160 18         148       last BLOCK unless defined ($refseq);
161              
162 18 50 66     291       if (ref($refseq) && ref($refseq) eq 'ARRAY') {
163 0         0 $arrayref = $refseq;
164 0         0 last BLOCK;
165                   }
166              
167 18 100 66     314       if (ref($refseq) && ($refseq->can('smapped'))) {
168 16 50       160 croak "Reference sequence has no common ancestor with sequence"
169             unless $self->parent eq $refseq->parent;
170 16         203 my ($a,$b,$c) = @{$refseq->{refseq}};
  16         175  
171             # $b += $refseq->offset;
172 16         160 $b += $refseq->offset;
173 16         202 $arrayref = [$refseq,$b,$refseq->strand];
174 16         178 last BLOCK;
175                   }
176              
177              
178             # look up reference sequence in database if we aren't given
179             # database object already
180 2 50       51       $refseq = $self->db->fetch('Sequence' => $refseq)
181             unless $refseq->isa('Ace::Object');
182 2 50       21       croak "Invalid reference sequence" unless $refseq;
183              
184             # find position of ref sequence in parent strand
185 2         25       my ($r_parent,$r_offset,$r_length,$r_strand) = find_parent($refseq);
186 2 50       22       croak "Reference sequence has no common ancestor with sequence"
187             unless $r_parent eq $self->{parent};
188              
189             # set to array reference containing this information
190 2         22       $arrayref = [$refseq,$r_offset,$r_strand];
191                 }
192 18         202     $self->{refseq} = $arrayref;
193               }
194 1018 50       11667   return unless $prev;
195 1018 50       9912   return $self->parent if $self->absolute;
196 1018 100       12546   return wantarray ? @{$prev} : $prev->[0];
  534         8019  
197             }
198              
199             # return strand
200 0     0 1 0 sub strand { return $_[0]->{strand} }
201              
202             # return reference strand
203             sub r_strand {
204 134     134 0 1125   my $self = shift;
205 134 50       1417   return "+1" if $self->absolute;
206 134 50       1332   if (my ($ref,$r_offset,$r_strand) = $self->refseq) {
207 134         2565     return $r_strand;
208               } else {
209 0         0     return $self->{strand}
210               }
211             }
212              
213 32     32 1 365 sub offset { $_[0]->{offset} }
214 0     0 0 0 sub p_offset { $_[0]->{p_offset} }
215              
216 0     0 0 0 sub smapped { 1; }
217 0     0 0 0 sub type { 'Sequence' }
218 0     0 0 0 sub subtype { }
219              
220             sub debug {
221 24     24 0 228   my $self = shift;
222 24         221   my $d = $self->{_debug};
223 24 50       218   $self->{_debug} = shift if @_;
224 24         292   $d;
225             }
226              
227             # return the database this sequence is associated with
228             sub db {
229 284   66 284 1 3818   return Ace->name2db($_[0]->{db} ||= $_[0]->source->db);
230             }
231              
232             sub start {
233 425     425 1 5482   my ($self,$abs) = @_;
234 425 100       4894   $abs = $self->absolute unless defined $abs;
235 425 100       4219   return $self->{p_offset} + $self->{offset} + 1 if $abs;
236              
237 389 50       3812   if ($self->refseq) {
238 389         4468     my ($ref,$r_offset,$r_strand) = $self->refseq;
239 389 100       5833     return $r_strand < 0 ? 1 + $r_offset - ($self->{p_offset} + $self->{offset})
240                                      : 1 + $self->{p_offset} + $self->{offset} - $r_offset;
241               }
242              
243               else {
244 0         0     return $self->{offset} +1;
245               }
246              
247             }
248              
249             sub end {
250 106     106 1 1196   my ($self,$abs) = @_;
251 106         1027   my $start = $self->start($abs);
252 106 100       1448   my $f = $self->{length} > 0 ? 1 : -1; # for stupid 1-based adjustments
253 106 100 100     1140   if ($abs && $self->refseq ne $self->parent) {
254 17         249     my $r_strand = $self->r_strand;
255 17 100 66     328     return $start - $self->{length} + $f
      100        
256                   if $r_strand < 0 or $self->{strand} < 0 or $self->{length} < 0;
257 8         97     return $start + $self->{length} - $f
258               }
259 89 100       1177   return $start + $self->{length} - $f if $self->r_strand eq $self->{strand};
260 7         112   return $start - $self->{length} + $f;
261             }
262              
263             # turn on absolute coordinates (relative to reference sequence)
264             sub absolute {
265 1611     1611 1 15964   my $self = shift;
266 1611         17364   my $prev = $self->{absolute};
267 1611 50       18331   $self->{absolute} = $_[0] if defined $_[0];
268 1611         17478   return $prev;
269             }
270              
271             # human readable string (for debugging)
272             sub asString {
273 59     59 1 598   my $self = shift;
274 59 50       1540   if ($self->absolute) {
    50          
275 0         0     return join '',$self->parent,'/',$self->start,',',$self->end;
276               } elsif (my $ref = $self->refseq){
277 59 50       944     my $label = $ref->isa('Ace::Sequence::Feature') ? $ref->info : "$ref";
278 59         789     return join '',$label,'/',$self->start,',',$self->end;
279              
280               } else {
281 0         0     join '',$self->source,'/',$self->start,',',$self->end;
282               }
283             }
284              
285             sub cmp {
286 3     3 0 33   my ($self,$arg,$reversed) = @_;
287 3 50 33     67   if (ref($arg) and $arg->isa('Ace::Sequence')) {
288 0   0     0     my $cmp = $self->parent cmp $arg->parent
289                   || $self->start <=> $arg->start;
290 0 0       0     return $reversed ? -$cmp : $cmp;
291               }
292 3         46   my $name = $self->asString;
293 3 50       34   return $reversed ? $arg cmp $name : $name cmp $arg;
294             }
295              
296             # Return the DNA
297             sub dna {
298 12     12 1 173   my $self = shift;
299 12 50       126   return $self->{dna} if $self->{dna};
300