4.7 Práctica 2: predicción de estructuras secundarias en proteínas¶

En esta parte de la práctica ajustaremos un modelo de Markov oculto para predecir las estructuras secundarias con la misma base de datos utilizada en la sección anterior.

# Import required Libraries
import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt

Ejercicio 1. Conjunto de datos y definiciones.¶

Carga el archivo prots-L30.txt.

Recuerda que las columnas del “data frame” asociadas a secuencias y estructuras han de ser de tipo carácter.

Define el alfabeto de aminoácidos (RKDEQSCHNTWYMAILFVPG) y calcula:

la longitud del alfabeto;
el número de instancias;
la longitud de cada secuencia.

dataFolder = './data/'
dataFile = 'prots-L30.txt'
dataFilepath = os.path.join(dataFolder,dataFile)
dataFilepath

df = pd.read_csv(dataFilepath)
df.head()

alfabeto=[ele for ele in "RKDEQSCHNTWYMAILFVPG"]

print("Longitud del alfabeto =", len(alfabeto))
print("Número de secuencias =", len(df['seq']))
print("Longitud de las secuencias =", len(df['seq'][0]))

Longitud del alfabeto = 20
Número de secuencias = 174
Longitud de las secuencias = 30

Ejemplo 1. Lista de estructuras.¶

Las estructuras que aparecen en el conjunto de datos son de tres tipos:

C (“coil”),
E (hojas),
H (hélices).

Define una lista con las tres estructuras, en orden alfabético, y almacena en la variable nsta el número de estructuras que tenemos en este conjunto de datos.

strlist = [ ele for ele in "CEH"]
nsta = len(strlist)
strlist, nsta

(['C', 'E', 'H'], 3)

Ejercicio 2. Conjunto de entrenamiento y validación.¶

Define dos “data frames”, llamados dftrain y dftest que contengan el conjunto de entrenamiento y el de validación.

Reserva un 70 % de las instancias para entrenamiento, y calcula los tamaños ntrain y ntest de cada conjunto.

from sklearn.model_selection import train_test_split

dftrain, dftest = train_test_split(df, test_size=0.3,random_state=2)
ntrain = dftrain.shape[0]
ntest = dftest.shape[0]
print('Número instancias en train: {}\nNúmero instancias en test: {}'.format(ntrain,ntest))
dftrain.head()

Número instancias en train: 121
Número instancias en test: 53

	pdb_id	seq	sst3	len
167	5ZNF	KTYQCQYCEYRSADSSNLKTHIKTKHSKEK	CCEECCCCCCEECCHHHHHHHHHHHCCCCC	30
98	2OT8	ERPAQNEKRKEKNIKRGGNRFEPYANPTKR	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC	30
145	4NIB	FVNQHLCGSHLVEALYLVCAERAFFYTPKT	CCCCCCCCHHHHHHHHHHHHHHCECCCCCC	30
166	5XZX	GSSPEGGEDSDREDGNYCPPVKRERTSSLT	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC	30
112	3CM8	GPLGSMDVNPTLLFLKVPAQNAISTTFPYT	CCCCCCEECHHHHHHCCCCCCCCCCCCCCC	30

Ejemplo 2. Obtención de los aminoácidos asociados a la estructura C en la primera secuencia.¶

Considera la primera secuencia del conjunto de entrenamiento, y su estructura secundaria asociada.

Encuentra la lista de aminoácidos que tienen la estructura C y almacena el resultado en la variable aasubset.

Calcula la longitud de esta lista (laasubset).

primeraInstanciaEntrenamiento = dftrain.iloc[0,:]
primeraInstanciaEntrenamiento

pdb_id                              5ZNF
seq       KTYQCQYCEYRSADSSNLKTHIKTKHSKEK
sst3      CCEECCCCCCEECCHHHHHHHHHHHCCCCC
len                                   30
Name: 167, dtype: object

primeraSecuencia = primeraInstanciaEntrenamiento['seq']
primeraSecuencia

