Generazione della sequenza di DNA sintetico con Subtitution Vota

Question

Attribuite questi ingressi:

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

voglio generare:

1. mille lunghezza-10 tag

2. tasso di sostituzione per ogni posizione in un tag è 0.003

uscita rendimento come:

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

C'è un modo compatto per farlo in Perl?

sono attaccato con la logica di questo script come nucleo:

#!/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

Come una piccola ottimizzazione, sostituire:

$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: Assumendo tasso di sostituzione è nel 0,001 a 1.000:

Oltre $roll, generare un'altra (pseudo) numero casuale nell'intervallo [1..1000], se è minore o uguale a (1000 * $ sub_rate) quindi eseguire la sostituzione, altrimenti non fare nulla (ad esempio, uscita 'A').

Tenere presente che è possibile introdurre un bias sottile a meno che non siano note le proprietà del generatore di numeri casuali.

Answer 3

Non è esattamente quello che stai cercando, ma vi consiglio di dare un'occhiata al modulo di Bioperl Bio::SeqEvolution::DNAPoint. Tuttavia, non richiede un tasso di mutazione come parametro. Piuttosto, chiede quale sia il limite inferiore dell'identità della sequenza con l'originale che si desidera.

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 } 

Tutti 1000 mutato sequenze saranno memorizzati nella matrice @mutated, le sequenze possono essere accessibili tramite il metodo seq.

Answer 4

In caso di una sostituzione, si desidera escludere la corrente di base dalle possibilità:

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

Inoltre si prega di vedere Mitch Wheat's answer per come implementare il tasso di sostituzione.

Answer 5

Non so se ho capito bene, ma mi piacerebbe fare qualcosa di simile a questo (pseudocodice):

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