← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/LocatableSeq.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/LocatableSeq.pm
1
 
# $Id: LocatableSeq.pm,v 1.39.4.4 2006/10/02 23:10:12 sendu Exp $
 
1
# $Id: LocatableSeq.pm 15084 2008-12-03 22:31:23Z cjfields $
2
2
#
3
3
# BioPerl module for Bio::LocatableSeq
4
4
#
12
12
 
13
13
=head1 NAME
14
14
 
15
 
Bio::LocatableSeq - A Sequence object with start/end points on it
 
15
Bio::LocatableSeq - A Bio::PrimarySeq object with start/end points on it
16
16
that can be projected into a MSA or have coordinates relative to
17
17
another seq.
18
18
 
19
19
=head1 SYNOPSIS
20
20
 
21
 
 
22
21
    use Bio::LocatableSeq;
23
 
    my $seq = new Bio::LocatableSeq(-seq => "CAGT-GGT",
24
 
                                    -id  => "seq1",
25
 
                                    -start => 1,
26
 
                                    -end   => 7);
27
 
 
28
 
 
29
 
=head1 DESCRIPTION
 
22
    my $seq = Bio::LocatableSeq->new(-seq => "CAGT-GGT",
 
23
                    -id  => "seq1",
 
24
                    -start => 1,
 
25
                    -end   => 7);
30
26
 
31
27
    # a normal sequence object
32
28
    $locseq->seq();
38
34
 
39
35
    # inherits off RangeI, so range operations possible
40
36
 
 
37
=head1 DESCRIPTION
 
38
 
 
39
The LocatableSeq sequence object was developed mainly because the SimpleAlign
 
40
object requires this functionality, and in the rewrite of the Sequence object we
 
41
had to decide what to do with this.
 
42
 
 
43
It is, to be honest, not well integrated with the rest of bioperl. For example,
 
44
the trunc() function does not return a LocatableSeq object, as some might have
 
45
thought. Also, the sequence is not a Bio::SeqI, so the location is simply
 
46
inherited from Bio::RangeI and is not stored in a Bio::Location. 
 
47
 
 
48
There are all sorts of nasty gotcha's about interactions between coordinate
 
49
systems when these sort of objects are used. Some mapping now occurs to deal
 
50
with HSP data, however it can probably be integrated in better and most methods
 
51
do not implement it correctly yet. Also, several PrimarySeqI methods (subseq(),
 
52
trunc(), etc.) do not behave as expected and must be used with care.
 
53
 
 
54
Due to this, LocatableSeq functionality is to be refactored in a future BioPerl
 
55
release. However, for alignment functionality it works adequately for the time
 
56
being
 
57
 
41
58
=head1 FEEDBACK
42
59
 
43
 
 
44
60
=head2 Mailing Lists
45
61
 
46
62
User feedback is an integral part of the evolution of this and other
50
66
  bioperl-l@bioperl.org                  - General discussion
51
67
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
52
68
 
53
 
The locatable sequence object was developed mainly because the
54
 
SimpleAlign object requires this functionality, and in the rewrite
55
 
of the Sequence object we had to decide what to do with this.
56
 
 
57
 
It is, to be honest, not well integrated with the rest of bioperl, for
58
 
example, the trunc() function does not return a LocatableSeq object,
59
 
as some might have thought. There are all sorts of nasty gotcha's
60
 
about interactions between coordinate systems when these sort of
61
 
objects are used.
62
 
 
63
 
 
64
69
=head2 Reporting Bugs
65
70
 
66
71
Report bugs to the Bioperl bug tracking system to help us keep track
69
74
 
70
75
  http://bugzilla.open-bio.org/
71
76
 
72
 
 
73
77
=head1 APPENDIX
74
78
 
75
79
The rest of the documentation details each of the object
85
89
 
86
90
use Bio::Location::Simple;
87
91
use Bio::Location::Fuzzy;
88
 
 
 
92
use vars qw($GAP_SYMBOLS $OTHER_SYMBOLS $FRAMESHIFT_SYMBOLS $RESIDUE_SYMBOLS $MATCHPATTERN);
 
93
 
 
94
# The following global variables contain symbols used to represent gaps,
 
95
# frameshifts, residues, and other valid symbols. These are set at compile-time;
 
96
# expect scoping errors when using 'local' and resetting $MATCHPATTERN (see
 