'KTYQCQYCEYRSADSSNLKTHIKTKHSKEK'

primeraEstSecundaria = primeraInstanciaEntrenamiento['sst3']
primeraEstSecundaria

'CCEECCCCCCEECCHHHHHHHHHHHCCCCC'

posicionesC = []
for posicion,letra in enumerate(primeraEstSecundaria):
    if letra =='C':
        posicionesC.append(posicion)
posicionesC

[0, 1, 4, 5, 6, 7, 8, 9, 12, 13, 25, 26, 27, 28, 29]

aasubset = []
for pos in posicionesC:
    AA = primeraSecuencia[pos]
    aasubset.append(AA)
laasubset = len(aasubset)
print(aasubset,laasubset)

['K', 'T', 'C', 'Q', 'Y', 'C', 'E', 'Y', 'A', 'D', 'H', 'S', 'K', 'E', 'K'] 15

# sin repeticiones
aasubset = list(set(aasubset))
laasubset = len(aasubset)
aasubset,laasubset

(['H', 'K', 'S', 'A', 'Y', 'D', 'Q', 'C', 'E', 'T'], 10)

Ejercicio 3. Frecuencias de aparición de los aminoácidos asociados a la estructura C en la primera secuencia.¶

Define una matriz frecs con una fila y un número de columnas igual al tamaño del alfabeto.

Etiqueta las columnas para mantener el mismo orden de aminoácidos que en la lista que define el alfabeto (RKDEQSCHNTWYMAILFVPG).

Almacena en este vector las frecuencias de cada aminoácido en la primera secuencia.

frecs = pd.DataFrame(0,index=['Frecuencias'],columns=alfabeto)
frecs

for posicion,letra in enumerate(primeraEstSecundaria):
    if letra =='C':
        # buscamos el AA
        AA = primeraSecuencia[posicion]
        # actualizamos la tabla
        frecs.loc['Frecuencias',AA]+=1

print('Suma de frecuencias: {}'.format(frecs.sum().sum()))
frecs

Suma de frecuencias: 15

	R	K	D	E	Q	S	C	H	N	T	W	Y	M	A	I	L	F	V	P	G
Frecuencias	0	3	1	2	1	1	2	1	0	1	0	2	0	1	0	0	0	0	0	0

Ejercicio 4. Frecuencias absolutas de los aminoácidos asociados a la estructura C.¶

Reutiliza el código del ejercicio anterior para calcular las frecuencias absolutas de los aminoácidos asociados a una estructura tipo C en el conjunto de entrenamiento.

numInstancias = dftrain.shape[0]

frecsAbsolutas = pd.DataFrame(0,index=['Frecuencias'],columns=alfabeto)

for ninst in range(numInstancias):
    instancia = dftrain.iloc[ninst,:]
    secuencia = instancia['seq']
    est2 = instancia['sst3']
    for posicion,letra in enumerate(est2):
        if letra =='C':
            # buscamos el AA
            AA = secuencia[posicion]
            # actualizamos la tabla
            frecsAbsolutas.loc['Frecuencias',AA]+=1
frecsAbsolutas

	R	K	D	E	Q	S	C	H	N	T	W	Y	M	A	I	L	F	V	P	G
Frecuencias	105	151	73	102	67	127	175	46	83	134	22	59	23	129	83	91	77	80	184	212

Ejercicio 5. Matriz de emisión.¶

Repite el ejercicio 4 para las otras dos estructuras.

Almacena las frecuencias en una matriz con 3 filas y 20 columnas, utilizando la ordenación de aminoácidos definida en el alfabeto y el orden alfabético en las estructuras.

todasFrecsAbsolutas = pd.DataFrame(0,index=strlist,columns=alfabeto)
todasFrecsAbsolutas

for stru in strlist:
    for ninst in range(numInstancias):
        instancia = dftrain.iloc[ninst,:]
        secuencia = instancia['seq']
        est2 = instancia['sst3']
        for posicion,letra in enumerate(est2):
            if letra ==stru:
                # buscamos el AA
                AA = secuencia[posicion]
                # actualizamos la tabla
                todasFrecsAbsolutas.loc[stru,AA]+=1
