Análisis de Datos Fenotípicos y Genómicos de Arroz
Juanita Torres
2025-07-26
Introducción
Este documento tiene como objetivo aplicar un flujo completo de análisis a un conjunto de líneas de arroz, combinando datos fenotípicos y genómicos. Se implementan métodos como análisis descriptivo, clustering no supervisado (K-means), PCA y clasificación supervisada. Finalmente, se evalúa el desempeño de los modelos.
1. Carga, limpieza y resumen de los datos
library(tidyverse)
library(glue)
library(caret)
library(FactoMineR)
library(factoextra)
library(rpart)
library(rpart.plot)
library(DT)
library(corrplot)
library(knitr)
# Lectura de datos
feno <- read.csv("quantitative_traits.csv")
geno <- read.csv("12k_ld_imputed.csv")
datos <- merge(feno, geno, by = "ID") %>% distinct() %>% drop_na()
feno <- datos[, 1:13] # ID + 12 rasgos
n_lineas <- nrow(feno) # número de líneas evaluadas
n_variables_feno <- ncol(feno) - 1 # número de rasgos fenotípicos
n_variables_gen <- ncol(geno) - 1 # SNPs
variables <- names(feno)[-1] # nombres sin el ID
# Mostrar resultados
knitr::asis_output(
glue("
### Resumen general de los datos
- Número de líneas evaluadas: {n_lineas}
- Número de variables fenotípicas: {n_variables_feno}
Variables fenotípicas disponibles:
{paste0(' - ', variables, collapse = '\n')}
- Número de variables genéticas (SNPs): {n_variables_gen}
")
)Resumen general de los datos
- Número de líneas evaluadas: 1865
- Número de variables fenotípicas: 12
Variables fenotípicas disponibles:- Culm.diameter
- Culm.length
- Culm.number
- Grain.length
- Grain.width
- Grain.weight
- Heading.date
- Ligule.length
- Leaf.length
- Leaf.width
- Panicle.length
- Seedling.height
- Número de variables genéticas (SNPs): 12486
2. Análisis exploratorio
2.1 Tabla descriptiva por línea
Se generó una tabla con medidas de tendencia central y rangos para cada rasgo
tabla <- feno %>%
group_by(ID) %>%
summarise(across(where(is.numeric), list(media = mean, min = min, max = max), .names = "{.col}_{.fn}"))
datatable(tabla, options = list(scrollX = TRUE))# Guardar archivo CSV para descarga
write.csv(tabla, "tabla_resumen.csv", row.names = FALSE)2.2 Histogramas
Los histogramas permitieron observar la distribución de las variables
fenotípicas. Algunas mostraron asimetrías y rangos amplios,
especialmente en variables como Heading.date y
Grain.weight
feno %>% select(-ID) %>%
pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
ggplot(aes(x = valor)) +
geom_histogram(bins = 30, fill = "#69b3a2", color = "white") +
facet_wrap(~variable, scales = "free") +
theme_minimal()
2.3 Matriz de correlación
Se calcularon y visualizaron correlaciones entre las variables fenotípicas. Se encontraron asociaciones
mat_cor <- cor(feno %>% select(-ID))
corrplot(mat_cor, method = "color", type = "upper", tl.col = "black", title = "Correlación entre rasgos", mar=c(0,0,1,0))
2.4 PCA
El análisis de componentes principales permitió visualizar la variabilidad entre líneas. Las dos primeras componentes explicaron gran parte de la varianza fenotípica
pca <- PCA(feno %>% select(-ID), graph = FALSE)
fviz_pca_ind(pca, geom = "point", title = "PCA - Individuos")
fviz_pca_var(pca, col.var = "contrib", gradient.cols = c("blue", "orange", "red"))
3. Agrupamiento no supervisado (K-means)
El método K-means permitió identificar 3 grupos distintos de líneas fenotípicas basados en sus coordenadas PCA. Estos grupos pueden reflejar diferentes perfiles agronómicos
datos_escalados <- scale(feno %>% select(-ID))
pca_coord <- as.data.frame(PCA(datos_escalados, graph = FALSE)$ind$coord[, 1:2])
set.seed(123)
kmeans_model <- kmeans(pca_coord, centers = 3, nstart = 25)
fviz_pca_ind(PCA(datos_escalados, graph = FALSE),
geom.ind = "point",
col.ind = as.factor(kmeans_model$cluster),
palette = "Set2",
addEllipses = TRUE,
legend.title = "Cluster")
4. Modelos supervisados de clasificación
4.1 Clasificación usando datos fenotípicos
4.1.1 Clasificación según Grain.weight (PCA de datos fenotípicos)
Se clasificaron líneas con base en si su peso de grano (‘Grain.weight’)
4.1.1.1 Árbol de decisión
pca_coords <- as.data.frame(pca$ind$coord[, 1:2])
pca_coords$grupo <- ifelse(feno$Grain.weight > median(feno$Grain.weight), "Alto", "Bajo") %>%
factor(levels = c("Bajo", "Alto"))
set.seed(123)
modelo1 <- train(grupo ~ ., data = pca_coords, method = "rpart")
pred1 <- predict(modelo1, pca_coords)
conf1 <- confusionMatrix(pred1, pca_coords$grupo)
data.frame(
Modelo = "Árbol - PCA fenotípico",
Accuracy = round(conf1$overall["Accuracy"], 3),
Sensibilidad = round(conf1$byClass["Sensitivity"], 3),
Especificidad = round(conf1$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Árbol - PCA fenotípico | 0.802 | 0.805 | 0.799 |
rpart.plot(modelo1$finalModel, main = "Árbol de clasificación - Grain.weight")
4.1.1.2 Regresión logística
# Modelo logístico con interacción entre las componentes PCA
modelo_inter <- glm(grupo ~ Dim.1 * Dim.2, data = pca_coords, family = "binomial")
# Predicción
pred_probs <- predict(modelo_inter, type = "response")
pred_class <- ifelse(pred_probs > 0.5, "Alto", "Bajo") %>% factor(levels = c("Bajo", "Alto"))
# Evaluación del modelo
conf_inter <- confusionMatrix(pred_class, pca_coords$grupo)
data.frame(
Modelo = "Logística - PCA fenotípico Grain.weight",
Accuracy = round(conf_inter$overall["Accuracy"], 3),
Sensibilidad = round(conf_inter$byClass["Sensitivity"], 3),
Especificidad = round(conf_inter$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Logística - PCA fenotípico Grain.weight | 0.814 | 0.859 | 0.755 |
4.1.2 Clasificación según Culm.length (PCA de datos fenotípicos)
Se clasificaron líneas en alto/bajo según Culm.length
usando árboles de decisión
4.1.2.1 Árbol de decisión
pca_coords2 <- as.data.frame(pca$ind$coord[, 1:2])
pca_coords2$grupo <- ifelse(feno$Culm.length > median(feno$Culm.length), "Alto", "Bajo")
pca_coords2$grupo <- as.factor(pca_coords2$grupo)
set.seed(123)
modelo2 <- train(grupo ~ ., data = pca_coords2, method = "rpart")
pred2 <- predict(modelo2, pca_coords2)
conf2 <- confusionMatrix(pred2, pca_coords2$grupo)
data.frame(
Modelo = "Árbol - PCA Culm.length",
Accuracy = round(conf2$overall["Accuracy"], 3),
Sensibilidad = round(conf2$byClass["Sensitivity"], 3),
Especificidad = round(conf2$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Árbol - PCA Culm.length | 0.827 | 0.801 | 0.853 |
rpart.plot(modelo2$finalModel, main = "Árbol de clasificación - Culm.length")
4.1.2.2 Regresión logística
# Modelo logístico para Culm.length
modelo_inter2 <- glm(grupo ~ Dim.1 * Dim.2, data = pca_coords2, family = "binomial")
# Predicción
pred_probs2 <- predict(modelo_inter2, type = "response")
pred_class2 <- ifelse(pred_probs2 > 0.5, "Alto", "Bajo") %>% factor(levels = c("Bajo", "Alto"))
# Evaluación del modelo
conf_inter2 <- confusionMatrix(pred_class2, pca_coords2$grupo)
data.frame(
Modelo = "Logística - PCA Culm.length",
Accuracy = round(conf_inter2$overall["Accuracy"], 3),
Sensibilidad = round(conf_inter2$byClass["Sensitivity"], 3),
Especificidad = round(conf_inter2$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Logística - PCA Culm.length | 0.169 | 0.164 | 0.175 |
4.1.3 Clasificación con Random Forest
Se entrenó un modelo de Random Forest utilizando todas las variables fenotípicas disponibles, excluyendo el ID y la variable objetivo (Grain.weight). El objetivo fue clasificar las líneas como de peso de grano “Alto” o “Bajo”, según la mediana. A continuación, se muestran los resultados del modelo, incluyendo la matriz de confusión y la importancia de las variables.
library(randomForest)
datos_rf <- feno %>%
mutate(grupo = ifelse(Grain.weight > median(Grain.weight), "Alto", "Bajo")) %>%
select(-ID, -Grain.weight)
datos_rf$grupo <- as.factor(datos_rf$grupo)
set.seed(123)
split <- sample(2, nrow(datos_rf), replace = TRUE, prob = c(0.7, 0.3))
train_rf <- datos_rf[split == 1, ]
test_rf <- datos_rf[split == 2, ]
# Entrenar Random Forest
modelo_rf <- randomForest(grupo ~ ., data = train_rf, importance = TRUE)
pred_rf <- predict(modelo_rf, test_rf)
conf_rf_feno <- confusionMatrix(pred_rf, test_rf$grupo)
data.frame(
Modelo = "Random Forest - Fenotípico",
Accuracy = round(conf_rf_feno$overall["Accuracy"], 3),
Sensibilidad = round(conf_rf_feno$byClass["Sensitivity"], 3),
Especificidad = round(conf_rf_feno$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Random Forest - Fenotípico | 0.758 | 0.663 | 0.835 |
plot(modelo_rf, main = "Error del modelo - Random Forest fenotípico")
varImpPlot(modelo_rf, main = "Importancia de variables fenotípicas")
5. Clasificación usando datos genotípicos
El PCA con datos de SNPs (12k) permitió visualizar relaciones
genéticas entre líneas. Se identificaron agrupamientos genéticos y se
usaron los datos para entrenar un modelo Random Forest, clasificando
líneas por Grain.weight
5.1 PCA con datos genotípicos
Se realizó un análisis de componentes principales sobre los SNPs para explorar la estructura genética de las líneas.
pca_gen <- PCA(geno[,-1], graph = FALSE)
fviz_pca_ind(pca_gen, geom = "point", title = "PCA con datos genotípicos")
snp_data <- geno %>%
select(-ID) %>%
select(where(is.numeric)) %>%
select(where(~ var(.) > 0))
grupo <- factor(ifelse(datos$Grain.weight > median(datos$Grain.weight, na.rm = TRUE), "Alto", "Bajo"))
ids_validos <- geno$ID %in% datos$ID
snp_data <- snp_data[ids_validos, ]
grupo <- grupo[match(geno$ID[ids_validos], datos$ID)]
datos_rf_solo_snps <- data.frame(ID = geno$ID[ids_validos], grupo, snp_data)
# Dividir y entrenar modelo
set.seed(123)
train_index <- createDataPartition(datos_rf_solo_snps$grupo, p = 0.7, list = FALSE)
train_rf <- datos_rf_solo_snps[train_index, ]
test_rf <- datos_rf_solo_snps[-train_index, ]
modelo_rf <- randomForest(grupo ~ ., data = train_rf %>% select(-ID), importance = TRUE)
pred_rf <- predict(modelo_rf, newdata = test_rf %>% select(-ID, -grupo))
conf_rf_solo_snps <- confusionMatrix(pred_rf, test_rf$grupo)
data.frame(
Modelo = "Random Forest - Solo SNPs",
Accuracy = round(conf_rf_solo_snps$overall["Accuracy"], 3),
Sensibilidad = round(conf_rf_solo_snps$byClass["Sensitivity"], 3),
Especificidad = round(conf_rf_solo_snps$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Random Forest - Solo SNPs | 0.673 | 0.479 | 0.818 |
plot(modelo_rf, main = "Error del modelo - Random Forest (solo SNPs)")
varImpPlot(modelo_rf, n.var = 15, main = "Top 15 SNPs más importantes")
5.2 Clasificación según Grain.weight (PCA con datos genotípicos)
Se utilizaron las dos primeras componentes principales obtenidas del PCA genotípico como predictores para clasificar las líneas según su peso de grano (Grain.weight), usando modelos de árbol de decisión y regresión logística.
ids_geno <- geno$ID
sub_feno <- datos %>%
filter(ID %in% ids_geno, !is.na(Grain.weight))
sub_feno <- sub_feno[match(ids_geno, sub_feno$ID), ]
sub_feno <- sub_feno[!is.na(sub_feno$Grain.weight), ]
pca_coords_raw <- as.data.frame(pca_gen$ind$coord[, 1:2])
pca_coords_raw <- pca_coords_raw[match(sub_feno$ID, ids_geno), ]
colnames(pca_coords_raw) <- c("Dim.1", "Dim.2")
stopifnot(nrow(sub_feno) == nrow(pca_coords_raw))
pca_coords_gen <- pca_coords_raw
pca_coords_gen$grupo <- factor(ifelse(sub_feno$Grain.weight > median(sub_feno$Grain.weight), "Alto", "Bajo"))
set.seed(123)
modelo_gen_tree <- train(grupo ~ ., data = pca_coords_gen, method = "rpart")
pred_gen_tree <- predict(modelo_gen_tree, pca_coords_gen)
conf_gen_tree <- confusionMatrix(pred_gen_tree, pca_coords_gen$grupo)
data.frame(
Modelo = "Árbol - PCA genotípico",
Accuracy = round(conf_gen_tree$overall["Accuracy"], 3),
Sensibilidad = round(conf_gen_tree$byClass["Sensitivity"], 3),
Especificidad = round(conf_gen_tree$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Árbol - PCA genotípico | 0.617 | 0.335 | 0.829 |
rpart.plot(modelo_gen_tree$finalModel, main = "Árbol de clasificación - PCA genotípico")
modelo_gen_log <- glm(grupo ~ Dim.1 * Dim.2, data = pca_coords_gen, family = "binomial")
pred_probs_gen <- predict(modelo_gen_log, type = "response")
pred_class_gen <- ifelse(pred_probs_gen > 0.5, "Alto", "Bajo") %>% factor(levels = c("Bajo", "Alto"))
conf_gen_log <- confusionMatrix(pred_class_gen, pca_coords_gen$grupo)
data.frame(
Modelo = "Logística - PCA genotípico",
Accuracy = round(conf_gen_log$overall["Accuracy"], 3),
Sensibilidad = round(conf_gen_log$byClass["Sensitivity"], 3),
Especificidad = round(conf_gen_log$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Logística - PCA genotípico | 0.396 | 0.649 | 0.206 |
5.3 Clasificación con Random Forest (SNPs completos)
En este modelo se utilizaron directamente todos los SNPs (tras eliminar aquellos con varianza cero), para predecir el peso de grano (Grain.weight) en categorías de Alto y Bajo.
datos_rf_geno <- datos %>%
mutate(grupo = factor(ifelse(Grain.weight > median(Grain.weight), "Alto", "Bajo"))) %>%
select(-Grain.weight)
datos_num_geno <- datos_rf_geno %>% select(where(is.numeric)) %>% select(where(~ var(.) > 0))
datos_rf_clean_geno <- bind_cols(datos_rf_geno %>% select(ID, grupo), datos_num_geno)
set.seed(123)
train_index_geno <- createDataPartition(datos_rf_clean_geno$grupo, p = 0.7, list = FALSE)
train_geno <- datos_rf_clean_geno[train_index_geno, ]
test_geno <- datos_rf_clean_geno[-train_index_geno, ]
modelo_rf_geno <- randomForest(grupo ~ ., data = train_geno %>% select(-ID), importance = TRUE)
pred_rf_geno <- predict(modelo_rf_geno, newdata = test_geno %>% select(-ID, -grupo))
conf_rf_geno <- confusionMatrix(pred_rf_geno, test_geno$grupo)
data.frame(
Modelo = "Random Forest - Genotípico",
Accuracy = round(conf_rf_geno$overall["Accuracy"], 3),
Sensibilidad = round(conf_rf_geno$byClass["Sensitivity"], 3),
Especificidad = round(conf_rf_geno$byClass["Specificity"], 3)
) %>% knitr::kable()| Modelo | Accuracy | Sensibilidad | Especificidad | |
|---|---|---|---|---|
| Accuracy | Random Forest - Genotípico | 0.662 | 0.392 | 0.865 |
varImpPlot(modelo_rf, n.var = 15, main = "Top 15 SNPs más importantes - Random Forest")
6. Comparación de modelos
6.1 Clasificación por variable objetivo
Aquí se comparan modelos para dos variables objetivo diferentes (Grain.weight y Culm.length) usando árboles de decisión.
comparacion_pca <- data.frame(
Modelo = c("Árbol - PCA Grain.weight","Árbol - PCA Culm.length",
"Logística - PCA Grain.weight","Logística - PCA Culm.length"),
Accuracy = c(conf1$overall["Accuracy"], conf_inter$overall["Accuracy"],
conf2$overall["Accuracy"], conf_inter2$overall["Accuracy"]),
Sensibilidad = c(conf1$byClass["Sensitivity"], conf_inter$byClass["Sensitivity"],
conf2$byClass["Sensitivity"], conf_inter2$byClass["Sensitivity"]),
Especificidad = c(conf1$byClass["Specificity"], conf_inter$byClass["Specificity"],
conf2$byClass["Specificity"], conf_inter2$byClass["Specificity"])
)
knitr::kable(comparacion_pca, digits = 3)| Modelo | Accuracy | Sensibilidad | Especificidad |
|---|---|---|---|
| Árbol - PCA Grain.weight | 0.802 | 0.805 | 0.799 |
| Árbol - PCA Culm.length | 0.814 | 0.859 | 0.755 |
| Logística - PCA Grain.weight | 0.827 | 0.801 | 0.853 |
| Logística - PCA Culm.length | 0.169 | 0.164 | 0.175 |
6.2 Comparación de clasificación usando PCA
Se comparan árboles de decisión y regresiones logísticas entrenados con las dos primeras componentes principales de los datos fenotípicos y genotípicos.
comparacion_pca <- data.frame(
Modelo = c("Árbol - PCA fenotípico","Árbol - PCA genotípico",
"Logística - PCA fenotípico","Logística - PCA genotípico"),
Accuracy = c(conf1$overall["Accuracy"], conf_inter$overall["Accuracy"],
conf_gen_tree$overall["Accuracy"], conf_gen_log$overall["Accuracy"]),
Sensibilidad = c(conf1$byClass["Sensitivity"], conf_inter$byClass["Sensitivity"],
conf_gen_tree$byClass["Sensitivity"], conf_gen_log$byClass["Sensitivity"]),
Especificidad = c(conf1$byClass["Specificity"], conf_inter$byClass["Specificity"],
conf_gen_tree$byClass["Specificity"], conf_gen_log$byClass["Specificity"])
)
knitr::kable(comparacion_pca, digits = 3)| Modelo | Accuracy | Sensibilidad | Especificidad |
|---|---|---|---|
| Árbol - PCA fenotípico | 0.802 | 0.805 | 0.799 |
| Árbol - PCA genotípico | 0.814 | 0.859 | 0.755 |
| Logística - PCA fenotípico | 0.617 | 0.335 | 0.829 |
| Logística - PCA genotípico | 0.396 | 0.649 | 0.206 |
6.3 Comparación de Random Forest (todos los datos)
Aquí comparamos los modelos Random Forest entrenados con: Todas las variables fenotípicas y todos los SNPs (genotípicos)
comparacion_rf <- data.frame(
Modelo = c("Random Forest - Fenotípico", "Random Forest - Solo SNPs", "Random Forest - SNPs completos"),
Accuracy = c(conf_rf_feno$overall["Accuracy"], conf_rf_solo_snps$overall["Accuracy"], conf_rf_geno$overall["Accuracy"]),
Sensibilidad = c(conf_rf_feno$byClass["Sensitivity"], conf_rf_solo_snps$byClass["Sensitivity"], conf_rf_geno$byClass["Sensitivity"]),
Especificidad = c(conf_rf_feno$byClass["Specificity"], conf_rf_solo_snps$byClass["Specificity"], conf_rf_geno$byClass["Specificity"])
)
knitr::kable(comparacion_rf, digits = 3)| Modelo | Accuracy | Sensibilidad | Especificidad |
|---|---|---|---|
| Random Forest - Fenotípico | 0.758 | 0.663 | 0.835 |
| Random Forest - Solo SNPs | 0.673 | 0.479 | 0.818 |
| Random Forest - SNPs completos | 0.662 | 0.392 | 0.865 |
7. Conclusiones
Desempeño con PCA
- PCA fenotípico ayudó a identificar patrones entre líneas.
- El modelo logístico con PCA fenotípico fue el mejor de su tipo (Accuracy: 0.814).
- Árbol de decisión con PCA para Culm.length tuvo el mejor desempeño general (Accuracy: 0.827).
Random Forest
- RF con datos fenotípicos fue el mejor entre todos los Random Forest (Accuracy: 0.758).
- RF con SNPs (solo o todos) tuvo menor desempeño (Accuracy: 0.662 - 0.673).
- Variables importantes fueron identificadas automáticamente por el modelo.
Comparación Genotípico vs. Fenotípico
- Modelos con datos fenotípicos superaron a los genotípicos en precisión.
- El PCA genotípico no logró separar claramente los grupos de peso de grano.
- RF con SNPs completos fue el mejor entre los genotípicos, aunque con menor desempeño.
Aplicabilidad en mejoramiento
- Los modelos permiten identificar líneas superiores y variables clave.
- Facilitan la selección indirecta usando características morfológicas o SNPs.
- Son herramientas prometedoras para acelerar procesos de selección en programas de mejoramiento.