La aplicación de la teoría QSPR-QSAR en la predicción de actividades biológicas (página 7)
Enviado por Eduardo Alberto Castro
9. Liu, R., Paxton, W. A., Choe, S., Ceradini, D., Martin, S. R., Horuk, R., MacDonald, M. E., Stuhlmann, H., Koup, R. A., Landau, N. R., Homozygous Defect in HIV-1 Coreceptor Accounts for Resistance of Some Multiply-Exposed Individuals to HIV-1 Infection. Cell 1996, 86, 367.
10. Samson, M., Libert, F., Doranz, B. J., Rucker, J., Liesnard, C., Farber, C.-M., Saragosti, S., Lapouméroulie, C., Cognaux, J., Forceille, C., Muyldermans, G., Verhofstede, C., Burtonboy, G., Georges, M., Imai, T., Rana, S., Yi, Y., Smyth, R. J., Collman, R. G., Doms, R. W., Vassart, G., Parmentier, M., Resistance to HIV-1 infection in Caucasian individuals bearing mutant alleles of the CCR-5 chemokine receptor gene. Nature 1996, 382, 722.
11. Michael, N. L., Chang, G., Louie, L. G., Mascola, J. R., Dondero,D., Birx, D. L., Sheppard, H. W., The role of viral phenotype and CCR-5 gene defects in HIV-1 transmission and disease progression. Nat. Med. 1997, 3, 338.
12. Imamura, S., Ichikawa, T., Nishikawa, Y., Kanzaki, N., Takashima, K., Niwa, S., Iizawa, Y., Baba, M., Sugihara, Y., Discovery of a Piperidine-4-carboxamide CCR5 Antagonist (TAK-220) with Highly Potent Anti-HIV-1 Activity. J. Med. Chem. 2006, 49, 2784.
13. Imamura, S., Ishihara, Y., Hattori, T., Kurasawa, O., Matsushita, Y., Sugihara, Y., Kanzaki, N., Iizawa, Y., Baba, M., Hashiguchi, S., CCR5 Antagonists as Anti-HIV-1 Agents. 1. Synthesis and Biological Evaluation of 5-Oxopyrrolidine-3-carboxamide Derivatives. Chem. Pharm. Bull. 2004, 52, (1), 63.
14. Imamura, S., Nishikawa, Y., Ichikawa, T., Hattori, T., Matsushita, Y., Hashiguachi, S., Kanzaki, N., Iizawa, Y., Baba, M., Susihara, Y., CCR5 antagonists as anti-HIV-1 agents. Part 3: Synthesis and biological evaluation of piperidine-4-carboxamide derivatives. Bioorganic & Medicinal Chemistry 2005, 13, 397.
15. Seto, M., Aramaki, Y., Imoto, H., Aikawa, K., Oda, T., Kanzaki, N., Iizawa, Y., Baba, M., Shiraishi, M., Orally Active CCR5 Antagonists as Anti-HIV-1 Agents 2: Synthesis and Biological Activities of Anilide Derivatives Containing a Pyridine N-Oxide Moiety. Chem. Pharm. Bull. 2004, 52, (7), 818.
16. Musha, D., Body water in man. I. Total body water in normal subjects and edematous patients. Tohoku J. Exp. Med. 1956, 63, 309.
17. Schuster, D., Laggner, C., Langer, T., Why drugs fail – a study on side effects in new chemical entities. Pharm. Des. 2005, 11, 3545.
18. Stegemann, S., Leveiller, F., Franchi, D., de Jong, H., Lindén, H., When poor solubility becomes an issue: from early stage to proof of concept. Eur J Pharm Sci. 2007, 31, 249.
19. Alsenz, J., Kansy, M., High throughput solubility measurement in drug discovery and development. Adv. Drug Deliv. Rev. 2007, 59, 546.
20. Bhattachar, S. N., Deschenes, L., Wesley, Solubility: it's not just for physical chemists. J. A. Drug Discov. Today 2006, 11, 1012.
21. Balakin, K. V., Savchuk, N. P., Tetko, I. V. , In Silico approaches to prediction of aqueous and DMSO Solubility of drug-like compounds: Trends, problems and solutions. Curr. Med. Chem. 2006, 13, 226.
22. Delaney, J. S., Predicting aqueous solubility from structure. Drug Discovery Today 2005, 10, (4), 289.
23. Duchowicz, P. R., Talevi, A., Bruno-Blanch, L. E., Castro, E. A., New QSPR study for the prediction of aqueous solubility of drug-like compounds. Bioorganic & Medicinal Chemistry 2008, 16, 7944.
24. Bradbury, S. P., Quantitative structure-activity relationships and ecological risk assessment: an overview of predictive aquatic toxicology research. Toxicology Letters 1995, 79, 229.
25. Auer, C. M., Nabholz, J.V., Baetcke, K.P., Mode of action and the assessment of chemical hazards in the presence of limited data: use of structure activity relationships (SARs) under TSCA. Environ. Health Perspect. 1990, 87, (5), 183.
26. Verhaar, H. J. M., van Leeuwen, C.J., Hermens, J.L.M., Classifying environmental pollutants. 1: Structure-activity relationships for prediction of aquatic toxicity. Chemosphere 1992, 25, 471.
27. Bradbury, S. P., Predicting modes of toxic action from chemical structure: an overview. SAR and QSAR Environ, Res. 1994, 2, 89.
Capítulo 4. Resultados
1. Introducción
En este capítulo se reportan los resultados más relevantes obtenidos a partir de la implementación de los distintos métodos de clasificación de objetos estudiados en el Capítulo 2, aplicados en el Trabajo de Tesina para el armado de conjuntos moleculares balanceados. Además, se describe cada uno de los algoritmos que fue necesario programar en Matlab para obtener los resultados que a continuación se presentan.
A la hora de armar un conjunto molecular balanceado para aplicaciones QSAR-QSPR, se busca que los errores cometidos por el modelo en la etapa de calibración sean similares a los encontrados durante la etapa de validación. Si esto se cumple, el modelo funciona con un carácter más general y predictivo sobre los datos, y se asigna igual preferencia al ajuste de los datos en los conjuntos de calibración y validación. De nada sirve, por ejemplo, ajustar muy bien el conjunto de calibración si luego las predicciones alcanzadas en el conjunto de validación presentan errores grandes, o viceversa. Sí resulta conveniente tratar de ajustar ambos conjuntos por igual, con error comparable.
Después de realizar la aplicación de alguno de los métodos de clasificación vistos, quedan definidos los conjuntos de calibración y validación. El error de cada conjunto está asociado a la propiedad predicha para las moléculas pertenecientes a dicho conjunto, resultante de la aplicación del modelo QSAR-QSPR consistente en una regresión lineal univariable. En este modelo, la variable dependiente es la propiedad experimental () y la variable independiente es el descriptor molecular que correlacione mejor con la propiedad:
(1)
donde a y b son los coeficientes de regresión, es el descriptor molecular, y la propiedad predicha. El error lo cuantificamos con el parámetro rrcm, la raíz cuadrada del residuo cuadrático medio:
(2)
(3)
Aquí, N es el número de moléculas ajustadas y es el residuo para la molécula i.
El problema que se presenta es el siguiente. Como se aprecia de la Ecs. (2) y (3), el valor del parámetro rrcm que mide el error en el conjunto depende de los valores numéricos que adopta la propiedad experimental modelada y de los valores que adopta el descriptor usado en la Ec. (1). Por tanto, si queremos que los errores de calibración y validación tengan magnitud comparable, debemos considerar estos dos factores en el diseño del conjunto balanceado. Obviamente, a la hora de armar un modelo QSAR-QSPR no es lícito considerar la propiedad experimental del conjunto de validación, pues este conjunto se utiliza solamente para probar la habilidad predictiva del modelo sin considerarlo durante la etapa de entrenamiento.
Lo anteriormente expuesto conduce al motivo principal por el cual decidimos usar dos conjuntos de calibración (cal1 y cal2) y un conjunto de validación (val) durante los análisis: usamos cal1 para calibrar el modelo QSAR-QSPR, comparamos el error rrcm de cal1 y cal2 para este modelo y así comprobar si se tienen conjuntos balanceados, mientras que con el conjunto val solamente verificamos el poder predictivo. De esta manera evitamos el inconveniente anteriormente planteado. Además, el conjunto cal2 sirve para pre-validar la relación cuantitativa obtenida, y constituye una suerte de transición menos abrupta entre los conjuntos de calibración y validación.
En cada caso, el número de moléculas incluidas en los conjuntos cal1 y cal2 representan el 70% del número total de moléculas, mientras que el 30% restante corresponde a moléculas de validación. Se utiliza igual número de moléculas en cal1 y cal2. Finalmente, cabe mencionar que el gran tamaño de los conjuntos moleculares correspondientes a las tres propiedades ensayadas en este trabajo, que poseen un número de moléculas superior a 100, hace posible usar dos conjuntos de calibración.
Una vez entendida la manera de obtener el error en los conjuntos cal1, cal2, y val, procedemos a describir los algoritmos basados en los distintos métodos de clasificación estudiados. Estos algoritmos permiten seleccionar los mejores descriptores para clasificar moléculas en cada método, a partir del análisis de 1497 descriptores provistos por Dragon para el conjunto molecular ensayado. En todos los casos, para la medida de distancia entre pares de moléculas se utiliza la distancia Euclídea, aunque también podría recurrirse a otras alternativas (ver Apéndice, sección II).
2. Algoritmos y Criterio Matemático Utilizados
2.1. Algoritmo clusterskmeans.m
Permite realizar las mejores agrupaciones de moléculas a partir de N consideradas, a través de la aplicación del método K-Medias y la exploración de D=1497 descriptores moleculares. Descarta los descriptores que conduzcan a valores negativos del parámetro silueta medio. A partir de los K grupos generados con cada descriptor clasificador y al considerar que los integrantes del grupo son equivalentes entre sí, se extraen moléculas representativas de cada uno de ellos y se arman los conjuntos cal1, cal2, y val que respeten las proporciones señaladas anteriormente (70% en cal1 y cal2, y 30% en val). Luego, se obtiene el parámetro rrcm en cada conjunto según la Ec. (2); este error se calcula con el descriptor molecular que correlacione mejor con la propiedad en el modelo de la Ec. (1). En consecuencia, el descriptor que se utiliza para realizar la clasificación molecular no necesariamente debe ser el mismo a utilizado en la ecuación QSAR-QSPR para calcular rrcm.
Es posible involucrar un mayor número de descriptores clasificadores en clusterskmeans.m (o en los algoritmos explicados en las siguientes secciones), que conduce a una clasificación más estricta de los datos, pero por razones de falta de tiempo no incluimos dichos resultados. A pesar de ello, los resultados presentados aquí no varían demasiado para descriptores clasificadores adicionales. Además, para simplificar el análisis tampoco consideramos modelos QSAR-QSPR que involucren un mayor número de descriptores en la Ec. (1).
El algoritmo se ejecuta con la sentencia siguiente:
[Resultkmeans]=clusterskmeans(p, tot, nclusters, percent); (4)
Aquí, "p" es la propiedad experimental objeto de estudio, "tot" es la matriz NxD, "nclusters" es K, y "percent" corresponde al porcentaje de moléculas de calibración (70%). El resultado "Resultkmeans" es una matriz; el formato que se utiliza para presentar los resultados del algoritmo clusterskmeans.m es el mismo al usado para los algoritmos explicados en las siguientes secciones.
Por ejemplo, para el conjunto de 166 solubilidades acuosas y la creación de 3 grupos, se muestran las primeras 29 filas de "Resultkmeans" en la Tabla 1.
Tabla 1. Resultado obtenido con clusterskmeans.m en 166 solubilidades acuosas. nclusters=3, percent=70. M es el número de moléculas de cada grupo.
M | |||||||||||||
808 | 557 | 1.182 | 1.244 | 1.643 | 0.062 | 4.98 | 87 | 65 | 14 | ||||
1377 | 557 | 1.185 | 1.269 | 1.620 | 0.085 | 6.67 | 30 | 1 | 135 | ||||
741 | 558 | 1.183 | 1.288 | 1.608 | 0.106 | 8.19 | 96 | 61 | 9 | ||||
895 | 1497 | 1.228 | 1.376 | 1.227 | 0.148 | 10.77 | 68 | 84 | 14 | ||||
646 | 1025 | 1.202 | 1.388 | 1.807 | 0.186 | 13.39 | 92 | 15 | 59 | ||||
1255 | 559 | 1.145 | 1.325 | 1.646 | 0.180 | 13.56 | 144 | 17 | 5 | ||||
766 | 1024 | 1.199 | 1.395 | 1.878 | 0.196 | 14.05 | 92 | 15 | 59 | ||||
1205 | 255 | 1.183 | 1.391 | 1.745 | 0.209 | 15.00 | 31 | 76 | 59 | ||||
887 | 1497 | 1.213 | 1.427 | 1.199 | 0.214 | 15.02 | 84 | 73 | 9 | ||||
951 | 1497 | 1.208 | 1.424 | 1.221 | 0.216 | 15.15 | 76 | 81 | 9 | ||||
245 | 1497 | 1.178 | 1.445 | 1.187 | 0.267 | 18.47 | 61 | 74 | 31 | ||||
130 | 1497 | 1.180 | 1.454 | 1.183 | 0.274 | 18.86 | 35 | 76 | 55 | ||||
861 | 1497 | 1.169 | 1.453 | 1.210 | 0.284 | 19.55 | 12 | 52 | 102 | ||||
440 | 1497 | 1.166 | 1.464 | 1.188 | 0.299 | 20.39 | 6 | 105 | 55 | ||||
439 | 1497 | 1.167 | 1.470 | 1.198 | 0.303 | 20.62 | 6 | 59 | 101 | ||||
904 | 1497 | 1.154 | 1.475 | 1.198 | 0.321 | 21.76 | 88 | 64 | 14 | ||||
1236 | 113 | 1.008 | 1.426 | 1.369 | 0.417 | 29.27 | 6 | 96 | 64 | ||||
1065 | 113 | 1.012 | 1.432 | 1.378 | 0.420 | 29.31 | 111 | 36 | 19 | ||||
55 | 113 | 1.021 | 1.445 | 1.335 | 0.424 | 29.32 | 64 | 16 | 86 | ||||
1144 | 113 | 1.017 | 1.445 | 1.367 | 0.428 | 29.60 | 105 | 21 | 40 | ||||
415 | 113 | 1.002 | 1.428 | 1.387 | 0.426 | 29.82 | 34 | 55 | 77 | ||||
262 | 113 | 1.006 | 1.434 | 1.356 | 0.429 | 29.89 | 139 | 7 | 20 | ||||
25 | 113 | 0.986 | 1.411 | 1.413 | 0.424 | 30.07 | 27 | 77 | 62 | ||||
826 | 113 | 1.007 | 1.445 | 1.360 | 0.438 | 30.30 | 76 | 76 | 14 | ||||
888 | 113 | 1.003 | 1.439 | 1.387 | 0.437 | 30.34 | 18 | 77 | 71 | ||||
431 | 113 | 1.011 | 1.451 | 1.341 | 0.441 | 30.38 | 56 | 105 | 5 | ||||
656 | 113 | 0.977 | 1.406 | 1.434 | 0.429 | 30.48 | 108 | 9 | 49 | ||||
639 | 1147 | 1.059 | 1.524 | 1.255 | 0.465 | 30.51 | 162 | 3 | 1 | ||||
580 | 113 | 1.009 | 1.456 | 1.357 | 0.447 | 30.68 | 73 | 73 | 20 |
En esta tabla se definen los parámetros y necesarios para el análisis.
(5)
(6)
Además, representa el número de moléculas presentes en el grupo
Los resultados de "Resultkmeans" se hallan ordenados según creciente. Este parámetro mide la diferencia porcentual en el error para los dos conjuntos de calibración: si los dos conjuntos están balanceados, este porcentaje debe ser bajo. Ahora, se aprecia que se tienen varias soluciones posibles (filas) en la matriz de resultados, consecuencia de explorar 1497 descriptores moleculares. Por tanto, debe especificarse un criterio matemático que permita rescatar una solución satisfactoria (balanceada) entre los varios resultados disponibles. La definición de este criterio será igualmente aplicable a los algoritmos presentados en las siguientes secciones, en vista que los resultados presentados por los diferentes métodos poseen el mismo formato al indicado en la Tabla 1.
2.2. Criterio Matemático
La especificación del criterio tiene su base en el estudio del comportamiento numérico de los datos, por lo que es aplicable en principio a cualquier conjunto molecular ensayado (cualquier propiedad experimental, independiente de la diversidad estructural molecular). La solución balanceada se busca entre los varios resultados posibles con el siguiente procedimiento:
a- se ordena la matriz de resultados según creciente.
b- la solución principal (prin) es la primer solución que tenga y sus parámetros son y
c- es posible encontrar una solución diferente (secundaria, sec) a la solución principal si para lo cual se define:
(7)
La solución secundaria es una situación de compromiso entre los siguientes requisitos, en ese orden:
i. con y alto
ii. y alto
iii. de la solución secundaria debe ser bajo.
iv. bajo
En algunos casos, si se cumple apreciablemente iii) para la solución secundaria y entonces igual se acepta este resultado como solución.
Por ejemplo, si se aplica el criterio a los datos de la Tabla 1, se encuentra que la solución principal se caracteriza con:
y
Es posible encontrar una solución secundaria (sombreada en gris) que cumple las condiciones i.-iv. y posee (2.52%):
y
Este criterio matemático adoptado para el análisis de las soluciones demuestra funcionar bastante bien y de carácter general para el armado de conjuntos moleculares balanceados, en las tres propiedades estudiadas en este trabajo. Como se observa de las condiciones del criterio, ninguna de ellas considera a parámetros derivados del conjunto de validación, por lo que estas reglas resultan válidas. El criterio tiene en cuenta que sea bajo entre los varios resultados. Además, el hecho de considerar que se atribuye a que esta elección evita soluciones en las que se ajusta extremadamente bien ambos conjuntos de calibración y se ajusta peor el conjunto de validación. Por ejemplo, este es el caso para la primera fila:
y
En definitiva, el criterio establecido no sólo permite rescatar una solución aceptable de la matriz de resultados, sino que también permite arribar a particiones moleculares en las que los modelos QSAR-QSPR resultantes de dichas particiones resultan más predictivos en el conjunto de validación.
2.3. Algoritmo clustersknn.m
La implementación del método K-Vecinos Más Cercanos se efectúa a través del desarrollo y posterior aplicación del algoritmo clusterknn.m, que funciona de manera similar a clusterskmeans.m y busca el descriptor clasificador entre los D disponibles que consiga las mejores agrupaciones. Sin embargo, a diferencia de K-Medias, la aplicación del método requiere conocer de antemano un conjunto de entrenamiento. Si se utiliza como conjunto de entrenamiento los centroides proporcionados por el método K-Medias, entonces los grupos formados por K-NN y K-Medias coinciden. Este resultado permitió definir el conjunto de entrenamiento a utilizar en la técnica K-NN, es decir, en vez de usar centroides de K-Medias se definieron nuevos centroides. Por ejemplo, para el caso de un descriptor clasificador, si se quieren formar dos grupos, los centros se ubican en los valores máximo y mínimo del descriptor. Si se busca generar tres grupos, entonces los centros se posicionan en los valores máximo, medio y mínimo del descriptor considerado. Para mayor cantidad de grupos, los centros se ubican en iguales intervalos del descriptor. La siguiente figura presenta la definición de los nuevos centros.
Figura 1. Esquema de selección de nuevos centroides de d1 en clustersknn
El algoritmo se ejecuta con la sentencia siguiente:
[Resultknn]=clustersknn(p, tot, nclusters, percent); (8)
El criterio matemático adoptado para la elección del resultado se mantiene. La descripción de las variables de entrada del algoritmo es la misma que para clusterskmeans.m.
2.4. Algoritmo clusterslda.m
La implementación del Análisis Discriminante Lineal a través del algoritmo clusterlda.m se lleva a cabo de la misma manera que para clusterknn.m. Al igual que en K-NN, el conjunto de entrenamiento está constituido por los nuevos centroides definidos en la Figura 1, y establecen las clases con las cuales el método discriminará a las moléculas.
El algoritmo se ejecuta con la sentencia siguiente:
[Resultlda]=clusterslda(p, tot, nclusters, percent); (9)
2.5. Algoritmo clusterspca.m
La aplicación del método PCA utiliza el algoritmo clusterspca.m. En este caso, se establecen las agrupaciones luego de analizar el signo de las coordenadas de los componentes principales. Por tanto, si L es el número de componentes principales, será el número de agrupamientos obtenidos. A modo de ejemplo, si se utiliza únicamente el primer componente principal se obtienen dos grupos, uno correspondiente a moléculas con coordenadas PC1>0 y otro correspondiente a moléculas con PC1<0. La Figura 2 ilustra esta idea.
Figura 2. Especificación de agrupaciones en clusterspca.m
Con objeto de calcular los componentes principales, se requiere disminuir la dimensión de la matriz de descriptores a tratar. En nuestro caso D=1497, y la aplicación directa de la técnica PCA no es posible. Por tanto, se considera solamente un conjunto de descriptores clasificadores tal que el coeficiente de correlación entre cada par de descriptores i y j () sea menor al valor límite 0.5. De esta manera, se evita utilizar descriptores que se encuentren muy correlacionados, pues varían de la misma forma y no contribuyen al proceso de clasificación.
Una vez que se tiene una matriz reducida de descriptores linealmente independientes, se buscan todas las combinaciones posibles de 2 descriptores, cada una de las cuales permite obtener al primer y segundo componente principal y, por tanto, a 2 y 4 grupos, respectivamente. Para obtener el tercer componente principal, se requiere la combinación de 3 descriptores. La búsqueda combinatorial exacta de 3 descriptores requiere una mayor demanda computacional, y el problema resulta más complejo aún si el número de descriptores aumenta. Por tanto, una vez que se elige la mejor solución de 2 descriptores clasificadores, se busca el tercer descriptor que mejor se combine con los 2 descriptores previamente elegidos.
El algoritmo se ejecuta con la sentencia siguiente:
[Resultpca]=clusterspca(p,tot,rlim,nclusters,percent); (10)
donde rlim es 0.5 en nuestro caso.
2.6. Método HCA
Al tratarse el Análisis de Agrupamiento Jerárquico de un método gráfico, no se lo utiliza en este trabajo por la dificultad que presenta a la hora de programar su algoritmo para la exploración de más de mil descriptores clasificadores. Sí puede aplicarse, en cambio, al tratamiento de unos pocos descriptores de unas pocas moléculas, lo cual no es ninguno de los casos aquí tratados.
3. Resultados
A continuación se presentan en tablas los principales resultados encontrados para cada propiedad ensayada, luego de aplicar los algoritmos basados en las distintas técnicas de clasificación. Cada uno de estos resultados se basó en el criterio matemático establecido previamente para seleccionar una solución balanceada.
Por otro lado, también se calcularon algunas soluciones en las que se aplican los métodos K-Medias, K-NN y LDA con el mismo descriptor clasificador. Esto se hizo a fines comparativos y de discusión de los resultados hallados si se considera la misma variable clasificadora.
El número de agrupamientos considerados en los cálculos para los distintos métodos fue K = 2, 3, 4, 10, 15, 20, 25, y 30, a excepción del método PCA, que admite 2, 4, 8, y 16.
Tabla 2. Actividad anti-VIH
Tabla 3. Solubilidades acuosas
Tabla 4. Toxicidades acuosas
4. Discusión
Se observa de las Tablas 2-4 que, para un número de agrupamientos determinado, los errores obtenidos en cal1 con los métodos de clasificación K-Medias, K-NN, LDA y PCA tienden a ser parecidos en las tres propiedades ensayadas, sin observar discrepancias apreciables. Lo mismo sucede para el caso de las comparaciones en cal2 y val. A su vez, los conjuntos de calibración y validación tienden a tener errores que no se diferencian demasiado entre sí, por lo que los conjuntos obtenidos para estos cuatro métodos tienden a ser balanceados.
Ahora, cómo se explica que diferentes métodos tiendan a generar soluciones parecidas, para un número de agrupamiento determinado? Para responder a este interrogante, debemos plantear lo que sucede con el desempeño de estos métodos si utilizan el mismo descriptor clasificador. Las Tablas 5-7 muestran discrepancias mayores en el error de los conjuntos de calibración y validación, dependiente del método que se use, por lo que no se obtienen conjuntos balanceados en estas condiciones. Esto sucede especialmente en los conjuntos moleculares de solubilidades acuosas y actividades anti-VIH-1, que al parecer son conjuntos más heterogéneos para el modelo QSAR-QSPR.
La conclusión a la que se llega es que la aplicación del criterio matemático para la selección de la mejor solución entre varias posibilidades, que surgen de explorar un gran número de descriptores, permite unificar el funcionamiento de los métodos clasificadores aquí estudiados. Esto surge de las tres propiedades estudiadas en el trabajo actual y podría extenderse al tratamiento de otras propiedades.
Capítulo 5. Conclusiones
El objetivo principal del Trabajo de Tesina consiste en establecer una correcta clasificación molecular que permita seleccionar de manera racional conjuntos moleculares balanceados para su posterior aplicación en la Teoría QSAR-QSPR.
En función de los resultados obtenidos, se puede concluir que:
A la hora de armar conjuntos moleculares balanceados, resulta necesario buscar el mejor resultado de clasificación de manera que tenga en cuenta tanto las características estructurales de las moléculas como a la propiedad experimental objeto de estudio.
El criterio matemático establecido resulta funcionar de manera general sobre los distintos datos, y permite unificar los métodos estudiados. Esto es valioso si para un determinado problema de clasificación uno no sabe a cual método recurrir.
Al tratarse HCA de un método gráfico, no se lo utilizó en este trabajo por la dificultad que presenta a la hora de programar su algoritmo para la exploración de más de mil descriptores.
Apéndice
I. Discriminación y Clasificación en LDA
Si consideramos una matriz de atributos o variables independientes para cada objeto que pertenece a una clase G1 determinada, este conjunto de muestras es el llamado conjunto de entrenamiento o calibración. El problema consiste entonces en hallar una buena predicción de la clase G1 de un objeto considerado, con el uso de la misma distribución del conjunto de entrenamiento, a través de los valores de las variables
La obtención de la función discriminante incluye una serie de aproximaciones, a saber:
i. las funciones densidad de probabilidad, y son ambas distribuciones normales, G1.
ii. si y son los parámetros media y covarianza de G1=0 y G1=1, respectivamente, se supone que las covarianzas son iguales.
Función Discriminante Lineal
Es posible encontrar una relación lineal que caracterice a cada objeto según los valores de sus atributos. De esta manera, si tenemos un conjunto de N objetos de los que se conocen D variables explicativas, y se observa que N1 de ellos pertenecen a la clase C1, y los N2 restantes a la clase C2, con N1+N2=N, es posible construir una función lineal en base a las D variables y puede usarse para predecir si pertenece a un grupo u otro con una probabilidad determinada. En la función lineal:
(AI.1)
es una variable clasificadora, es la i-ésima variable atributo, y es su coeficiente. El objetivo principal de tal función lineal desde el punto de vista de la varianza consiste en responder a la pregunta de si dos o más grupos son significativamente distintos uno a otro respecto a la medida de una variable en particular. Debe tenerse presente que si la media de una variable es significativamente diferente en varios grupos, puede decirse que esta variable discrimina bien entre grupos.
En caso de que sea posible identificar más de dos grupos en los datos, pueden estimarse funciones discriminantes múltiples, cada una de ellas similares a la presentada en la Ec. (AI.1). Por ejemplo, cuando se tienen tres grupos, puede estimarse: 1) una función para discriminar entre el grupo 1 y los grupos 2 y 3 combinados; y 2) otra función para discriminar entre el grupo 2 y el grupo 3. Además, se pueden considerar sólo las funciones discriminantes múltiples que resulten más significativas: si se observan los coeficientes estandarizados de las variables de cada una de las funciones escogidas, cuanto mayor sean estos coeficientes, más alta es la contribución a la discriminación especificada. Finalmente, pueden considerarse las medias de las funciones discriminantes significativas para analizar entre cuales grupos éstas discriminan.
II. Tipos de Medida de Distancia
Sea una matriz X (mxn) que es tratada como m vectores fila x1, x2, …, xm. Varios tipos de medida de distancia que son posibles definir para el par de objetos r y s, con xr y xs, se incluyen a continuación:
i. Distancia Euclídea: es la distancia entre dos puntos que se mide en el espacio euclídeo, se define como
(AII.1)
ii. Distancia Euclídea Estandarizada: cada coordenada en la suma cuadrática se pesa inversamente por la varianza muestral de esa coordenada.
(AII.2)
donde es la matriz con elementos diagonales dados por vj2, que se refiere a la varianza de la variable sobre los m objetos.
iii. Distancia Mahalanobis: es una forma de determinar la similitud entre dos variables aleatorias multidimensionales. A diferencia de la distancia Euclídea, esta medida tiene en cuenta la correlación de las variables.
(AII.3)
donde es la matriz de covarianza muestral.
iv. Distancia Manhattan: aquí la distancia entre dos puntos es la suma de las diferencias (absolutas) de sus coordenadas.
(AII.4)
v. Distancia Minkowski
(AII.5)
En el caso especial p=1, la distancia Minkowski coincide con la distancia Manhattan, y para el caso especial p=2, la distancia Minkowski coincide con la distancia Euclídea.
vi. Distancia Coseno: uno menos del ángulo incluido entre los puntos (en forma de vector).
(AII.6)
vii. Distancia de Correlación: uno menos la correlación entre los puntos (tratado como secuencias de valores).
(AII.7)
donde y
viii. Distancia Hamming: es el porcentaje de coordenadas que difieren.
(AII.8)
III. Métodos de Enlace o Vinculación
En el método HCA, se utiliza una función vinculante que crea un árbol de agrupamiento jerárquico a partir de las distancias entre pares de objetos previamente obtenidas. La función puede utilizar diversos métodos de vinculación, los cuales difieren en la forma que se calculan las distancias entre los agrupamientos.
La solución que se obtiene en HCA es una matriz (m-1)x3 llamada Q, donde m es el número de observaciones en el conjunto original de datos. Las primera y segunda columnas de Q contienen a los índices de los agrupamientos vinculados de a pares, para formar el árbol binario. La tercera columna contiene la distancia de vinculación entre los agrupamientos formados.
La siguiente notación se utiliza para describir los distintos métodos de vinculación:
Un agrupamiento r es formado a partir de los agrupamientos p y q
nr es el número de objetos en el agrupamiento r.
xir es el i-ésimo objeto en el agrupamiento r.
a. Vinculación individual, también llamado vecino más cercano, utiliza la menor distancia entre dos objetos en los dos agrupamientos.
(AIII.1)
b. Vinculación completa, también llamado vecino más lejano, utiliza la mayor distancia entre dos objetos en los dos agrupamientos.
(AIII.2)
c. Vinculación promedio, utiliza la distancia promedio entre todos los pares de objetos en cualquiera de los dos agrupamientos.
(AIII.3)
d. Vinculación centroide, utiliza la distancia Euclídea entre los centroides de los dos agrupamientos.
(AIII.4)
donde y se refiere a la distancia Euclídea.
e. Vinculación media, utiliza la distancia Euclídea entre los centroides ponderados de los dos agrupamientos,
(AIII.5)
donde y son los centroides pesados para los agrupamientos r y s. Si el agrupamiento r fue creado por la combinación de los agrupamientos p y q, es definido recursivamente como
f. Vinculación de Ward, utiliza la suma incremental de los cuadrados, es decir, el incremento en la suma total de los cuadrados dentro del agrupamiento, como resultado de la unión de dos grupos. La suma de los cuadrados dentro del agrupamiento es definida como la suma del cuadrado de las distancias entre todos los objetos en el agrupamiento y el centroide del agrupamiento. La distancia equivalente es:
(AIII.6)
g. Promedio ponderado de vinculación, utiliza una definición recursiva para la distancia entre dos agrupamientos. Si el agrupamiento r fue creado mediante la combinación de los agrupamientos p y q, la distancia entre r y otro agrupamiento s se define como el promedio de las distancias entre p y s y la distancia entre q y s:
(AIII.7)
IV. Eliminación de Mínimos Locales y Descripción del Cálculo Iterativo en el Método K-Medias.
Eliminación de mínimos locales
Al igual que sucede en muchos otros problemas de optimización numérica, la solución que se alcanza con el método K-Medias depende a menudo del punto de partida, en este caso la posición inicial del centroide de cada agrupamiento. Es posible así alcanzar un mínimo local, donde la reasignación de cualquier punto a un nuevo agrupamiento debería incrementar la suma total de distancias centroide-punto, pero donde puede existir realmente una mejor solución. Para solucionar este problema, es posible especificar en el método el número de "réplicas", es decir, el número de veces en que se repetirá el proceso de agrupación, cada uno con un nuevo conjunto de posiciones iniciales del centroide del agrupamiento. Por supuesto, la mejor solución será aquella para la cual la suma de las distancias centroide-punto para cada uno de los agrupamientos sea mínima.
Descripción del algoritmo
El algoritmo consta de dos partes:
Primera fase. Cada iteración consiste en la reasignación colectiva de elementos al centroide del agrupamiento más cercano, todos a la vez, seguida de un nuevo cálculo de las posiciones de los centroides. La primera fase ocasionalmente converge a soluciones que son un mínimo local; es más probable alcanzar un mínimo global si se trabaja con pequeños grupos de datos. La fase de actualización colectiva es rápida, pero posiblemente sólo aproxime una solución que sea el punto de partida de la segunda fase.
Segunda fase. Los elementos son reasignados individualmente si con ello se reduce la suma de distancias, y en cada reasignación se calcula la ubicación del centroide del agrupamiento. Cada iteración consiste en la ubicación de todos los elementos. En esta fase la solución converge a un mínimo local, aunque puede haber otro mínimo local con menor suma total de distancias. Generalmente el problema de hallar un mínimo global puede ser resulto únicamente por medio de una selección exhaustiva de los puntos de partida, aunque la utilización de varias réplicas con puntos de partida aleatorios generalmente converge a una solución que es un mínimo global.
V. Definición del Parámetro Silueta
La definición de los valores silueta es la siguiente: dado un grupo G1 y un objeto i asignado a éste, la disimilitud promedio de i para todos los objetos j en G1 está dada por:
número de objetos en G1 (AV.1)
donde es la distancia de cada objeto i a cada objeto j en G1. La menor disimilitud correspondiente a i respecto a cualquier otro agrupamiento, b(i), también es calculada. Si i es más similar a los objetos en un grupo G2 que a los del grupo G1, entonces:
número de objetos en G2 (AV.2)
Por lo tanto, el valor silueta s(i) definido para un objeto i es:
(AV.3)
donde se refiere al mayor valor entre y
VI. El Clasificador del Método K-Vecinos Más Cercanos
Si se quiere conocer la clase a la que pertenece un objeto dado, entre varias clases posibles, se introduce el concepto de clasificador. Un clasificador es una función que asigna un objeto a una clase determinada, para lo cual se basa en el conocimiento de sus variables atributo. Existen dos tipos de clasificadores:
Paramétricos: asumen que la distribución estadística que sigue el conjunto de variables es conocido, y trata de estimar los parámetros de dicha distribución.
No-paramétricos: no asume ninguna distribución en particular. El clasificador se construye únicamente con los datos del conjunto de entrenamiento.
Entre los clasificadores no-paramétricos, el más conocido es el basado en el método K-vecinos más cercanos: si Ki es el número de objetos que pertenecen a la clase Gi entre los K vecinos más cercanos al objeto considerado x, la probabilidad a posteriori (la probabilidad de que la clase sea Gi cuando x se describe con el conjunto de variables d ) se estima como De esta manera, el clasificador asigna x a la clase más frecuente entre sus K vecinos más cercanos, según una cierta medida de distancia.
VII. Algunos Algoritmos Utilizados en Matlab
Algoritmo clusterskmeans.m
function [Resultkmeans] = clusterskmeans(p,tot,nclusters,percent)
% INPUT:
% p = experimental property
% tot = data matrix with pool of descriptors
% nclusters = number of clusters to be created
% percent = percent of total compounds to be used as part of train and
% val sets
% OUTPUT:
% Resultkmeans = returns in each row in the following order:
% the descriptor number, mean silhouette value, minimum silhouette
% value, rmse(train), rmse(val), rmse(test), rmse(train)-mse(val),
% the within-cluster sums of point-to-centroid distances, the
% sum of within-cluster sums of point-to-centroid distances, and
% the number of objects in each cluster
% Pablo R. Duchowicz
% INIFTA, La Plata, Argentina
% Created: 2nd March 2011
warning off
resj=[];
[d1, d2]=size(tot);
for j=1:d2
apt=0;
vecto=[j];
[IDX0, C, sumd, D, silh3, msilh3, ncomp] = kmca(tot(:,vecto), nclusters, 'distance', 'sqEuclidean', 'off', 50);
if msilh3==-999
apt=1;
else
if msilh3==1
apt=1;
else
if silh3>0
apt=0;
else
apt=1;
end
end
end
if apt==0
res=sumd';
[Idmol Ncomp] = idmol(IDX0, nclusters);
vari=0;
for k=1:nclusters
if Ncomp(k)==0
vari=1;
end
end
if vari==0
[train, val, test] = analysis2(Idmol,percent);
[Res] = rmsr(p(train), [1], tot0(train,:), 1);
[rsstrain rssval rsstest] = trainvaltestrss(train,val,test,p,tot0,Res(2));
resj=[resj;[vecto,msilh3,min(silh3),Res(2),rsstrain,rssval,rsstest,abs(rsstrain-rssval), (1-rsstrain/rssval)*100, res,sum(res),ncomp]];
end
end
end
if size(resj,1)>0
resj=sortrows(resj,(dv+9));
end
Resultkmeans=resj;
Algoritmo clustersknn.m
function [Resultknn] = clustersknn(p,tot,nclusters,percent)
% INPUT:
% p = experimental property
% tot = data matrix with pool of descriptors
% nclusters = number of clusters to be created
% percent = percentage of total compounds to be used as part of train and
% val sets
% OUTPUT:
% Resultknn = returns in each row the results in the following order:
% the classifying descriptor's number, the best descriptor of
% linear model, rmse(train), rmse(val), rmse(test), rmse(train)-
% rmse(val), the percentage difference between train and val, and
% the number of objects in each cluster
% Pablo R. Duchowicz
% INIFTA, La Plata, Argentina
% Created: 15th March of 2011
warning off
resj = [];
[d1, d2] = size(tot);
for j=1:d2
vecto = [j];
[CN] = centers(tot(:,vecto),nclusters);
c = knnclassify(tot(:,vecto),CN,(1:nclusters)');
[Idmol Ncomp] = idmol(c, nclusters);
vari=0;
for k = 1:nclusters
if Ncomp(k)==0
vari=1;
end
end
if vari == 0
[train, val, test] = analysis2(Idmol,percent);
[Res] = rmsr(p(train), [1], tot(train,:), 1);
[rsstrain rssval rsstest] = trainvaltestrss(train,val,test,p,tot,Res(2));
resj=[resj;[vecto,Res(2),rsstrain,rssval,rsstest,abs(rsstrain-rssval),(1-rsstrain/rssval)*100,Ncomp]];
end
end
if size(resj,1)>0
resj=sortrows(resj,(dv+7));
end
Resultknn=res
Página anterior | Volver al principio del trabajo | Página siguiente |