Erzeugung synthetischen DNA-Sequenz mit Subtitution Rate

Question

Bei diesen Eingängen:

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

Ich mag generieren:

1. Eintausend Länge-10-Tags für jede Position

2. Substitutionsrate in ein Tag ist 0.003

Nachgeben Ausgabe wie:

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

Gibt es eine kompakte Art und Weise es in Perl zu tun?

Ich bin mit der Logik dieses Skripts als Kern stuck:

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

Als kleine Optimierung ersetzen:

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

mit

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

Answer 2

EDIT: Rate Substitution Unter der Annahme, liegt im Bereich von 0,001 bis 1,000:

Neben $roll, erzeugt ein anderes (pseudo) Zufallszahl im Bereich [1..1000], wenn sie kleiner oder gleich (1000 * $ sub_rate) ist, dann führe die Substitution durch, ansonsten tue nichts (dh gib 'A' aus).

Beachten Sie, dass Sie geringfügige Verzerrungen einführen können, wenn die Eigenschaften Ihres Zufallsgenerators nicht bekannt sind.

Answer 3

nicht genau das, was Sie suchen, aber ich schlage vor, Sie Bio::SeqEvolution::DNAPoint Modul einen Blick auf BioPerl des nehmen. Es erfordert jedoch keine Mutationsrate als Parameter. Es fragt vielmehr nach der Untergrenze der Sequenzidentität mit dem Original, das Sie möchten.

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 } 

1000 Alle mutierten Sequenzen in der @mutated Array gespeichert werden sollen, können ihre Sequenzen über die seq Verfahren zugegriffen werden.

Answer 4

Im Falle einer Substitution, möchten Sie die aktuelle Basis von den Möglichkeiten auszuschließen:

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

Bitte beachten Sie auch Mitch Wheat's answer dafür, wie die Substitutionsrate zu implementieren.

Answer 5

Ich weiß nicht, ob ich so (Pseudo-Code) korrekt, aber ich würde etwas tun, verstehen:

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