Génération d'une séquence d'ADN synthétique avec taux de substitution

Question

Étant donné 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 balises de longueur 10 -

2. Le taux de substitution pour chaque position dans une balise est de 0,003

Production de résultats tels que:

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

Existe-t-il un moyen compact de le faire en Perl?

Je suis coincé avec la logique de ce script en tant que 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

Comme petite optimisation, remplacez:

    $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 soit compris entre 0,001 et 1 000:

En plus de $ roll , générez un autre nombre (pseudo) aléatoire compris entre [1..1000], s'il est inférieur ou égal à (1000 * $ sub_rate), puis exécutez la procédure suivante: substitution, sinon ne faites rien (c.-à-d. la sortie 'A').

Sachez que vous pouvez introduire un biais subtil si les propriétés de votre générateur de nombres aléatoires ne sont pas connues.

Answer 3

Ce n'est pas exactement ce que vous recherchez, mais je vous suggère de consulter le module Bio :: SeqEvolution :: DNAPoint . Cependant, le taux de mutation n’est pas pris comme paramètre. Il demande plutôt quelle est la limite inférieure de l’identité de séquence avec l’original de votre choix.

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 seront stockées dans le tableau @mutated. Vous pouvez accéder à leurs séquences via la méthode seq .

Answer 4

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

my @other_bases = grep { 
 En cas de substitution, vous souhaitez  exclure la base actuelle  des possibilités suivantes: 

<*>

 Voir également la  réponse de Mitch Wheat  pour comment mettre en œuvre le taux de substitution. 
 ne substr($init_seq, $i, 1) } @dna;
$base = @other_bases[int(rand 3)];

Voir également la réponse de Mitch Wheat pour comment mettre en œuvre le taux de substitution.

Answer 5

Je ne sais pas si j'ai bien compris, mais je ferais quelque chose comme ça (pseudocode):

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