view release on metacpan or  search on metacpan

lib/Bio/GeneDesign/Random.pm  view on Meta::CPAN

#
# GeneDesign libraries for Random DNA generation
#

=head1 NAME

GeneDesign::Random

=head1 VERSION

Version 5.56

=head1 DESCRIPTION

Random DNA Generators

=head1 AUTHOR

Sarah Richardson <SMRichardson@lbl.gov>

=cut

package Bio::GeneDesign::Random;
require Exporter;

use Bio::GeneDesign::Codons qw(_find_in_frame);
use Bio::GeneDesign::Basic qw(@BASES %NTIDES);
use List::Util qw(shuffle);

use strict;
use warnings;

our $VERSION = 5.56;

use base qw(Exporter);
our @EXPORT_OK = qw(
  _randomDNA
  _randombase
  _randombase_weighted
  _replace_ambiguous_bases
  _weighted_rand
  _random_index
);
our %EXPORT_TAGS =  (GD => \@EXPORT_OK);

=head1 Functions

=head2 _randomDNA()

  takes a target length and a GC percentage and generates a random nucleotide
  sequence, with or without stops in the first frame
  in: nucleotide sequence length (scalar),
      GC percentage (0 <= scalar <= 100),
      stop codon prevention(0 stops allowed, else no stops),
      codon table (hash reference)
  out: nucleotide sequence (string)

=cut

sub _randomDNA
{
  my ($len, $GCperc, $stopswit, $codon_t) = @_;

  return q{} if ($len == 0);
  return _randombase_weighted($GCperc) if ($len == 1);

  #GC
  my $GCtotal = sprintf "%.0f",  $GCperc * $len / 100;
  my $Gcount  = sprintf "%.0f", rand( $GCtotal );
  my $Gstr = 'G' x $Gcount;
  my $Ccount  = $GCtotal - $Gcount;
  my $Cstr = 'C' x $Ccount;

  #AT
  my $ATtotal = $len - $GCtotal;
  my $Acount  = sprintf "%.0f", rand( $ATtotal );
  my $Astr = 'A' x $Acount;
  my $Tcount  = $ATtotal - $Acount;
  my $Tstr = 'T' x $Tcount;

  my @randomarray = shuffle( split( '', $Gstr . $Cstr . $Astr . $Tstr) );
  my $DNA = join('', @randomarray);

  if ($stopswit)
  {
    my $stophsh = _find_in_frame($DNA, "*", $codon_t);
    while (scalar keys %{$stophsh})
    {
      foreach my $pos (keys %{$stophsh})
      {
        my $bit = substr $DNA, $pos, 3;
        substr $DNA, $pos, 3, scalar reverse $bit;
        if (int(rand(1)+.5) == 1)
        {
          my $bat = substr $DNA, $pos, 2;
          substr $DNA, $pos, 2, scalar reverse$bat;
        }
      }
      $stophsh = _find_in_frame($DNA, "*", $codon_t);
    }
  }
  return $DNA;
}

=head2 _randombase()

  when you just want one random base

=cut

sub _randombase
{
  my $int = _random_index(4);