todasFrecsAbsolutas

	R	K	D	E	Q	S	C	H	N	T	W	Y	M	A	I	L	F	V	P	G
C	105	151	73	102	67	127	175	46	83	134	22	59	23	129	83	91	77	80	184	212
E	41	12	6	8	12	17	84	5	4	36	18	58	2	34	28	16	36	22	9	18
H	51	53	36	95	37	64	49	44	27	32	20	50	18	126	55	183	32	100	12	57

Comprueba que la suma de frecuencias coincide con los datos que estamos usando.

print('Suma frecuencias absolutas: {}'.format(todasFrecsAbsolutas.sum().sum()))
print('Suma elementos en datos: {}'.format(dftrain.shape[0]*len(dftrain['seq'].iloc[0])))

Suma frecuencias absolutas: 3630
Suma elementos en datos: 3630

Transforma frecuencias absolutas en relativas para obtener una matriz estocástica.

MatrizEmision = todasFrecsAbsolutas.copy()

for stru in strlist:
    sumRow = todasFrecsAbsolutas.loc[stru,:].sum()
    MatrizEmision.loc[stru,:]=MatrizEmision.loc[stru,:]/sumRow
MatrizEmision

	R	K	D	E	Q	S	C	H	N	T	W	Y	M	A	I	L	F	V	P	G
C	0.051903	0.074642	0.036085	0.050420	0.033119	0.062778	0.086505	0.022739	0.041028	0.066238	0.010875	0.029165	0.011369	0.063767	0.041028	0.044983	0.038062	0.039545	0.090954	0.104795
E	0.087983	0.025751	0.012876	0.017167	0.025751	0.036481	0.180258	0.010730	0.008584	0.077253	0.038627	0.124464	0.004292	0.072961	0.060086	0.034335	0.077253	0.047210	0.019313	0.038627
H	0.044698	0.046450	0.031551	0.083260	0.032428	0.056091	0.042945	0.038563	0.023663	0.028046	0.017528	0.043821	0.015776	0.110429	0.048203	0.160386	0.028046	0.087642	0.010517	0.049956

Suma las filas y las columnas de la matriz de emisión para comprobar que es estocástica.

MatrizEmision.sum(axis=0)

R    0.184584
K    0.146843
D    0.080512
E    0.150848
Q    0.091298
S    0.155350
C    0.309707
H    0.072031
N    0.073275
T    0.171537
W    0.067030
Y    0.197449
M    0.031437
A    0.247158
I    0.149317
L    0.239703
F    0.143361
V    0.174398
P    0.120784
G    0.193378
dtype: float64

MatrizEmision.sum(axis=1)

C    1.0
E    1.0
H    1.0
dtype: float64

Ejercicio 6. Matriz de transiciones entre estados y vector de probabilidades inicial.¶

Al igual que hicimos en los ejercicios 1, 2 y 3 de la primera parte de la práctica (cadenas de Markov sin variables ocultas), tenemos que calcular las probabilidades de transición entre estados y las probabilidades iniciales asociadas al primer estado las estructuras.

En este caso hay que obtener una matriz \(3 \times 3\) y un vector de 3 componentes.

Utiliza el conjunto de entrenamiento para estimar las probabilidades. Nombra las filas y columnas de la matriz (y las columnas del vector inicial) con las etiquetas correspondientes a los tres estados: \(C\), \(E\) y \(H\).

Una vez hecho este ejercicio habrás finalizado la fase de entrenamiento (aprendizaje) del modelo de Markov oculto.

Empieza calculando el vector de probabilidad inicial.

# vector de probabilidades inicial

vectorProbInicial = pd.DataFrame(0.0,columns=strlist,index=['FrecuenciasAbsolutas','Probabilidades'])
vectorProbInicial

