Bio/SeqEvolution/Factory.pm

   1 #
   2 # BioPerl module for Bio::SeqEvolution::Factory
   3 #
   4 # Please direct questions and support issues to <bioperl-l@bioperl.org>
   5 #
   6 # Cared for by Heikki Lehvaslaiho <heikki at bioperl dot org>
   7 #
   8 # Copyright Heikki Lehvaslaiho
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs before the code
  13
  14 =head1 NAME
  15
  16 Bio::SeqEvolution::Factory - Factory object to instantiate sequence evolving classes
  17
  18 =head1 SYNOPSIS
  19
  20     # not an instantiable class
  21
  22 =head1 DESCRIPTION
  23
  24 This is the factory class that can be used to call for a specific
  25 model to mutate a sequence.
  26
  27 Bio::SeqEvolution::DNAPoint is the default for nucleotide sequences
  28 and the only implementation at this point.
  29
  30 =head1 FEEDBACK
  31
  32 =head2 Mailing Lists
  33
  34 User feedback is an integral part of the evolution of this and other
  35 Bioperl modules. Send your comments and suggestions preferably to
  36 the Bioperl mailing list.  Your participation is much appreciated.
  37
  38   bioperl-l@bioperl.org                  - General discussion
  39   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  40
  41 =head2 Support
  42
  43 Please direct usage questions or support issues to the mailing list:
  44
  45 I<bioperl-l@bioperl.org>
  46
  47 rather than to the module maintainer directly. Many experienced and
  48 reponsive experts will be able look at the problem and quickly
  49 address it. Please include a thorough description of the problem
  50 with code and data examples if at all possible.
  51
  52 =head2 Reporting Bugs
  53
  54 Report bugs to the Bioperl bug tracking system to help us keep track
  55 of the bugs and their resolution. Bug reports can be submitted via the
  56 web:
  57
  58   https://github.com/bioperl/bioperl-live/issues
  59
  60 =head1 AUTHOR
  61
  62   Heikki Lehvaslaiho E<lt>heikki at bioperl dot orgE<gt>
  63
  64 =head1 CONTRIBUTORS
  65
  66 Additional contributor's names and emails here
  67
  68 =head1 APPENDIX
  69
  70 The rest of the documentation details each of the object methods.
  71 Internal methods are usually preceded with a _
  72
  73 =cut
  74
  75
  76 # Let the code begin...
  77
  78
  79 package Bio::SeqEvolution::Factory;
  80 use strict;
  81 use Bio::Root::Root;
  82 use Bio::SeqEvolution::EvolutionI;
  83 use base qw(Bio::Root::Root Bio::SeqEvolution::EvolutionI);
  84
  85 =head2 new
  86
  87   Title   : new
  88   Usage   : my $obj = Bio::SeqEvolution::Factory->new();
  89   Function: Builds a new Bio:SeqEvolution::EvolutionI object
  90   Returns : Bio:SeqEvolution::EvolutionI object
  91   Args    : -type           => class name
  92
  93 See L<Bio:SeqEvolution::EvolutionI>
  94
  95 =cut
  96
  97 sub new {
  98     my($caller,@args) = @_;
  99     my $class = ref($caller) || $caller;
 100
 101     my %param = @args;
 102     @param{ map { lc $_ } keys %param } = values %param; # lowercase keys
 103
 104     if ( $class eq 'Bio::SeqEvolution::Factory') {
 105         #my %param = @args;
 106         #@param{ map { lc $_ } keys %param } = values %param; # lowercase keys
 107
 108         if (exists $param{'-type'}) {
 109 #            $self->type($param{'-type'});
 110         } else {
 111             $param{'-type'} = 'Bio::SeqEvolution::DNAPoint';
 112             #$self->type('Bio::SeqEvolution::DNAPoint'} unless $seq->alphabet == 'protein'
 113         }
 114         my $type = $param{'-type'};
 115         return unless( $class->_load_format_module($param{'-type'}) );
 116         return $type->new(%param);
 117     } else {
 118         my ($self) = $class->SUPER::new(%param);
 119         $self->_initialize(%param);
 120         return $self;
 121     }
 122 }
 123
 124 sub _initialize {
 125     my($self, @args) = @_;
 126
 127     $self->SUPER::_initialize(@args);
 128     my %param = @args;
 129     @param{ map { lc $_ } keys %param } = values %param; # lowercase keys
 130
 131     exists $param{'-seq'} && $self->seq($param{'-seq'});
 132     exists $param{'-set_mutated_seq'} && $self->set_mutated_seq($param{'-set_mutated_seq'});
 133     exists $param{'-identity'} && $self->identity($param{'-identity'});
 134     exists $param{'-pam'} && $self->pam($param{'-pam'});
 135     exists $param{'-mutation_count'} && $self->mutation_count($param{'-mutation_count'});
 136
 137 }
 138
 139
 140 =head2 _load_format_module
 141
 142  Title   : _load_format_module
 143  Usage   : *INTERNAL SeqIO stuff*
 144  Function: Loads up (like use) a module at run time on demand
 145  Example :
 146  Returns :
 147  Args    :
 148
 149 =cut
 150
 151 sub _load_format_module {
 152         my ($self, $format) = @_;
 153         my $module =  $format;
 154         my $ok;
 155
 156         eval {
 157                 $ok = $self->_load_module($module);
 158         };
 159         if ( $@ ) {
 160                 print STDERR <<END;
 161 $self: $format cannot be found
 162 Exception $@
 163 END
 164                 ;
 165         }
 166         return $ok;
 167 }
 168
 169
 170 =head2 type
 171
 172  Title   : type
 173  Usage   : $obj->type($newval)
 174  Function: Set used evolution model. It is set by giving a
 175            valid Bio::SeqEvolution::* class name
 176  Returns : value of type
 177  Args    : newvalue (optional)
 178
 179 Defaults to Bio::SeqEvolution::DNAPoint.
 180
 181 =cut
 182
 183 sub type{
 184    my $self = shift;
 185    if (@_) {
 186        $self->{'_type'} = shift @_;
 187        $self->_load_module($self->{'_type'});
 188    }
 189    return $self->{'_type'} || 'Bio::SeqEvolution::DNAPoint';
 190 }
 191
 192 =head1 mutation counters
 193
 194 The next three methods set a value to limit the number of mutations
 195 introduced the the input sequence.
 196
 197 =cut
 198
 199 =head2 identity
 200
 201  Title   : identity
 202  Usage   : $obj->identity($newval)
 203  Function: Set the desired identity between original and mutated sequence
 204  Returns : value of identity
 205  Args    : newvalue (optional)
 206
 207 =cut
 208
 209 sub identity{
 210    my $self = shift;
 211    $self->{'_identity'} = shift @_ if @_;
 212    return $self->{'_identity'};
 213 }
 214
 215
 216 =head2 pam
 217
 218  Title   : pam
 219  Usage   : $obj->pam($newval)
 220  Function: Set the wanted Percentage of Accepted Mutations, PAM
 221  Returns : value of PAM
 222  Args    : newvalue (optional)
 223
 224 When you are measuring sequence divergence, PAM needs to be
 225 estimated. When you are generating sequences, PAM is simply the count
 226 of mutations introduced to the reference sequence normalised to the
 227 original sequence length.
 228
 229 =cut
 230
 231 sub pam{
 232    my $self = shift;
 233    $self->{'_pam'} = shift @_ if @_;
 234    return $self->{'_pam'};
 235 }
 236
 237 =head2 mutation_count
 238
 239  Title   : mutation_count
 240  Usage   : $obj->mutation_count($newval)
 241  Function: Set the number of wanted mutations to the sequence
 242  Returns : value of mutation_count
 243  Args    : newvalue (optional)
 244
 245 =cut
 246
 247 sub mutation_count{
 248    my $self = shift;
 249    $self->{'_mutation_count'} = shift @_ if @_;
 250    return $self->{'_mutation_count'};
 251 }
 252
 253
 254
 255 =head2 seq
 256
 257  Title   : seq
 258  Usage   : $obj->seq($newval)
 259  Function: Set the sequence object for the original sequence
 260  Returns : The sequence object
 261  Args    : newvalue (optional)
 262
 263 Setting this will reset mutation and generated mutation counters.
 264
 265 =cut
 266
 267 sub seq {
 268    my $self = shift;
 269    if (@_) {
 270        $self->{'_seq'} = shift @_ ;
 271        return $self->{'_seq'};
 272        $self->reset_mutation_counter;
 273        $self->reset_sequence_counter;
 274    }
 275    return $self->{'_seq'};
 276 }
 277
 278 =head2 seq_type
 279
 280  Title   : seq_type
 281  Usage   : $obj->seq_type($newval)
 282  Function: Set the returned seq_type to one needed
 283  Returns : value of seq_type
 284  Args    : newvalue (optional)
 285
 286 Defaults to Bio::PrimarySeq.
 287
 288 =cut
 289
 290 sub seq_type{
 291    my $self = shift;
 292    if (@_) {
 293        $self->{'_seq_type'} = shift @_;
 294        $self->_load_module($self->{'_seq_type'});
 295    }
 296    return $self->{'_seq_type'} || 'Bio::PrimarySeq';
 297 }
 298
 299
 300 =head2 get_mutation_counter
 301
 302  Title   : get_mutation_counter
 303  Usage   : $obj->get_mutation_counter()
 304  Function: Get the count of sequences created
 305  Returns : value of counter
 306  Args    : -
 307
 308 =cut
 309
 310 sub get_mutation_counter{
 311    return shift->{'_mutation_counter'};
 312 }
 313
 314
 315 =head2 reset_mutation_counter
 316
 317  Title   : reset_mutation_counter
 318  Usage   : $obj->reset_mutation_counter()
 319  Function: Resert the counter of mutations
 320  Returns : value of counter
 321  Args    : -
 322
 323 =cut
 324
 325 sub reset_mutation_counter{
 326    shift->{'_mutation_counter'} = 0;
 327    return 1;
 328 }
 329
 330
 331 =head2 get_sequence_counter
 332
 333  Title   : get_sequence_counter
 334  Usage   : $obj->get_sequence_counter()
 335  Function: Get the count of sequences created
 336  Returns : value of counter
 337  Args    : -
 338
 339 =cut
 340
 341 sub get_sequence_counter{
 342    return shift->{'_sequence_counter'};
 343 }
 344
 345 =head2 reset_sequence_counter
 346
 347  Title   : reset_sequence_counter
 348  Usage   : $obj->reset_sequence_counter()
 349  Function: Resert the counter of sequences created
 350  Returns : value of counter
 351  Args    : -
 352
 353 This is called when ever mutated sequences are reassigned new values
 354 using methods seq() and mutated_seq().  As a side affect, this method
 355 also recreates the intermal alignment that is used to calculate the
 356 sequence identity.
 357
 358 =cut
 359
 360 sub reset_sequence_counter{
 361    my $self = shift;
 362    $self->{'_sequence_counter'} = 0;
 363    $self->_init_alignment;
 364    return 1;
 365 }
 366
 367
 368
 369 =head2 each_seq
 370
 371  Title   : each_seq
 372  Usage   : $obj->each_seq($int)
 373  Function:
 374  Returns : an array of sequences mutated from the reference sequence
 375            according to evolutionary parameters given
 376  Args    : -
 377
 378 =cut
 379
 380 sub each_seq{
 381    my $self = shift;
 382    my $number = shift;
 383
 384    $self->throw("[$number] ". ' should be a positive integer')
 385        unless $number =~ /^[+\d]+$/;
 386
 387    my @array;
 388    for (my $count=1; $count<$number; $count++) {
 389        push @array, $self->next_seq();
 390
 391    }
 392    return @array;
 393 }
 394
 395
 396
 397 =head2 each_mutation
 398
 399   Title   : each_mutation
 400   Usage   : $obj->each_mutation
 401   Function: return the mutations leading to the last generated
 402             sequence in objects
 403   Returns : an array of Bio::Variation::DNAMutation objects
 404   Args    : optional argument to return an array of  stringified names
 405
 406 =cut
 407
 408 sub each_mutation {
 409     my $self = shift;
 410     my $string = shift;
 411
 412     return @{$self->{'_mutations'}} if $string;
 413
 414     return map {
 415         /(\d+)(\w*)>(\w*)/;
 416 #        print;
 417         my $dnamut = Bio::Variation::DNAMutation->new
 418             ('-start'         => $1,
 419              '-end'           => $1,
 420              '-length'        => 1,
 421              '-isMutation'    => 1
 422             );
 423         $dnamut->allele_ori( Bio::Variation::Allele->new(-seq => $2,
 424                                                          -alphabet => 'dna') );
 425         $dnamut->add_Allele( Bio::Variation::Allele->new(-seq => $3,
 426                                                          -alphabet => 'dna') );
 427         $dnamut;
 428     } @{$self->{'_mutations'}}
 429 }
 430
 431
 432 sub get_alignment_identity  {
 433     my $self = shift;
 434     return $self->{'_align'}->overall_percentage_identity;
 435 }
 436
 437
 438 sub get_alignmet {
 439    my $self = shift;
 440    return $self->{'_align'}->remove_gaps('-', 'all-gaps');
 441 }
 442
 443
 444 =head1 Internal methods
 445
 446 =cut
 447
 448
 449 =head2 _increase_mutation_counter
 450
 451  Title   : _increase_mutation_counter
 452  Usage   : $obj->_increase_mutation_counter()
 453  Function: Internal method to increase the counter of mutations performed
 454  Returns : value of counter
 455  Args    : -
 456
 457 =cut
 458
 459 sub _increase_mutation_counter{
 460    return shift->{'_mutation_counter'}++;
 461 }
 462
 463
 464
 465 =head2 _increase_sequence_counter
 466
 467  Title   : _increase_sequence_counter
 468  Usage   : $obj->_increase_sequence_counter()
 469  Function: Internal method to increase the counter of sequences created
 470  Returns : value of counter
 471  Args    : -
 472
 473 =cut
 474
 475 sub _increase_sequence_counter{
 476    return shift->{'_sequence_counter'}++;
 477 }
 478
 479
 480 1;
 481