置換率を使用した合成DNAシーケンスの生成
-
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 );
生成したい:
-
1000個の長さ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;};
with
$base = $dna[int(rand 4)];
他のヒント
編集:置換率が0.001から1.000の範囲であると仮定します:
$ roll
と同様に、[1..1000]の範囲で別の(擬似)乱数を生成し、それが(1000 * $ sub_rate)以下の場合は、置換、それ以外は何もしません(つまり、 'A'を出力します)。
乱数ジェネレーターのプロパティがわからない限り、微妙なバイアスを導入する可能性があることに注意してください。
探しているものとは異なりますが、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
メソッドを介してアクセスできます。
置換の場合、現在のベースを除外する可能性から除外します:
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
所属していません StackOverflow