Comment modifier aléatoirement la matrice de transition en utilisant une chaîne de Markov?

Question

J'utilise une matrice de transition de 1ère étape pour générer les séquences d'ADN. Maintenant, je dois donner une probabilité à la matrice de transition de changer tous les 1000 pas. Disons que tous les 1000 pas, il y a 40% de probabilité que la matrice de transition change. Chaque ligne doit être ajoutée à 1 après le changement. Maintenant, je ne sais pas comment accéder à la valeur dans les données de dictionnaire imbriquées dans python, et comment implémenter le changement de probabilité de 40%. Je joins mon code ici, toute suggestion est appréciée ~

#!/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]) 

Answer 1

EDIT: Ajout d'une mise en œuvre rapide de code qui génère de nouvelles fréquences de transition. Vous devrez peut-être jouer pour trouver quel générateur de nombres aléatoires fonctionne le mieux pour votre cas, et pour voir si vous pouvez obtenir une génération plus sensible en utilisant certaines contraintes sur les probabilités aléatoires.

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])