97
# LocatableSeq.t)
 
98
 
 
99
$GAP_SYMBOLS = '\-\.=~';
 
100
$FRAMESHIFT_SYMBOLS = '\\\/';
 
101
$OTHER_SYMBOLS = '\?';
 
102
$RESIDUE_SYMBOLS = '0-9A-Za-z\*';
 
103
$MATCHPATTERN = $RESIDUE_SYMBOLS.$GAP_SYMBOLS.$FRAMESHIFT_SYMBOLS.$OTHER_SYMBOLS;
89
104
 
90
105
use base qw(Bio::PrimarySeq Bio::RangeI);
91
106
 
93
108
    my ($class, @args) = @_;
94
109
    my $self = $class->SUPER::new(@args);
95
110
 
96
 
    my ($start,$end,$strand) =
97
 
        $self->_rearrange( [qw(START END STRAND)],
98
 
                           @args);
99
 
 
 
111
    my ($start,$end,$strand, $mapping, $fs, $nse) =
 
112
    $self->_rearrange( [qw(START
 
113
                        END
 
114
                        STRAND
 
115
                        MAPPING
 
116
                        FRAMESHIFTS
 
117
                        FORCE_NSE
 
118
                        )],
 
119
               @args);
 
120
    
 
121
    $mapping ||= [1,1];
 
122
    $self->mapping($mapping);
 
123
    $nse || 0;
 
124
    $self->force_nse($nse);
 
125
    defined $fs    && $self->frameshifts($fs);
100
126
    defined $start && $self->start($start);
101
127
    defined $end   && $self->end($end);
102
128
    defined $strand && $self->strand($strand);
108
134
 
109
135
 Title   : start
110
136
 Usage   : $obj->start($newval)
111
 
 Function:
 
137
 Function: Get/set the 1-based start position of this sequence in the original
 
138
           sequence. '0' means before the original sequence starts.
112
139
 Returns : value of start
113
140
 Args    : newvalue (optional)
114
141
 
115
142
=cut
116
143
 
117
144
sub start{
118
 
   my $self = shift;
119
 
   if( @_ ) {
120
 
      my $value = shift;
121
 
      $self->{'start'} = $value;
122
 
  }
123
 
   return $self->{'start'} if defined $self->{'start'};
124
 
   return 1                if $self->seq;
125
 
   return;
 
145
    my $self = shift;
 
146
    if( @_ ) {
 
147
        my $value = shift;
 
148
        $self->{'start'} = $value;
 
149
    }
 
150
    return $self->{'start'} if defined $self->{'start'};
 
151
    return 1                if $self->seq;
 
152
    return;
126
153
}
127
154
 
128
155
=head2 end
129
156
 
130
157
 Title   : end
131
158
 Usage   : $obj->end($newval)
132
 
 Function:
 
159
 Function: Get/set the 1-based end position of this sequence in the original
 
160
           sequence. '0' means before the original sequence starts.
133
161
 Returns : value of end
134
162
 Args    : newvalue (optional)
 
163
 Note    : although this is a get/set, it checks passed values against the
 
164
           calculated end point ( derived from the sequence and based on
 
165
           $GAP_SYMBOLS and possible frameshifts() ).  If there is no match,
 
166
           it will warn and set the proper value.  Probably best used for
 
167
           debugging proper sequence calculations.
135
168
 
136
169
=cut
137
170
 
