Generazione di sequenze di DNA sintetico con tasso di sostituzione

Question

Dati questi input:

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 tag di lunghezza 10

2. Il tasso di sostituzione per ogni posizione in un tag è 0,003

Produrre output come:

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

Esiste un modo compatto per farlo in Perl?

Sono bloccato con la logica di questo script come core:

#!/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 piccola ottimizzazione, sostituisci:

    $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: supponendo che il tasso di sostituzione sia compreso nell'intervallo da 0,001 a 1.000:

Oltre a $ roll , genera un altro numero (pseudo) casuale nell'intervallo [1..1000], se è inferiore o uguale a (1000 * $ sub_rate), quindi esegui il sostituzione, altrimenti non fare nulla (ad es. output 'A').

Tieni presente che puoi introdurre una distorsione sottile a meno che non siano note le proprietà del tuo generatore di numeri casuali.

Answer 3

Non esattamente quello che stai cercando, ma ti suggerisco di dare un'occhiata a modulo Bio :: SeqEvolution :: DNAPoint . Tuttavia, non utilizza il tasso di mutazione come parametro. Piuttosto, chiede quale sia il limite inferiore dell'identità della sequenza con l'originale desiderato.

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 }

Tutte le 1000 sequenze mutate verranno archiviate nell'array @mutated, le cui sequenze sono accessibili tramite il metodo seq .

Answer 4

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

my @other_bases = grep { 
 In caso di sostituzione, si desidera  escludere la base corrente  dalle possibilità: 

<*>

 Vedi anche  Mitch Wheat's answer  per come applicare il tasso di sostituzione. 
 ne substr($init_seq, $i, 1) } @dna;
$base = @other_bases[int(rand 3)];

Vedi anche Mitch Wheat's answer per come applicare il tasso di sostituzione.

Answer 5

Non so se ho capito bene ma farei una cosa del genere (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