6

これらの入力が与えられた場合：

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

生成したい：

  1. 1000の長さ-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++;
    }
4

5 に答える 5

5

小さな最適化として、以下を置き換えます。

    $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)];
于 2009-03-02T10:21:58.370 に答える
2

あなたが探しているものとは正確には異なりますが、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メソッドを介してアクセスできます。

于 2009-03-02T20:14:57.137 に答える