https://doi.org/10.22319/rmcp.v14i1.6182

Artículo

Regresión cuantil para predicción de caracteres complejos en bovinos Suizo Europeo usando marcadores SNP y pedigrí

Jonathan Emanuel Valerio-Hernández a

Paulino Pérez-Rodríguez b*

Agustín Ruíz-Flores a

a Universidad Autónoma Chapingo. Posgrado en Producción Animal. Carretera Federal México-Texcoco Km 38.5, 56227, Texcoco, Estado de México, México.

b Colegio de Postgraduados. Socio Economía Estadística e Informática. Carretera Federal México-Texcoco Km 36.5, 56230, Texcoco, Estado de México.

*Autor para correspondencia: perpdgo@colpos.mx

Resumen:

Los modelos de predicción genómica generalmente suponen que los errores se distribuyen como variables aleatorias normales, independientes e idénticamente distribuidas con media cero e igual varianza. Esto no siempre se cumple, además puede haber fenotipos distantes de la media poblacional, los que usualmente se eliminan al realizar la predicción. La regresión cuantil (QR) afronta aspectos estadísticos como alta dimensionalidad, multicolinealidad y distribución fenotípica diferente de la normal. El objetivo de este trabajo fue comparar QR utilizando información de marcadores y pedigrí con los métodos alternativos tales como mejor predicción lineal insesgada genómica (GBLUP) y mejor predicción lineal insesgada genómica en un solo paso (ssGBLUP) para analizar los pesos al nacimiento (PN), destete (PD) y al año (PA) de bovinos Suizo Europeo y datos simulados con diferente grado de asimetría y proporción de datos atípicos. La capacidad predictiva de los modelos se evaluó mediante validación cruzada. El desempeño predictivo de QR tanto sólo con información de marcadores como con marcadores más pedigrí, con el conjunto de datos reales, fue mejor que las metodologías GBLUP y ssGBLUP para PD y PA. Para PN GBLUP y ssGBLUP fueron mejores, sin embargo, solo se utilizaron los cuantiles 0.25, 0.50 y 0.75, y la distribución de PN no fue asimétrica. 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. Las ventajas de QR fueron más notorias con distribución asimétrica de los fenotipos y con mayor proporción de datos atípicos, como fue el caso del conjunto de datos simulados.

Palabras clave: Regresión cuantil, GBLUP, ssGBLUP.

Recibido: 30/03/2022

Aceptado: 04/08/2022

Introducción

La motivación principal del método de regresión cuantil (QR) reside en que la mayoría de los modelos para evaluación genética asumen normalidad, lo que no siempre se cumple. Otro problema es que en ocasiones registros fenotípicos muy alejados del promedio de la población se consideran como errores de registro o datos atípicos y por consiguiente se eliminan de los análisis, visto desde el punto de vista genómico se está perdiendo información valiosa de marcadores asociados a ciertas regiones de ADN con fuerte influencia en características de interés.

Con el método QR se obtienen resultados robustos y una visión amplia de las variables explicativas sobre las dependientes(1). Los datos generados a partir de experimentos ómicos frecuentemente son complejos y de gran dimensión, por consiguiente existe un desafío estadístico para extraer información biológica relevante del gran volumen de datos(2,3). El uso de un enfoque sólido como QR hace que la inferencia sea menos sesgada y esté menos sujeta a falsos positivos(2). En estudios recientes que utilizan QR, se describen aplicaciones diversas como estudios de asociación genética(4), genética de poblaciones(5), expresión génica(6,7) y selección genómica(8–10).