138
171
sub end {
139
 
   my $self = shift;
140
 
   if( @_ ) {
141
 
      my $value = shift;
142
 
      my $string = $self->seq;
143
 
      if ($self->seq) {
144
 
          my $len = $self->_ungapped_len;
145
 
          my $id = $self->id;
146
 
          $self->warn("In sequence $id residue count gives end value $len.
147
 
Overriding value [$value] with value $len for Bio::LocatableSeq::end().")
148
 
              and $value = $len if $len != $value and $self->verbose > 0;
149
 
      }
150
 
 
151
 
      $self->{'end'} = $value;
152
 
    }
153
 
 
154
 
   return $self->{'end'} || $self->_ungapped_len;
 
172
    my $self = shift;
 
173
    if( @_ ) {
 
174
        my $value = shift;
 
175
        my $st = $self->start;
 
176
        # start of 0 usually means the sequence is all gaps but maps to
 
177
        # other sequences in an alignment
 
178
        if ($self->seq && $st != 0 ) {
 
179
            my $len = $self->_ungapped_len;
 
180
            my $calend = $st + $len - 1;
 
181
            my $id = $self->id || 'unknown';
 
182
            if ($calend != $value) {
 
183
                $self->warn("In sequence $id residue count gives end value ".
 
184
                "$calend.  \nOverriding value [$value] with value $calend for ".
 
185
                "Bio::LocatableSeq::end().\n".$self->seq);
 
186
                $value = $calend;
 
187
            }
 
188
        }
 
189
        $self->{'end'} = $value;
 
190
    }
 
191
 
 
192
    if (defined $self->{'end'}) {
 
193
        return $self->{'end'}
 
194
    } elsif ( my $len = $self->_ungapped_len) {
 
195
        return $len + $self->start - 1;
 
196
    } else {
 
197
        return;
 
198
    }
155
199
}
156
200
 
 
201
# changed 08.10.26 to return ungapped length, not the calculated end
 
202
# of the sequence
157
203
sub _ungapped_len {
158
204
    my $self = shift;
159
 
    my $string = $self->seq || '';
160
 
    $string =~ s/[\.\-]+//g;
161
 
    $self->seq ? (return $self->start + CORE::length($string) - 1 ) : undef;
 
205
    return unless my $string = $self->seq;
 
206
    my ($map_res, $map_coord) = $self->mapping;
 
207
    my $offset = 0;
 
208
    if (my %data = $self->frameshifts) {
 
209
        map {$offset += $_} values %data;
 
210
    }
 
211
    $string =~ s{[$GAP_SYMBOLS$FRAMESHIFT_SYMBOLS]+}{}g;
 
212
    return CORE::length($string)/($map_res/$map_coord) + $offset/($map_coord/$map_res);
162
213
}
163
214
 
164
215
=head2 strand
165
216
 
166
217
 Title   : strand
167
218
 Usage   : $obj->strand($newval)
168
 
 Function:
169
 
 Returns : value of strand
170
 
 Args    : newvalue (optional)
 
219
 Function: return or set the strandedness
 
220
 Returns : the value of the strandedness (-1, 0 or 1)
 
221
 Args    : the value of the strandedness (-1, 0 or 1)
171
222
 
172
223
=cut
173
224
 
174
225
sub strand{
175
226
   my $self = shift;
176
227
   if( @_ ) {
177
 
      my $value = shift;
178
 
      $self->{'strand'} = $value;
 
228
        my $value = shift;
 
229
        $self->{'strand'} = $value;
179
230
    }
180
231
    return $self->{'strand'};
181
232
}
182
233
 
 
234
=head2 mapping
 
235
 
 
236
 Title   : mapping
 
237
 Usage   : $obj->mapping($newval)
 
238
 Function: return or set the mapping indices (indicates # symbols/positions in
 
239
           the source string mapping to # of coordinate positions)
 
240
 Returns : two-element array (# symbols => # coordinate pos)
 
241
 Args    : two elements (# symbols => # coordinate pos); this can also be
 
242
           passed in as an array reference of the two elements (as might be
 
243
           passed upon Bio::LocatableSeq instantiation, for instance).
 
244
 
 
245
=cut
 
246
 
 
247
sub mapping {
 
248
    my $self = shift;
 
249
    if( @_ ) {
 
250
        my @mapping = (ref($_[0]) eq 'ARRAY') ? @{$_[0]} : @_;
 
251
        $self->throw("Must pass two values (# residues mapped to # positions)")
 
252
            if @mapping != 2;
 
253
        if ((grep {$_ != 1 && $_ != 3} @mapping) || ($mapping[0] == 3 && $mapping[1] == 3)) {
 
254
            $self->throw("Mapping values other than 1 or 3 are not currently supported")
 
255
        }
 
256
        $self->{'_mapping'} = \@mapping;
 
257
    }
 
258
    $self->throw('Mapping for LocatableSeq not set') if !exists $self->{'_mapping'};
 
259
    return @{ $self->{'_mapping'} };
 
260
}
 
261
 
 
262
=head2 frameshifts
 
263
 
 
264
 Title   : frameshifts
 
265
 Usage   : $obj->frameshifts($newval)
 
266
 Function: get/set the frameshift hash, which contains sequence positions as
 
267
           keys and the shift (-2, -1, 1, 2) as the value
 
268
 Returns : hash
 
269
 Args    : hash or hash reference
 
270
 
 
271
=cut
 
272
 
 
273
sub frameshifts {
 
274
    my $self = shift;
 
275
    if( @_ ) {
 
276
        if (ref $_[0] eq 'HASH') {
 
277
            $self->{_frameshifts} = $_[0];
 
278
        } else {
 
279
            # assume this is a full list to be converted to a hash
 
280
            $self->{_frameshifts} = \%{@_} # coerce into hash ref
 
281
        }
 
282
    }
 
283
    (defined $self->{_frameshifts} && ref $self->{_frameshifts} eq 'HASH') ?
 
284
        return %{$self->{_frameshifts}} : return ();
 
285
}
 
286
 
183
287
=head2 get_nse
184
288
 
185
289
 Title   : get_nse
196
300
 
197
301
   $char1 ||= "/";
198
302
   $char2 ||= "-";
199
 
 
200
 
   $self->throw("Attribute id not set") unless defined($self->id());
201
 
   $self->throw("Attribute start not set") unless defined($self->start());
202
 
   $self->throw("Attribute end not set") unless defined($self->end());
203
 
 
204
 
   return $self->id() . $char1 . $self->start . $char2 . $self->end ;
205
 
 
206
 
}
207
 
 
 
303
   
 
304
   my ($id, $st, $end)  = ($self->id(), $self->start(), $self->end());
 
305
   
 
306
   if ($self->force_nse) {
 
307
        $id  ||= '';
 
308
        $st  ||= 0;
 
309
        $end ||= 0;
 
310
   }
 
311
   
 
312
   $self->throw("Attribute id not set") unless defined($id);
 
313
   $self->throw("Attribute start not set") unless defined($st);
 
314
   $self->throw("Attribute end not set") unless defined($end);
 
315
   
 
316
   #Stockholm Rfam includes version if present so it is optional
 
317
   my $v = $self->version ? '.'.$self->version : ''; 
 
318
   return $id . $v. $char1 . $st . $char2 . $end ;
 
319
}
 
320
 
 
321
=head2 force_nse
 
322
 
 
323
 Title   : force_nse
 
324
 Usage   : $ls->force_nse()
 
325
 Function: Boolean which forces get_nse() to build an NSE, regardless
 
326
           of whether id(), start(), or end() is set
 
327
 Returns : Boolean value
 
328
 Args    : (optional) Boolean (1 or 0)
 
329
 Note    : This will convert any passed value evaluating as TRUE/FALSE to 1/0
 
330
           respectively
 
331
 
 
332
=cut
 
333
 
 
334
sub force_nse {
 
335
    my ($self, $flag) = @_;
 
336
    if (defined $flag) {
 
337
        $flag ? (return $self->{'_force_nse'} = 1) : (return $self->{'_force_nse'} = 0);
 
338
    }
 
339
    return $self->{'_force_nse'};
 
340
}
208
341
 
209
342
=head2 no_gap
210
343
 
211
344
 Title   : no_gaps
212
345
 Usage   :$self->no_gaps('.')
213
 
 Function:
214
 
 
215
 
           Gets number of gaps in the sequence. The count excludes
 
346
 Function:Gets number of gaps in the sequence. The count excludes
216
347
           leading or trailing gap characters.
217
348
 
218
349
           Valid bioperl sequence characters are [A-Za-z\-\.\*]. Of
219
350
           these, '.' and '-' are counted as gap characters unless an
220
351
           optional argument specifies one of them.
221
352
 
222
 
 Returns : number of internal gaps in the sequnce.
 
353
 Returns : number of internal gaps in the sequence.
223
354
 Args    : a gap character (optional)
224
355
 
225
356
=cut
229
360
    my ($seq, $count) = (undef, 0);
230
361
 
231
362
    # default gap characters
232
 
    $char ||= '-.';
 
363
    $char ||= $GAP_SYMBOLS;
233
364
 
234
365
    $self->warn("I hope you know what you are doing setting gap to [$char]")
235
 
        unless $char =~ /[-.]/;
 
366
        unless $char =~ /[$GAP_SYMBOLS]/;
236
367
 
237
368
    $seq = $self->seq;
238
369
    return 0 unless $seq; # empty sequence does not have gaps
239
370
 
240
371
    $seq =~ s/^([$char]+)//;
241
372
    $seq =~ s/([$char]+)$//;
242
 
    $count++ while $seq =~ /[$char]+/g;
 
373
    while ( $seq =~ /[$char]+/g ) {
 
374
        $count++;
 
375
    }
243
376
 
244
377
    return $count;
245
 
 
246
378
}
247
379
 
248
380
 
256
388
           (i.e. column number) of the given residue number in the
257
389
           sequence. For example, for the sequence
258
390
 
259
 
             Seq1/91-97 AC..DEF.GH
 
391
         Seq1/91-97 AC..DEF.GH
260
392
 
261
 
           column_from_residue_number(94) returns 5.
 
393
           column_from_residue_number(94) returns 6.
262
394
 
263
395
           An exception is thrown if the residue number would lie
264
396
           outside the length of the aligment
275
407
    my ($self, $resnumber) = @_;
276
408
 
277
409
    $self->throw("Residue number has to be a positive integer, not [$resnumber]")
278
 
        unless $resnumber =~ /^\d+$/ and $resnumber > 0;
 
410
    unless $resnumber =~ /^\d+$/ and $resnumber > 0;
279
411
 
280
412
    if ($resnumber >= $self->start() and $resnumber <= $self->end()) {
281
 
        my @residues = split //, $self->seq;
282
 
        my $count = $self->start();
283
 
        my $i;
284
 
        my ($start,$end,$inc,$test);
 
413
    my @residues = split //, $self->seq;
 
414
    my $count = $self->start();
 
415
    my $i;
 
416
    my ($start,$end,$inc,$test);
285
417
        my $strand = $self->strand || 0;
286
 
        # the following bit of "magic" allows the main loop logic to be the
287
 
        # same regardless of the strand of the sequence
288
 
        ($start,$end,$inc,$test)= ($strand == -1)?
 
418
    # the following bit of "magic" allows the main loop logic to be the
 
419
    # same regardless of the strand of the sequence
 
420
    ($start,$end,$inc,$test)= ($strand == -1)?
289
421
            (scalar(@residues-1),0,-1,sub{$i >= $end}) :
290
422
                (0,scalar(@residues-1),1,sub{$i <= $end});
291
423
 
292
 
        for ($i=$start; $test->(); $i+= $inc) {
293
 
            if ($residues[$i] ne '.' and $residues[$i] ne '-') {
294
 
                $count == $resnumber and last;
295
 
                $count++;
296
 
            }
297
 
        }
298
 
        # $i now holds the index of the column.
 
424
    for ($i=$start; $test->(); $i+= $inc) {
 
425
        if ($residues[$i] ne '.' and $residues[$i] ne '-') {
 
426
        $count == $resnumber and last;
 
427
        $count++;
 
428
        }
 
429
    }
 
430
    # $i now holds the index of the column.
299
431
        # The actual column number is this index + 1
300
432
 
301
 
        return $i+1;
 
433
    return $i+1;
302
434
    }
303
435
 
304
436
    $self->throw("Could not find residue number $resnumber");
305
437
 
306
438
}
307
439
 
