生成具有替换率的合成 DNA 序列答案

【问题标题】：Generating Synthetic DNA Sequence with Substitution Rate生成具有替换率的合成 DNA 序列
【发布时间】：2009-03-02 09:19:34
【问题描述】：

鉴于这些输入：

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

我要生成：

千长10个标签
标签中每个位置的替换率为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++;
    }

【问题讨论】：

这是作业，对吧？：birg.cs.wright.edu/resources/perl/hw3.shtml
不，米奇，这不是家庭作业。真的。
您应该检查是否有重复项。
其他人已经编写了程序来生成发生突变的序列。您可能想看看这些现有程序是否满足您的需求，或者您是否需要重新实现轮子。

标签： perl algorithm bioinformatics dna-sequence

【解决方案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)];

【讨论】：

+0。这是一个很好的优化，但它允许从 G 到 G 的“突变”。
G->G“自我突变”实际上是计算生物学中的替换矩阵考虑的真正突变。有两种理由，一种是生化的，一种是统计的。在生化方面，碱基被化学修饰但被 DNA 修复酶修复的可能性是有限的。从统计学上讲，大多数变异矩阵都描述了马尔可夫过程，因此必须考虑自转移或保持相同状态的概率。

【解决方案2】：

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

与$roll一样，在[1..1000]范围内生成另一个（伪）随机数，如果小于或等于（1000 * $sub_rate）则执行替换，否则不执行任何操作（即输出'A'）。

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

【讨论】：

rand() 返回一个 [0,1) 范围内的数字，因此可以直接与 $sub_rate 比较，无需任何 1000 *。

【解决方案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 方法访问。

【讨论】：

【解决方案4】：

如果发生替换，您希望排除当前基础的可能性：

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

另请参阅Mitch Wheat's 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

【讨论】：