Hi guys.

I am trying to learn and understand the ideas of markov models. I tried to develop a code to accomplish my understanding of this.

The idea at the moment is to build or create a random sequence with the model using the code i tried to develop with python. The final idea is to use it to compare the k-mer from a genome against the random sequence.

If you have time, patience to check if the code is correct in it ideas I really appreciate or if you have any other site or I could get help I would appreciate.

My code is this:

import random


random.seed(31)


def random_weighted_base(dist_dict):
    '''This function use a specified frequency distribution and choose and returns 
    a random symbol according to the given distribution '''
    # generate a random number in [0,1)
    bases = list(dist_dict.keys())
    freqs = get_count_frequencies(dist_dict)
    p = list(freqs.values())
    return ''.join(random.choices(bases, weights=p))


def transition_matrix(base_freq, mer_freq):
    """Returns a transition matrix from a base_count dictionary and a
    k-mer frequncy dictionary as inputs."""
    tm = defaultdict(dict)
    for char1, cnt1 in mer_freq.items():
        tm[char1] = tm.get(char1, {})
        for char2, cnt2 in base_freq.items():
            tm[char1][char2] = tm[char1].get(char2, cnt1 * cnt2)
    return tm


def kmer_count(sequence, alphabet=set(['A', 'C', 'G', 'T']), k=1):
    """Return a dictionary with the k-mers as keys and it counts
    as values."""                          
    seq = sequence.upper()
    seq_len = len(seq)                                            
    kmers = [seq[i:i + k] for i in range(0, seq_len - k + 1)]     
    filterd_kmers = [kmer for kmer in kmers if                    
                     all(base in set(alphabet) for base in kmer)] 
    return Counter(filterd_kmers)


def get_count_frequencies(counts):
    """Receive a dictionary as key:counts."""
    total = sum(counts.values())
    return {key: (cnt/total) for key, cnt in counts.items()}


def markov_higher_order(symbol_weights, mer_count, transition_matrix, length):
    """Receive a symbol_weights(dictionary wiht the count of the symbols), 
    a function, that converts counts in frequencies, a transition matrix and 
    the length of the random sequence that is the final output.
    """
    seq = []

    mer = get_random_weighted_mer(mer_count)
    for _ in range(length - 1):
        base, _, _ = get_random_weighted_base(symbol_weights)
        nuc = transition_matrix[mer]
        for char, cnt in nuc.items():
            if char == base:
                seq.append(char)
    return ''.join(seq)

I make the modifications you suggest and I make something different making the sequence adding the mer+base, once the markov order is dependent in this case of base[i-2] and base[i-1] to chose the base[i] right?

def markov_second_order_(symbol_weights, mer_count, transition_matrix, length):
    """Receive a symbol_weights, that is a dictionary with symbol counts, a function
    that converts counts in frequencies, a transition matrix and the length of the random
    sequence that is the output final of the markov_first_order.
    """
    seq = []
    for _ in range(length - 1):
        mer = get_random_weighted_mer(mer_count)
        base, _, _ = get_random_weighted_base(symbol_weights)
        nuc = transition_matrix[mer]
        if base in nuc:
            seq.append(mer+base)
    return ''.join(seq)

If you have any nice reference about the subject, I would really appreciate.

Thank you for your time.

Paulo