308
 
 
309
440
=head2 location_from_column
310
441
 
311
442
 Title   : location_from_column
318
449
           L<Bio::Range> where values can be undefined. For example,
319
450
           for the sequence:
320
451
 
321
 
             Seq/91-97 .AC..DEF.G.
 
452
         Seq/91-96 .AC..DEF.G.
322
453
 
323
 
           location_from_column( 3 ) position 93
324
 
           location_from_column( 2 ) position 92^93
325
 
           location_from_column(10 ) position 97^98
 
454
           location_from_column( 3 ) position 92
 
455
           location_from_column( 4 ) position 92^93
 
456
           location_from_column( 9 ) position 95^96
326
457
           location_from_column( 1 ) position undef
327
458
 
328
459
           An exact position returns a Bio::Location::Simple object
348
479
    my ($self, $column) = @_;
349
480
 
350
481
    $self->throw("Column number has to be a positive integer, not [$column]")
351
 
        unless $column =~ /^\d+$/ and $column > 0;
 
482
    unless $column =~ /^\d+$/ and $column > 0;
352
483
    $self->throw("Column number [$column] is larger than".
353
 
                 " sequence length [". $self->length. "]")
354
 
        unless $column <= $self->length;
 
484
         " sequence length [". $self->length. "]")
 
