Bio/Tools/IUPAC.pm

   1 # $Id$
   2 #
   3 # BioPerl module for IUPAC
   4 #
   5 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   6 #
   7 # Cared for by Aaron Mackey <amackey@virginia.edu>
   8 #
   9 # Copyright Aaron Mackey
  10 #
  11 # You may distribute this module under the same terms as perl itself
  12
  13 # POD documentation - main docs before the code
  14
  15 =head1 NAME
  16
  17 Bio::Tools::IUPAC - Generates unique Seq objects from an ambiguous Seq object
  18
  19 =head1 SYNOPSIS
  20
  21  use Bio::Seq;
  22  use Bio::Tools::IUPAC;
  23
  24  my $ambiseq = Bio::Seq->new(-seq => 'ARTCGUTGR', -alphabet => 'dna');
  25  my $stream  = Bio::Tools::IUPAC->new(-seq => $ambiseq);
  26
  27  while ($uniqueseq = $stream->next_seq()) {
  28      # process the unique Seq object.
  29  }
  30
  31 =head1 DESCRIPTION
  32
  33 IUPAC is a tool that produces a stream of unique, "strict"-satisfying Seq
  34 objects from an ambiquous Seq object (containing non-standard characters given
  35 the meaning shown below)
  36
  37         Extended DNA / RNA alphabet :
  38         (includes symbols for nucleotide ambiguity)
  39         ------------------------------------------
  40         Symbol       Meaning      Nucleic Acid
  41         ------------------------------------------
  42          A            A           Adenine
  43          C            C           Cytosine
  44          G            G           Guanine
  45          T            T           Thymine
  46          U            U           Uracil
  47          M          A or C
  48          R          A or G
  49          W          A or T
  50          S          C or G
  51          Y          C or T
  52          K          G or T
  53          V        A or C or G
  54          H        A or C or T
  55          D        A or G or T
  56          B        C or G or T
  57          X      G or A or T or C
  58          N      G or A or T or C
  59
  60         IUPAC-IUB SYMBOLS FOR NUCLEOTIDE NOMENCLATURE:
  61           Cornish-Bowden (1985) Nucl. Acids Res. 13: 3021-3030.
  62
  63 -----------------------------------
  64
  65        Amino Acid alphabet:
  66         ------------------------------------------
  67         Symbol           Meaning
  68         ------------------------------------------
  69         A        Alanine
  70         B        Aspartic Acid, Asparagine
  71         C        Cystine
  72         D        Aspartic Acid
  73         E        Glutamic Acid
  74         F        Phenylalanine
  75         G        Glycine
  76         H        Histidine
  77         I        Isoleucine
  78         J        Isoleucine/Leucine
  79         K        Lysine
  80         L        Leucine
  81         M        Methionine
  82         N        Asparagine
  83         O        Pyrrolysine
  84         P        Proline
  85         Q        Glutamine
  86         R        Arginine
  87         S        Serine
  88         T        Threonine
  89         U        Selenocysteine
  90         V        Valine
  91         W        Tryptophan
  92         X        Unknown
  93         Y        Tyrosine
  94         Z        Glutamic Acid, Glutamine
  95         *        Terminator
  96
  97
  98         IUPAC-IUP AMINO ACID SYMBOLS:
  99           Biochem J. 1984 Apr 15; 219(2): 345-373
 100           Eur J Biochem. 1993 Apr 1; 213(1): 2
 101
 102 =head1 FEEDBACK
 103
 104 =head2 Mailing Lists
 105
 106 User feedback is an integral part of the evolution of this and other
 107 Bioperl modules. Send your comments and suggestions preferably to one
 108 of the Bioperl mailing lists.  Your participation is much appreciated.
 109
 110   bioperl-l@bioperl.org                  - General discussion
 111   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 112
 113 =head2 Support
 114
 115 Please direct usage questions or support issues to the mailing list:
 116
 117 L<bioperl-l@bioperl.org>
 118
 119 rather than to the module maintainer directly. Many experienced and
 120 reponsive experts will be able look at the problem and quickly
 121 address it. Please include a thorough description of the problem
 122 with code and data examples if at all possible.
 123
 124 =head2 Reporting Bugs
 125
 126 Report bugs to the Bioperl bug tracking system to help us keep track
 127 the bugs and their resolution.  Bug reports can be submitted via the
 128 web:
 129
 130   http://bugzilla.open-bio.org/
 131
 132 =head1 AUTHOR - Aaron Mackey
 133
 134 Email amackey-at-virginia.edu
 135
 136 =head1 APPENDIX
 137
 138 The rest of the documentation details each of the object
 139 methods. Internal methods are usually preceded with a _
 140
 141 =cut
 142
 143
 144 # Let the code begin...
 145
 146 package Bio::Tools::IUPAC;
 147
 148 use strict;
 149 use vars qw(%IUP %IUB %REV_IUB $AUTOLOAD);
 150
 151 BEGIN {
 152     %IUB = ( A => [qw(A)],
 153              C => [qw(C)],
 154              G => [qw(G)],
 155              T => [qw(T)],
 156              U => [qw(U)],
 157              M => [qw(A C)],
 158              R => [qw(A G)],
 159              W => [qw(A T)],
 160              S => [qw(C G)],
 161              Y => [qw(C T)],
 162              K => [qw(G T)],
 163              V => [qw(A C G)],
 164              H => [qw(A C T)],
 165              D => [qw(A G T)],
 166              B => [qw(C G T)],
 167              X => [qw(G A T C)],
 168              N => [qw(G A T C)]
 169              );
 170         %REV_IUB = (A   => 'A',
 171                                 T       => 'T',
 172                                 C       => 'C',
 173                                 G       => 'G',
 174                                 AC      => 'M',
 175                                 AG      => 'R',
 176                                 AT      => 'W',
 177                                 CG      => 'S',
 178                                 CT      => 'Y',
 179                                 'GT'    => 'K',
 180                                 ACG     => 'V',
 181                                 ACT     => 'H',
 182                                 AGT     => 'D',
 183                                 CGT     => 'B',
 184                                 ACGT=> 'N',
 185                                 N       => 'N'
 186                                 );
 187
 188
 189     %IUP = (A => [qw(A)],
 190             B => [qw(D N)],
 191             C => [qw(C)],
 192             D => [qw(D)],
 193             E => [qw(E)],
 194             F => [qw(F)],
 195             G => [qw(G)],
 196             H => [qw(H)],
 197             I => [qw(I)],
 198         J => [qw(I L)],
 199             K => [qw(K)],
 200             L => [qw(L)],
 201             M => [qw(M)],
 202             N => [qw(N)],
 203         O => [qw(O)],
 204             P => [qw(P)],
 205             Q => [qw(Q)],
 206             R => [qw(R)],
 207             S => [qw(S)],
 208             T => [qw(T)],
 209             U => [qw(U)],
 210             V => [qw(V)],
 211             W => [qw(W)],
 212             X => [qw(X)],
 213             Y => [qw(Y)],
 214             Z => [qw(E Q)],
 215             '*' => ['*']
 216             );
 217
 218 }
 219 use base qw(Bio::Root::Root);
 220
 221 =head2 new
 222
 223  Title   : new
 224  Usage   : Bio::Tools::IUPAC->new( $seq)
 225  Function: returns a new seq stream (akin to SeqIO)
 226  Returns : a Bio::Tools::IUPAC stream object that will produce unique
 227            Seq objects on demand.
 228  Args    : an ambiguously coded Seq.pm object that has a specified 'alphabet'
 229
 230
 231 =cut
 232
 233
 234 sub new {
 235     my($class,@args) = @_;
 236     my $self = $class->SUPER::new(@args);
 237
 238     my ($seq) = $self->_rearrange([qw(SEQ)],@args);
 239     if((! defined($seq)) && @args && ref($args[0])) {
 240         # parameter not passed as named parameter?
 241         $seq = $args[0];
 242     }
 243     $seq->isa('Bio::Seq') or
 244         $self->throw("Must supply a Seq.pm object to IUPAC!");
 245     $self->{'_SeqObj'} = $seq;
 246     if ($self->{'_SeqObj'}->alphabet() =~ m/^[dr]na$/i ) {
 247         # nucleotide seq object
 248         $self->{'_alpha'} = [ map { $IUB{uc($_)} }
 249                               split('', $self->{'_SeqObj'}->seq()) ];
 250     } elsif ($self->{'_SeqObj'}->alphabet() =~ m/^protein$/i ) {
 251         # amino acid seq object
 252         $self->{'_alpha'} = [ map { $IUP{uc($_)} }
 253                                split('', $self->{'_SeqObj'}->seq()) ];
 254     } else { # unknown type: we could make a guess, but let's not.
 255         $self->throw("You must specify the 'type' of sequence provided to IUPAC");
 256     }
 257     $self->{'_string'} = [(0) x length($self->{'_SeqObj'}->seq())];
 258     scalar @{$self->{'_string'}} or $self->throw("Sequence has zero-length!");
 259     $self->{'_string'}->[0] = -1;
 260     return $self;
 261 }
 262
 263 =head2 next_seq
 264
 265  Title   : next_seq
 266  Usage   : $iupac->next_seq()
 267  Function: returns the next unique Seq object
 268  Returns : a Seq.pm object
 269  Args    : none.
 270
 271
 272 =cut
 273
 274 sub next_seq{
 275     my ($self) = @_;
 276
 277     for my $i ( 0 .. $#{$self->{'_string'}} ) {
 278         next unless $self->{'_string'}->[$i] || @{$self->{'_alpha'}->[$i]} > 1;
 279         if ( $self->{'_string'}->[$i] == $#{$self->{'_alpha'}->[$i]} ) { # rollover
 280             if ( $i == $#{$self->{'_string'}} ) { # end of possibilities
 281                 return;
 282             } else {
 283                 $self->{'_string'}->[$i] = 0;
 284                 next;
 285             }
 286         } else {
 287             $self->{'_string'}->[$i]++;
 288             my $j = -1;
 289             $self->{'_SeqObj'}->seq(join('', map { $j++; $self->{'_alpha'}->[$j]->[$_]; } @{$self->{'_string'}}));
 290             my $desc = $self->{'_SeqObj'}->desc();
 291             if ( !defined $desc ) { $desc = ""; }
 292
 293             $self->{'_num'}++;
 294             1 while $self->{'_num'} =~ s/(\d)(\d\d\d)(?!\d)/$1,$2/;
 295             $desc =~ s/( \[Bio::Tools::IUPAC-generated\sunique sequence # [^\]]*\])|$/ \[Bio::Tools::IUPAC-generated unique sequence # $self->{'_num'}\]/;
 296             $self->{'_SeqObj'}->desc($desc);
 297             $self->{'_num'} =~ s/,//g;
 298             return $self->{'_SeqObj'};
 299         }
 300     }
 301 }
 302
 303 =head2 iupac_iup
 304
 305  Title   : iupac_iup
 306  Usage   : my %aasymbols = $iupac->iupac_iup
 307  Function: Returns a hash of PROTEIN symbols -> symbol components
 308  Returns : Hash
 309  Args    : none
 310
 311 =cut
 312
 313 sub iupac_iup{
 314    return %IUP;
 315
 316 }
 317
 318 =head2 iupac_iub
 319
 320  Title   : iupac_iub
 321  Usage   : my %dnasymbols = $iupac->iupac_iub
 322  Function: Returns a hash of DNA symbols -> symbol components
 323  Returns : Hash
 324  Args    : none
 325
 326 =cut
 327
 328 sub iupac_iub{
 329    return %IUB;
 330 }
 331
 332 =head2 iupac_rev_iub
 333
 334  Title   : iupac_rev_iub
 335  Usage   : my %dnasymbols = $iupac->iupac_rev_iub
 336  Function: Returns a hash of nucleotide combinations -> IUPAC code
 337            (a reverse of the iupac_iub hash).
 338  Returns : Hash
 339  Args    : none
 340
 341 =cut
 342
 343 sub iupac_rev_iub{
 344    return %REV_IUB;
 345 }
 346
 347 =head2 count
 348
 349  Title   : count
 350  Usage   : my $total = $iupac->count();
 351  Function: Calculates the number of unique, unambiguous sequences that
 352            this ambiguous sequence could generate
 353  Return  : int
 354  Args    : none
 355
 356 =cut
 357
 358 sub count {
 359     my ($self) = @_;
 360
 361     my $count = 1;
 362     $count *= scalar(@$_) for (@{$self->{'_alpha'}});
 363     return $count;
 364 }
 365
 366
 367 sub AUTOLOAD {
 368
 369     my $self = shift @_;
 370     my $method = $AUTOLOAD;
 371     $method =~ s/.*:://;
 372     return $self->{'_SeqObj'}->$method(@_)
 373         unless $method eq 'DESTROY';
 374 }
 375
 376 1;
 377