Generación de secuencia de ADN sintético con Subtitution Anote

Question

Dadas 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); 

I desea generar:

1. Uno mil longitud 10-tags

2. tasa de sustitución para cada posición en una etiqueta es 0.003

salida de rendimiento como:

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

¿Hay una manera compacta de hacerlo en Perl?

estoy atascado con la lógica de este script como núcleo:

#!/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 una pequeña optimización, reemplace:

$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;}; 

con

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

Answer 2

EDIT: Suponiendo tasa de sustitución está en el intervalo de 0,001 a 1,000:

Así como $roll, generar otra (pseudo) número aleatorio en el rango [1..1000], si es menor o igual que (1000 * $ sub_rate) luego realice la sustitución; de lo contrario, no haga nada (es decir, haga salir 'A').

Tenga en cuenta que puede introducir sesgos sutiles a menos que se conozcan las propiedades de su generador de números aleatorios.

Answer 3

No es exactamente lo que está buscando, pero yo sugiero que eche un vistazo a módulo de BioPerl Bio::SeqEvolution::DNAPoint. Sin embargo, no toma la tasa de mutación como un parámetro. Más bien, pregunta cuál es el límite inferior de la identidad de secuencia con el original que desea.

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 1000 mutado secuencias se almacenan en la matriz @mutated, sus secuencias se puede acceder a través del método seq.

Answer 4

En el caso de una sustitución, que desea excluye la base actual de las posibilidades:

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

También por favor ver Mitch Wheat's answer de la forma de aplicar la tasa de sustitución.

Answer 5

No sé si he entendido bien, pero me gustaría hacer algo como esto (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