문제

이러한 입력이 주어지면 :

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++;
    }
도움이 되었습니까?

해결책

작은 최적화로 : 대체 :

    $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
라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top