Генерация синтетической последовательности ДНК со скоростью замещения
-
03-07-2019 - |
Вопрос
Учитывая эти входные данные:
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++;
}
Решение
В качестве небольшой оптимизации замените:
$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,001 до 1,000:
А также $roll
, сгенерируйте другое (псевдо)случайное число в диапазоне [1..1000], если оно меньше или равно (1000 * $sub_rate), то выполните замену, в противном случае ничего не делайте (т.е.выход «А»).
Имейте в виду, что вы можете внести незначительную предвзятость, если не известны свойства вашего генератора случайных чисел.
Это не совсем то, что вы ищете, но я предлагаю вам взглянуть на BioPerl. Био::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
метод.
В случае замены вы хотите исключить текущую базу из возможностей:
my @other_bases = grep { $_ ne substr($init_seq, $i, 1) } @dna;
$base = @other_bases[int(rand 3)];
Также, пожалуйста, посмотрите Ответ Митча Уита о том, как реализовать коэффициент замещения.
Не знаю, правильно ли я понял, но я бы сделал что-то вроде этого (псевдокод):
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