Genotipos

Se utilizó información de 300 animales (236 hembras, 64 machos) nacidos de 2001 a 2016 en ocho hatos ubicados en el Este, Centro y Oeste de México. Se recolectaron muestras de pelo para su genotipado por la compañía GeneSeek (Lincoln, https://www.neogen.com/, NE, EE. UU.), utilizando el panel GeneSeek® Genomic Profiler Bovine LDv.4, con 30,000 y 50,000 marcadores SNP, 150 animales con cada Chip. La genotipación se realizó en dos ocasiones distintas, inicialmente 150 individuos con el Chip de 30K y posteriormente otro grupo de 150 individuos con el Chip de 50K ya que el Chip de 30K no estuvo disponible en ese momento. Se utilizaron los SNP en común entre los chips de 30K y 50K (12,835 SNP). Se calcularon las proporciones de valores faltantes para cada marcador y para cada individuo. El promedio de valores faltantes por individuo fue 2.09 % con una desviación estándar de 7.50 %. La tasa de llamada promedio (proporción no faltante para cada marcador) fue 97.90 % con una desviación estándar del 4.66 %. Se eliminaron los marcadores con una tasa de llamada inferior al 95 %. Los genotipos se recodificaron como AA= 0, AB= 1 y BB= 2, de donde se obtuvo una matriz con 300 filas (individuos) y 12,835 columnas (marcadores) cuyas celdas toman valores en el conjunto {0,1,2,-}, donde “-” denota un valor faltante. Para los 12,835 marcadores en común de los dos chips, los valores faltantes se imputaron al azar, generando muestras de la distribución Binomial(2,p^) donde p^ es la frecuencia del alelo mayor, calculado a partir de los genotipos de marcadores observados. Se eliminaron los marcadores monomórficos o aquellos con frecuencia de alelos menores (MAF) menor que 0.04. Después del control de calidad, 9,628 de los 12,835 SNP en común entre los dos chips estuvieron disponibles para análisis adicionales.

Fenotipos

La información fenotípica y de pedigrí de la población de ganado Suizo Europeo se obtuvo de la base de datos de la Asociación Mexicana de Criadores de Ganado Suizo de Registro. Los registros de pesos al nacer (PN), al destete (PD) y al año (PA) se utilizaron para el análisis. La edición de los fenotipos fue similar para PN, PD y PA, se descartaron los registros de animales no relacionados genéticamente con aquellos genotipados o con información faltante para hato, edad de la madre y manejo. Los grupos contemporáneos (GC) se definieron al eliminar los animales en GC de 2 individuos o con varianza igual a cero. Para PN los GC se definieron combinando los efectos del hato (8 hatos), año (1998 a 2016) y temporada de nacimiento; las temporadas de nacimiento se definieron considerando el calendario juliano, de 80 a 171 días, primavera; de 172 a 264 días, verano; de 265 a 354 días, otoño; de 355 a 366 días y de 1 a 79 días, invierno. Después de la edición de datos para PN se obtuvieron 330 registros. Para PD y PA los GC se definieron combinando los efectos del hato (6 hatos), año (de 1998 a 2016), temporada de nacimiento (igual que PN) y manejo. En el caso de PD los grupos de manejo se definieron de tres formas: terneros alimentados con leche de su madre; leche de su madre más alimento balanceado; y leche de su madre y nodriza más alimentación balanceada. Para PA los grupos de manejo se definieron de tres maneras: animales en pastoreo; animales en pastoreo con concentrado alimenticio; y animales estabulados con alimentación equilibrada. La edición de datos de PD y PA finalizó con 267 y 232 registros para análisis posteriores. En el Cuadro 1 se presenta un resumen del número de animales genotipados, y fenotipados para PN, PD y PA. La Figura 1 muestra las gráficas de violín para PN, PD y PA, la media muestral está representada por el punto rojo y la mediana muestral por la línea horizontal dentro del recuadro, de la gráfica es claro que las variables respuestas presentan una distribución asimétrica y los círculos con relleno sólido en la misma sugieren la presencia de datos atípicos.

La media muestral está representada por el punto rojo y la mediana muestral por la línea horizontal dentro del recuadro

Figura 1 Gráficas de violín de peso al nacimiento (PN), al destete (PD) y al año (PA) en una población de bovinos Suizo Europeo 

Modelos

Modelo de regresión cuantil con marcadores (QRM)

El modelo para regresión cuantil es:

donde yi es el valor del fenotipo del i-ésimo animal; μ es un intercepto; xit=(xi1,…,xip) representa la i-ésima fila de la matriz de marcadores, β=β1,…,βpt es el vector de coeficientes de regresión asociados con marcadores y wi son variables aleatorias independientes tales que su cuantil θ∈(0,1) es cero. La estimación de los coeficientes de regresión para un cuantil de interés θ fijo se obtiene al resolver el problema de minimización siguiente:

donde ∑j=1pβj es la suma de los valores absolutos de los coeficientes de regresión; λ es el parámetro de penalización que controla la intensidad de la regularización; y ρθ(⋅) es la función definida como¹:

donde ti=yi-μ-xitβ. Después de estimar los parámetros del modelo, los valores de cría estimados mediante marcadores (GEBV) se obtienen mediante la siguiente expresión:

donde β^j(τ) es el efecto del j-ésimo marcador, definido por la relación funcional obtenida para el cuantil de interés.

El modelo QR puede extenderse para incluir otros términos, en particular para las características de crecimiento, se utiliza el siguiente modelo:

donde yi es el valor del fenotipo de la característica analizada (PN, PD o PA) del i-ésimo animal, μ es un intercepto; sit=(si1,...,sif) la i-ésima fila de la matriz de incidencia para los efectos fijos (sexo, edad de la madre, manejo), ς=(ς1,...,ςf)t los coeficientes de regresión para los efectos fijos, cit=(ci1,...,cit) la i-ésima fila de la matriz de incidencia para efectos aleatorios de grupo contemporáneo (54, 43 y 37 para PN, PD y PA), ϱ=(ϱ1,...,ϱt)t efectos aleatorios de grupo contemporáneo, el resto de los términos como se describieron anteriormente.

GBLUP

El modelo está dado por:

donde zit=(zi1,..,zin) es la i-ésima fila de la matriz que conecta fenotipos con genotipos, u=(u1,...,un)t es el vector de efectos aleatorios para animales. Las varianzas genéticas aditiva, de grupo contemporáneo y residual se asumen Varu=Gσu2, Varc=Iσcg2, y Vare=Iσe2, respectivamente. La matriz de relaciones genómicas, G, se calcula como lo describen Lopez-Cruz et al¹⁴ y Pérez-Rodríguez et al¹⁵; brevemente, G = WW’/p, donde W es la matriz de marcadores estandarizada y centrada (cada marcador centrado restando la media de frecuencia de alelo y estandarizado, dividiendo por la desviación estándar de la muestra de la frecuencia de alelos), p es el número total de marcadores, ei variables aleatorias normales e independientes con distribución normal con media 0 y varianza σe2.

Modelo de regresión cuantil en un solo paso (QRH)

Este método se considera una extensión del modelo de cuantiles para una matriz de relaciones construida utilizando matrices de relaciones para animales genotipados y no genotipados y de los cuales se tiene disponible un pedigrí. La matriz que resulta se conoce en la literatura como matriz H^16,17, esta matriz está dada por:

donde, A _gg es una submatriz de A para animales genotipados, G _a = β G + α;β y α se obtienen al resolver el sistema de ecuaciones:

El modelo QRH está dado por:

donde Varu=σH2H, el resto de los términos como se describieron previamente.

Modelo de regresión GBLUP en un solo paso (ssGBLUP)

El modelo ssGBLUP es equivalente al modelo GBLUP descrito anteriormente con la diferencia de que se reemplaza la matriz de relaciones genómicas G por la matriz de relaciones genéticas extendida H, se asume que Varu=HσH2.

Validación cruzada

La capacidad predictiva de los modelos se evaluó mediante validación cruzada, la cual se realizó de la siguiente manera. El conjunto de datos se dividió en cinco grupos del mismo tamaño {S1,S2,…,S5}, 80 % de los datos se usaron para entrenamiento del modelo, el 20 % restante para validación. Por ejemplo, {S2} se usa como grupo de validación y el conjunto {S1,S3,…,S5} para entrenamiento del modelo. Los modelos se ajustaron utilizando el conjunto de entrenamiento, y el modelo ajustado se utilizó para obtener predicciones para el conjunto de validación. Este procedimiento se repitió cinco veces y se obtuvieron predicciones para cada grupo. Las correlaciones entre fenotipos observados y predichos se calcularon y promediaron para los conjuntos de prueba¹⁸. Note que debido a que se trata de valores reales, no se conocen los valores de crías verdaderos, sino solamente se tienen los fenotipos observados, el modelo ajustado proporciona predicciones para valores de cría y las predicciones de otros efectos fijos y ambientales, con lo cual se obtiene una predicción del fenotipo mismo que se contrasta con el valor verdadero del fenotipo.

Simulación

Con la finalidad de evaluar el poder predictivo del modelo QR frente a GBLUP se realizó además una simulación de datos asimétricos con presencia de datos atípicos; la simulación del presente trabajo es análoga a la utilizada por Pérez-Rodríguez et al¹⁰. La idea principal es resaltar que el modelo de regresión por cuantiles funciona de manera adecuada en la presencia de observaciones atípicas, varianzas no homogéneas y variables respuesta con respuestas con distribución asimétrica. En el contexto de selección, por ejemplo, no es poco usual contar con distribuciones asimétricas para los fenotipos debido al proceso mismo, ya que, si se selecciona para alguna característica Y, y si existe además de esto otra característica de interés O, entonces la distribución condicional de Y |O>o¹⁹ es asimétrica. En el contexto de selección genómica es también usual encontrar subconjuntos de observaciones que difieren significativamente del resto y estas observaciones podrían considerarse como atípicas. Montesinos-López et al²⁰ propusieron un modelo con errores Laplace y mostraron que predice de manera adecuada aun en presencia de datos atípicos, el modelo propuesto es un caso especial del modelo de regresión cuantil que se estudia en el presente trabajo. Se consideraron los 9,628 SNP resultantes del control de calidad descrito anteriormente para 300 animales, la simulación de los datos se realizó considerando el modelo lineal:

donde i=1,…, 300, con μ=39 para PN, se supuso que los errores vienen de una distribución normal sesgada (SNc) con media 0, varianza σ2 (parámetro de escala σ) e índice de asimetría γ1, es decir ei~SNC(0,σ,γ1), con σ=1-h2, se consideró h2 con un valor de 0.35, γ1=2πρ34π-11-2ρ2π-3/2, ρ∈0.950, 0.975, 0.999 lo que lleva a diferentes grados de sesgo positivo. Solo se consideraron valores positivos de γ1, ya que el sesgo negativo se obtiene simplemente cambiando el signo de las ei's y por lo tanto las conclusiones obtenidas para el caso de sesgo positivo serán también válidas para el caso negativo^21,22. Se fijaron 50 marcadores con efecto diferente de cero simulándolos de una distribución normal con media 0 y varianza 1-h2/50, el resto de los marcadores se establecieron en 0; las posiciones de los marcadores muestreados se tomaron al azar. Para introducir datos atípicos en los fenotipos, se generó cierta proporción de los residuos de ei~SNC(0,3,γ1) al azar, se consideraron dos proporciones, 5 y 10 %, por lo que se tomaron muestras de una mezcla de dos componentes de distribuciones normales sesgadas. Se generaron seis conjuntos de datos, tres diferentes coeficientes de asimetría 0.950, 0.975, 0.999 con sus dos alternativas de proporción de datos atípicos 5 % y 10 %. La distribución normal asimétrica se ha utilizado en predicción genómica²² y también se ha sugerido su uso en selección canalizada²³. Una vez generados los datos se ajustó el modelo QRM con θ=0.25, 0.50, 0.75 para compararlo con GBLUP. La selección de los cuantiles se realizó de acuerdo a Nascimento et al¹¹ quienes consideran estas tres posibilidades cuando la distribución de los fenotipos es asimétrica θ∈{0.25,0.75} o bien cuando la distribución es simétrica θ 0.50, pues nuestro interés fundamental en este trabajo se centró en la modelación de datos posiblemente asimétricos y con la presencia de valores atípicos. La selección de los parámetros también se realizó por conveniencia computacional ya que el ajuste del modelo se realiza mediante el uso de técnicas de cómputo intensivas basadas en cadenas de Markov Monte Carlo como se menciona en la sección de software y ajuste de los modelos. Para cada análisis se calculó la correlación entre los β's verdaderos y estimados, la correlación entre las señales verdaderas Xβ y estimadas Xβ^ y el componente de varianza asociado con los residuos para cada modelo, el cual es una forma de evaluar la bondad de ajuste de los modelos. También se consideró el Criterio de Información de Desviación (DIC), éste se puede usar para seleccionar modelos candidatos; los modelos con menor DIC se prefieren contra modelos de mayor DIC²⁴.

Software y ajuste de los modelos

Los modelos de regresión por cuantiles se ajustaron usando una estrategia computacional similar a la descrita por Pérez-Rodríguez et al¹⁰. Las adaptaciones de los algoritmos para incluir efectos fijos y aleatorios no presentan gran dificultad computacional. Los códigos para el ajuste de los modelos se desarrollaron en los lenguajes de programación R²⁵ y C. Los códigos para el ajuste de los modelos fueron organizados de tal forma que puedan ser ejecutados fácilmente desde el software estadístico R y están disponibles solicitándolos al primer autor del presente estudio. En todos los casos se seleccionaron tres cuantiles, θ={0.25, 0.50, 0.75}. Los modelos GBLUP y ssGBLUP se ajustaron con la librería de R BGLR²⁶.

Datos reales

Los Cuadros 2, 3 y 4 muestran los resultados del experimento realizado con datos de PN, PD y PA de una población de bovinos Suizo Europeo, evaluados bajo dos escenarios 1) solo con información de marcadores, y 2) información de marcadores y pedigrí. De manera general las correlaciones entre valores observados y predichos más altas se obtuvieron con QR, excepto para PN donde las correlaciones de GBLUP y ssGBLUP fueron más altas que las obtenidas con QRM y QRH, sin embargo, las correlaciones de QRM θ=0.75 y QRH θ=0.75 fueron cercanas a las obtenidas con GBLUP y ssGBLUP (0.7902 vs 0.7924), (0.6981 vs 0.7055), respectivamente. Los valores de MSE más bajos se obtuvieron con QRM θ=0.75 y QRH θ=0.75 en la característica de PD, mientras que en las características PN y PA los valores más bajos se obtuvieron con GBLUP y ssGBLUP. Los componentes de varianza asociados con el error obtenidos con QRM y QRH fueron menores que los obtenidos con GBLUP y ssGBLUP. Los valores de DIC más bajos de manera general se obtuvieron con QRM θ=0.75 y QRH θ=0.75, excepto para PN con el escenario de solo marcadores, donde el DIC más bajo se obtuvo con QRM θ=0.25.

