치환 속도로 합성 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 );
생성하고 싶습니다 :
천 길이 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
, (1000 * $ sub_rate)보다 작거나 같은 경우 범위에서 다른 (Pseudo) 랜덤 숫자를 생성 한 다음 대체를 수행하고 그렇지 않으면 아무것도하지 않습니다 (예 : 출력 'A').
임의의 숫자 생성기의 속성이 알려지지 않으면 미묘한 편향을 도입 할 수 있습니다.
정확히 당신이 찾고있는 것이 아니지만 Bioperl 's를 살펴 보는 것이 좋습니다. 바이오 :: 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)];
또한 참조하십시오 미치 밀의 대답 대체율을 구현하는 방법.
내가 올바르게 이해하는지 모르겠지만 (pseudocode)와 같은 일을 할 것입니다.
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