Uno de los primeros estudios donde se utilizó QR para predecir el mérito genético individual lo presentaron Nascimento et al(11), quienes utilizaron datos simulados encontrando ventajas al usar QR frente a metodologías convencionales. El mismo año(12) se publicaron resultados usando QR para ajustar curvas de crecimiento con datos de cerdos y marcadores moleculares; no solo ajustaron con éxito las curvas de crecimiento sino que identificaron marcadores de importancia asociados con la característica estudiada. Otro trabajo similar del mismo equipo de investigadores fue presentado por Nascimento et al(13), pero con datos de frijol. Recientemente Pérez-Rodríguez et al(10) extendieron el modelo de regresión cuantil para incluir información de pedigrí a través del uso de la matriz de relaciones genéticas aditivas, mejorando aún más la habilidad predictiva de los modelos y al mismo tiempo identificando las proporciones de las varianzas atribuidas a marcadores, relaciones entre individuos y el residual, lo cual permite obtener de manera precisa una partición de la varianza fenotípica.

El objetivo del presente estudio fue estudiar el poder predictivo del modelo de regresión cuantil utilizando datos simulados y datos reales (pesos al nacimiento, destete y al año) de bovinos Suizo Europeo y se consideraron los modelos: 1) QR con información de marcadores moleculares SNP (QRM), 2) QR incluyendo simultáneamente información de marcadores moleculares y genealógica derivada del pedigrí (QRH); 3) GBLUP el cual al igual que QRM solo incluyó información de marcadores moleculares, y 4) evaluación genómica en un solo paso (ssGBLUP) que incluyó información de marcadores y pedigrí.

Material y métodos

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 , 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  donde  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.

Cuadro 1: Número de animales genotipados y fenotipados para el análisis de los pesos al nacer, al destete y al año de una población de ganado Suizo Europeo

QRM=Regresión cuantil usando información de marcadores, QRH=Regresión cuantil usando información de marcadores y pedigrí, GBLUP=Mejor predictor lineal insesgado genómico, ssGBLUP=Evaluación genómica en un solo paso.

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

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

Modelos

Modelo de regresión cuantil con marcadores (QRM)

El modelo para regresión cuantil es:

donde  es el valor del fenotipo del i-ésimo animal;  es un intercepto;  representa la i-ésima fila de la matriz de marcadores,  es el vector de coeficientes de regresión asociados con marcadores y  son variables aleatorias independientes tales que su cuantil  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  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(1):

donde . 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  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  es el valor del fenotipo de la característica analizada (PN, PD o PA) del i-ésimo animal,  es un intercepto;  la i-ésima fila de la matriz de incidencia para los efectos fijos (sexo, edad de la madre, manejo),  los coeficientes de regresión para los efectos fijos,  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),  efectos aleatorios de grupo contemporáneo, el resto de los términos como se describieron anteriormente.

GBLUP

El modelo está dado por:

,

donde  es la i-ésima fila de la matriz que conecta fenotipos con genotipos,  es el vector de efectos aleatorios para animales. Las varianzas genéticas aditiva, de grupo contemporáneo y residual se asumen , , y , respectivamente. La matriz de relaciones genómicas, , se calcula como lo describen Lopez-Cruz et al(14) y Pérez-Rodríguez et al(15); 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,  variables aleatorias normales e independientes con distribución normal con media 0 y varianza .

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, Agg es una submatriz de A para animales genotipados, Ga = βG + α;  y  se obtienen al resolver el sistema de ecuaciones:

El modelo QRH está dado por:

,

donde , 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 .

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 , 80 % de los datos se usaron para entrenamiento del modelo, el 20 % restante para validación. Por ejemplo,  se usa como grupo de validación y el conjunto  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(18). 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(10). 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(19) 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(20) 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 , con  para PN, se supuso que los errores vienen de una distribución normal sesgada () con media 0, varianza  (parámetro de escala ) e índice de asimetría , es decir , con , se consideró  con un valor de 0.35, ,  lo que lleva a diferentes grados de sesgo positivo. Solo se consideraron valores positivos de , ya que el sesgo negativo se obtiene simplemente cambiando el signo de las  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 , 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  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(22) y también se ha sugerido su uso en selección canalizada(23). Una vez generados los datos se ajustó el modelo QRM con  para compararlo con GBLUP. La selección de los cuantiles se realizó de acuerdo a Nascimento et al(11) quienes consideran estas tres posibilidades cuando la distribución de los fenotipos es asimétrica  o bien cuando la distribución es simétrica , 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  verdaderos y estimados, la correlación entre las señales verdaderas  y estimadas  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(24).

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(10). 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(25) 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, . Los modelos GBLUP y ssGBLUP se ajustaron con la librería de R BGLR(26).