Datos simulados

Los resultados del ejercicio de simulación donde se compara QR con GBLUP bajo diferentes grados de asimetría y proporciones de datos atípicos se muestran en el Cuadro 5. En la columna 2 se registran las correlaciones entre los efectos de marcador “verdaderos” y los efectos de marcador estimados, las correlaciones obtenidas con QR fueron más altas que las obtenidas con GBLUP. La columna 3 muestra las correlaciones entre las “señales verdaderas” y las estimadas, las correlaciones más altas se obtuvieron con QR. La columna 4 registra la estimación de los componentes de varianza asociados al error y la columna 5 los valores de DIC, los valores más bajos en ambas columnas se obtuvieron con QR θ=0.75.

Datos reales

Las correlaciones de fenotipos observados y predichos obtenidos de la validación cruzada con datos reales fueron más altas al usar QRM y QRH en las características de PD y PA. Para PN las correlaciones más altas se obtuvieron con GBLUP y ssGBLUP; sin embargo, en este estudio solo se probaron tres cuantiles 0.25, 0.50 y 0.75, hay evidencia en otros estudios donde QR es mejor que GBLUP como en el caso del trabajo de Nascimento et al⁴, quienes compararon QR con modelos como BLASSO, BayesB y RR-BLUP. Estos autores encontraron una ganancia en la capacidad predictiva de QR frente a RR-BLUP de 15.15 %, cabe señalar que matemáticamente RR-es equivalente a GBLUP, además de que los conjuntos de datos usados en este experimento presentaron asimetría.

