#-------------------------------------------------------------------------------------------------
# OBJETIVO DEL TRABAJO
#-------------------------------------------------------------------------------------------------
#Este proyecto aborda el análisis de diversidad genética en un panel de 275 accesiones de sandía (Citrullus spp.), utilizando datos de genotipado de alta densidad (5K SNP), medidas fenotípicas de firmeza (Pheno), y datos de origen geográfico (Country). El objetivo es caracterizar la variabilidad genética y relacionarla con diferencias fenotípicas, tanto en el conjunto completo de especies como en el subconjunto cultivado C. lanatus.
#-------------------------------------------------------------------------------------------------
# 1) CARGA DE DATOS
#-------------------------------------------------------------------------------------------------
### 1.1) Establecer directorio de trabajo y cargar paquetes necesarios
setwd("/home/IACA")
# Paquetes para manejo de datos y gráficos
library(dplyr) # Manipulación de datos
library(ggplot2) # Gráficos
# Paquetes para tablas y resúmenes en PDF
library(gridExtra) # Presentación de tablas en PDF (grid.table)
# Paquetes para dendrogramas y representaciones circulares
library(dendextend) # Personalización de dendrogramas (colores y tamaño de etiquetas)
library(circlize) # Dendrogramas circulares (circlize_dendrogram, circos.par, circos.clear)
# Paquetes para heatmaps y paletas de color
library(RColorBrewer) # Paletas de color predefinidas (brewer.pal para el gradiente del heatmap)
library(pheatmap) # Heatmaps de matrices de distancias con dendrogramas y anotaciones
### 1.2) Leer los dos archivos CSV
#download.file("https://cursoronline.usal.es/WM_GBK_5k_Rinput.csv","WM_GBK_5k_Rinput.csv")
geno <- read.csv(file="WM_GBK_5k_Rinput.csv",sep=";", stringsAsFactors = FALSE)
#download.file("https://cursoronline.usal.es/WM_5k_metadata.csv","WM_5k_metadata.csv")
meta <- read.csv(file="WM_5k_metadata.csv",sep=";", stringsAsFactors = FALSE)
### 1.3) Generar tibbles por facilidad de visualización de datos y manejo.
geno_tibble <- as_tibble(geno)
meta_tibble <- as_tibble(meta)
# Ver un resumen de los datos disponibles de genotipado y los metadatos. Así visualizamos fácilmente que ambas matrices no tienen el mismo tamaño
geno_tibble
## # A tibble: 275 × 5,766
## ID M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 M12
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 S1 A U B U U B B B B B A A
## 2 S10 A A A A A A A A A B B A
## 3 S100 A B A A A B B B B A A A
## 4 S101 A H H H H H H H H H H A
## 5 S102 A B A A A B B B B A A A
## 6 S103 A A A A A A A A A H H A
## 7 S104 A A A A A A A A A A A A
## 8 S105 A H H H H H H H H H H A
## 9 S106 A A A A A A A A A B B A
## 10 S107 A A A A A A A A A B B A
## # … with 265 more rows, and 5,753 more variables: M13 <chr>, M14 <chr>,
## # M15 <chr>, M16 <chr>, M17 <chr>, M18 <chr>, M19 <chr>, M20 <chr>,
## # M21 <chr>, M22 <chr>, M23 <chr>, M24 <chr>, M25 <chr>, M26 <chr>,
## # M27 <chr>, M28 <chr>, M29 <chr>, M30 <chr>, M31 <chr>, M32 <chr>,
## # M33 <chr>, M34 <chr>, M35 <chr>, M36 <chr>, M37 <chr>, M38 <chr>,
## # M39 <chr>, M40 <chr>, M41 <chr>, M42 <chr>, M43 <chr>, M44 <chr>,
## # M45 <chr>, M46 <chr>, M47 <chr>, M48 <chr>, M49 <chr>, M50 <chr>, …
meta_tibble
## # A tibble: 294 × 4
## ID Species Country Pheno
## <chr> <chr> <chr> <dbl>
## 1 S1 colocynthis Turkey 8.8
## 2 S2 lanatus Spain 0.8
## 3 S3 amarus Turkey 7
## 4 S4 mucosospermus China 9
## 5 S5 rehmii China 9.7
## 6 S6 rehmii China 10.9
## 7 S7 mucosospermus China 8.4
## 8 S8 lanatus Spain 1.3
## 9 S9 lanatus Spain 1.3
## 10 S10 lanatus Spain 2
## # … with 284 more rows
### 1.4) Unir genotipos + resto de datos por columna comun 'ID', conservando solo las muestras que tienen datos de genotipado
data <- as_tibble(inner_join(geno_tibble, meta_tibble, by="ID"))
### 1.5) Reordenar columnas, por preferencia de estructura y orden
meta_cols <- c("ID", "Species", "Country", "Pheno")
# El resto de columnas son marcadores
marker_cols <- setdiff(names(data), meta_cols)
# Reordenamos: datos primero, luego marcadores
data_ordered <- data %>%
select(all_of(meta_cols), all_of(marker_cols))
### 1.6) Guardar dataset fusionado y ordenado
readr::write_csv(data_ordered, "WM_GBK_5K_alldata.csv")
### 1.7) Sustituir U por NA y guardamos un objeto nuevo: 'data_clean'
data_clean <- data_ordered %>%
mutate(
across(all_of(marker_cols), ~ ifelse(.x == "U", NA, .x)))
#---------------------------------------------------------
# 2) FILTRADO DE DATOS POR PARÁMETROS DE CALIDAD (QC)
#---------------------------------------------------------
### 2.1) QC1: Eliminar marcadores con > 10% NA
# Calculamos proporción de NA por marcador y nos quedamos con los marcadores que tienen <= 10% NA
marker_cols_clean <- setdiff(names(data_clean), meta_cols)
prop_na_markers <- colMeans(is.na(data_clean[, marker_cols_clean]))
marker_cols_qc1 <- names(prop_na_markers)[prop_na_markers <= 0.10]
#Creamos un nuevo tibble: metadatos + marcadores filtrados
data_qc1 <- data_clean %>%
select(all_of(meta_cols), all_of(marker_cols_qc1))
### 2.2) QC2: Eliminar muestras con > 10% NA
# Calculamos proporción de NA por muestra y conservamos muestras con <= 10% de NA
marker_cols_qc1 <- setdiff(names(data_qc1), meta_cols)
prop_na_samples <- rowMeans(is.na(data_qc1[, marker_cols_qc1]))
data_qc2 <- data_qc1[prop_na_samples <= 0.10, ]
# Comprobaciones: Marcadores que han pasado el control de calidad
n_markers_before <- length(setdiff(names(data_clean), meta_cols))
n_markers_after <- length(setdiff(names(data_qc1), meta_cols))
n_markers_before
## [1] 5765
n_markers_after
## [1] 5108
# Confirmamos que ningún marcador de data_qc1 supera el 10% de NA (el valor devuelto debería ser ≤ 0.10)
max(colMeans(is.na(data_qc1[marker_cols_qc1])))
## [1] 0.09818182
# Comprobaciones: Muestras que han pasado el control de calidad
n_samples_before <- nrow(data_qc1)
n_samples_after <- nrow(data_qc2)
n_samples_before
## [1] 275
n_samples_after
## [1] 254
# Confirmamos que todas las muestras cumplen el criterio <= 10% NA
max(rowMeans(is.na(data_qc2[marker_cols_qc1])))
## [1] 0.08574785
### 2.3) QC3: Eliminar marcadores con MAF < 0.10 (poco polimórficos)
# Identificar columnas de marcadores
marker_cols_qc2 <- setdiff(names(data_qc2), meta_cols)
# Creamos función para calcular MAF de un vector de genotipos
calc_maf <- function(g) {
g <- g[!is.na(g)] # quitar NAs
A <- sum(g == "A")
H <- sum(g == "H")
B <- sum(g == "B")
total_genotipos <- A + H + B
if (total_genotipos == 0) return(NA_real_)
total_alelos <- 2 * total_genotipos
B_alleles <- 2 * B + H # Alelos B (H cuenta como 1, B cuenta como 2)
pB <- B_alleles / total_alelos # frecuencia del alelo B
maf <- min(pB, 1 - pB) # frecuencia menor
return(maf)
}
# Calculamos MAF para todos los marcadores y comprobamos cuántos marcadores tienen MAF < 0.10
maf_values <- sapply(data_qc2[marker_cols_qc2], calc_maf)
sum(maf_values < 0.10, na.rm = TRUE)
## [1] 351
# Seleccionamos los marcadores que pasan el umbral MAF >= 0.10 y construimos un nuevo tibble: datos + marcadores filtrados
marker_cols_maf <- names(maf_values)[maf_values >= 0.10]
data_filtered <- data_qc2 %>%
select(all_of(meta_cols), all_of(marker_cols_maf))
# Guardamos este dataset nuevo
readr::write_csv(data_filtered, "WM_GBK_5K_alldata_filtered.csv")
## Resumen QC y lo guardamos en PDF y csv
qc_summary <- data.frame(
Etapa = c(
"Inicial",
"QC1: Excluir marcadores > 10% NA",
"QC2: Excluir muestras > 10% NA",
"QC3: Excluir marcadores con MAF < 0.10"
),
Muestras = c(
nrow(data_ordered),
nrow(data_qc1),
nrow(data_qc2),
nrow(data_filtered)
),
Marcadores = c(
length(setdiff(names(data_ordered), meta_cols)),
length(setdiff(names(data_qc1), meta_cols)),
length(setdiff(names(data_qc2), meta_cols)),
length(setdiff(names(data_filtered), meta_cols))
)
)
qc_summary
## Etapa Muestras Marcadores
## 1 Inicial 275 5765
## 2 QC1: Excluir marcadores > 10% NA 275 5108
## 3 QC2: Excluir muestras > 10% NA 254 5108
## 4 QC3: Excluir marcadores con MAF < 0.10 254 4757
pdf("WM_QC_resumen_limpieza_datos.pdf", width = 10, height = 4)
grid.table(qc_summary, rows = NULL) # No mostrar números de fila
dev.off()
## png
## 2
write.csv(qc_summary, "WM_QC_resumen_limpieza_datos.csv", row.names = FALSE)
#------------------------------------------
# 3) ANÁLISIS – TODAS LAS ESPECIES
#------------------------------------------
### 3.1) Barplot muestras por especie
conteo_df <- data_filtered %>% count(Species)
barplot_species <- ggplot(
conteo_df,
aes(x = n, y = reorder(Species, -n), fill = Species)) +
geom_col() +
geom_text(aes(label = n), hjust = -0.1, size = 3) +
labs(
title = "Número de muestras por especie",
x = "Número de muestras",
y = "Especie"
) +
theme_minimal() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5)) +
xlim(0, max(conteo_df$n) * 1.1)
barplot_species
ggsave("1_barplot_muestras_especie.png", barplot_species, width = 8, height = 6, dpi = 300)
### 3.2) Pie chart muestras por países
percent_country <- data_filtered %>%
count(Country) %>%
mutate(pct = n / sum(n) * 100)
gg_pie_country <- ggplot(
percent_country,
aes(x = "", y = n, fill = Country)
) +
geom_col(color = "white") +
coord_polar(theta = "y") +
theme_void() +
labs(title = "Distribución de muestras por país",
fill = "País") +
theme(plot.title = element_text(hjust = 0.5)) +
geom_text(
aes(label = paste0(round(pct, 1), "%")),
position = position_stack(vjust = 0.5),
size = 5
)
gg_pie_country
ggsave("2_piechart_muestras_pais.png", gg_pie_country, width = 8, height = 6, dpi = 300)
### 3.3) Barplot muestras por especie y país
conteo_especie_pais <- data_filtered %>%
count(Country, Species)
barplot_sp_country <- ggplot(
conteo_especie_pais,
aes(x = Country, y = n, fill = Species)) +
geom_col(position = "stack") +
labs(
title = "Número de muestras por país y especie",
x = "País",
y = "Número de muestras",
fill = "Especie"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 40, hjust = 1)
)
barplot_sp_country
ggsave("3_barplot_muestras_especie_pais.png", barplot_sp_country, width = 8, height = 6, dpi = 300)
### 3.4) Boxplot firmeza por especie
ggplot(
data_filtered,
aes(x = Species, y = Pheno, fill = Species)) +
geom_boxplot() +
labs(title = "Distribución de firmeza por especie",
x = "Especie",
y = "Firmeza (kg/cm²)") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5),
legend.position = "none",
axis.text.x = element_text(angle = 30, hjust = 1))
ggsave("4_boxplot_pheno_especie.png", width = 8, height = 6, dpi = 300)
## 3.4.1) ANOVA firmeza por especie: ¿hay diferencias significativas de firmeza entre especies?
# Filtramos especies con n >= 2 para poder aplicar estadística
df_sp <- data_filtered %>%
group_by(Species) %>%
filter(n() >= 2) %>%
ungroup()
modelo_sp <- aov(Pheno ~ Species, data = df_sp)
summary(modelo_sp)
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 3 1342.6 447.5 175.5 <2e-16 ***
## Residuals 249 635.1 2.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comprobamos normalidad de los residuos y la homogeneidad de las varianzas. Como el resultado es negativo, hacemos un Welch ANOVA.
shapiro.test(residuals(modelo_sp))
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_sp)
## W = 0.66475, p-value < 2.2e-16
bartlett.test(Pheno ~ Species, data = df_sp)
##
## Bartlett test of homogeneity of variances
##
## data: Pheno by Species
## Bartlett's K-squared = 23.385, df = 3, p-value = 3.357e-05
welch_sp <- oneway.test(Pheno ~ Species, data = df_sp)
welch_sp
##
## One-way analysis of means (not assuming equal variances)
##
## data: Pheno and Species
## F = 664.81, num df = 3.000, denom df = 13.051, p-value = 1.627e-14
# Comparaciones entre especies (Welch). En este dataset, hay diferencias significativas de firmeza entre todas las especies. La que más firmeza tiene es C. amarus.
pairwise.t.test(df_sp$Pheno, df_sp$Species,
p.adjust.method = "BH",
pool.sd = FALSE)
##
## Pairwise comparisons using t tests with non-pooled SD
##
## data: df_sp$Pheno and df_sp$Species
##
## amarus colocynthis lanatus
## colocynthis 9.1e-09 - -
## lanatus < 2e-16 9.1e-09 -
## mucosospermus 3.6e-07 0.0022 < 2e-16
##
## P value adjustment method: BH
### 3.5) PCA con prcomp() – todas las especies
# Separar metadatos y genotipos (A/H/B/NA): meta_all y geno_all
meta_all <- data_filtered[, meta_cols]
marker_cols_filtered <- setdiff(names(data_filtered), meta_cols)
geno_all <- data_filtered[, marker_cols_filtered]
# Convertir A/H/B -> 0/1/2, manteniendo NA
geno_num <- apply(
geno_all, 2,
function(x) ifelse(x == "A", 0,
ifelse(x == "H", 1,
ifelse(x == "B", 2, NA)))
)
rownames(geno_num) <- meta_all$ID
# Imputación de NA por la media del marcador
geno_num_pca <- geno_num
for (i in 1:ncol(geno_num_pca)) {
nas <- is.na(geno_num_pca[, i])
if (any(nas)) {
geno_num_pca[nas, i] <- mean(geno_num_pca[, i], na.rm = TRUE)
}
}
# calcular varianza de cada marcador y mantener solo los de varianza > 0 (evitamos errores en prcomp)
var_cols_all <- apply(geno_num_pca, 2, var)
geno_num_pca <- geno_num_pca[, var_cols_all > 0, drop = FALSE]
# Ejecutar PCA
pca_res <- prcomp(
geno_num_pca,
center = TRUE,
scale. = TRUE # center = TRUE y scale. = TRUE para trabajar con marcadores en misma escala
)
# Calcular % de varianza explicada por cada componente
var_exp <- 100 * (pca_res$sdev^2 / sum(pca_res$sdev^2))
# Miramos los dos primeros para comprobar
var_exp[1:2]
## [1] 16.009428 8.196635
# var_exp[1] = % varianza explicada por PC1
# var_exp[2] = % varianza explicada por PC2
# Construir data.frame con scores del PCA + metadatos para análisis y gráficos posteriores
pca_df <- as.data.frame(pca_res$x) %>%
mutate(
ID = meta_all$ID,
Species = meta_all$Species,
Country = meta_all$Country,
Pheno = meta_all$Pheno
)
## 3.5.1) Gráfico PCA todas las especies (color = especie)
gg_all_species <- ggplot(
pca_df,
aes(x = PC1, y = PC2, color = Species)
) +
geom_point(size = 2) +
coord_equal() +
theme_minimal() +
labs(
title = "PCA - Todas las especies",
x = paste0("PC1 (", round(var_exp[1], 2), "%)"),
y = paste0("PC2 (", round(var_exp[2], 2), "%)")
) +
theme(plot.title = element_text(hjust = 0.5))
gg_all_species
ggsave("5_PCA_especies.png", gg_all_species, width = 8, height = 6, dpi = 300)
## 3.5.2) Gráfico PCA todas las especies (color = firmeza)
gg_all_pheno <- ggplot(
pca_df,
aes(x = PC1, y = PC2, color = Pheno)
) +
geom_point(size = 2) +
scale_color_viridis_c(option = "plasma") +
coord_equal() +
theme_minimal() +
labs(
title = "PCA - Todas las especies (firmeza)",
x = paste0("PC1 (", round(var_exp[1], 2), "%)"),
y = paste0("PC2 (", round(var_exp[2], 2), "%)"),
color = "Firmeza (kg/cm²)"
) +
theme(plot.title = element_text(hjust = 0.5))
gg_all_pheno
ggsave("6_PCA_especies_firmeza.png", gg_all_pheno, width = 8, height = 6, dpi = 300)
## 3.5.3) Gráfico PCA todas las especies (color = Pais)
gg_all_country <- ggplot(
pca_df,
aes(x = PC1, y = PC2, color = Country)
) +
geom_point(size = 2) +
coord_equal() +
theme_minimal() +
labs(
title = "PCA - Todas las especies (país)",
x = paste0("PC1 (", round(var_exp[1], 2), "%)"),
y = paste0("PC2 (", round(var_exp[2], 2), "%)"),
color = "Pais"
) +
theme(plot.title = element_text(hjust = 0.5))
gg_all_country
ggsave("7_PCA_especies_pais.png", gg_all_country, width = 8, height = 6, dpi = 300)
### 3.6) Estudio de correlación entre componentes principales y firmeza
## 3.6.1) Test de normalidad para PC1 y firmeza
# Comprobamos si las variables siguen una distribución normal. Elegimos el test Shapiro-Wilk, en lugar de Kolmogorov-Smirnof, más adecuado en este caso
shapiro.test(pca_df$PC1) # normalidad de PC1
##
## Shapiro-Wilk normality test
##
## data: pca_df$PC1
## W = 0.77444, p-value < 2.2e-16
shapiro.test(pca_df$Pheno) # normalidad de firmeza
##
## Shapiro-Wilk normality test
##
## data: pca_df$Pheno
## W = 0.73767, p-value < 2.2e-16
# Ambas variables muestran p < 0.05 -> se rechaza la normalidad.
# Por tanto, no es adecuado usar Pearson para el test de correlación posterior. Aunque Spearman funciona, produce avisos por valores repetidos ("ties"). Kendall es la opción no paramétrica más adecuada en este caso.
## 3.6.2) Correlación no paramétrica (Kendall) PC1 ~ firmeza
cor.test(pca_df$PC1, pca_df$Pheno, method = "kendall")
##
## Kendall's rank correlation tau
##
## data: pca_df$PC1 and pca_df$Pheno
## z = 5.6433, p-value = 1.669e-08
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.2426514
# Resultado: asociación positiva moderada y significativa, indicando que la estructura genética principal de PC1 está asociada con diferencias en firmeza entre accesiones.
## 3.6.3) Scatterplot PC1–firmeza coloreado por especie
# En el siguiente scatterplot vemos cómo las especies silvestres, especialmente amarus (mayor firmeza), se sitúa a la derecha en el eje del PC1 respecto a lanatus (menor firmeza)
gg_PC1_firm <- ggplot(
pca_df,
aes(x = PC1, y = Pheno)
) +
geom_point(aes(color = Species), size = 2, alpha = 0.8) +
theme_minimal(base_size = 13) +
labs(
title = "Relación PC1–firmeza (todas las especies)",
x = paste0("PC1 (", round(var_exp[1], 2), "% var. explicada)"),
y = "Firmeza (kg/cm²)",
color = "Especie"
) +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "right"
)
gg_PC1_firm
ggsave("8_PC1_firmeza_scatterplot.png", gg_PC1_firm, width = 8, height = 6, dpi = 300)
## 3.6.4) Test de normalidad para PC2 y correlación con firmeza
# Repetimos el análisis para PC2.
shapiro.test(pca_df$PC2) # Tampoco es normal
##
## Shapiro-Wilk normality test
##
## data: pca_df$PC2
## W = 0.935, p-value = 3.798e-09
# 3.6.5) Correlación no paramétrica (Kendall) PC2 ~ firmeza
cor.test(pca_df$PC2, pca_df$Pheno, method = "kendall")
##
## Kendall's rank correlation tau
##
## data: pca_df$PC2 and pca_df$Pheno
## z = 1.1189, p-value = 0.2632
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.04811026
# Resultado: En este caso no hay correlación significativa entre PC2 y firmeza
### 3.7) Dendrograma circular (todas las especies)
# Cálculo de distancias genéticas a partir de matriz 0-2 (geno_num_pca) calculada previamente para el PCA
dist_ind_all <- dist(geno_num_pca, method = "euclidean")
hc_all <- hclust(dist_ind_all, method = "average")
dend_all <- as.dendrogram(hc_all)
# Colorear etiquetas por especie
species_factor <- factor(meta_all$Species)
species_levels <- levels(species_factor)
species_colors <- setNames(rainbow(length(species_levels)), species_levels)
cols_by_id <- species_colors[species_factor]
names(cols_by_id) <- meta_all$ID
cols_ordered <- cols_by_id[labels(dend_all)]
labels_colors(dend_all) <- cols_ordered
labels_cex(dend_all) <- 0.6
#Primero mostramos el dendrograma en el documento
# Parámetros de la circunferencia
par(mar = c(1, 1, 4, 1))
circos.par(start.degree = 90, gap.after = 0, track.margin = c(0.01, 0.01))
circlize_dendrogram(
dend_all,
facing = "outside",
labels_track_height = 0.35, # espacio para etiquetas
dend_track_height = 0.6 # espacio para el árbol
)
title(
main = "Dendrograma circular de todas las especies", cex.main = 1.6, line = -1
)
legend(
x = 0.65, y = 0.65,
legend = species_levels,
col = species_colors,
pch = 19,
bty = "n",
cex = 1.4
)
circos.clear()
#Guardamos el dendrograma
png("9_Dendrograma_circular_especies.png", width = 2200, height = 2200, res = 200)
par(mar = c(1, 1, 4, 1))
circos.par(start.degree = 90, gap.after = 0, track.margin = c(0.01, 0.01))
circlize_dendrogram(
dend_all,
facing = "outside",
labels_track_height = 0.35, # espacio para etiquetas
dend_track_height = 0.6 # espacio para el árbol
)
title(
main = "Dendrograma circular de todas las especies", cex.main = 1.6, line = -1
)
legend(
x = 0.65, y = 0.65,
legend = species_levels,
col = species_colors,
pch = 19,
bty = "n",
cex = 1.4
)
circos.clear()
dev.off()
## png
## 2
### 3.8) Heatmap todas las especies
# Matriz de distancias normalizada 0–1. Primero convertimos el objeto 'dist' a matriz.
dist_mat_all <- as.matrix(dist_ind_all)
# Aseguramos nombres de filas/columnas = IDs
rownames(dist_mat_all) <- meta_all$ID
colnames(dist_mat_all) <- meta_all$ID
# Normalización: distancia / distancia máxima → escala 0–1
max_dist_all <- max(dist_mat_all)
dist_mat_all_scaled <- dist_mat_all / max_dist_all
# Anotaciones por especie (barras de color)
annotation_all <- data.frame(
Species = factor(meta_all$Species)
)
rownames(annotation_all) <- meta_all$ID
species_levels_all <- levels(annotation_all$Species)
species_colors_all <- setNames(rainbow(length(species_levels_all)), species_levels_all)
ann_colors_all <- list(Species = species_colors_all)
#Mostrar heatmap en el documento
pheatmap(
mat = dist_mat_all_scaled,
cluster_rows = hc_all,
cluster_cols = hc_all,
annotation_row = annotation_all,
annotation_col = annotation_all,
annotation_colors = ann_colors_all,
color = colorRampPalette(brewer.pal(9, "YlOrRd"))(200),
show_rownames = FALSE,
show_colnames = FALSE,
treeheight_row = 80,
treeheight_col = 80,
main = "Heatmap de distancias geneticas (todas las especies)",
filename = NA
)
#Guardar el heatmap
pheatmap(
mat = dist_mat_all_scaled,
cluster_rows = hc_all,
cluster_cols = hc_all,
annotation_row = annotation_all,
annotation_col = annotation_all,
annotation_colors = ann_colors_all,
color = colorRampPalette(brewer.pal(9, "YlOrRd"))(200),
show_rownames = FALSE,
show_colnames = FALSE,
treeheight_row = 80,
treeheight_col = 80,
main = "Heatmap de distancias geneticas (todas las especies)",
filename = "10_Heatmap_especies_normalizado.png"
)
#------------------------------------------
# 4) ANÁLISIS – SOLO C. LANATUS
#------------------------------------------
### 4.1) Boxplot firmeza por país (lanatus)
data_lan <- data_filtered %>% filter(Species == "lanatus")
gg_box_lan_country <- ggplot(
data_lan,
aes(x = Country, y = Pheno, fill = Country)
) +
geom_boxplot() +
labs(
title = "Distribución de firmeza por país (C. lanatus)",
x = "Pais",
y = "Firmeza (kg/cm²)"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(angle = 40, hjust = 1),
legend.position = "none"
)
gg_box_lan_country
ggsave("11_boxplot_firmeza_lanatus_pais.png", gg_box_lan_country, width = 8, height = 6, dpi = 300)
## 4.1.1) ANOVA firmeza por pais: ¿hay diferencias significativas de firmeza entre países para C. lanatus?
# Hacemos este analisis solo en C. lanatus porque es la que tiene un número más representativo de muestras de diferentes paises.
modelo <- aov(Pheno ~ Country, data = data_lan)
summary(modelo)
## Df Sum Sq Mean Sq F value Pr(>F)
## Country 4 189.9 47.48 23.77 5.56e-16 ***
## Residuals 192 383.5 2.00
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comprobamos normalidad de los residuos y la homogeneidad de las varianzas. Como el resultado es negativo, hacemos un Welch ANOVA.
shapiro.test(residuals(modelo))
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo)
## W = 0.84921, p-value = 5.109e-13
bartlett.test(Pheno ~ Country, data = data_lan)
##
## Bartlett test of homogeneity of variances
##
## data: Pheno by Country
## Bartlett's K-squared = 303.15, df = 4, p-value < 2.2e-16
welch <- oneway.test(Pheno ~ Country, data = data_lan)
welch
##
## One-way analysis of means (not assuming equal variances)
##
## data: Pheno and Country
## F = 7.7351, num df = 4.000, denom df = 93.541, p-value = 1.981e-05
# Comparaciones entre paises (Welch). Se detectan únicamente diferencias significativas entre Peru y el resto de países.
pairwise.t.test(data_lan$Pheno, data_lan$Country,
p.adjust.method = "BH",
pool.sd = FALSE)
##
## Pairwise comparisons using t tests with non-pooled SD
##
## data: data_lan$Pheno and data_lan$Country
##
## China Peru Spain Turkey
## Peru 5.1e-06 - - -
## Spain 0.74 5.0e-06 - -
## Turkey 0.74 5.0e-06 0.74 -
## USA 0.74 5.0e-06 1.00 0.74
##
## P value adjustment method: BH
### 4.2) PCA para lanatus
lan_idx <- meta_all$Species == "lanatus" # Nos aseguramos de que el valor en la columna Species es "lanatus"
meta_lan <- meta_all[lan_idx, ]
geno_num_pca_lan <- geno_num_pca[lan_idx, , drop = FALSE]
# Calcular varianza de cada marcador
var_cols_lan <- apply(geno_num_pca_lan, 2, var)
# Mantener solo con marcadores de varianza > 0
geno_num_pca_lan2 <- geno_num_pca_lan[, var_cols_lan > 0, drop = FALSE]
pca_res_lan <- prcomp(
geno_num_pca_lan2,
center = TRUE,
scale. = TRUE
)
var_exp_lan <- 100 * (pca_res_lan$sdev^2 / sum(pca_res_lan$sdev^2))
var_exp_lan[1:2]
## [1] 9.295461 8.265650
pca_df_lan <- as.data.frame(pca_res_lan$x) %>%
mutate(
ID = meta_lan$ID,
Country = meta_lan$Country,
Pheno = meta_lan$Pheno
)
## 4.2.1) Gráfico PCA lanatus (color = País)
gg_pca_lan_country <- ggplot(
pca_df_lan,
aes(x = PC1, y = PC2, color = Country)
) +
geom_point(size = 2.5) +
coord_equal() +
theme_minimal() +
labs(
title = "PCA - C.lanatus (color por país)",
x = paste0("PC1 (", round(var_exp_lan[1], 2), "%)"),
y = paste0("PC2 (", round(var_exp_lan[2], 2), "%)")
) +
theme(plot.title = element_text(hjust = 0.5))
gg_pca_lan_country
ggsave("12_PCA_lanatus_pais.png", gg_pca_lan_country, width = 8, height = 6, dpi = 300)
## 4.2.2) Gráfico PCA lanatus (color = Firmeza)
gg_pca_lan_pheno <- ggplot(
pca_df_lan,
aes(x = PC1, y = PC2, color = Pheno)
) +
geom_point(size = 2) +
scale_color_viridis_c(option = "plasma") +
coord_equal() +
theme_minimal() +
labs(
title = "PCA - C. lanatus (color por firmeza)",
x = paste0("PC1 (", round(var_exp_lan[1], 2), "%)"),
y = paste0("PC2 (", round(var_exp_lan[2], 2), "%)"),
color = "Firmeza (kg/cm²)"
) +
theme(plot.title = element_text(hjust = 0.5))
gg_pca_lan_pheno
ggsave("13_PCA_lanatus_firmeza.png", gg_pca_lan_pheno, width = 8, height = 6, dpi = 300)
### 4.3) Estudio de correlación entre componentes principales y firmeza en C.lanatus
## 4.3.1) Test de normalidad para PC1 y firmeza (C.lanatus)
# Comprobamos si las variables siguen una distribución normal. Elegimos nuevamente el test Shapiro-Wilk.
shapiro.test(pca_df_lan$PC1)
##
## Shapiro-Wilk normality test
##
## data: pca_df_lan$PC1
## W = 0.84469, p-value = 3.111e-13
shapiro.test(pca_df_lan$Pheno)
##
## Shapiro-Wilk normality test
##
## data: pca_df_lan$Pheno
## W = 0.56868, p-value < 2.2e-16
# Ambas variables muestran p < 0.05 -> se rechaza la normalidad. Usaremos Kendall.
## 4.3.2) Correlación no paramétrica (Kendall) PC1 ~ firmeza
cor.test(pca_df_lan$PC1, pca_df_lan$Pheno, method = "kendall")
##
## Kendall's rank correlation tau
##
## data: pca_df_lan$PC1 and pca_df_lan$Pheno
## z = -0.6905, p-value = 0.4899
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.03406607
# Resultado: En lanatus, la variación genética capturada por PC1 no está asociada a la firmeza.
## 4.3.3) Test de normalidad para PC2 y correlación con firmeza
# Repetimos el análisis para PC2.
shapiro.test(pca_df_lan$PC2) # Tampoco es normal
##
## Shapiro-Wilk normality test
##
## data: pca_df_lan$PC2
## W = 0.92084, p-value = 8.214e-09
## 4.3.4) Correlación no paramétrica (Kendall) PC2 ~ firmeza
cor.test(pca_df_lan$PC2, pca_df_lan$Pheno, method = "kendall")
##
## Kendall's rank correlation tau
##
## data: pca_df_lan$PC2 and pca_df_lan$Pheno
## z = 0.44335, p-value = 0.6575
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## 0.02187288
# Resultado: En este caso tampoco no hay correlación significativa entre PC2 y firmeza.
### 4.4) Dendrograma lanatus (color por país)
# Distancias genéticas (solo lanatus)
dist_lan <- dist(geno_num_pca_lan2, method = "euclidean")
hc_lan <- hclust(dist_lan, method = "average")
dend_lan <- as.dendrogram(hc_lan)
# Colorear etiquetas por país
country_factor <- factor(meta_lan$Country)
country_levels <- levels(country_factor)
country_colors <- setNames(rainbow(length(country_levels)), country_levels)
cols_by_id_lan <- country_colors[country_factor]
names(cols_by_id_lan) <- meta_lan$ID
# Reordenar según labels del dendrograma
cols_ordered_lan <- cols_by_id_lan[labels(dend_lan)]
labels_colors(dend_lan) <- cols_ordered_lan
labels_cex(dend_lan) <- 0.6
#Primero mostramos en el documento
par(mar = c(1, 1, 4, 1))
circos.par(start.degree = 90, gap.after = 0, track.margin = c(0.01, 0.01))
circlize_dendrogram(
dend_lan,
facing = "outside",
labels_track_height = 0.35,
dend_track_height = 0.6
)
title(
main = "Dendrograma circular lanatus (color por país)", cex.main = 1.6, line = -1
)
legend(
x = 0.9, y = 0.5,
legend = country_levels,
col = country_colors,
pch = 19,
bty = "n",
cex = 1.4,
)
circos.clear()
#Guardamos el dendrograma
png("14_dendrograma_lanatus_pais.png", width = 2200, height = 2200, res = 200)
par(mar = c(1, 1, 4, 1))
circos.par(start.degree = 90, gap.after = 0, track.margin = c(0.01, 0.01))
circlize_dendrogram(
dend_lan,
facing = "outside",
labels_track_height = 0.35,
dend_track_height = 0.6
)
title(
main = "Dendrograma circular lanatus (color por país)", cex.main = 1.6, line = -1
)
legend(
x = 0.9, y = 0.5,
legend = country_levels,
col = country_colors,
pch = 19,
bty = "n",
cex = 1.4,
)
circos.clear()
dev.off()
## png
## 2
### 4.5) Heatmap lanatus
# Matriz de distancias
dist_mat_lan <- as.matrix(dist_lan)
# Nombres de filas/columnas = ID
rownames(dist_mat_lan) <- meta_lan$ID
colnames(dist_mat_lan) <- meta_lan$ID
# Normalización 0–1
max_dist_lan <- max(dist_mat_lan)
dist_mat_lan_scaled <- dist_mat_lan / max_dist_lan
# Anotación por país
annotation_lan <- data.frame(
Country = factor(meta_lan$Country)
)
rownames(annotation_lan) <- meta_lan$ID
country_levels_hm <- levels(annotation_lan$Country)
country_colors_hm <- setNames(
rainbow(length(country_levels_hm)),
country_levels_hm
)
ann_colors_lan <- list(Country = country_colors_hm)
pheatmap(
mat = dist_mat_lan_scaled,
cluster_rows = hc_lan,
cluster_cols = hc_lan,
annotation_row = annotation_lan,
annotation_col = annotation_lan,
annotation_colors = ann_colors_lan,
color = colorRampPalette(brewer.pal(9, "YlOrRd"))(200),
show_rownames = FALSE,
show_colnames = FALSE,
treeheight_row = 80,
treeheight_col = 80,
main = "Heatmap distancias geneticas lanatus",
filename = NA
)
pheatmap(
mat = dist_mat_lan_scaled,
cluster_rows = hc_lan,
cluster_cols = hc_lan,
annotation_row = annotation_lan,
annotation_col = annotation_lan,
annotation_colors = ann_colors_lan,
color = colorRampPalette(brewer.pal(9, "YlOrRd"))(200),
show_rownames = FALSE,
show_colnames = FALSE,
treeheight_row = 80,
treeheight_col = 80,
main = "Heatmap distancias geneticas lanatus",
filename = "15_heatmap_lanatus_normalizado.png"
)
# De estos datos concluimos que en el caso de C.lanatus no hay una estructura genética clara asociada al país de origen.
# Encontramos diferencias significativas en firmeza entre países, destacando que varias accesiones procedentes de Perú presentan valores superiores al resto.