Clustering: Agrupación de muestras en función de su similitud
1 julio 2021
Abstract
En general, el objetivo de la exploración de datos a través del aprendizaje no supervisado se centra en comprender cómo se relacionan entre sí las variables de una base de datos y cómo se relacionan entre sí las muestras definidas por esas variables. Estos puntos de información se pueden utilizar para generar una hipótesis, encontrar valores atípicos en las muestras o identificar grupos de muestras que necesitan más puntos de datos.
Contents
1 Motivación
Con mucha frecuencia se va a querer evaluar cómo las muestras se relacionan entre sí. Por ejemplo, en las aplicaciones, se pueden abordar las siguientes preguntas:
¿Las muestras son similares entre sí?
¿Las muestras del mismo grupo de tratamiento tienen señales similares en todo el genoma?
Para paciente con enfermedades similares, ¿se tiene un perfil de expresión génica similar?
Para responder algunas de estas interrogantes se necesitan definir métricas de distancia o similitud entre perfiles de expresión de cada paciente y usar esa métrica para encontrar agrupaciones de pacientes que sean más similares entre sí. Dada una métrica de distancia, se emplea una metodología para encontrar agrupaciones auto-similares. La agrupación es un procedimiento presente en cualquier campo que se ocupe de datos de alta dimensión.
2 Métricas de distancia
La métrica de distancia es simplemente una medida de cuán similares son las expresiones entre sí. Hay muchas opciones para las métricas de distancia y la elección de la métrica es bastante importante para la agrupación. Por ejemplo, se cuenta con muestras de cuatro pacientes y la expresión de tres genes se mide en la tabla 1 y se desean encontrar similitudes en función de sus perfiles de expresión genética.
| Muestra | IRX4 | OCT4 | PAX6 |
|---|---|---|---|
| paciente 1 | 11 | 10 | 1 |
| paciente 2 | 13 | 13 | 3 |
| paciente 3 | 12 | 4 | 10 |
| paciente 4 | 1 | 3 | 9 |
Al graficar el perfil de expresión génica de cada paciente (que se muestra en la Figura 1), se observa que los perfiles de expresión de la paciente \(1\) y de la paciente \(2\) son más similares entre sí que la paciente \(3\) o paciente \(4\).
Figure 1: Valores de expresión génica para diferentes pacientes
Paciente 1 y paciente 2 tienen valores de expresión genética similares.
Para cuantificar la visualización, una métrica simple para la distancia entre los vectores de expresión génica entre una pareja de pacientes dada es la suma de la diferencia absoluta entre los valores de expresión génica. Esto se puede formular de la siguiente manera: \(d_{AB}=\sum_{i=1}^n|e_{Ai}-e_{Bi}|\), donde \(d_{AB}\) es la distancia entre las pacientes \(A\) y \(B\), y \(e_{Ai}\) y \(e_{Bi}\) son valores de expresión del gen \(i-\)ésimo para las pacientes \(A\) y \(B\). Esta métrica de distancia se llama distancia de Manhattan o norma L1.
Otra métrica de distancia usa la suma de distancias al cuadrado y toma la raíz cuadrada del valor resultante; esta métrica se puede formular como: \(d_{AB}=\sqrt{\sum_{i=1}^n(e_{Ai}-e_{Bi})^2}\). Esta distancia se llama Distancia euclidiana o norma L2. Esta suele ser la métrica de distancia predeterminada para muchos algoritmos de agrupamiento. Debido a la operación de cuadratura, los valores que son muy diferentes obtienen una mayor contribución a la distancia. Debido a esto, en comparación con la distancia de Manhattan, puede verse más afectada por valores atípicos. Pero, en general, si los valores atípicos son raros, esta métrica de distancia funciona bien.
Una tercera métrica es la distancia de correlación. Esta es simplemente \(d_{AB}=1-{\rho}\), donde \({\rho}\) es el coeficiente de correlación de Pearson entre dos vectores; en el ejemplo anterior, esos vectores son perfiles de expresión génica de pacientes. Usando esta distancia, los vectores de expresión génica que tienen un patrón similar tendrán una distancia pequeña, mientras que cuando los vectores tienen patrones diferentes tendrán una distancia grande. En esta métrica, la correlación lineal entre vectores importa, aunque la escala de los vectores puede ser diferente.
De la tabla de expresión génica por paciente, se pueden ir calculandoe estas distancias.
gexp## # A tibble: 4 x 4
## Muestra IRX4 OCT4 PAX6
## <chr> <dbl> <dbl> <dbl>
## 1 paciente 1 11 10 1
## 2 paciente 2 13 13 3
## 3 paciente 3 12 4 10
## 4 paciente 4 1 3 9# distancia de Manhattan
dist(gexp[,-1],method="manhattan")## 1 2 3
## 2 7
## 3 16 17
## 4 25 28 13# distancia euclidena
dist(gexp[,-1],method="euclidean")## 1 2 3
## 2 4.123106
## 3 10.862780 11.445523
## 4 14.594520 16.733201 11.090537# distancia de correlación
as.dist(1-cor(t(gexp[,-1]))) ## 1 2 3
## 2 0.004129405
## 3 1.188982237 1.277350098
## 4 1.988522468 1.970725343 0.961538462La escala de los vectores en la matriz de datos puede afectar el cálculo de la distancia. Si los valores de expresión de un gen están en una escala mucho más alta que los otros genes, ese gen afectará la distancia más que otros cuando se use la distancia euclidiana o de Manhattan. La forma tradicional de escalar variables es estandarizando. Si esto se hace en todos los genes, cada gen tendrá el mismo efecto en las medidas de distancia. Si los valores de expresión génica han sido normalizados previamente entre pacientes, tener genes que dominen la métrica de distancia podría tener un significado biológico y, por lo tanto, puede que no sea deseable escalar más las variables.
scale(gexp[,-1])## IRX4 OCT4 PAX6
## [1,] 0.3147326 0.5212860 -1.0733721
## [2,] 0.6744270 1.1468293 -0.6214260
## [3,] 0.4945798 -0.7298004 0.9603856
## [4,] -1.4837394 -0.9383149 0.7344125
## attr(,"scaled:center")
## IRX4 OCT4 PAX6
## 9.25 7.50 5.75
## attr(,"scaled:scale")
## IRX4 OCT4 PAX6
## 5.560276 4.795832 4.4253063 Agrupación jerárquica
Esta es una de las algorítmicas de agrupamiento más recurrentes. Con este algoritmo, se puede ver la relación de los puntos de datos individuales y las relaciones de los clústeres o conglomerados. Esto se logra uniendo sucesivamente pequeñas agrupaciones entre sí en función de la distancia entre grupos. Finalmente, se obtiene una estructura de árbol o un dendrograma que muestra la relación entre los puntos de datos individuales y cada agrupación. La altura del dendrograma es la distancia entre agrupaciones. La función base para hacer agrupaciones jerárquicas es hclust(). El árbol de agrupamiento o dendrograma resultante para el ejemplo de \(4\) pacientes se muestra en la Figura 2.
# se reconstruye un data frame, pues una tibble rechaza las rownames
df <- as.data.frame(gexp[,-1])
rownames(df) <- Muestra
d=dist(df) # por defecto se toma la norma L2
hc=hclust(d,method = "complete")
plot(hc) # la gráfica retiene las rownames
Figure 2: Dendrograma de matriz de distancias
En el de código anterior, se ha utilizado el argumento method = "complete", el cual define los criterios que dirigen cómo se fusionan las subagrupaciones. Durante la agrupación, comenzando con grupos de un solo miembro, los grupos se fusionan en función de la distancia entre ellos. El argumento method controla eso, que puede tomar las siguientes opciones:
completesignifica Complete Linkage y la distancia entre dos conglomerados se define como la mayor distancia entre los miembros de las dos agrupaciones.singlesignifica Single Linkage y la distancia entre dos conglomerados se define como la distancia más pequeña entre los miembros de las dos agrupaciones.averagesignifica Average Linkage o más precisamente el método UPGMA (Unweighted Pair Group Method with Arithmetic Mean). En este caso, la distancia entre dos conglomerados se define como la distancia media entre cualquier miembro de las dos agrupaciones.ward.D2yward.Drepresentan diferentes implementaciones del método de variación mínima de Ward. Este método tiene como objetivo encontrar conglomerados esféricos y compactos mediante la selección de agrupaciones para fusionar en función del cambio en las variaciones de la agrupación. Los conglomerados se fusionan si el aumento de la varianza combinada sobre la suma de las varianzas específicas del conglomerado es el mínimo en comparación con las operaciones de fusión alternativas.
Con tidymodels se puede usar una recipe para codificar de manera más limpia los pasos anteriores. Con factoextra se visualizan las distinatas formas de enlace en las Figuras 3, 4 y 5
library(tidymodels)
# se crea la recipe sin la columna tipo <chr> y estandarizando
gexp_scaled <- recipe(~ ., data = gexp) %>%
step_rm(Muestra) %>%
step_normalize(all_predictors()) %>%
prep() %>%
bake(new_data = NULL) # aplica operaciones a los datos y crear matriz de diseño.
# se pasa la función dist bajo varios argumentos a la recipe
gexp_complete <- gexp_scaled %>%
dist() %>%
hclust(method = "complete")
gexp_average <- gexp_scaled %>%
dist() %>%
hclust(method = "average")
gexp_single <- gexp_scaled %>%
dist() %>%
hclust(method = "single")
# paquetería que permite visualizar resultados de análisis multivariados
library(factoextra)
# se grafican los conglomerados bajo cada argumento
fviz_dend(gexp_complete, main = "Complete")
Figure 3: Dendrograma de matriz de distancias para expresiones génicas
Enlace completo
fviz_dend(gexp_average, main = "Average")
Figure 4: Dendrograma de matriz de distancias para expresiones génicas
Enlace promedio
fviz_dend(gexp_single, main = "Single")
Figure 5: Dendrograma de matriz de distancias para expresiones génicas
Enlace simple
En un ejemplo más robusto, se cuenta con perfiles de expresión de miles de genes y normalmente se tendrán muchas más pacientes que en el ejemplo anterior. Se cuenta ahora con la expresión génica de \(60\) muestras de médula ósea de pacientes con uno de los cuatro tipos principales de leucemia (ALL, AML, CLL, CML) o controles sin leucemia. Se recortan en los \(1000\) genes más variables, ya que los genes que no son muy variables no contribuyen mucho a las distancias entre pacientes y se agrupaa a las pacientes y se mustran los valores en forma de un mapa de calor y de un dendrograma. Un mapa de calor muestra los valores de expresión de genes en las pacientes de una manera codificada por colores. La función de mapa de calor, pheatmap(), se usa también para realizar la agrupación. La matriz que contiene expresiones de genes tiene los genes en las filas y las pacientes en las columnas. Por lo tanto, también se usa un código de color en el lado de la columna para marcar a las pacientes según su tipo de leucemia. Para el agrupamiento jerárquico, se usa el método de Ward designado por el argumento clustering_method de la funcióin pheatmap(). El mapa de calor resultante se muestra en la Figura 6.
# se cargan las librerías necesarias
library(compGenomRData)
library(pheatmap)
# función para obtener y generar la ruta a los archivos que están instalados en la librería.
expFile <- system.file("extdata","leukemiaExpressionSubset.rds",
package="compGenomRData")
mat <- readRDS(expFile)
# establecer la anotación del tipo de leucemia para cada muestra
annotation_col <- data.frame(
LeukemiaType =substr(colnames(mat),1,3))
rownames(annotation_col)=colnames(mat)
# creación del mapa de calor
pheatmap(mat,show_rownames=FALSE,show_colnames=FALSE,
annotation_col=annotation_col,
scale = "none",clustering_method="ward.D2",
clustering_distance_cols="euclidean")
Figure 6: Mapa de calor de los valores de expresión génica de pacientes con leucemia
Cada columna representa a una paciente. Las columnas se agrupan utilizando la expresión génica y se codificadan por colores según el tipo de enfermedad: ALL, AML, CLL, CML o no leucemia
Como se observa del mapa de calor, cada conglomerado tiene un conjunto distinto de valores de expresión. Estos conglomerados principales distinguen casi perfectamente los tipos de leucemia. Solo una paciente con leucemia mieloide crónica se agrupa como muestra sin leucemia. Esto podría significar que los perfiles de expresión génica son suficientes para clasificar el tipo de leucemia.
Para un segundo ejemplo, se cuenta con la base de datos NCI60. Este es un conjunto de datos genómicos, que contiene datos de microarreglos de líneas de células cancerosas, que consta de \(6830\) mediciones de expresión génica en \(64\) líneas de células cancerosas. Los datos vienen como una lista que contiene una matriz y sus etiquetas, la cual se convierte en una tibble.
# se cargan las librerías necesarias
library(ISLR)
library(magrittr)
data(NCI60, package = "ISLR")
nci60 <- NCI60$data %>%
as_tibble() %>%
set_colnames(., paste0("V_", 1:ncol(.))) %>%
mutate(label = factor(NCI60$labs)) %>%
relocate(label)Viendo ahora qué sucede al realizar la agrupación en conglomerados en la base de datos nci60, primero, se crea una versión escalada usando una recipe y se elimina la columna de etiquetas para poder realizar la estandarización.
nci60_scaled <- recipe(~ ., data = nci60) %>%
step_rm(label) %>%
step_normalize(all_predictors()) %>%
prep() %>%
bake(new_data = NULL)Y se ajustan múltiples modelos de agrupamiento jerárquico usando diferentes métodos de aglomeración.
nci60_complete <- nci60_scaled %>%
dist() %>%
hclust(method = "complete")
nci60_average <- nci60_scaled %>%
dist() %>%
hclust(method = "average")
nci60_single <- nci60_scaled %>%
dist() %>%
hclust(method = "single")Luego se visualizan para checar si se generan algunas buenas separaciones naturales (Figuras 7, 8 y 9).
fviz_dend(nci60_complete, main = "Complete")
Figure 7: Dendrograma de matriz de distancias para líneas cancerígenas NCI60
fviz_dend(nci60_average, main = "Average")
Figure 8: Dendrograma de matriz de distancias para líneas cancerígenas NCI60
fviz_dend(nci60_single, main = "Single")
Figure 9: Dendrograma de matriz de distancias para líneas cancerígenas NCI60
3.1 ¿Dónde cortar el árbol?
Del ejemplo de tipos de leucemia parece claro que se puedan seleccionar agrupaciones del dendrograma a simple vista. Esto se debe principalmente al método de Ward, donde se prefieren los conglomerados compactos. Sin embargo, como suele ser el caso, no se tienen etiquetas de pacientes y sería difícil saber qué hojas (pacientes) del dendrograma se deben considerar como parte del mismo grupo. En otras palabras, qué tan profundo se debe cortar el dendrograma para que cada muestra de pacientes aún conectada a través de los subdendrogramas restantes constituya agrupaciones.
Con la codificación más limpia, también se puede colorer de acuerdo con el argumento k = 3 y se obtienen las siguientes separaciones (Figura 10 y Figura 11).
mat2 <- read_rds(expFile) %>%
as_tibble() %>%
set_colnames(., paste0("V_", 1:ncol(.)))
mat_scaled <- recipe(~ ., data = mat2) %>%
step_normalize(all_predictors()) %>%
prep() %>%
bake(new_data = NULL)
mat_complete <- mat_scaled %>%
dist() %>%
hclust(method = "complete")
fviz_dend(mat_complete, main = "Complete")
Figure 10: Dendrograma de matriz de distancias para pacientes con leucemia
mat_complete %>%
fviz_dend(k = 3, main = "hclust(complete) para tipos de leucemia")
Figure 11: Dendrograma de matriz de distancias para pacientes con leucemia
División por conglomerados
La función cutree() proporciona la funcionalidad para generar el número deseado de conglomerados o conglomerados obtenidos al cortar el dendrograma con una determinada cantidad de racimos.
clu.k3=cutree(mat_complete,k=3) # cortar el árbol para que haya 5 racimos
table(clu.k3) # número de muestras para cada conglomerado ## clu.k3
## 1 2 3
## 617 300 83A continuación, se agrupan a las expresiones génicas de NCI60 con agrupamiento jerárquico utilizando el método predeterminado complete linkage (Figura 12).
nci60_complete %>%
fviz_dend(k = 4, main = "hclust(complete) para nci60")
Figure 12: Dendrograma de matriz de distancias con 4 divisoines
Líneas cancerígenas NCI60
Con esto se toma la etiqueta de agrupamiento extrayendo con cutree y se calcula qué etiqueta es la más común dentro de cada grupo.
tibble(
label = nci60$label,
cluster_id = cutree(nci60_complete, k = 4)
) %>%
count(label, cluster_id) %>%
group_by(cluster_id) %>%
mutate(prop = n / sum(n)) %>%
slice_max(n = 1, order_by = prop) %>%
ungroup() ## # A tibble: 6 x 4
## label cluster_id n prop
## <fct> <int> <int> <dbl>
## 1 MELANOMA 1 8 0.2
## 2 NSCLC 1 8 0.2
## 3 RENAL 1 8 0.2
## 4 BREAST 2 3 0.429
## 5 LEUKEMIA 3 6 0.75
## 6 COLON 4 5 0.5564 Agrupación de k-medias
Otro algoritmo de agrupación en conglomerados muy común es k-medias. Este método divide o particiona los puntos de datos, las bases de datos de pacientes por ejemplo, en un número \(k\) predeterminado de agrupaciones. Por lo tanto, estos tipos de métodos generalmente se denominan métodos de partición. El algoritmo se inicializa eligiendo al azar \(k\) centros o centroides. En cierto sentido, un centroide es un punto de datos con múltiples valores. En la base de datos de leucemia, un centroide es una paciente hipotética con valores de expresión génica. Pero en la fase de inicialización, esos valores de expresión génica se eligen al azar dentro de los límites de las distribuciones de expresión génica de pacientes reales. Como siguiente paso en el algoritmo, a cada paciente se le asigna el centroide más cercano y, en la siguiente iteración, los centroides se establecen en la media de los valores de los genes en la agrupación. Este proceso de establecer centroides y asignar pacientes a las agrupaciones se repite hasta que se minimiza la suma de las distancias cuadradas a los centroides de la agrupación.
Este algoritmo de conglomerado comienza con centroides iniciales aleatorios. Esta característica puede producir resultados diferentes para cada ejecución del algoritmo.
Primero se aplica la algorítmica para la base de datos de expresión génica de pacientes con leucemia.
set.seed(101)
# se transpone la matrix para calcualr distancias
leu_clu=kmeans(t(mat),centers=5) Este modelo ajustado tiene muchos tipos diferentes de información, por lo que se pueden usar funciones de la librería broom para extraer información en formatos tidy. La función tidy() devuelve información para cada agrupación, incluida su posición, tamaño y suma de cuadrados dentro del grupo.
tidy(leu_clu)## # A tibble: 5 x 1,003
## ENSG00000224137 ENSG00000153253 ENSG00000096006 ENSG00000229807
## <dbl> <dbl> <dbl> <dbl>
## 1 3.64 7.20 5.41 5.43
## 2 4.07 3.35 11.3 4.51
## 3 3.57 3.12 11.9 4.24
## 4 10.3 3.19 3.48 5.36
## 5 3.54 4.25 5.17 6.45
## # … with 999 more variables: ENSG00000138772 <dbl>, ENSG00000169575 <dbl>,
## # ENSG00000102935 <dbl>, ENSG00000102837 <dbl>, ENSG00000173391 <dbl>,
## # ENSG00000154262 <dbl>, ENSG00000121742 <dbl>, ENSG00000107447 <dbl>,
## # ENSG00000124469 <dbl>, ENSG00000164330 <dbl>, ENSG00000196735 <dbl>,
## # ENSG00000132185 <dbl>, ENSG00000241911 <dbl>, ENSG00000196415 <dbl>,
## # ENSG00000118113 <dbl>, ENSG00000118523 <dbl>, ENSG00000170549 <dbl>,
## # ENSG00000156575 <dbl>, ENSG00000166349 <dbl>, ENSG00000133048 <dbl>,
## # ENSG00000172232 <dbl>, ENSG00000229314 <dbl>, ENSG00000090269 <dbl>,
## # ENSG00000095585 <dbl>, ENSG00000179869 <dbl>, ENSG00000231439 <dbl>,
## # ENSG00000242550 <dbl>, ENSG00000134827 <dbl>, ENSG00000104043 <dbl>,
## # ENSG00000118785 <dbl>, ENSG00000163751 <dbl>, ENSG00000111863 <dbl>,
## # ENSG00000164932 <dbl>, ENSG00000118520 <dbl>, ENSG00000163661 <dbl>,
## # ENSG00000198795 <dbl>, ENSG00000078399 <dbl>, ENSG00000130508 <dbl>,
## # ENSG00000125845 <dbl>, ENSG00000139329 <dbl>, ENSG00000012223 <dbl>,
## # ENSG00000007062 <dbl>, ENSG00000197632 <dbl>, ENSG00000124882 <dbl>,
## # ENSG00000197561 <dbl>, ENSG00000148346 <dbl>, ENSG00000117318 <dbl>,
## # ENSG00000115828 <dbl>, ENSG00000122585 <dbl>, ENSG00000198692 <dbl>,
## # ENSG00000145416 <dbl>, ENSG00000136929 <dbl>, ENSG00000143297 <dbl>,
## # ENSG00000162692 <dbl>, ENSG00000150625 <dbl>, ENSG00000086548 <dbl>,
## # ENSG00000099882 <dbl>, ENSG00000174944 <dbl>, ENSG00000139211 <dbl>,
## # ENSG00000228278 <dbl>, ENSG00000149516 <dbl>, ENSG00000047597 <dbl>,
## # ENSG00000133063 <dbl>, ENSG00000073756 <dbl>, ENSG00000164047 <dbl>,
## # ENSG00000198947 <dbl>, ENSG00000179087 <dbl>, ENSG00000106004 <dbl>,
## # ENSG00000169397 <dbl>, ENSG00000235688 <dbl>, ENSG00000177508 <dbl>,
## # ENSG00000163710 <dbl>, ENSG00000005381 <dbl>, ENSG00000170365 <dbl>,
## # ENSG00000163421 <dbl>, ENSG00000162747 <dbl>, ENSG00000138449 <dbl>,
## # ENSG00000100448 <dbl>, ENSG00000164821 <dbl>, ENSG00000134198 <dbl>,
## # ENSG00000110777 <dbl>, ENSG00000133818 <dbl>, ENSG00000143416 <dbl>,
## # ENSG00000124490 <dbl>, ENSG00000159339 <dbl>, ENSG00000127955 <dbl>,
## # ENSG00000135218 <dbl>, ENSG00000117586 <dbl>, ENSG00000100721 <dbl>,
## # ENSG00000111261 <dbl>, ENSG00000118946 <dbl>, ENSG00000127990 <dbl>,
## # ENSG00000181631 <dbl>, ENSG00000120833 <dbl>, ENSG00000156738 <dbl>,
## # ENSG00000104921 <dbl>, ENSG00000186431 <dbl>, ENSG00000129682 <dbl>,
## # ENSG00000136573 <dbl>, ENSG00000119888 <dbl>, …La función glance() devuelve métricas basadas en modelos. Uno de ellos es tot.withinss, que es la suma de cuadrados total dentro del conglomerado que se busca minimizar al realizar el conglomerado de k-medias.
glance(leu_clu)## # A tibble: 1 x 4
## totss tot.withinss betweenss iter
## <dbl> <dbl> <dbl> <int>
## 1 168850. 64666. 104184. 3Por último, se puede ver a qué agrupación pertenece cada observación mediante el uso de augment() que predice a qué conglomerado pertenece una determinada observación.
augment(leu_clu, data = t(mat)) %>%
relocate(.cluster)## # A tibble: 60 x 1,002
## .cluster .rownames ENSG00000224137 ENSG00000153253 ENSG00000096006
## <fct> <chr> <dbl> <dbl> <dbl>
## 1 1 ALL_GSM330151.CEL 5.33 6.44 6.94
## 2 1 ALL_GSM330153.CEL 3.51 9.56 8.84
## 3 1 ALL_GSM330154.CEL 3.46 7.19 4.60
## 4 1 ALL_GSM330157.CEL 3.47 2.95 4.74
## 5 1 ALL_GSM330171.CEL 3.64 6.95 4.18
## 6 1 ALL_GSM330174.CEL 3.39 9.10 3.05
## 7 1 ALL_GSM330178.CEL 3.28 3.80 9.74
## 8 1 ALL_GSM330182.CEL 3.68 12.0 3.34
## 9 1 ALL_GSM330185.CEL 3.41 3.22 7.30
## 10 1 ALL_GSM330186.CEL 3.60 6.74 3.70
## # … with 50 more rows, and 997 more variables: ENSG00000229807 <dbl>,
## # ENSG00000138772 <dbl>, ENSG00000169575 <dbl>, ENSG00000102935 <dbl>,
## # ENSG00000102837 <dbl>, ENSG00000173391 <dbl>, ENSG00000154262 <dbl>,
## # ENSG00000121742 <dbl>, ENSG00000107447 <dbl>, ENSG00000124469 <dbl>,
## # ENSG00000164330 <dbl>, ENSG00000196735 <dbl>, ENSG00000132185 <dbl>,
## # ENSG00000241911 <dbl>, ENSG00000196415 <dbl>, ENSG00000118113 <dbl>,
## # ENSG00000118523 <dbl>, ENSG00000170549 <dbl>, ENSG00000156575 <dbl>,
## # ENSG00000166349 <dbl>, ENSG00000133048 <dbl>, ENSG00000172232 <dbl>,
## # ENSG00000229314 <dbl>, ENSG00000090269 <dbl>, ENSG00000095585 <dbl>,
## # ENSG00000179869 <dbl>, ENSG00000231439 <dbl>, ENSG00000242550 <dbl>,
## # ENSG00000134827 <dbl>, ENSG00000104043 <dbl>, ENSG00000118785 <dbl>,
## # ENSG00000163751 <dbl>, ENSG00000111863 <dbl>, ENSG00000164932 <dbl>,
## # ENSG00000118520 <dbl>, ENSG00000163661 <dbl>, ENSG00000198795 <dbl>,
## # ENSG00000078399 <dbl>, ENSG00000130508 <dbl>, ENSG00000125845 <dbl>,
## # ENSG00000139329 <dbl>, ENSG00000012223 <dbl>, ENSG00000007062 <dbl>,
## # ENSG00000197632 <dbl>, ENSG00000124882 <dbl>, ENSG00000197561 <dbl>,
## # ENSG00000148346 <dbl>, ENSG00000117318 <dbl>, ENSG00000115828 <dbl>,
## # ENSG00000122585 <dbl>, ENSG00000198692 <dbl>, ENSG00000145416 <dbl>,
## # ENSG00000136929 <dbl>, ENSG00000143297 <dbl>, ENSG00000162692 <dbl>,
## # ENSG00000150625 <dbl>, ENSG00000086548 <dbl>, ENSG00000099882 <dbl>,
## # ENSG00000174944 <dbl>, ENSG00000139211 <dbl>, ENSG00000228278 <dbl>,
## # ENSG00000149516 <dbl>, ENSG00000047597 <dbl>, ENSG00000133063 <dbl>,
## # ENSG00000073756 <dbl>, ENSG00000164047 <dbl>, ENSG00000198947 <dbl>,
## # ENSG00000179087 <dbl>, ENSG00000106004 <dbl>, ENSG00000169397 <dbl>,
## # ENSG00000235688 <dbl>, ENSG00000177508 <dbl>, ENSG00000163710 <dbl>,
## # ENSG00000005381 <dbl>, ENSG00000170365 <dbl>, ENSG00000163421 <dbl>,
## # ENSG00000162747 <dbl>, ENSG00000138449 <dbl>, ENSG00000100448 <dbl>,
## # ENSG00000164821 <dbl>, ENSG00000134198 <dbl>, ENSG00000110777 <dbl>,
## # ENSG00000133818 <dbl>, ENSG00000143416 <dbl>, ENSG00000124490 <dbl>,
## # ENSG00000159339 <dbl>, ENSG00000127955 <dbl>, ENSG00000135218 <dbl>,
## # ENSG00000117586 <dbl>, ENSG00000100721 <dbl>, ENSG00000111261 <dbl>,
## # ENSG00000118946 <dbl>, ENSG00000127990 <dbl>, ENSG00000181631 <dbl>,
## # ENSG00000120833 <dbl>, ENSG00000156738 <dbl>, ENSG00000104921 <dbl>,
## # ENSG00000186431 <dbl>, ENSG00000129682 <dbl>, ENSG00000136573 <dbl>, …Y se puede visualizar el resultado de augment() para ver qué tan bien funcionó la agrupación en las dos primeras expresiones (Figura 13).
augment(leu_clu, data = t(mat)) %>%
ggplot(aes(ENSG00000224137, ENSG00000153253, color = .cluster)) +
geom_point()
Figure 13: Gráfica de dispersión de expresiones génicas ENSG00000224137 y ENSG00000153253
Se pueden distinguir dos conglomerados claramente, pero sería bueno si se pudieran probar varias agrupaciones diferentes y luego encontrar la mejor. Con un combo de mutate() y map() se ajustan múltiples modelos y se extrae información de ellos.
set.seed(1234)
multi_kmeans <- tibble(k = 1:10) %>%
mutate(
model = purrr::map(k, ~ kmeans(t(mat), centers = .x, nstart = 20)),
tot.withinss = purrr::map_dbl(model, ~ glance(.x)$tot.withinss)
)
multi_kmeans## # A tibble: 10 x 3
## k model tot.withinss
## <int> <list> <dbl>
## 1 1 <kmeans> 168850.
## 2 2 <kmeans> 105826.
## 3 3 <kmeans> 82408.
## 4 4 <kmeans> 69888.
## 5 5 <kmeans> 64666.
## 6 6 <kmeans> 59921.
## 7 7 <kmeans> 56262.
## 8 8 <kmeans> 52923.
## 9 9 <kmeans> 50549.
## 10 10 <kmeans> 48784.Ahora que se cuenta con la suma de cuadrados total dentro del conglomerado, se grafican contra \(k\) en la Figura 14, de modo que se pueda usar el método del codo para encontrar el número óptimo de conglomerados.
multi_kmeans %>%
ggplot(aes(k, tot.withinss)) +
geom_point() +
geom_line()
Figure 14: Gráfica de codo para expresiones génicas ENSG00000224137 y ENSG00000153253
Se observa un codo en \(k = 3\), lo que indica la creación específica para tener \(3\) agrupaciones. Ahora se extrae el modelo donde \(k = 3\) de multi_kmeans.
final_kmeans <- multi_kmeans %>%
filter(k == 3) %>%
pull(model) %>%
pluck(1)Y se finaliza visualizando en la Figura 15 los conglomerados encontrados.
augment(final_kmeans, data = t(mat)) %>%
ggplot(aes(ENSG00000224137, ENSG00000153253, color = .cluster)) +
geom_point()
Figure 15: Gráfica de dispersión para expresiones génicas ENSG00000224137 y ENSG00000153253
Con 3 conglomerados óptimos
También se ajusta un clúster de k-medias para las expresiones de NCI60. Se seguirán empleando las \(4\) agrupaciones encontradas por aglomeracioón jerárquica.
set.seed(2)
res_kmeans_scaled <- kmeans(nci60_scaled, centers = 4, nstart = 50)Y se observa la información usando los formatos tidy de broom.
tidy(res_kmeans_scaled) %>%
select(cluster, size, withinss)## # A tibble: 4 x 3
## cluster size withinss
## <fct> <int> <dbl>
## 1 1 20 108801.
## 2 2 27 154545.
## 3 3 9 37150.
## 4 4 8 44071.Se pueden comparar k_medias y agrupación jerárquica usando conf_mat() y un mapa de calor (Figura 16).
cluster_kmeans <- res_kmeans_scaled$cluster
cluster_hclust <- cutree(nci60_complete, k = 4)
tibble(
kmeans = factor(cluster_kmeans),
hclust = factor(cluster_hclust)
) %>%
conf_mat(kmeans, hclust) %>%
autoplot(type = "heatmap")
Figure 16: Mapa de calor para matriz de confusión
No hay mucha congruencia entre las etiquetas, lo cual tiene sentido, ya que las etiquetas en sí se agregan arbitrariamente. Lo importante es que tienden a estar bastante de acuerdo al ser matriz de confusión dispersa.
5 Cómo elegir “k”, el número de conglomerados
Hasta este punto, se ha evitado la cuestión de seleccionar conglomerados de números óptimos. ¿Cómo se puede saber dónde cortar un dendrograma o qué \(k\) elegir? En primer lugar, esta es una pregunta difícil. Por lo general, los conglomerados tienen una granularidad diferente. Algunas agrupaciones son estrechas y compactas y algunas son anchas, y ambos tipos de agrupaciones pueden estar en la misma base de datos. Cuando se visualizan, algunas agrupaciones grandes pueden parecer que pueden tener subagrupaciones. Entonces, ¿se debe considerar la agrupación grande como un conglomerado o se deben considerar las subagrupaciones como conglomerados individuales? Hay algunas métricas que pueden ayudar, pero no hay una respuesta definitiva.
5.1 Silueta
Una forma de determinar la calidad de la agrupación es medir la naturaleza auto-similar esperada de los puntos en un conjunto de agrupaciones. El valor de silueta hace precisamente eso y es una medida de cuán similar es un punto de datos a su propia agrupación en comparación con otras agrupaciones. El valor de la silueta varía de \(-1\) a \(+1\), donde los valores que son positivos indican que el punto de datos está bien emparejado con su propio conglomerado, si el valor es cero es un caso límite, y si el valor es negativo significa que el el punto de datos puede estar mal agrupado porque es más similar a un conglomerado vecino. Si la mayoría de los puntos de datos tienen un valor alto, entonces el agrupamiento es apropiado. Idealmente, se pueden crear muchas agrupaciones diferentes, cada una con una \(k\) diferente que indica el número de conglomerados, y evaluar su idoneidad utilizando los valores medios de silueta. En las implementaciones en R, los valores de silueta se denominan silhouette widths.
Se calcula un valor de silueta para cada punto de la base de datos. En el ejemplo de pacientes con leucemia, cada paciente obtendrá valores de silueta que muestran qué tan bien se corresponden con sus agrupaciones asignadas. Formalmente, esto se calcula de la siguiente manera. Para cada punto de datos \(i\), se calcula \(a(i)\), que denota la distancia promedio entre \(i\) y todos los demás puntos de datos dentro del mismo conglomerado. Esto muestra qué tan bien encaja el punto en esa agrupación. Para el mismo punto de datos, también se calcula \(b(i)\), que denota la distancia promedio más baja de \(i\) a todos los puntos en cualquier otra conglomerado, del cual \(i\) no es miembro. El conglomerado con este promedio más bajo \(b(i)\) es el conglomerado vecino del punto de datos \(i\), ya que es el siguiente conglomerado que mejor se ajusta a ese punto de datos. Entonces, el valor de silueta para un punto de datos dado es \(s(i) = \frac{b(i) − a(i)}{max\{a(i), b(i)\}}\).
Como se describe, esta cantidad es positiva cuando \(b(i)\) es alto y \(a(i)\) es bajo, lo que significa que el punto de datos \(i\) es auto-similar a su conglomerado. Y el valor de silueta, \(s(i)\), es negativo si es más similar a sus vecinos que a su conglomerado asignado.
En R, se calculan valores de silueta usando la función cluster::silhouette(). El siguiente código gráfica (Figura 17) y calcula este enfoque para \(1\) a \(15\) conglomerados en la base de datos para pacientes con leucemia.
# function para calcualr siluetas promedio para k clusters
library(cluster)
avg_sil <- function(k) {
km.res <- kmeans(t(mat), centers = k, nstart = 25)
ss <- cluster::silhouette(km.res$cluster, dist(t(mat)))
mean(ss[, 3])
}
# calcular y graficar wss para k = 2 a k = 15
k.values <- 2:15
# extraer silueta promedio para 2-15 clusters
avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Número de conglomerados k",
ylab = "Siluetas Promedio")
Figure 17: Gráfica de codo para el método de silueta promedio
Similar al método del codo, la Figura 18 muestra este proceso para calcular el método de silueta promedio, envuelto en una sola función fviz_nbclust().
fviz_nbclust(t(mat), kmeans, method = "silhouette")
Figure 18: Gráfica de codo para el método de silueta promedio
Con 2 conglomerados óptimos.
En este caso, parece que el mejor valor para \(k\) es \(2\).
6 Estadística de brecha (gap)
Este enfoque se puede aplicar a cualquier método de agrupación (i.e., k-medias, agrupación jerárquica). La estadística de brecha compara la variación intragrupo total para diferentes valores de \(k\) con sus valores esperados bajo una distribución de referencia nula de los datos (es decir, una distribución sin agrupamiento obvio). El conjunto de datos de referencia se genera utilizando simulaciones de Monte Carlo del proceso de muestreo. Es decir, para cada variable \(x_i\) en el conjunto de datos se calcula su rango \([\min(x_i), \max(x_j)]\) y se generan valores para los \(n\) puntos uniformemente desde el intervalo mínimo al máximo.
Para los datos observados y los datos de referencia, la variación intragrupo total se calcula utilizando diferentes valores de \(k\). La estadística de brecha para una \(k\) dada se define de la siguiente manera:
\[\begin{equation} Gap_n(k)=E_n^*\log(W_k)−\log(W_k). \tag{1} \end{equation}\]Donde \(E_n^*\) denota la esperanza bajo un tamaño de muestra \(n\) de la distribución de referencia. \(E_n^*\) se define mediante bootstrapping (\(B\)) generando \(B\) copias de los conjuntos de datos de referencia y calculando el logaritmo promedio \(\log W_k^*\). La estadística de brecha mide la desviación del valor \(W_k\) observado de su valor esperado bajo la hipótesis nula. La estimación de los conglomerados óptimos \(\widehat{k}\) será el valor que maximice \(Gap_n(k)\). Esto significa que la estructura de agrupamiento está lejos de la distribución uniforme de puntos.
En resumen, el algoritmo implica los siguientes pasos:
Agrupar los datos observados, variando el número de grupos de \(k = 1,..., k_{max}\), y calcular la \(W_k\) correspondiente.
Generar \(B\) conjuntos de datos de referencia y agrupar cada uno de ellos con un número variable de grupos \(k = 1,..., k_{max}\). Calcular las estadísticas de brecha estimadas que se presentan en la ecuación (1).
Sea \(\overline{w} = (1 / B) ∑_b\log(W_{k_b}^*)\), calcular la desviación estándar \(sd (k) = \sqrt{(1 / b) ∑_b [\log (W^∗_{k_b}) −\overline{w}]^2}\) y se define \(s_k = sd_k \times \sqrt{1 + 1 / B}\).
Elegir el número de conglomerados como la más pequeña \(k\) tal que \(Gap (k) ≥Gap (k + 1) − s_{k + 1}\).
Para calcular la estadística de brecha, se puede usar la función clusGap() que proporciona la brecha y el error estándar para una salida.
# calcular estadística gap
set.seed(123)
gap_stat <- clusGap(t(mat), FUN = kmeans, nstart = 25,
K.max = 10, B = 50)
# imprimir resultados
print(gap_stat, method = "firstmax")## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = t(mat), FUNcluster = kmeans, K.max = 10, B = 50, nstart = 25)
## B=50 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstmax'): 10
## logW E.logW gap SE.sim
## [1,] 6.979442 7.204019 0.2245772 0.01665346
## [2,] 6.727785 7.107497 0.3797121 0.01275230
## [3,] 6.600435 7.056767 0.4563317 0.01166583
## [4,] 6.517912 7.015651 0.4977393 0.01183251
## [5,] 6.460466 6.979839 0.5193728 0.01116938
## [6,] 6.415094 6.946476 0.5313826 0.01057627
## [7,] 6.377243 6.915169 0.5379258 0.01051512
## [8,] 6.339501 6.884625 0.5451243 0.01033359
## [9,] 6.305496 6.855133 0.5496373 0.01058301
## [10,] 6.269824 6.825824 0.5560001 0.01066139Se pueden visualizar los resultados con fviz_gap_stat() (Figura 19).
fviz_gap_stat(gap_stat)
Figure 19: Gráfica de codo para conglomerados óptimos en la estadística de brecha
7 Comentarios adicionales
La agrupación en conglomerados de k-medias es un algoritmo muy simple y rápido. Además, puede tratar de manera eficiente conjuntos de datos muy grandes. Sin embargo, existen algunas debilidades del enfoque de k-medias.
Una posible desventaja de la agrupación de k-medias es que requiere especificar previamente el número de conglomerados. La agrupación jerárquica es un enfoque alternativo que no requiere una elección particular de agrupaciones. El agrupamiento jerárquico tiene una ventaja adicional sobre el agrupamiento de k-medias en que da como resultado una representación atractiva basada en árboles de las observaciones, llamado dendrograma.
Una desventaja adicional de k-medias es que es sensible a valores atípicos y pueden producirse resultados diferentes si cambia el orden de los datos.
Además de estos enfoques de uso común, el paquete NbClust, proporciona \(30\) índices para determinar el número relevante de conglomerados y propone a los usuarios el mejor esquema de agrupamiento a partir de los diferentes resultados obtenidos al variar todas las combinaciones del número de conglomerados, medidas de distancia y métodos de agrupamiento.
Sin embargo, se debe tener en cuenta que la agrupación en conglomerados es una técnica exploratoria. Si se tienen etiquetas sólidas para los puntos de datos, tal vez la agrupación por conglomerados sea sólo una herramienta de verificación, y en su lugar, se debería hacer un modelo predictivo. Para las aplicaiones en biología rara vez hay etiquetas sólidas y las bases de datos tienen diferente granularidad. Como en el caso de las pacientes con leucemia que se estuvo usando, se sabe que los tipos de leucemia tienen subtipos y que los subtipos ienen diferentes perfiles de mutación y diferentes firmas moleculares. Debido a esto, algunas técnicas óptimas de números de conglomerados se encuentran con que más conglomerados son apropiados. Por otro lado, la leucemia mieloide crónica es una enfermedad de progresión lenta y tal vez sus firmas moleculares estén más cerca de los pacientes sin leucemia, por lo que los algoritmos de agrupación en conglomerados pueden confundirse según la granularidad con la que estén operando. Siempre es bueno mirar los mapas de calor después de la agrupación, si se tienen puntos de datos auto-similares significativos, incluso si las etiquetas que se tienen no están de acuerdo en que puede haber diferentes agrupaciones, se pueden realizar análisis posteriores para comprender mejor las subagrupaciones. Como se ha visto, se puede estimar el número óptimo de conglomerados, pero no se puede tomar esa estimación como la verdad absoluta. Dados más puntos de datos o un conjunto diferente de firmas de expresión, es posible que se produzcan diferentes agrupaciones óptimas, o la supuesta agrupación óptima podría pasar por alto subagrupaciones previamente conocidas en las bases de datos originales.
8 Referencias
Akalin A. Computational Genomics with R. CRC Press, 2021: https://compgenomr.github.io/book/
ISLR tidymodels Labs: https://emilhvitfeldt.github.io/ISLR-tidymodels-labs/index.html
UC Business Analytics R Programming Guide, K-means Cluster Analysis: https://uc-r.github.io/kmeans_clustering
Hartigan, and Wong. 1979. “Algorithm as 136: A K-Means Clustering Algorithm.” Journal of the Royal Statistical Society. Series C (Applied Statistics) 28 (1): 100–108.
Tidymodels. K-means clustering with tidy data principles: https://www.tidymodels.org/learn/statistics/k-means/
Tibshirani, R. et al. 2001. “Estimating the number of clusters ina data set via the gap statistic.” Journal of the Royal Statistical Society. Series B (Statistical Methodology) 63 (2): 411-423.