Programa RStudio

La interfaz gráfica del programa R, llamado RStudio, puede ser instalado en los sistemas operativos Windows, MacOS y Linux (https://www.rstudio.com/products/rstudio/) (Figura 1).

Figura 1 Interfaz gráfica de RStudio 

Paqueterías de R

Las paqueterías o librerías disponibles para análisis de regresión son numerosas, las empleadas en este estudio QSAR son:

Modelo QSAR/QSPR

Las herramientas matemáticas empleadas en QSAR incluyen: análisis de regresión lineal y no lineal, análisis de componentes principales (PCA), redes neuronales artificiales (ANN), análisis comparativo de campos moleculares (CoMFA), y análisis comparativo de índices de similitud molecular (CoMSIA), etc. (^{Dearden & Cronin, 2006}). La metodología de este trabajo se enfoca en el empleo de la regresión lineal múltiple, donde la actividad o propiedad es la variable dependiente, y los descriptores son las variables independientes.

Las instrucciones o comandos empleados para llevar a cabo el estudio QSAR/QSPR en R se especifican en el Anexo 1 y los ocho pasos centrales del diagrama de flujo para la generación del modelo se describen a continuación.

Paso 1. Librerías. Cargar las librerías de R requeridas para el estudio.

Paso 2. Documentos. Cargar la información de actividad o propiedad de las moléculas a analizar en un archivo de texto arreglado en columnas (ver Figura 2A) y el archivo de coordenadas de los compuestos en un formato utilizado en bases de datos, por ejemplo el formato “.sdf” o “.smi” (ver Figura 2B).

Figura 2 Ejemplo del contenido de los archivos (A) En columnas separadas por comas “.cvs” y (B) de coordenadas “.sdf”. 

Paso 3. Cálculo de descriptores. La librería rcdk de R puede calcular 265 descriptores disponibles de cinco categorías distintas: [1] híbridos, [2] constitucionales, [3] topológicos, [4] electrónicos y [5] geométricos. La instrucción permite generar los descriptores deseados, omitiendo aquellos cuyos valores no puedan ser asignados. La información obtenida se almacena en una tabla con la instrucción write (ver Nota 1 del Anexo 1).

Paso 4. Descriptores colineales. La colinealidad de los descriptores en un mismo modelo produce una variación excesiva cuando se cambian las i-ésimas observaciones. La colinealidad implica una dependencia casi lineal entre los descriptores x_j, e impacta significativamente en la estimación de los coeficientes (^{Montgomery C.D., Peck, & Vining, 2002}). La colinealidad entre descriptores se puede evaluar mediante la matriz triangular simétrica de correlación, la cual es de tamaño j×j y contiene los valores de los coeficientes de correlación de Pearson r para cada par de los j-ésimos descriptores.

El coeficiente de correlación de Pearson r es la covarianza normalizada de dos variables, por lo que toma valores en el intervalo de -1 a 1. Cuando el valor absoluto de r tiende a uno, existe una alta correlación lineal entre las variables analizadas (^{Montgomery, Peck, & Vining, 2002}). Cuando r→0 se dice que no hay correlación. Uno de los criterios empíricos utilizados para considerar que dos descriptores son colineales es la obtención de un valor de r > 0.6 (^{Maran et al., 2007}), aunque estrictamente el valor depende de la naturaleza de los datos. Los descriptores colineales no deben combinarse en el modelo.

Es deseable que los descriptores empleados para generar un modelo QSAR/QSPR no provengan de propiedades relacionadas. Por tanto, se evalúa la colinealidad de los descriptores antes de construir el modelo, y se excluye a los pares de descriptores identificados con un coeficiente de correlación r > 0.6.

Paso 5. Grupo de entrenamiento y de prueba. Para el ajuste del modelo QSAR/QSPR es necesario dividir el grupo de moléculas en dos subgrupos, el grupo de estimación o de entrenamiento y el grupo de predicción o de prueba. El grupo de entrenamiento de n estructuras permitirá generar el modelo de QSAR/QSPR con los descriptores moleculares y la actividad/propiedad, mientras que el grupo de prueba permitirá establecer si el modelo es suficientemente predictivo. Para este estudio, se seleccionó el 80% de los datos para el grupo de entrenamiento y el 20% restante para el grupo de prueba.

Paso 6. Modelo QSPR. Se seleccionan los descriptores x_j no colineales, con los cuales se genera la regresión lineal múltiple (ecuación 1) mediante el método de cuadrados mínimos para el grupo de entrenamiento.

Donde β_j son los coeficientes de los descriptores y ε es el error asociado a la regresión.

El método de cuadrados mínimos busca que la suma de los cuadrados de las diferencias entre los valores observado y predicho sea mínima. El informe de la regresión lineal múltiple resume los valores de: (1) coeficiente de determinación r², el cual indica la variabilidad de la variable dependiente que el modelo puede explicar; (2) el coeficiente de determinación ajustado adj r², que es la variabilidad de la variable dependiente y que es sensible al número de descriptores; (3) el estadístico F de Fisher, que indica la proximidad relativa del análisis de regresión al valor F de la función de densidad de Fisher correspondiente; y (4) el valor de probabilidad p, el cual es el área bajo la curva de la densidad de Fisher y se interpreta como la probabilidad de que los estadísticos presenten los valores obtenidos, o más extremos (^{Wasserstein, Lazar, Wasserstein, Lazar, & Asa, 2016}). Además, se genera el gráfico de actividad observada vs actividad calculada. El sobreajuste de un modelo de regresión puede aumentar la variabilidad explicada al aumentar r², pero pierde considerablemente la capacidad predictiva ante nuevas observaciones, por lo que el estadístico adj r² es de importancia fundamental.

Paso 7. Valores atípicos. Para obtener un mejor ajuste del modelo, los valores atípicos deben detectarse y eliminarse empleando un método estadístico. La eliminación de valores atípicos para cualquier modelo puede realizarse mediante un análisis de residuales estandarizados o de residuales estudentizados, debiéndose evitar la identificación y selección visual de estos valores, ya que pueden afectar notablemente la significancia del modelo. Una vez descartados los valores atípicos, se calcula un nuevo modelo con las mismas variables iniciales.

Paso 8. Validación. Finalmente, se analiza la capacidad predictiva del modelo generado mediante una validación interna y externa. Se consideran los estadísticos r², r₀ ², r’₀ ², que son los coeficientes de determinación para las variables dependientes observada y calculada; y las pendientes del gráfico de actividad observada contra actividad calculada ajustado al origen y viceversa.

Base de datos

Para este caso de estudio se siguieron los ocho pasos presentados en la metodología, empleando los valores experimentales de logBB de 202 compuestos recopilados por ^{Abraham y colaboradores (Abraham, Ibrahim, Zhao, & Acree, 2006}). Las estructuras químicas de los compuestos se descargaron en formato “.sdf” a partir de los SMILES de la base de datos de PubChem.

Cálculo de descriptores

El programa permitió calcular 278 descriptores para las 202 moléculas. Para la predicción de logBB se emplearon como descriptores del modelo: (1) el logaritmo de la constante de reparto 1-octanol-agua calculado (XLogP) y (2) el área superficial polar topológica (TopoPSA). Análogamente, la ecuación de ^{Clark (Clark, 1999}) emplea el logP calculado (clogP) y el área superficial polar (PSA). La sustitución de PSA por TopoPSA está justificada por su alto coeficiente de correlación (r = 0.99 con n = 34810) (^{Ertl, Rohde, & Selzer, 2000}).

Descriptores colineales

En este caso, los descriptores XLogP y TopoPSA no presentaron una relación colineal (r = 0.10), y son los más significativos en un modelo lineal múltiple.

Modelo QSAR/QSPR

Se generó el modelo QSAR/QSPR con el grupo de entrenamiento, el cual se obtuvo de manera aleatoria. Al escoger los datos de esta manera, es posible que se seleccionen aquellas que disminuyen la correlación y el ajuste en el modelo. Es necesario mencionar que existen algoritmos específicos para la división de datos (^{Montgomery, Peck, & Geoffrey, 2001}), además de la selección empírica.

La instrucción lm(logBB ~ XLogP + TopoPSA, entrenamiento) efectuó la regresión lineal múltiple a partir de los datos de 161 moléculas en el grupo de entrenamiento, el 80% del total. En la Figura 3 se muestra el resumen proveniente del análisis de regresión lineal múltiple, mostrando los valores de los coeficientes (Estimate), sus errores estándar (Std. Error), su significancia (Pr), y los estadísticos de la regresión. También se reporta el valor del residual mínimo y máximo y los cuartiles (1Q, mediana o 2Q y 3Q). Estos últimos son útiles para el análisis de normalidad de los residuales.

Figura 3 Resumen de la regresión lineal múltiple 

Grados de libertad

Es conocido que el cociente de dos varianzas s² ₁ y s² ₂ de dos muestras independientes tiene una distribución F de Fisher (Weimer, 2002), con ν₁ y ν₂ grados de libertad respectivamente. En el análisis de regresión lineal, la variación dentro del grupo de descriptores y la variación entre éste y el error constituyen las varianzas del análisis, y se conocen como raíces de los cuadrados medios.

En nuestro caso, los grados de libertad del análisis de varianza son ν₁ = 2 (j) para la regresión y ν = 158 (n-j-1) para el error; donde n es el número de componentes (compuestos) y j el número de variables independientes (descriptores).

Cuadrado medio residual (CMres)

El valor esperado de CMres es un estimador insesgado de la varianza del modelo, lo que significa que es un estimado adecuado de ésta. Si el valor de CMres es considerablemente mayor al cuadrado medio de la regresión, disminuye la significancia global del modelo.

En nuestro modelo, el valor del cuadrado medio residual fue de CMres = 0.3866, que no influyó de manera notable en la significancia.

Significancia de los descriptores

Se realizó una prueba de hipótesis para la regresión que busca examinar la relación entre la actividad o característica biológica y los descriptores. La hipótesis nula H₀ (ecuación 2) asume que los coeficientes no tienen significancia en el modelo.

El rechazo de H₀ implica que al menos uno de los j-ésimos descriptores contribuye al modelo en forma significativa (^{Mendenhall & Sincich, 1997}a; ^{Montgomery C.D. et al., 2002}).

En una prueba de hipótesis, se calcula un estadístico de prueba que sigue una distribución conocida y se contrasta con el valor crítico de esta distribución para el cual el área bajo la curva en su dominio es el nivel de confianza, generalmente 0.95 o 95%. El resto del área bajo la curva se conoce como área de significancia. Un estadístico de prueba que se encuentre en el área de significancia indica que existe una diferencia estadísticamente significativa en el contraste que se está realizando. El área bajo la curva en el intervalo del estadístico de prueba y el límite superior del dominio es el valor p. Éste valor indica la probabilidad de obtener un resultado tan extremo como el obtenido considerando que el modelo no es significativo (^{Wasserstein, 2016}).

El cociente del coeficiente del descriptor y su error estándar siguen una distribución t de Student. Para los descriptores del modelo lineal, los valores p fueron p_XLogP = 3.35×10^-8 y p_TopoPSA < 2×10^-16, que son menores al área de significancia de 0.05, por lo que son significativos para el modelo. Este hallazgo es de esperarse ya que el descriptor XLogP describe el predominio de la solubilidad en medio lipídico o acuoso, siendo liposoluble cuando es mayor a 0.3 (^{Clark, 1999}). El descriptor TopoPSA es una medida del área polar de la molécula, y un valor bajo favorece el transporte a través de la barrera hematoencefálica (^{Ertl et al., 2000}).

Estadístico F de Fisher

En este caso de estudio, el valor del estadístico F del modelo con 2 y 158 grados de libertad es F = 97.68, y el valor p < 2.2×10^-16. Lo cual indica en la prueba de hipótesis global de regresión que el modelo es significativo.

Coeficiente de determinación

El coeficiente de determinación r² es la proporción de la varianza de la observación explicada por los descriptores (^{Montgomery et al., 2002}). En general, aumenta al incluir un descriptor más al modelo, sin considerar su contribución. El estadístico r² ajustado (adj r²) sólo aumenta cuando el descriptor reduce el valor del CMres. Al igual que el estadístico F de Fisher, r² ajustado evita el sobreajuste del modelo, esto es, la adición de términos que son innecesarios.

Para este modelo lineal, el coeficiente de determinación es r² = 0.5529 y el coeficiente ajustado es adj r² = 0.5472. Los cuales no concuerdan con los valores esperados dada la significancia de los descriptores en el modelo, por lo que se debe a la dispersión de los valores respecto a la línea recta de ajuste (r = 0.7435).

Al graficar con la opción xyplot, se puede observar que la relación de los valores calculados y los valores observados (Figura 4) presentan una tendencia aproximadamente lineal, con múltiples pares ordenados dispersos.

Figura 4 Gráfico de valores calculados contra valores observados. 

Se espera que la relación entre las variables sea aproximadamente lineal, y la dispersión excesiva de los valores de los descriptores x_j afecta el modelo QSAR al disminuir el valor de r² y aumentar la varianza del j-ésimo descriptor.

Valores atípicos

Residuales estandarizados y estudentizados

Los residuales (e_i) son la diferencia de los valores observados y_i con los valores calculados ŷ_i para cada observación.

Los residuales estandarizados son el resultado de dividir cada residual por la raíz del cuadrado medio residual. Los residuales estudentizados no consideran esta aproximación y emplean la desviación exacta del i-ésimo residual, y pueden tener mayor confiabilidad en la detección de valores atípicos (^{Montgomery, Peck, & Vining, 2002}). Se puede considerar el intervalo de -2 a 2 de los residuales distribuidos normalmente como región de confianza (aproximadamente el 95.4%) para su aceptación en el modelo.

Es importante mencionar que existen otros métodos para la detección de valores atípicos, como las distancias de Cook (^{Weisberg & Cook, 1985}) y los intervalos de confianza (^{Mendenhall & Sincich, 1997}b).

Reconstrucción del modelo de QSPR

Al excluir los valores atípicos de los datos iniciales, se vuelve a crear un modelo lineal con las mismas variables, del cual se espera que aumenten los coeficientes de determinación y el estadístico F. Así, una alta determinación y un valor p mínimo justifica el modelo propuesto. Al eliminar los valores atípicos, tomando como base la información de los residuales estudentizados, se obtuvo el siguiente modelo:

logBB=0.1453±0.0131XLogP

-0.0108±0.0006TopoPSA

Con n = 129, CMres = 0.2280, r = 0.8989, r² = 0.8048, adj r² = 0.8017, F = 259.8, p < 2×10^-16. Los coeficientes β_j siguen siendo significativos para el modelo (p_Descriptor < 2×10^-16).

El nuevo modelo (3) presenta una mejora significativa con el inicial y con el modelo propuesto por Clark y colaboradores (^{Clark, 1999}), el cual se representó de la siguiente forma:

logBB=0.152±0.036clogP

+0.0148(±0.001)PSA

Con n = 55, r = 0.887, r² = 0.787, s = 0.354, F = 95.8

Validación del modelo

La validación permite explicar la capacidad predictiva del modelo, y puede ser interna o externa. La primera emplea la información de los integrantes del modelo para establecer la robustez del mismo, mientras que para la segunda se emplea un grupo de compuestos externos de los cuales se conoce su actividad/propiedad pero no fueron incluidos durante la construcción del modelo, es decir, el grupo de prueba.

Validación interna

La validación interna cruzada consiste en el cálculo del estadístico q², que mide en forma aproximada cuánto cabe esperar que el modelo explique la variabilidad de las nuevas observaciones (^{Montgomery et al. 2002}). Para el cálculo del estadístico se puede utilizar el método de “dejar uno fuera” (LOO, por sus siglas en inglés), que calcula el valor de la i-ésima observación a partir de un modelo generado excluyendo esa misma i-ésima observación (ŷ_i,LOO) (^{Dearden & Cronin, 2006}; ^{Kiralj & Ferreira, 2009}; Montgomery et al., 2001) y también se utilizan los métodos “dejar un subconjunto fuera” y mediante remuestreos.

El valor de los estadísticos de la validación interna son q² = 0.7988 y r² _int = 0.8048. Se espera que el valor de q² sea mayor a 0.5 para considerar una buena predicción (^{Veerasamy et al., 2011}).

Validación externa

Cuando no es posible reunir nuevos datos para fines de validación, un procedimiento razonable es dividir los datos disponibles en dos grupos, los datos de estimación o entrenamiento y los datos de predicción o de prueba. El grupo de entrenamiento se usa para formar el modelo lineal, y el grupo de predicción se usa para analizar la capacidad predictiva del modelo (^{Montgomery et al., 2001}). Si en el grupo de prueba hay presencia de valores atípicos que impacten en los estadísticos de validación externa, se debe realizar el procedimiento de generación del modelo para este grupo, así como detección y eliminación de valores atípicos y generación de un nuevo modelo, de manera análoga a como se realizó en el grupo de entrenamiento.

La validación externa puede ser cruzada o no. Para la validación externa cruzada se calcula el estadístico r² _ext para el grupo de prueba generado. Para la validación externa no cruzada se predice la actividad o característica (ŷ_i,pred) de un grupo de moléculas no considerado inicialmente (^{Kiralj & Ferreira, 2009}).