Hi there!!
I really need your knowledge and experience in Python. I am trying to write a code that allows me to get the Markov Model of a sequence. So, my input is the sequence and the alphabet this sequence has and the output I would like to get is a matrix where collect the probabilities of each element.
My code is composed by two for loops: the first counts the elements and collects the result in a matrix (mat_rec) and the second loop calculate the probabilities and put the result in the transition matrix (matTra). Everything seems to be okay, but that is my wish. The problem (I think) is located in the second loop. I want to replace a row by the result but I don't get that. For this code I am using pandas (where I am a novice, as well as I am in Python) and probably I don't know how to use it correctly. Here it's my code and what I get:
from Bio.Alphabet import IUPAC
alf = list(IUPAC.protein.letters)
def estimaModeloMarkov(sec, alf):
import pandas as pd
mat_rec = pd.DataFrame(0,columns = alf, index=alf)
matTra = pd.DataFrame(0, columns = alf, index = alf)
for i in range(len(sec)-1):
mat_rec.loc[sec[i],sec[i+1]] = mat_rec.loc[sec[i],sec[i+1]] +1
for i in range(len(alf)):
letra = alf[i]
row = mat_rec.loc[letra,]
if sum(mat_rec.loc[letra,]) == 0:
matTra.loc[letra, alf[2]] = 1
frec = row.divide(sum(row))
fila_T = matTra.loc[letra,]
matTra = matTra.replace(fila_T, frec)
return(mat_rec, matTra)
Result
estimaModeloMarkov(sec, alf)
Out[2]:
( A C D E F G H I K L M N P Q R S T V W Y
A 4 1 0 0 0 2 0 0 2 1 0 0 0 0 0 0 1 1 0 0
C 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
D 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
E 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0
F 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0
H 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
I 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
K 1 0 2 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0
L 3 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0
M 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
N 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
P 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Q 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
R 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
S 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
T 2 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0
V 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0
W 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,
A C D E F G H I K ... N P Q R S T V W Y
A NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
C NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
D NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
E NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
F NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
G NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
H NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
I NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
K NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
L NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
M NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
N NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
P NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
Q NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
R NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
S NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
T NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
V NaN NaN 0.0 NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
W NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
Y NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN
[20 rows x 20 columns])
What I am doing wrong?? Thank you so much!!!
general recommendation - don't operate with probabilities, but with their logarithms. instead of dividing by e.g. 10^-100 you will just subtract log2(10^-100) - which is ~ -332 and you wont suffer from floating point arithmetic so much.