pomegranate icon indicating copy to clipboard operation
pomegranate copied to clipboard

Empty viterbi_path

Open malonzm1 opened this issue 2 years ago • 0 comments

Hi!

viterbi_likelihood, viterbi_path = model.viterbi(sequence) returns -inf for viterbi_likelihood and an empty viterbi_path. Below is the code. For a different dataset it worked. total_reads_Y.txt.tar.gz methylated_reads_Y.txt.tar.gz

import numpy as np
from pomegranate import *
import os

def hidden_states(bsTot_list,bsMeth_list):
    ratio = np.divide(bsMeth_list,bsTot_list)
    control = np.nanmean(ratio[:,:2],axis=1)
    case = np.nanmean(ratio[:,2:4],axis=1)
    diff = case - control
    return diff

chroms = ['Y']
sequences = []
for chrom in chroms:
    bsTot_file = 'total_reads_%s.txt'%(chrom)
    bsTot_list = np.loadtxt(bsTot_file,delimiter='\t',skiprows=0,dtype='float', usecols=range(1, 5))
    bsTot_list[bsTot_list < 5] = np.nan
    bsMeth_file = 'methylated_reads_%s.txt'%(chrom)
    bsMeth_list = np.loadtxt(bsMeth_file,delimiter='\t',skiprows=0,dtype='float', usecols=range(1, 5))
    sequence=hidden_states(bsTot_list,bsMeth_list)
    sequences.append(sequence)
s0 = State(NormalDistribution(0, 0.08), name = 's0')
s1 = State(NormalDistribution(0.3, 0.06), name = 's1')
s2 = State(NormalDistribution(-0.3, 0.06), name = 's2')
model = HiddenMarkovModel()
model.add_states(s0, s1, s2)
model.add_transition(model.start, s0, 0.333)
model.add_transition(model.start, s1, 0.333)
model.add_transition(model.start, s2, 0.333)
model.add_transition(s0, s0, 0.5)
model.add_transition(s0, s1, 0.25)
model.add_transition(s0, s2, 0.25)
model.add_transition(s1, s0, 0.25)
model.add_transition(s1, s1, 0.5)
model.add_transition(s1, s2, 0.25)
model.add_transition(s2, s0, 0.25)
model.add_transition(s2, s1, 0.25)
model.add_transition(s2, s2, 0.5)
model.add_transition(s0, model.end, 0.333)
model.add_transition(s1, model.end, 0.333)
model.add_transition(s2, model.end, 0.333)
model.bake()
model.fit(sequences, distribution_inertia=1.0)
viterbi_likelihood, viterbi_path = model.viterbi(sequences[0])

Pls. advise.

malonzm1 avatar Jun 09 '22 06:06 malonzm1