生成具有取代率的合成 DNA 序列

Question

给定这些输入：

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

我想生成：

1. 一千个长度为 10 的标签

2. 标签中每个位置的替换率为 0.003

产生如下输出：

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

是否有一种紧凑的方法可以在 Perl 中做到这一点？

我坚持这个脚本的核心逻辑：

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

作为一个小的优化，将：替换

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

为

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

Answer 2

编辑：假设替换率在 0.001 到 1.000 范围内：

以及 $roll ，生成 [1..1000] 范围内的另一个（伪）随机数，如果它小于或等于到 (1000 * $sub_rate) 然后执行替换，否则不执行任何操作（即输出“A”）。

请注意，除非随机数生成器的属性已知，否则您可能会引入微妙的偏差。

Answer 3

不完全是你要找的，但我建议你看看 BioPerl 的 Bio::SeqEvolution::DNAPoint 模块。 但它不将突变率作为参数。 相反，它询问与您想要的原始序列同一性的下限是多少。

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 个突变序列将存储在 @mutated 数组中，它们的序列可以通过 seq 方法访问。

Answer 4

如果发生替换，您希望从可能性中排除当前碱基：

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

另请参阅Mitch Wheat 对于如何实现替代率的回答。

Answer 5

我不知道我是否理解正确，但我会做这样的事情（伪代码）：

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