← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/Tools/IUPAC.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Tools/IUPAC.pm
13
13
 
14
14
=head1 NAME
15
15
 
16
 
Bio::Tools::IUPAC - Generates unique Seq objects from an ambiguous Seq object
 
16
Bio::Tools::IUPAC - Generates unique sequence objects or regular expressions from
 
17
an ambiguous IUPAC sequence
17
18
 
18
19
=head1 SYNOPSIS
19
20
 
20
 
 use Bio::Seq;
 
21
 use Bio::PrimarySeq;
21
22
 use Bio::Tools::IUPAC;
22
23
 
23
 
 my $ambiseq = Bio::Seq->new(-seq => 'ARTCGUTGR', -alphabet => 'dna');
24
 
 my $stream  = Bio::Tools::IUPAC->new(-seq => $ambiseq);
25
 
 
26
 
 while ($uniqueseq = $stream->next_seq()) {
27
 
     # process the unique Seq object.
 
24
 # Get the IUPAC code for proteins
 
25
 my %iupac_prot = Bio::Tools::IUPAC->new->iupac_iup;
 
26
 
 
27
 # Create a sequence with degenerate residues
 
28
 my $ambiseq = Bio::PrimarySeq->new(-seq => 'ARTCGUTGN', -alphabet => 'dna');
 
29
 
 
30
 # Create all possible non-degenerate sequences
 
31
 my $iupac = Bio::Tools::IUPAC->new(-seq => $ambiseq);
 
32
 while ($uniqueseq = $iupac->next_seq()) {
 
33
     # process the unique Bio::Seq object.
28
34
 }
29
35
 
 
36
 # Get a regular expression that matches all possible sequences
 
37
 my $regexp = $iupac->regexp();
 
38
 
30
39
=head1 DESCRIPTION
31
40
 
32
 
IUPAC is a tool that produces a stream of unique, "strict"-satisfying Seq
33
 
objects from an ambiquous Seq object (containing non-standard characters given
34
 
the meaning shown below)
35
 
 
36
 
        Extended DNA / RNA alphabet :
37
 
        (includes symbols for nucleotide ambiguity)
38
 
        ------------------------------------------
39
 
        Symbol       Meaning      Nucleic Acid
40
 
        ------------------------------------------
41
 
         A            A           Adenine
42
 
         C            C           Cytosine
43
 
         G            G           Guanine
44
 
         T            T           Thymine
45
 
         U            U           Uracil
46
 
         M          A or C
47
 
         R          A or G
48
 
         W          A or T
49
 
         S          C or G
50
 
         Y          C or T
51
 
         K          G or T
52
 
         V        A or C or G
53
 
         H        A or C or T
54
 
         D        A or G or T
55
 
         B        C or G or T
56
 
         X      G or A or T or C
57
 
         N      G or A or T or C
58
 
 
59
 
        IUPAC-IUB SYMBOLS FOR NUCLEOTIDE NOMENCLATURE:
60
 
          Cornish-Bowden (1985) Nucl. Acids Res. 13: 3021-3030.
61
 
 
62
 
-----------------------------------
63
 
 
64
 
       Amino Acid alphabet:
65
 
        ------------------------------------------
66
 
        Symbol           Meaning
67
 
        ------------------------------------------
68
 
        A        Alanine
69
 
        B        Aspartic Acid, Asparagine
70
 
        C        Cystine
71
 
        D        Aspartic Acid
72
 
        E        Glutamic Acid
73
 
        F        Phenylalanine
74
 
        G        Glycine
75
 
        H        Histidine
76
 
        I        Isoleucine
77
 
        J        Isoleucine/Leucine
78
 
        K        Lysine
79
 
        L        Leucine
80
 
        M        Methionine
81
 
        N        Asparagine
82
 
        O        Pyrrolysine
83
 
        P        Proline
84
 
        Q        Glutamine
85
 
        R        Arginine
86
 
        S        Serine
87
 
        T        Threonine
88
 
        U        Selenocysteine
89
 
        V        Valine
90
 
        W        Tryptophan
91
 
        X        Unknown
92
 
        Y        Tyrosine
93
 
        Z        Glutamic Acid, Glutamine
94
 
        *        Terminator
95
 
 
96
 
 
97
 
        IUPAC-IUP AMINO ACID SYMBOLS:
98
 
          Biochem J. 1984 Apr 15; 219(2): 345-373
99
 
          Eur J Biochem. 1993 Apr 1; 213(1): 2
 
41
Bio::Tools::IUPAC is a tool that manipulates sequences with ambiguous residues
 
42
following the IUPAC conventions. Non-standard characters have the meaning 
 
43
described below:
 
44
 
 
45
    IUPAC-IUB SYMBOLS FOR NUCLEOTIDE (DNA OR RNA) NOMENCLATURE:
 
46
      Cornish-Bowden (1985) Nucl. Acids Res. 13: 3021-3030
 
47
 
 
48
    ------------------------------------------
 
49
    Symbol       Meaning      Nucleic Acid
 
50
    ------------------------------------------
 
51
     A            A           Adenine
 
52
     C            C           Cytosine
 
53
     G            G           Guanine
 
54
     T            T           Thymine
 
55
     U            U           Uracil
 
56
     M          A or C
 
57
     R          A or G
 
58
     W          A or T
 
59
     S          C or G
 
60
     Y          C or T
 
61
     K          G or T
 
62
     V        A or C or G
 
63
     H        A or C or T
 
64
     D        A or G or T
 
65
     B        C or G or T
 
66
     X      G or A or T or C
 
67
     N      G or A or T or C
 
68
 
 
69
 
 
70
    IUPAC-IUP AMINO ACID SYMBOLS:
 
71
      Biochem J. 1984 Apr 15; 219(2): 345-373
 
72
      Eur J Biochem. 1993 Apr 1; 213(1): 2
 
73
 
 
74
    ------------------------------------------
 
75
    Symbol           Meaning
 
76
    ------------------------------------------
 
77
    A        Alanine
 
78
    B        Aspartic Acid, Asparagine
 
79
    C        Cysteine
 
80
    D        Aspartic Acid
 
81
    E        Glutamic Acid
 
82
    F        Phenylalanine
 
83
    G        Glycine
 
84
    H        Histidine
 
85
    I        Isoleucine
 
86
    J        Isoleucine/Leucine
 
87
    K        Lysine
 
88
    L        Leucine
 
89
    M        Methionine
 
90
    N        Asparagine
 
91
    O        Pyrrolysine
 
92
    P        Proline
 
93
    Q        Glutamine
 
94
    R        Arginine
 
95
    S        Serine
 
96
    T        Threonine
 
97
    U        Selenocysteine
 
98
    V        Valine
 
99
    W        Tryptophan
 
100
    X        Unknown
 
101
    Y        Tyrosine
 
102
    Z        Glutamic Acid, Glutamine
 
103
    *        Terminator
 
104
 
 
105
There are a few things Bio::Tools::IUPAC can do for you:
 
106
 
 
107
=over
 
108
 
 
109
=item *
 
110
 
 
111
report the IUPAC mapping between ambiguous and non-ambiguous residues
 
112
 
 
113
=item *
 
114
 
 
115
produce a stream of all possible corresponding unambiguous Bio::Seq objects given
 
116
an ambiguous sequence object
 
117
 
 
118
=item *
 
119
 
 
120
convert an ambiguous sequence object to a corresponding regular expression
 
121
 
 
122
=back
100
123
 
101
124
=head1 FEEDBACK
102
125
 
140
163
=cut
141
164
 
142
165
 
143
 
# Let the code begin...
144
 
 
145
166
package Bio::Tools::IUPAC;
146
167
 
147
168
use strict;
148
 
use vars qw(%IUP %IUB %REV_IUB $AUTOLOAD);
 
169
use base qw(Bio::Root::Root);
 
170
use vars qw(%IUB %IUB_AMB %REV_IUB %IUP %IUP_AMB $AUTOLOAD);
149
171
 
150
172
BEGIN {
151
 
    %IUB = ( A => [qw(A)],
152
 
             C => [qw(C)],
153
 
             G => [qw(G)],
154
 
             T => [qw(T)],
155
 
             U => [qw(U)],
156
 
             M => [qw(A C)],
157
 
             R => [qw(A G)],
158
 
             W => [qw(A T)],
159
 
             S => [qw(C G)],
160
 
             Y => [qw(C T)],
161
 
             K => [qw(G T)],
162
 
             V => [qw(A C G)],
163
 
             H => [qw(A C T)],
164
 
             D => [qw(A G T)],
165
 
             B => [qw(C G T)],
166
 
             X => [qw(G A T C)],
167
 
             N => [qw(G A T C)]
168
 
             );
169
 
        %REV_IUB = (A   => 'A',
170
 
                                T       => 'T',
171
 
                                C       => 'C',
172
 
                                G       => 'G',
173
 
                                AC      => 'M',
174
 
                                AG      => 'R',
175
 
                                AT      => 'W',
176
 
                                CG      => 'S',
177
 
                                CT      => 'Y',
178
 
                                'GT'    => 'K',
179
 
                                ACG     => 'V',
180
 
                                ACT     => 'H',
181
 
                                AGT     => 'D',
182
 
                                CGT     => 'B',
183
 
                                ACGT=> 'N',
184
 
                                N       => 'N'
185
 
                                );
186
 
 
187
 
 
188
 
    %IUP = (A => [qw(A)],
189
 
            B => [qw(D N)],
190
 
            C => [qw(C)],
191
 
            D => [qw(D)],
192
 
            E => [qw(E)],
193
 
            F => [qw(F)],
194
 
            G => [qw(G)],
195
 
            H => [qw(H)],
196
 
            I => [qw(I)],
 
173
    # Ambiguous nucleic residues are matched to unambiguous residues
 
174
    %IUB = (
 
175
        A => [qw(A)],
 
176
        C => [qw(C)],
 
177
        G => [qw(G)],
 
178
        T => [qw(T)],
 
179
        U => [qw(U)],
 
180
        M => [qw(A C)],
 
181
        R => [qw(A G)],
 
182
        S => [qw(C G)],
 
183
        W => [qw(A T)],
 
184
        Y => [qw(C T)],
 
185
        K => [qw(G T)],
 
186
        V => [qw(A C G)],
 
187
        H => [qw(A C T)],
 
188
        D => [qw(A G T)],
 
189
        B => [qw(C G T)],
 
190
        N => [qw(A C G T)],
 
191
        X => [qw(A C G T)],
 
192
    );
 
193
 
 
194
    # Same as %IUB but ambiguous residues are matched to ambiguous residues only
 
195
    %IUB_AMB = (
 
196
        M => [qw(M)],
 
197
        R => [qw(R)],
 
198
        W => [qw(W)],
 
199
        S => [qw(S)],
 
200
        Y => [qw(Y)],
 
201
        K => [qw(K)],
 
202
        V => [qw(M R S V)],
 
203
        H => [qw(H M W Y)],
 
204
        D => [qw(D K R W)],
 
205
        B => [qw(B K S Y)],
 
206
        N => [qw(B D H K M N R S V W Y)],
 
207
    );
 
208
 
 
209
    # The inverse of %IUB
 
210
    %REV_IUB = (
 
211
        A    => 'A',
 
212
        T    => 'T',
 
213
        U    => 'U',
 
214
        C    => 'C',
 
215
        G    => 'G',
 
216
        AC   => 'M',
 
217
        AG   => 'R',
 
218
        AT   => 'W',
 
219
        CG   => 'S',
 
220
        CT   => 'Y',
 
221
        GT   => 'K',
 
222
        ACG  => 'V',
 
223
        ACT  => 'H',
 
224
        AGT  => 'D',
 
225
        CGT  => 'B',
 
226
        ACGT => 'N',
 
227
        N    => 'N'
 
228
    );
 
229
 
 
230
    # Same thing with proteins now
 
231
    %IUP = (
 
232
        A => [qw(A)],
 
233
        B => [qw(D N)],
 
234
        C => [qw(C)],
 
235
        D => [qw(D)],
 
236
        E => [qw(E)],
 
237
        F => [qw(F)],
 
238
        G => [qw(G)],
 
239
        H => [qw(H)],
 
240
        I => [qw(I)],
197
241
        J => [qw(I L)],
198
 
            K => [qw(K)],
199
 
            L => [qw(L)],
200
 
            M => [qw(M)],
201
 
            N => [qw(N)],
 
242
        K => [qw(K)],
 
243
        L => [qw(L)],
 
244
        M => [qw(M)],
 
245
        N => [qw(N)],
202
246
        O => [qw(O)],
203
 
            P => [qw(P)],
204
 
            Q => [qw(Q)],
205
 
            R => [qw(R)],
206
 
            S => [qw(S)],
207
 
            T => [qw(T)],
208
 
            U => [qw(U)],
209
 
            V => [qw(V)],
210
 
            W => [qw(W)],
211
 
            X => [qw(X)],
212
 
            Y => [qw(Y)],
213
 
            Z => [qw(E Q)],
214
 
            '*' => ['*']
215
 
            );
 
247
        P => [qw(P)],
 
248
        Q => [qw(Q)],
 
249
        R => [qw(R)],
 
250
        S => [qw(S)],
 
251
        T => [qw(T)],
 
252
        U => [qw(U)],
 
253
        V => [qw(V)],
 
254
        W => [qw(W)],
 
255
        X => [qw(X)],
 
256
        Y => [qw(Y)],
 
257
        Z => [qw(E Q)],
 
258
        '*' => [qw(*)],
 
259
    );
 
260
 
 
261
    %IUP_AMB = (
 
262
        B => [qw(B)],
 
263
        J => [qw(J)],
 
264
        Z => [qw(Z)],
 
265
    );
216
266
 
217
267
}
218
 
use base qw(Bio::Root::Root);
 
268
 
219
269
 
220
270
=head2 new
221
271
 
222
272
 Title   : new
223
 
 Usage   : Bio::Tools::IUPAC->new( $seq)
224
 
 Function: returns a new seq stream (akin to SeqIO)
225
 
 Returns : a Bio::Tools::IUPAC stream object that will produce unique
226
 
           Seq objects on demand.
227
 
 Args    : an ambiguously coded Seq.pm object that has a specified 'alphabet'
228
 
 
 
273
 Usage   : Bio::Tools::IUPAC->new($seq);
 
274
 Function: Create a new IUPAC object, which acts as a sequence stream (akin to
 
275
           SeqIO)
 
276
 Args    : an ambiguously coded sequence object that has a specified 'alphabet'
 
277
 Returns : a Bio::Tools::IUPAC object.
229
278
 
230
279
=cut
231
280
 
232
 
 
233
281
sub new {
234
 
    my($class,@args) = @_;
 
282
    my ($class,@args) = @_;
235
283
    my $self = $class->SUPER::new(@args);
236
 
 
237
284
    my ($seq) = $self->_rearrange([qw(SEQ)],@args);
238
 
    if((! defined($seq)) && @args && ref($args[0])) {
239
 
        # parameter not passed as named parameter?
240
 
        $seq = $args[0];
241
 
    }
242
 
    $seq->isa('Bio::Seq') or
243
 
        $self->throw("Must supply a Seq.pm object to IUPAC!");
244
 
    $self->{'_SeqObj'} = $seq;
245
 
    if ($self->{'_SeqObj'}->alphabet() =~ m/^[dr]na$/i ) {
246
 
        # nucleotide seq object
247
 
        $self->{'_alpha'} = [ map { $IUB{uc($_)} }
248
 
                              split('', $self->{'_SeqObj'}->seq()) ];
249
 
    } elsif ($self->{'_SeqObj'}->alphabet() =~ m/^protein$/i ) {
250
 
        # amino acid seq object
251
 
        $self->{'_alpha'} = [ map { $IUP{uc($_)} }
252
 
                               split('', $self->{'_SeqObj'}->seq()) ];
253
 
    } else { # unknown type: we could make a guess, but let's not.
254
 
        $self->throw("You must specify the 'type' of sequence provided to IUPAC");
255
 
    }
256
 
    $self->{'_string'} = [(0) x length($self->{'_SeqObj'}->seq())];
257
 
    scalar @{$self->{'_string'}} or $self->throw("Sequence has zero-length!");
 
285
 
 
286
    if ( (not defined $seq) && @args && ref($args[0]) ) {
 
287
        # parameter not passed as named parameter?
 
288
        $seq = $args[0];
 
289
    }
 
290
 
 
291
    if (defined $seq) {
 
292
        if (not $seq->isa('Bio::PrimarySeqI')) {
 
293
            $self->throw('Must supply a sequence object');
 
294
        }
 
295
        if (length $seq->seq == 0) {
 
296
            $self->throw('Sequence had zero-length');
 
297
        }
 
298
        $self->{'_seq'} = $seq;
 
299
    }
 
300
 
 
301
    return $self;
 
302
}
 
303
 
 
304
 
 
305
sub _initialize {
 
306
    my ($self) = @_;
 
307
    my %iupac = $self->iupac;
 
308
    $self->{'_alpha'} = [ map { $iupac{uc $_} } split('', $self->{'_seq'}->seq) ];
 
309
    $self->{'_string'} = [(0) x length($self->{'_seq'}->seq())];
258
310
    $self->{'_string'}->[0] = -1;
259
 
    return $self;
260
311
}
261
312
 
 
313
 
262
314
=head2 next_seq
263
315
 
264
316
 Title   : next_seq
265
 
 Usage   : $iupac->next_seq()
266
 
 Function: returns the next unique Seq object
267
 
 Returns : a Seq.pm object
 
317
 Usage   : $iupac->next_seq();
 
318
 Function: returns the next unique sequence object
268
319
 Args    : none.
269
 
 
 
320
 Returns : a Bio::Seq object
270
321
 
271
322
=cut
272
323
 
273
 
sub next_seq{
 
324
sub next_seq {
274
325
    my ($self) = @_;
275
326
 
 
327
    if (not exists $self->{'_string'}) {
 
328
        $self->_initialize();
 
329
    }
 
330
 
276
331
    for my $i ( 0 .. $#{$self->{'_string'}} ) {
277
 
        next unless $self->{'_string'}->[$i] || @{$self->{'_alpha'}->[$i]} > 1;
278
 
        if ( $self->{'_string'}->[$i] == $#{$self->{'_alpha'}->[$i]} ) { # rollover
279
 
            if ( $i == $#{$self->{'_string'}} ) { # end of possibilities
280
 
                return;
281
 
            } else {
282
 
                $self->{'_string'}->[$i] = 0;
283
 
                next;
284
 
            }
285
 
        } else {
286
 
            $self->{'_string'}->[$i]++;
287
 
            my $j = -1;
288
 
            $self->{'_SeqObj'}->seq(join('', map { $j++; $self->{'_alpha'}->[$j]->[$_]; } @{$self->{'_string'}}));
289
 
            my $desc = $self->{'_SeqObj'}->desc();
290
 
            if ( !defined $desc ) { $desc = ""; }
291
 
 
292
 
            $self->{'_num'}++;
293
 
            1 while $self->{'_num'} =~ s/(\d)(\d\d\d)(?!\d)/$1,$2/;
294
 
            $desc =~ s/( \[Bio::Tools::IUPAC-generated\sunique sequence # [^\]]*\])|$/ \[Bio::Tools::IUPAC-generated unique sequence # $self->{'_num'}\]/;
295
 
            $self->{'_SeqObj'}->desc($desc);
296
 
            $self->{'_num'} =~ s/,//g;
297
 
            return $self->{'_SeqObj'};
298
 
        }
299
 
    }
300
 
}
 
332
        next unless $self->{'_string'}->[$i] || @{$self->{'_alpha'}->[$i]} > 1;
 
333
        if ( $self->{'_string'}->[$i] == $#{$self->{'_alpha'}->[$i]} ) { # rollover
 
334
            if ( $i == $#{$self->{'_string'}} ) { # end of possibilities
 
335
                return;
 
336
            } else {
 
337
                $self->{'_string'}->[$i] = 0;
 
338
                next;
 
339
            }
 
340
        } else {
 
341
            $self->{'_string'}->[$i]++;
 
342
            my $j = -1;
 
343
            my $seqstr = join('', map { $j++; $self->{'_alpha'}->[$j]->[$_]; } @{$self->{'_string'}});
 
344
            my $desc   = $self->{'_seq'}->desc() || '';
 
345
            $self->{'_num'}++;
 
346
            1 while $self->{'_num'} =~ s/(\d)(\d\d\d)(?!\d)/$1,$2/;
 
347
            $desc =~ s/( \[Bio::Tools::IUPAC-generated\sunique sequence # [^\]]*\])|$/ \[Bio::Tools::IUPAC-generated unique sequence # $self->{'_num'}\]/;
 
348
            $self->{'_num'} =~ s/,//g;
 
349
 
 
350
            # Return a fresh sequence object
 
351
            return Bio::PrimarySeq->new(-seq  => $seqstr, -desc => $desc);
 
352
        }
 
353
    }
 
354
}
 
355
 
 
356
 
 
357
=head2 iupac
 
358
 
 
359
 Title   : iupac
 
360
 Usage   : my %symbols = $iupac->iupac;
 
361
 Function: Returns a hash of symbols -> symbol components of the right type
 
362
           for the given sequence, i.e. it is the same as iupac_iup() if
 
363
           Bio::Tools::IUPAC was given a proteic sequence, or iupac_iub() if the 
 
364
           sequence was nucleic. For example, the key 'M' has the value ['A', 'C'].
 
365
 Args    : none
 
366
 Returns : Hash
 
367
 
 
368
=cut
 
369
 
 
370
sub iupac {
 
371
    my ($self) = @_;
 
372
    my $alphabet = lc( $self->{'_seq'}->alphabet() );
 
373
    if ( ($alphabet eq 'dna') or ($alphabet eq 'rna') ) {
 
374
        return %IUB; # nucleic
 
375
    } elsif ( $alphabet eq 'protein' ) {
 
376
        return %IUP; # proteic
 
377
    } else {
 
378
        $self->throw("The input sequence had the unknown alphabet '$alphabet'\n");
 
379
    }
 
380
}
 
381
 
 
382
 
 
383
 
 
384
=head2 iupac_amb
 
385
 
 
386
 Title   : iupac_amb
 
387
 Usage   : my %symbols = $iupac->iupac_amb;
 
388
 Function: Same as iupac() but only contains a mapping between ambiguous residues
 
389
           and the ambiguous residues they map to. For example, the key 'N' has
 
390
           the value ['R', 'Y', 'K', 'M', 'S', 'W', 'B', 'D', 'H', 'V', 'N'],
 
391
           i.e. it matches all other ambiguous residues.
 
392
 Args    : none
 
393
 Returns : Hash
 
394
 
 
395
=cut
 
396
 
 
397
sub iupac_amb {
 
398
    my ($self) = @_;
 
399
    my $alphabet = lc( $self->{'_seq'}->alphabet() );
 
400
    if ( ($alphabet eq 'dna') or ($alphabet eq 'rna') ) {
 
401
        return %IUB_AMB; # nucleic
 
402
    } elsif ( $alphabet eq 'protein' ) {
 
403
        return %IUP_AMB; # proteic
 
404
    } else {
 
405
        $self->throw("The input sequence had the unknown alphabet '$alphabet'\n");
 
406
    }
 
407
}
 
408
 
301
409
 
302
410
=head2 iupac_iup
303
411
 
304
412
 Title   : iupac_iup
305
 
 Usage   : my %aasymbols = $iupac->iupac_iup
306
 
 Function: Returns a hash of PROTEIN symbols -> symbol components
 
413
 Usage   : my %aasymbols = $iupac->iupac_iup;
 
414
 Function: Returns a hash of PROTEIN symbols -> non-ambiguous symbol components
 
415
 Args    : none
307
416
 Returns : Hash
308
 
 Args    : none
309
417
 
310
418
=cut
311
419
 
312
 
sub iupac_iup{
 
420
sub iupac_iup {
313
421
   return %IUP;
314
 
 
315
 
}
 
422
}
 
423
 
 
424
 
 
425
=head2 iupac_iup_amb
 
426
 
 
427
 Title   : iupac_iup_amb
 
428
 Usage   : my %aasymbols = $iupac->iupac_iup_amb;
 
429
 Function: Returns a hash of PROTEIN symbols -> ambiguous symbol components
 
430
 Args    : none
 
431
 Returns : Hash
 
432
 
 
433
=cut
 
434
 
 
435
sub iupac_iup_amb {
 
436
   return %IUP_AMB;
 
437
}
 
438
 
316
439
 
317
440
=head2 iupac_iub
318
441
 
319
442
 Title   : iupac_iub
320
 
 Usage   : my %dnasymbols = $iupac->iupac_iub
321
 
 Function: Returns a hash of DNA symbols -> symbol components
 
443
 Usage   : my %dnasymbols = $iupac->iupac_iub;
 
444
 Function: Returns a hash of DNA symbols -> non-ambiguous symbol components
 
445
 Args    : none
322
446
 Returns : Hash
323
 
 Args    : none
324
447
 
325
448
=cut
326
449
 
327
 
sub iupac_iub{
 
450
sub iupac_iub {
328
451
   return %IUB;
329
452
}
330
453
 
 
454
 
 
455
=head2 iupac_iub_amb
 
456
 
 
457
 Title   : iupac_iub_amb
 
458
 Usage   : my %dnasymbols = $iupac->iupac_iub;
 
459
 Function: Returns a hash of DNA symbols -> ambiguous symbol components
 
460
 Args    : none
 
461
 Returns : Hash
 
462
 
 
463
=cut
 
464
 
 
465
sub iupac_iub_amb {
 
466
   return %IUB_AMB;
 
467
}
 
468
 
 
469
 
331
470
=head2 iupac_rev_iub
332
471
 
333
472
 Title   : iupac_rev_iub
334
 
 Usage   : my %dnasymbols = $iupac->iupac_rev_iub
 
473
 Usage   : my %dnasymbols = $iupac->iupac_rev_iub;
335
474
 Function: Returns a hash of nucleotide combinations -> IUPAC code
336
475
           (a reverse of the iupac_iub hash).
 
476
 Args    : none
337
477
 Returns : Hash
338
 
 Args    : none
339
478
 
340
479
=cut
341
480
 
342
 
sub iupac_rev_iub{
 
481
sub iupac_rev_iub {
343
482
   return %REV_IUB;
344
483
}
345
484
 
 
485
 
346
486
=head2 count
347
487
 
348
488
 Title   : count
349
489
 Usage   : my $total = $iupac->count();
350
490
 Function: Calculates the number of unique, unambiguous sequences that
351
491
           this ambiguous sequence could generate
 
492
 Args    : none
352
493
 Return  : int
353
 
 Args    : none
354
494
 
355
495
=cut
356
496
 
357
497
sub count {
358
498
    my ($self) = @_;
359
 
 
 
499
    if (not exists $self->{'_string'}) {
 
500
        $self->_initialize();
 
501
    }
360
502
    my $count = 1;
361
503
    $count *= scalar(@$_) for (@{$self->{'_alpha'}});
362
504
    return $count;
363
505
}
364
506
 
365
507
 
 
508
=head2 regexp
 
509
 
 
510
 Title   : regexp
 
511
 Usage   : my $re = $iupac->regexp();
 
512
 Function: Converts the ambiguous sequence into a regular expression that
 
513
           matches all of the corresponding ambiguous and non-ambiguous sequences.
 
514
           You can further manipulate the resulting regular expression with the
 
515
           Bio::Tools::SeqPattern module. After you are done building your
 
516
           regular expression, you might want to compile it and make it case-
 
517
           insensitive:
 
518
              $re = qr/$re/i;
 
519
 Args    : 1 to match RNA: T and U characters will match interchangeably
 
520
 Return  : regular expression
 
521
 
 
522
=cut
 
523
 
 
524
sub regexp {
 
525
    my ($self, $match_rna) = @_;
 
526
    my $re;
 
527
    my $seq = $self->{'_seq'}->seq;
 
528
    my %iupac = $self->iupac;
 
529
    my %iupac_amb = $self->iupac_amb;
 
530
    for my $pos (0 .. length($seq)-1) {
 
531
        my $res = substr $seq, $pos, 1;
 
532
        my $iupacs = $iupac{$res};
 
533
        my $iupacs_amb = $iupac_amb{$res} || [];
 
534
        if (not defined $iupacs) {
 
535
            $self->throw("Primer sequence '$seq' is not a valid IUPAC sequence.".
 
536
                         " Offending character was '$res'.\n");
 
537
        }
 
538
        my $part = join '', (@$iupacs, @$iupacs_amb);
 
539
        if ($match_rna) {
 
540
            $part =~ s/T/TU/i || $part =~ s/U/TU/i;
 
541
        }
 
542
        if (length $part > 1) {
 
543
           $part = '['.$part.']';
 
544
        }
 
545
        $re .= $part;
 
546
    }
 
547
    return $re;
 
548
}
 
549
 
 
550
 
366
551
sub AUTOLOAD {
367
 
 
368
552
    my $self = shift @_;
369
553
    my $method = $AUTOLOAD;
370
554
    $method =~ s/.*:://;
371
 
    return $self->{'_SeqObj'}->$method(@_)
372
 
        unless $method eq 'DESTROY';
 
555
    return $self->{'_seq'}->$method(@_)
 
556
        unless $method eq 'DESTROY';
373
557
}
374
558
 
375
559
1;