Gerando sequência de ADN sintético com Substituição Classificação

Question

Tendo em conta estas entradas:

my $init_seq = "AAAAAAAAAA" #length 10 bp 
my $sub_rate = 0.003;
my $nof_tags = 1000;
my @dna = qw( A C G T );

Eu quero gerar:

1. Mil Comprimento-10 etiquetas

2. Taxa de Substituição para cada posição em que um marcador é 0,003

Cedendo uma saída como:

AAAAAAAAAA
AATAACAAAA
.....
AAGGAAAAGA # 1000th tags

Existe uma maneira compacta para fazê-lo em Perl?

Eu estou preso com a lógica deste script como core:

#!/usr/bin/perl

my $init_seq = "AAAAAAAAAA" #length 10 bp 
my $sub_rate = 0.003;
my $nof_tags = 1000;
my @dna = qw( A C G T );

    $i = 0;
    while ($i < length($init_seq)) {
        $roll = int(rand 4) + 1;       # $roll is now an integer between 1 and 4

        if ($roll == 1) {$base = A;}
        elsif ($roll == 2) {$base = T;}
        elsif ($roll == 3) {$base = C;}
        elsif ($roll == 4) {$base = G;};

        print $base;
    }
    continue {
        $i++;
    }

Answer 1

Como uma pequena otimização, substitua:

    $roll = int(rand 4) + 1;       # $roll is now an integer between 1 and 4

    if ($roll == 1) {$base = A;}
    elsif ($roll == 2) {$base = T;}
    elsif ($roll == 3) {$base = C;}
    elsif ($roll == 4) {$base = G;};

com

    $base = $dna[int(rand 4)];

Answer 2

EDIT: Assumindo que a taxa de substituição é na gama 0,001-1,000:

Como bem como $roll, gerar outro (pseudo) de número aleatório no intervalo [1..1000], se é menor do que ou igual a (1000 * $ sub_rate), em seguida, realizar a substituição, caso contrário, não fazer nada (ou seja, de saída 'A').

Esteja ciente de que você pode introduzir viés sutil menos que sejam conhecidas as propriedades do seu gerador de números aleatórios.

Answer 3

Não exatamente que você está procurando, mas eu sugiro que você dê uma olhada BioPerl de Bio :: SeqEvolution :: DNAPoint módulo. Não é preciso taxa de mutação como parâmetro embora. Pelo contrário, ela pergunta o que o limite inferior de identidade de sequência com o original que deseja.

use strict;
use warnings;
use Bio::Seq;
use Bio::SeqEvolution::Factory;

my $seq = Bio::Seq->new(-seq => 'AAAAAAAAAA', -alphabet => 'dna');

my $evolve = Bio::SeqEvolution::Factory->new (
   -rate     => 2,      # transition/transversion rate
   -seq      => $seq
   -identity => 50      # At least 50% identity with the original
);


my @mutated;
for (1..1000) { push @mutated, $evolve->next_seq }

Todos os 1000 sequências mutadas serão armazenados na matriz @mutated, as suas sequências podem ser acedidas através do método seq.

Answer 4

No caso de uma substituição, você quer excluir a base atual a partir das possibilidades:

my @other_bases = grep { $_ ne substr($init_seq, $i, 1) } @dna;
$base = @other_bases[int(rand 3)];

Além disso consulte o Mitch Wheat resposta para como implementar a taxa de substituição.

Answer 5

Eu não sei se eu entendi corretamente, mas eu faria algo assim (pseudocódigo):

digits = 'ATCG'
base = 'AAAAAAAAAA'
MAX = 1000
for i = 1 to len(base)
  # check if we have to mutate
  mutate = 1+rand(MAX) <= rate*MAX
  if mutate then
    # find current A:0 T:1 C:2 G:3
    current = digits.find(base[i])
    # get a new position 
    # but ensure that it is not current
    new = (j+1+rand(3)) mod 4        
    base[i] = digits[new]
  end if
end for