Los valores del cuadrado medio del error (MSE) miden el promedio de los errores al cuadrado, es decir, la diferencia entre el estimador y lo que se estima, por lo que se prefieren valores bajos; los promedios de MSE de QRM y QRH fueron menores que los obtenidos con GBLUP y ssGBLUP únicamente para PD. El estimador de la varianza residual es un indicativo de que tan bien o mal se ajusta el modelo a los datos observados, se prefieren valores bajos; los componentes de varianza del error menores se obtuvieron con QRM y QRH para las tres características analizadas. Finalmente, el valor de DIC se usa para seleccionar modelos candidatos y al igual que MSE y los componentes de varianza del error se prefieren valores bajos. Los valores de DIC más bajos se obtuvieron con QRM θ=0.75 y QRH θ=0.75, excepto en el escenario de solo marcadores y PN, donde el DIC más bajo se obtuvo con QRM θ=0.25. El cuadrado medio del error, la varianza residual y el DIC son valores que ayuda a escoger el modelo de mejor ajuste. Al examinar estos valores de manera conjunta se observa que QRM y QRH son mejores en algunos de ellos mientras que en otros no, es decir, QR presenta un desempeño igual o superior a GBLUP y ssGBLUP; aunque debe resaltarse que solo se probaron tres cuantiles y que QR presenta ventajas al usarse en datos asimétricos y datos atípicos, para este caso solo había datos atípicos y las distribuciones no fueron asimétricas. Mendes et al²⁷ compararon QR con el método bayesiano del LASSO (BLASSO), estos autores reportaron un aumento de 6.7 y 20.0 % en la precisión y consideraron los cuantiles 0.15 y 0.45 en la evaluación del rendimiento de la canal y el grosor del tocino, respectivamente, sin embargo, las características evaluadas en su estudio fueron asimétricas.

