Génération ADN synthétique Séquence avec Subtitution taux

Question

Compte tenu de ces entrées:

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

Je veux générer:

1. Mille longueur-10 balises

2. Taux de substitution pour chaque position une étiquette est 0,003

sortie Cédant comme:

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

est-il un moyen compact de le faire en Perl?

Je suis coincé avec la logique de ce script comme noyau:

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

En tant que petite optimisation, Replace:

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

avec

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

Answer 2

EDIT: En supposant que le taux de substitution est comprise entre 0,001 à 1,000:

En plus $roll, générer une autre (pseudo) nombre aléatoire dans la plage [1..1000], s'il est inférieur ou égal à (1000 * $ sub_rate), alors effectuez la substitution, sinon ne faites rien (ex: sortie 'A'). Sachez que vous pouvez introduire un biais subtil à moins que les propriétés de votre générateur de nombres aléatoires soient connues.

Answer 3

Pas exactement ce que vous cherchez, mais je vous suggère de jeter un coup d'œil sur le module Bio::SeqEvolution::DNAPoint de BioPerl. Cependant, il ne prend pas le taux de mutation comme paramètre. Au contraire, il demande quelle est la limite inférieure de l'identité de séquence avec l'original que vous voulez.

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 } 

Toutes les 1000 séquences mutées sont stockées dans la matrice @mutated, leurs séquences sont accessibles par la méthode seq.

Answer 4

En cas de substitution, vous voulez exclure la base actuelle des possibilités:

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

Aussi s'il vous plaît voir Mitch Wheat's answer pour savoir comment mettre en œuvre le taux de substitution.

Answer 5

Je ne sais pas si je comprends bien, mais je ferais quelque chose comme ça (pseudo-code):

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