for ninst in range(numInstancias):
    instancia = dftrain.iloc[ninst,:]
    #secuencia = instancia['seq'] # no lo observamos
    est2 = instancia['sst3']
    primeraEstructura = est2[0]
    vectorProbInicial.loc['FrecuenciasAbsolutas',primeraEstructura]+=1
vectorProbInicial.loc['Probabilidades'] = vectorProbInicial.loc['FrecuenciasAbsolutas']/sum(vectorProbInicial.loc['FrecuenciasAbsolutas'])
vectorProbInicial

	C	E	H
FrecuenciasAbsolutas	121.0	0.0	0.0
Probabilidades	1.0	0.0	0.0

Calcula la matriz que contenga las transiciones entre estados en frecuencias absolutas.

matrizTransicion = pd.DataFrame(0.0,columns=strlist,index=strlist)
matrizTransicion

for ninst in range(numInstancias):
    instancia = dftrain.iloc[ninst,:]
    #secuencia = instancia['seq'] # no lo observamos
    est2 = instancia['sst3']
    for pos in range(len(est2)-1):
        pos1 = est2[pos]
        pos2 = est2[pos+1]
        matrizTransicion.loc[pos1,pos2]+=1
matrizTransicion

	C	E	H
C	1655.0	143.0	104.0
E	144.0	322.0	0.0
H	103.0	1.0	1037.0

Comprueba que la suma de los elementos de la matriz coincide con el conjunto de datos.

print('Suma elementos matriz transicion: {}'.format(matrizTransicion.sum().sum()))
print('Suma elementos en datos: {}'.format(dftrain.shape[0]*(len(dftrain['sst3'].iloc[0])-1)))

Suma elementos matriz transicion: 3509.0
Suma elementos en datos: 3509

Normaliza la matriz de transición para obtener una matriz estocástica por filas.

#Falta normalizar por filas para que sea una matriz estocástica (por filas).
from sklearn.preprocessing import Normalizer
matrizTransicionNormalizada = pd.DataFrame().reindex_like(matrizTransicion)
matrizTransicionNormalizada.iloc[:,:] = Normalizer(norm='l1').fit_transform(matrizTransicion)
matrizTransicionNormalizada

	C	E	H
C	0.870137	0.075184	0.054679
E	0.309013	0.690987	0.000000
H	0.090272	0.000876	0.908852

Calcula la suma de las filas y las columnas de la matriz obtenida para comprobar que es estocástica por filas.

matrizTransicionNormalizada.sum(axis=0)

C    1.269421
E    0.767048
H    0.963531
dtype: float64

matrizTransicionNormalizada.sum(axis=1)

C    1.0
E    1.0
H    1.0
dtype: float64

Ejercicio 7. Predicción de una estructura secundaria.¶

Elige una secuencia del conjunto de validación de forma aleatoria. Vamos a predecir su estructura secundaria de acuerdo con el modelo de Markov oculto que acabamos de ajustar.

A la secuencia le asignaremos la estructura secundaria que maximiza la verosimilitud, de acuerdo con el algoritmo de Viterbi.

Para ello usaremos la librería hidden_markov. La documentación está aquí.

En esta web puedes ver la documentación del algoritmo de Viterbi.

Una vez obtenida la estructura predicha, la compararemos con la estructura anotada en los datos.

Crea un modelo de markov oculto usando la librería hidden_markov

from hidden_markov import hmm

# parámetros del modelo
states = tuple(strlist)
possible_observations = tuple(alfabeto)
start_prob = np.asmatrix(vectorProbInicial.loc['Probabilidades',:].to_numpy())
trans_prob = np.asmatrix(matrizTransicionNormalizada.to_numpy())
em_prob = np.asmatrix(MatrizEmision.to_numpy())

# modelo
mod_HMM = hmm(states=states,
             observations=possible_observations,
             start_prob=start_prob,
             trans_prob=trans_prob,
             em_prob=em_prob)
mod_HMM

<hidden_markov.hmm_class.hmm at 0x7f82179399a0>