Resultados

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  y QRH  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  y QRH  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  y QRH , excepto para PN con el escenario de solo marcadores, donde el DIC más bajo se obtuvo con QRM .

Cuadro 2: Promedios de la correlación Pearson y la desviación estándar (entre paréntesis) entre los valores fenotípicos observados () y valores fenotípicos predichos (), cuadrado medio del error, componentes de varianza asociados al error (, ) y criterio de información de desviación (DIC) para peso al nacimiento

Cor()= correlación entre fenotipos observados y predichos, MSE=cuadrados medios del error,  o =componentes de varianza asociados al error, DIC=criterio de información de desviación.

Cuadro 3: Promedios de la correlación Pearson y la desviación estándar (entre paréntesis) entre los valores fenotípicos observados () y valores fenotípicos predichos (), cuadrado medio del error, componentes de varianza asociados al error (, ) y criterio de información de desviación (DIC) para peso al destete

Cor()= correlación entre fenotipos observados y predichos, MSE=cuadrados medios del error,  o =componentes de varianza asociados al error, DIC=criterio de información de desviación.

Cuadro 4: Promedios de la correlación Pearson y la desviación estándar (entre paréntesis) entre los valores fenotípicos observados () y valores fenotípicos predichos (), cuadrado medio del error, componentes de varianza asociados al error (, ) y criterio de información de desviación (DIC) para peso al año

Cor()= correlación entre fenotipos observados y predichos, MSE=cuadrados medios del error,  o =componentes de varianza asociados al error, DIC=criterio de información de desviación.

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 .

Cuadro 5: Promedios de la correlación Pearson y desviación estándar (entre paréntesis) entre efectos de marcador “verdaderos y estimados, señales “verdaderas” y estimadas, componentes de varianza asociados al error y valores de DIC para datos simulados con diferentes grados de asimetría y proporción de valores atípicos

Cor()= correlación entre efectos de marcador “verdaderos” y estimados, Cor()= correlación entre señales “verdaderas” y estimadas,  o =componentes de varianza asociados al error, DIC=criterio de información de desviación.

Discusión

En este estudio se compararon las metodologías de análisis QR con GBLUP y ssGBLUP. Esta comparación se hizo con fenotipos simulados con diferentes grados de asimetría y proporciones de datos atípicos y datos reales de pesos al nacimiento, al destete y al año.

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(4), 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  y QRH , excepto en el escenario de solo marcadores y PN, donde el DIC más bajo se obtuvo con QRM . 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(27) 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 , 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(10), 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(4), 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(24). Los estimadores de varianza residual y los valores de DIC más bajos se obtuvieron con QR , 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.

Conclusiones e implicaciones

El desempeño predictivo de QR tanto sólo con información de marcadores como con marcadores más pedigrí, con el conjunto de datos reales, fue mejor que las metodologías GBLUP y ssGBLUP para PD y PA. Para PN GBLUP y ssGBLUP fueron mejores; sin embargo, solo se utilizaron los cuantiles 0.25, 0.50 y 0.75, y la distribución de PN no fue asimétrica. 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. Las ventajas de QR fueron más notorias con distribución asimétrica de los fenotipos y con mayor proporción de datos atípicos, como fue el caso del conjunto de datos simulados.

Agradecimientos

Al Consejo Nacional de Ciencia y Tecnología, México, por el apoyo financiero para el primer autor durante sus estudios de doctorado. Los autores también agradecen a la Asociación Mexicana de Criadores de Ganado Suizo de Registro por permitir el uso de sus bases de datos, a los criadores cooperantes por su amable cooperación en este estudio.