485
    unless $column <= $self->length;
355
486
 
356
487
    my ($loc);
357
488
    my $s = $self->subseq(1,$column);
363
494
    my $strand = $self->strand() || 1;
364
495
    my $relative_pos = ($strand == -1)
365
496
        ? ($self->end - $pos + 1)
366
 
        : ($pos + $start - 1);
 
497
    : ($pos + $start - 1);
367
498
    if ($self->subseq($column, $column) =~ /[a-zA-Z\*]/ ) {
368
 
        $loc = new Bio::Location::Simple
369
 
            (-start  => $relative_pos,
370
 
             -end    => $relative_pos,
371
 
             -strand => 1,
372
 
             );
 
499
    $loc = Bio::Location::Simple->new
 
500
        (-start  => $relative_pos,
 
501
         -end    => $relative_pos,
 
502
         -strand => 1,
 
503
         );
373
504
    } elsif ($pos == 0 and $self->start == 1) {
374
505
    } else {
375
506
      my ($start,$end) = ($relative_pos, $relative_pos + $strand);
376
507
      if ($strand == -1) {
377
 
        ($start,$end) = ($end,$start);
 
508
    ($start,$end) = ($end,$start);
378
509
      }
379
 
        $loc = new Bio::Location::Simple
380
 
            (-start         => $start,
381
 
             -end           => $end,
382
 
             -strand        => 1,
383
 
             -location_type => 'IN-BETWEEN'
384
 
             );
 
510
    $loc = Bio::Location::Simple->new
 
511
        (-start         => $start,
 
512
         -end           => $end,
 
513
         -strand        => 1,
 
514
         -location_type => 'IN-BETWEEN'
 
515
         );
385
516
    }