Elige una secuencia del conjunto de validación de forma aleatoria.

ntest = dftest.shape[0]
elementoAleatorio = np.random.randint(ntest)
data = dftest.iloc[elementoAleatorio]
observations = data.seq
observations

'SGEADCGLRPLFEKKSLEDKTERELLESYI'

Usa el modelo de markov oculto para predecir la secuencia de estados ocultos asociada a la secuencia de AA observada. Imprime el resultado. Fíjate qué tipo de dato devuelve la librería hidden_markov.

print(mod_HMM.viterbi(observations))

['C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H']

Convierte la predicción anterior en un dato tipo str (cadena).

predicted = ''
for x in mod_HMM.viterbi(observations):
    predicted+=x
predicted

'CCCCCCCCCCHHHHHHHHHHHHHHHHHHHH'

Compara la predicción con la secuencia de estados ocultos anotada: llamamos actual.

actual = data.sst3
print('Actual: '+actual)
print('Predicted:'+predicted)

Actual: CCCCCCCCCCCCHHHCCCCCCHHHHHHHCC
Predicted:CCCCCCCCCCHHHHHHHHHHHHHHHHHHHH

Crea la matriz de confusión para esta observación. Las filas son los valores de actual y las columnas los de predicted.

actualPredicted = pd.DataFrame(0,index = strlist,columns=strlist)
actualPredicted

for i,act in enumerate(actual):
    pre = predicted[i]
    actualPredicted.loc[act,pre]+=1
actualPredicted

	C	H
C	10	10
E	0	0
H	0	10

Ejercicio 8. Matriz de confusión y tasa de error.¶

Haz la predicción de las estructuras secundarias en todo el conjunto de validación y, para cada secuencia, calcula la matriz de confusión asociada a los tres estados \(C\), \(E\) y \(H\) comparando la estructura secundaria predicha y la verdadera.

Acumula las frecuencias de acierto o fallo en cada uno de los 9 elementos de la matriz de confusión al recorrer el conjunto de validación.

confusionTable = pd.DataFrame(0,index = strlist,columns=strlist)
confusionTable

for i in range(ntest):
    data = dftest.iloc[i]
    observations = data.seq
    actual = data.sst3
    predictedHMM = mod_HMM.viterbi(observations)
    predicted = "".join(predictedHMM)
    
    # actualizamos la matriz de confusion
    for i,act in enumerate(actual):
        pre = predicted[i]
        confusionTable.loc[act,pre]+=1
    
confusionTable

	C	E	H
C	666	10	217
E	118	34	12
H	178	0	355

Con esta matriz de confusión agregada, calcula la tasa de error del método.

CFT = confusionTable.to_numpy()
CFTnonDiagonal = CFT.copy()
np.fill_diagonal(CFTnonDiagonal,0)

CFTnonDiagonal

sumaNonDiagonal = CFTnonDiagonal.sum()
sumaNonDiagonal

sumaTotal = CFT.sum()
sumaTotal

tasaError = sumaNonDiagonal/sumaTotal
tasaError

0.33647798742138363

	R	K	D	E	Q	S	C	H	N	T	W	Y	M	A	I	L	F	V	P	G
C	105	151	73	102	67	127	175	46	83	134	22	59	23	129	83	91	77	80	184	212
E	41	12	6	8	12	17	84	5	4	36	18	58	2	34	28	16	36	22	9	18
H	51	53	36	95	37	64	49	44	27	32	20	50	18	126	55	183	32	100	12	57

	R	K	D	E	Q	S	C	H	N	T	W	Y	M	A	I	L	F	V	P	G
C	105	151	73	102	67	127	175	46	83	134	22	59	23	129	83	91	77	80	184	212
E	41	12	6	8	12	17	84	5	4	36	18	58	2	34	28	16	36	22	9	18
H	51	53	36	95	37	64	49	44	27	32	20	50	18	126	55	183	32	100	12	57

Introducción al Aprendizaje Automático