Tạo trình tự ADN tổng hợp với Subtitution Rate

Question

Với những nguyên liệu đầu vào:

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

Tôi muốn tạo:

1. Một nghìn chiều dài-10 thẻ

2. tỷ lệ Thay cho mỗi vị trí trong thẻ là 0,003

Sản lượng đầu ra như:

AAAAAAAAAA 
AATAACAAAA 
..... 
AAGGAAAAGA # 1000th tags 

Có cách nào nhỏ gọn để làm điều đó trong Perl?

Tôi bị mắc kẹt với logic của kịch bản này là cốt lõi:

#!/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++; 
    } 

Answer 1

Là một tối ưu hóa nhỏ, thay thế:

$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;}; 

với

$base = $dna[int(rand 4)]; 

Answer 2

EDIT: Giả sử tỷ lệ thay thế là trong khoảng 0,001-1,000:

Cũng như $roll, tạo khác (giả) số ngẫu nhiên trong phạm vi [1..1000], nếu nó nhỏ hơn hoặc bằng (1000 * $ sub_rate) thì thực hiện thay thế, nếu không thì không có gì (tức là đầu ra 'A').

Lưu ý rằng bạn có thể giới thiệu độ lệch tinh tế trừ khi các thuộc tính của trình tạo số ngẫu nhiên của bạn được biết.

Answer 3

Không chính xác những gì bạn đang tìm kiếm, nhưng tôi đề nghị bạn hãy xem môđun Bio::SeqEvolution::DNAPoint BioPerl của. Nó không có tỷ lệ đột biến như một tham số mặc dù. Thay vào đó, nó sẽ hỏi những gì thấp hơn ràng buộc của bản sắc trình tự với bản gốc bạn muốn.

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 } 

Tất cả 1000 chuỗi đột biến sẽ được lưu trữ trong mảng @mutated, trình tự của chúng có thể được truy cập thông qua phương pháp seq.

Answer 4

Trong trường hợp thay người, bạn muốn loại trừ các cơ sở hiện từ các khả năng:

my @other_bases = grep { $_ ne substr($init_seq, $i, 1) } @dna; 
$base = @other_bases[int(rand 3)]; 

Cũng xin xem Mitch Wheat's answer cho làm thế nào để thực hiện các tỷ lệ thay thế.

Answer 5

Tôi không biết nếu tôi hiểu đúng nhưng tôi muốn làm điều gì đó như thế này (giả):

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