Conflictos de interés

Los autores declaran que no existen conflictos de interés.

1. Koenker R, Bassett G. Regression quantiles. Econometrica 1978;46(1):33. https://doi.org/10.2307/1913643.
2. Briollais L, Durrieu G. Application of quantile regression to recent genetic and -omic studies. Hum Genet 2014;133(8):951–966. https://doi.org/10.1007/s00439-014-1440-6.
3. Wang L, Wu Y, Li R. Quantile regression for analyzing heterogeneity in ultra-high dimension. J Am Stat Assoc 2012;107(497):214-222. https://doi.org/10.1080/01621459.2012.656014.
4. Nascimento M, Nascimento ACC, Silva FF, Barili LD, Do Vale NM, Carneiro JE, Cruz CD, Carneiro PCS, Serão NVL. Quantile regression for genome-wide association study of flowering time-related traits in common bean. PLoS One 2018;13(1):1-14. https://doi.org/10.1371/journal.pone.0190303.
5. Fisher E, Schweiger R, Rosset S. Efficient construction of test inversion confidence intervals using quantile regression. J Comput Graph Stat 2016;29:140-148, http://arxiv.org/abs/1612.02300.
6. Logan JAR, Petrill SA, Hart SA, Schatschneider C, Thompson LA, Deater-Deckard K, de Thorne LS, Bartlett C. Heritability across the distribution: An application of quantile regression. Behav Genet 2012;42(2):256–267. https://doi.org/10.1007/s10519-011-9497-7.
7. Vinciotti V, Yu K. M-quantile regression analysis of temporal gene expression data. Stat Appl Genet Mol Biol 2009;8(1). https://doi.org/10.2202/1544-6115.1452.
8. Gianola D, Cecchinato A, Naya H, Schön CC. Prediction of complex traits: Robust alternatives to best linear unbiased prediction. Front Genet 2018;9:195. https://doi.org/10.3389/fgene.2018.00195.
9. Oliveira GF, Nascimento ACC, Nascimento M, Sant’Anna IdeC, Romero JV, Azevedo CF, Bhering LL, Moura ETC. Quantile regression in genomic selection for oligogenic traits in autogamous plants: A simulation study. PLoS One 2021;16(1):1-12. https://doi.org/10.1371/journal.pone.0243666.
10. Pérez-Rodríguez P, Montesinos-López OA, Montesinos-López A, Crossa J. Bayesian regularized quantile regression: A robust alternative for genome-based prediction of skewed data. Crop J 2020;8(5):713-722. https://doi.org/10.1016/j.cj.2020.04.009.
11. Nascimento M, e Silva FF, de Resende MDV, Cruz CD, Nascimento ACC, Viana JMS, Azebedo CF, Barroso LMA. Regularized quantile regression applied to genome-enabled prediction of quantitative traits. Genet Mol Res 2017;16(1). https://doi.org/10.4238/GMR16019538.
12. Barroso LMA, Nascimento M, Nascimento ACC, Silva FF, Serão NVL, Cruz, CD, et al. Regularized quantile regression for SNP marker estimation of pig growth curves, J. Anim Sci Biotechnol 2017;8:59. https://doi.org/10.1186/s40104-017-0187-z.
13. Nascimento AC, Nascimento M, Azevedo C, Silva F, Barili L, Vale N, et al. Quantile regression applied to genome-enabled prediction of traits related to flowering time in the common bean. Agronomy 2019;9(12):1-11. https://doi.org/10.3390/agronomy9120796.
14. López-Cruz M, Crossa J, Bonnett D, Dreisigacker S, Poland J, Jannink JL, Singh RP, Autrique E, de los Campos G. Increased prediction accuracy in wheat breeding trials using a Marker × Environment interaction genomic selection model. G3 Genes Genom Genet 2015;5(4):569–582. https://doi.org/10.1534/g3.114.016097.
15. Pérez-Rodríguez P, Crossa J, Rutkoski J, Poland J, Singh R, Legarra A, et al. Single-step genomic and Pedigree Genotype × Environment interaction models for predicting wheat lines in international environments. Plant Genome 2017;10(2). https://doi.org/10.3835/plantgenome2016.09.0089.
16. Christensen O, Lund M. Genomic relationship matrix when some animals are not genotyped. Genet Sel Evol 2010;42(2):1-8. http://www.gsejournal.org/content/42/1/2.
17. Aguilar I, Misztal I, Johnson DL, Legarra A, Tsuruta S, Lawlor TJ. Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score. J Dairy Sci 2010;93(2):743-752. https://doi.org/10.3168/jds.2009-2730.
18. Crossa J, Pérez P, de los Campos G, Mahuku G, Dreisigacker S, Magorokosho C. Genomic selection and prediction in plant breeding. J Crop Improv 2011;25(3):239-261. https://doi.org/10.1080/15427528.2011.558767.
19. Arnold BC, Beaver RJ. Hidden truncation models. Shankhya. Indian J Stat 2000;62(1):23–35. http://www.jstor.org/stable/25051286. Accessed  Jul 6, 2022.
20. Montesinos-López A, Montesinos-López OA, Villa-Diharce ER, Gianola D, Crossa J. A robust Bayesian genome-based median regression model. Theor Appl Genet 2019;132(5):1587-1606. https://doi.org/10.1007/s00122-019-03303-6.
21. Pérez-Rodríguez P, Villaseñor-Alva JA. On testing the skew normal hypothesis. J Stat Plan Inference 2010;140(11):3148-3159. https://doi.org/10.1016/j.jspi.2010.04.013.
22. Pérez-Rodríguez P, Acosta-Pech R, Pérez-Elizalde S, Cruz CV, Espinosa JS, Crossa J. A Bayesian genomic regression model with skew normal random errors. G3 Genes|Genom|Genet 2018;8(5):1771–1785. https://doi.org/10.1534/g3.117.300406.
23. Domínguez-Viveros J. Parámetros genéticos en la varianza residual de variables de comportamiento en toros de lidia. Arch Zoot 2020;69(267):354–358. https://doi.org/10.21071/az.v69i267.5354.
24. Spiegelhalter DJ, Best NG, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc: Series B Stat Methodol 2002;64(4):583–639. https://doi.org/10.1111/1467-9868.00353.
25. R Core Team. R: A language and environment for statistical computing. R Foundation for statistical computing 2021. Vienna, Austria. https://www.R-project.org/.
26. Pérez P, de los Campos G. Genome-wide regression and prediction with the BGLR statistical package. Genetics 2014;198(2):483-495. https://doi.org/10.1534/genetics.114.164442.
27. Mendes dos Santos P, Nascimento ACC, Nascimento M, Fonseca e Silva F, Azevedo CF, Mota RR, et al. Use of regularized quantile regression to predict the genetic merit of pigs for asymmetric carcass traits. Pesqui Agropecu Bras 2018;53(9):1011–1017. https://doi.org/10.1590/S0100-204X2018000900004.
28. Gianola D. Priors in whole-genome regression: The Bayesian alphabet returns. Genetics 2013;194(3):573-596. https://doi.org/10.1534/genetics.113.151753.
29. de los Campos G, Hickey JM, Pong-Wong R, Daetwyler HD, Calus MPL. Whole-genome regression and prediction methods applied to plant and animal breeding. Genetics 2013;193(2):327-345. https://doi.org/10.1534/genetics.112.143313.
30. Ferragina A, de los Campos G, Vazquez AI, Cecchinato A, Bittante G. Bayesian regression models outperform partial least squares methods for predicting milk components and technological properties using infrared spectral data. J Dairy Sci 2015;98(11):8133-8151. https://doi.org/10.3168/jds.2014-9143.