Generowanie sekwencji DNA syntetyczny z Subtitution Rate

Question

Biorąc pod uwagę te wejść:

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

Chcę, aby wygenerować:

1. tysiąc długości 10 Tagi

2. stopy Zmiana dla każdej pozycji w tag to 0,003

Plonowanie wyjście jak:

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

Czy istnieje zwarty sposób to zrobić w Perl?

Ja zablokowany z logiką skryptu jako rdzeń:

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

Jako mały optymalizacji, należy wymienić:

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

z

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

Answer 2

Edycja: Zakładając, że szybkość podstawienia wynosi w zakresie od 0,001 do 1,000:

jak również $roll generują inny (pseudo) liczba losowa z przedziału [1..1000], jeśli jest mniejsza lub równa (1000 * $ sub_rate), należy wykonać podstawienie, w przeciwnym razie nic nie robić (tzn. wynik "A").

Pamiętaj, że możesz wprowadzić subtelne odchylenie, chyba że znane są właściwości generatora liczb losowych.

Answer 3

Nie dokładnie to, czego szukasz, ale proponuję spojrzeć na module BioPerl za Bio::SeqEvolution::DNAPoint. Nie bierze jednak współczynnika mutacji jako parametru. Zamiast tego prosi o to, jaką dolną granicę identyczności sekwencji z oryginałem chcesz.

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 } 

Wszystkie 1000 zmutowane sekwencje będą przechowywane w @mutated tablicy ich sekwencje można uzyskać metodą seq.

Answer 4

W przypadku zmiany, chcesz wykluczyć aktualną bazę z możliwości:

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

Także proszę zobaczyć Mitch Wheat's answer dla sposobu realizacji stopy substytucji.

Answer 5

ja nie wiem, czy dobrze rozumiem, ale chciałbym zrobić coś takiego (Pseudokod):

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