386
517
    return $loc;
387
518
}
402
533
 
403
534
sub revcom {
404
535
    my ($self) = @_;
405
 
 
 
536
    # since we don't know whether sequences without 1 => 1 correlation can be
 
537
    # revcom'd, kick back
 
538
    if (grep {$_ != 1} $self->mapping) {
 
539
        $self->warn('revcom() not supported for sequences with mapped values of > 1');
 
540
        return;
 
541
    }
406
542
    my $new = $self->SUPER::revcom;
407
 
    $new->strand($self->strand * -1);
 
543
    $new->strand($self->strand * -1) if $self->strand;
408
544
    $new->start($self->start) if $self->start;
409
545
    $new->end($self->end) if $self->end;
410
546
    return $new;
411
547
}
412
548
 
413
 
 
414
549
=head2 trunc
415
550
 
416
551
 Title   : trunc
440
575
    return $new;
441
576
}
442
577
 
 
578
=head2 validate_seq
 
579
 
 
580
 Title   : validate_seq
 
581
 Usage   : if(! $seq->validate_seq($seq_str) ) {
 
582
                print "sequence $seq_str is not valid for an object of
 
583
                alphabet ",$seq->alphabet, "\n";
 
584
            }
 
585
 Function: Validates a given sequence string. A validating sequence string
 
586
           must be accepted by seq(). A string that does not validate will
 
587
           lead to an exception if passed to seq().
 
588
 
 
589
           The implementation provided here does not take alphabet() into
 
590
           account. Allowed are all letters (A-Z), numbers [0-9] 
 
591
           and common symbols used for gaps, stop codons, unknown residues,
 
592
           and frameshifts, including '-','.','*','?','=',and '~'.
 
593
 
 
594
 Example :
 
595
 Returns : 1 if the supplied sequence string is valid for the object, and
 
596
           0 otherwise.
 
597
 Args    : The sequence string to be validated.
 
598
 
 
599
=cut
 
600
 
 
601
sub validate_seq {
 
602
    my ($self,$seqstr) = @_;
 
603
    if( ! defined $seqstr ){ $seqstr = $self->seq(); }
 
604
    return 0 unless( defined $seqstr);
 
605
    
 
606
    if((CORE::length($seqstr) > 0) &&
 
607
       ($seqstr !~ /^([$MATCHPATTERN]+)$/)) {
 
608
        $self->warn("seq doesn't validate with [$MATCHPATTERN], mismatch is " .
 
609
            join(",",($seqstr =~ /([^$MATCHPATTERN]+)/g)));
 
610
        return 0;
 
611
    }
 
612
    return 1;
 
613
}
 
614
 
443
615
1;