0

DNA配列を生成するために第1ステップの遷移行列を使用しています。ここで、1000ステップごとに変化する確率を遷移行列に与える必要があります。たとえば、1000ステップごとに、遷移行列が変化する確率は40％です。変更後、すべての行が1に追加されます。Pythonでネストされた辞書データの値にアクセスする方法と、40％の確率変更を実装する方法がわかりません。ここにコードを添付しました。提案をいただければ幸いです。

#!/usr/bin/env python

import sys, random


length = 10000

tran_matrix = {'a': {'a':0.495,'c':0.113,'g':0.129,'t':0.263},
               'c': {'a':0.129,'c':0.063,'g':0.413,'t':0.395},
               't': {'a':0.213,'c':0.495,'g':0.263,'t':0.029},
               'g': {'a':0.263,'c':0.129,'g':0.295,'t':0.313}}

initial_p = {'a':0.25,'c':0.25,'t':0.25,'g':0.25}             

def choose(dist):
    r = random.random()
    sum = 0.0
    keys = dist.keys()
    for k in keys:
        sum += dist[k]
        if sum > r:
        return k
    return keys[-1]
c = choose(initial_p)
for i in range(length):
    sys.stdout.write(c) 
    c = choose(tran_matrix[c])
4

1 に答える 1

0

編集: 新しい遷移周波数を生成するいくつかのコードの簡単な実装を追加しました。ケースに最適な乱数ジェネレーターを見つけ、ランダム確率にいくつかの制約を使用してより賢明な生成ができるかどうかを確認するために、いろいろ試してみる必要があるかもしれません。

import sys, random


LENGTH = 10000
CHANGE_EVERY = 1000
CHANGE_PROB = 0.4

tran_matrix = {'a': {'a':0.495,'c':0.113,'g':0.129,'t':0.263},
               'c': {'a':0.129,'c':0.063,'g':0.413,'t':0.395},
               't': {'a':0.213,'c':0.495,'g':0.263,'t':0.029},
               'g': {'a':0.263,'c':0.129,'g':0.295,'t':0.313}}

initial_p = {'a':0.25,'c':0.25,'t':0.25,'g':0.25}             


def choose(dist):
    r = random.random()
    sum = 0.0
    keys = dist.keys()
    for k in keys:
        sum += dist[k]
        if sum > r:
            return k
    return keys[-1]


def new_probs(precision=2):
    """
    Generate a dictionary of random transition frequencies, of the form
    {'a':0.495,'c':0.113,'g':0.129,'t':0.263}
    """
    probs = []
    total_prob = 0
    # Choose a random probability p1 from a uniform distribution in
    # the range (0, 1), then choose p2 in the range (0, 1 - p1), etc.
    for i in range(3):
        up_to = 1 - total_prob
        p = round(random.uniform(0, up_to), precision)
        probs.append(p)
        total_prob += p
    # Final probability is 1 - (sum of first 3 probabilities)
    probs.append(1 - total_prob)
    # Assign randomly to bases
    # If you don't shuffle the order of the bases each time, 't'
    # would end up with consistently lower probabilities
    bases = ['a', 'c', 'g', 't']
    random.shuffle(bases)
    new_prob_dict = {}
    for base, prob in zip(bases, probs):
        new_prob_dict[base] = prob
    return new_prob_dict

c = choose(initial_p)
for i in range(LENGTH):
    if i % CHANGE_EVERY == 0:
        dice_roll = random.random()
        if dice_roll < CHANGE_PROB:
            for base in tran_matrix:
                # Generate a new probability dictionary for each
                # base in the transition matrix
                tran_matrix[base] = new_probs()
    sys.stdout.write(c) 
    c = choose(tran_matrix[c])
于 2012-12-10T03:36:26.777 に答える