En el análisis de datos reales, una limitación del presente estudio es el tamaño de muestra, lo cual puede impactar la variabilidad de los parámetros estimados con los modelos y en consecuencia la variabilidad de las predicciones, sin embargo, todos los modelos se ajustaron utilizando la misma información y por tanto la comparación de la capacidad predictiva de los modelos se considera razonable, lo ideal sería contar con tamaños de muestra grandes, pero por cuestión de limitaciones económicas esto no siempre es posible. Por otro lado, actualmente es muy común utilizar modelos de predicción en el que el número de registros fenotípicos es más pequeño que el número de predictores (SNPS), es decir n≪p, aun bajo este contexto numerosos estudios han mostrado que los métodos bayesianos proveen herramientas sofisticadas que permiten derivar predicciones razonables siempre y cuando los parámetros de regularización sean seleccionados de manera adecuada por ejemplo utilizando métodos de validación cruzada^28-30.

Datos simulados

En el experimento de datos simulados las correlaciones entre efectos de marcador “verdadero” y efectos estimados, así como las correlaciones de señales “verdaderas” y estimadas fueron más altas cuando se usó QR comparado con GBLUP. Estos resultados son similares a los obtenidos por otros investigadores¹⁰, quienes simularon datos con tres diferentes coeficientes de asimetría 0.75, 0.95, 0.999 con 5 % y 10 % de datos atípicos y encontraron que las correlaciones obtenidas con QR fueron más altas que las obtenidas con la regresión de cresta bayesiana (BRR), además, este patrón fue más evidente con una mayor asimetría y proporción de datos atípicos. En este estudio se realizaron simulaciones con coeficientes de asimetría de 0.950, 0.975, 0.999 y los cuantiles con los que se obtuvieron correlaciones más altas variaron entre 0.25 y 0.50; la ventaja de QR es que se puede probar diferentes cuantiles obteniendo mejores resultados dependiendo del cuantil utilizado, esta ventaja en la capacidad de predicción de efectos de marcadores y señales ha sido aprovechada por otros investigadores⁴, quienes no encontraron ninguna asociación de rasgo usando el modelo tradicional GWAS de SNP único, pero, cuando usaron QR con cuantiles extremos como 0.1 el modelo fue capaz de encontrar hasta 7 SNP asociados con las características estudiadas.

Los coeficientes de varianza del error indican qué tan bien se ajusta el modelo propuesto a los datos estudiados, entre más pequeño sea este valor, el ajuste es mejor, el DIC es otro valor que se usa para comparar modelos candidatos. Los modelos con un DIC más pequeño se prefieren a los modelos con un DIC mayor²⁴. Los estimadores de varianza residual y los valores de DIC más bajos se obtuvieron con QR θ=0.75, quizá esto se deba a que en la simulación se usaron coeficientes de asimetría altos 0.950, 0.975, 0.999, por consiguiente, se espera que un cuantil que se ajuste mejor sea el más alto, en este caso 0.75.

QR tuvo un desempeño igual o mejor que GBLUP y ssGBLUP para predecir características de crecimiento PN, PD y PA, las ventajas de este método son más notorias cuando los datos están más sesgados y presentan mayor proporción de datos atípicos como en el caso del experimento de simulación.