Informe 4
Daniel Celis
2024-05-17
4.K-Means Clustering
La metodología de agrupación mediante K-medias, formulada por MacQueen en 1967, sobresale como el principal algoritmo de aprendizaje no supervisado empleado para dividir un conjunto específico de datos en k grupos (o clústeres), donde k es establecido de antemano por el analista. Su función primordial consiste en la categorización de elementos en diversos grupos, con el objetivo de maximizar la semejanza intraclase (alta concordancia entre elementos dentro del mismo grupo) y minimizar la semejanza interclase (baja concordancia entre elementos de grupos diferentes). Durante este procedimiento, cada grupo se caracteriza a través de su centroide, que representa la media de los puntos asignados a ese grupo.
4.3 Computing k-means clustering in R
Para llevar a cabo la agrupación de k-medias en R, se empleará el conjunto de datos de muestra “USArrests”. Es esencial preparar los datos según las indicaciones del capítulo 2, asegurándose de que solo incluyan variables continuas, ya que el algoritmo k-medias se basa en promedios variables.
data("USArrests") #loading the data set
df<- scale(USArrests) # Scaling the data
#view the firt 3 rows of the data
head(df,n=3)## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388#Required R packages and functions
kmeans(10, 1, iter.max=10, nstart=1)## K-means clustering with 1 clusters of sizes 1
##
## Cluster means:
## [,1]
## 1 10
##
## Clustering vector:
## [1] 1
##
## Within cluster sum of squares by cluster:
## [1] 0
## (between_SS / total_SS = NaN %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"library(factoextra)## Cargando paquete requerido: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBafviz_nbclust(df,kmeans,method="wss") +
geom_vline(xintercept=4,linetype=2)
4.3.4 Computing k-means clustering
Como el algoritmo de agrupamiento de k-medias comienza con k centroides seleccionados aleatoriamente, siempre se recomienda usar la función set.seed() para establecer una semilla para el algoritmo aleatorio de R.
generador de números. El objetivo es hacer reproducibles los resultados, de modo que el lector de este artículo obtenga exactamente los mismos resultados que los que se muestran a continuación.
#Compute k-means with k=4
set.seed(123)
km.res<- kmeans(df, 4, nstart=25)
#Print the results
print(km.res)## K-means clustering with 4 clusters of sizes 8, 13, 16, 13
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 1.4118898 0.8743346 -0.8145211 0.01927104
## 2 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 3 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 4 0.6950701 1.0394414 0.7226370 1.27693964
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 4 4 1 4
## Colorado Connecticut Delaware Florida Georgia
## 4 3 3 4 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 2 4 3 2
## Kansas Kentucky Louisiana Maine Maryland
## 3 2 1 2 4
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 4 2 1 4
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 4 2 3
## New Mexico New York North Carolina North Dakota Ohio
## 4 4 1 2 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 4 3 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 2 2 3
##
## Within cluster sum of squares by cluster:
## [1] 8.316061 11.952463 16.212213 19.922437
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"aggregate(USArrests,by=list(cluster=km.res$cluster),mean)## cluster Murder Assault UrbanPop Rape
## 1 1 13.93750 243.62500 53.75000 21.41250
## 2 2 3.60000 78.53846 52.07692 12.17692
## 3 3 5.65625 138.87500 73.87500 18.78125
## 4 4 10.81538 257.38462 76.00000 33.19231dd<- cbind(USArrests, cluster=km.res$cluster)
head(dd)## Murder Assault UrbanPop Rape cluster
## Alabama 13.2 236 58 21.2 1
## Alaska 10.0 263 48 44.5 4
## Arizona 8.1 294 80 31.0 4
## Arkansas 8.8 190 50 19.5 1
## California 9.0 276 91 40.6 4
## Colorado 7.9 204 78 38.7 44.3.5 Accesing to the results of kmeans() function
La función kmeans() devuelve una lista de componentes, que incluye:
• cluster: Un vector de números enteros (de 1:k) que indica el cluster al que se asigna cada punto
• Centers: una matriz de centros de conglomerados (medios de conglomerados)
• totss: La suma total de cuadrados (TSS), es decir (x,-)². TSS mide la varianza total de los datos. • insidess: Vector de suma de cuadrados dentro del grupo, un componente por grupo
• tot.withinss: suma total de cuadrados dentro del grupo, es decir, sum(withinss)
•Betweenss: suma de cuadrados entre grupos, es decir, totss-tot.withinss
•size: el número de observaciones en cada grupo
#Cluster number for each of the observations
km.res$cluster## Alabama Alaska Arizona Arkansas California
## 1 4 4 1 4
## Colorado Connecticut Delaware Florida Georgia
## 4 3 3 4 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 2 4 3 2
## Kansas Kentucky Louisiana Maine Maryland
## 3 2 1 2 4
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 4 2 1 4
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 4 2 3
## New Mexico New York North Carolina North Dakota Ohio
## 4 4 1 2 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 4 3 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 2 2 3#Cluster means
km.res$centers## Murder Assault UrbanPop Rape
## 1 1.4118898 0.8743346 -0.8145211 0.01927104
## 2 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 3 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 4 0.6950701 1.0394414 0.7226370 1.276939644.3.6 Visualizing k-means clusters
Visualizar los resultados de la agrupación de k-medias es esencial para evaluar el número óptimo de conglomerados y comparar análisis de agrupación. La representación gráfica se realiza mediante un diagrama de dispersión, coloreando cada punto según su asignación de grupo.
fviz_cluster(km.res, data=df,
palette=c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
ellipse.type = "euclid", #Concentration ellipse
star.plot=TRUE, #Add segments form centroids to items
repel=TRUE, #Avoid label overplotting(slow)
ggtheme=theme_minimal()
)
CHAPTER 5
5.3 Computing PAM in R
LEl algoritmo k-medoids, relacionado con la agrupación k-medias (véase capítulo 4), divide un conjunto de datos en k grupos o clusters. En este enfoque, cada grupo está representado por un punto de datos específico, denominado medoide del cúmulo. El término “medoide” se refiere a un objeto dentro del grupo cuya disimilitud promedio con los demás miembros del grupo es mínima, siendo este el punto más céntrico del cluster. A diferencia de la agrupación k-medias, donde el centro del grupo es el valor medio de todos los puntos, los medoides son objetos representativos de cada grupo.
5.3.1 Data
Usaremos los conjuntos de datos de demostración “USArrests”, que comenzamos escalando (Capítulo 2) usando la función R scale() de la siguiente manera:
data("USArrests") #Load the data set
df<- scale(USArrests) #Scale the data
head(df,n=3) #View the first 3 rows of the data## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.0428783885.3.2 Required R packages and functions
library(cluster)
pam(df,4 , metric="euclidean",stand=FALSE)## Medoids:
## ID Murder Assault UrbanPop Rape
## Alabama 1 1.2425641 0.7828393 -0.5209066 -0.003416473
## Michigan 22 0.9900104 1.0108275 0.5844655 1.480613993
## Oklahoma 36 -0.2727580 -0.2371077 0.1699510 -0.131534211
## New Hampshire 29 -1.3059321 -1.3650491 -0.6590781 -1.252564419
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 2 2 1 2
## Colorado Connecticut Delaware Florida Georgia
## 2 3 3 2 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 4 2 3 4
## Kansas Kentucky Louisiana Maine Maryland
## 3 3 1 4 2
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 2 4 1 3
## Montana Nebraska Nevada New Hampshire New Jersey
## 3 3 2 4 3
## New Mexico New York North Carolina North Dakota Ohio
## 2 2 1 4 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 4 1 2 3 4
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 4 4 3
## Objective function:
## build swap
## 1.035116 1.027102
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"5.3.3 Estimating the optimal number of clusters
Para determinar el número óptimo de grupos, empleamos el método de silueta promedio en la estimación del algoritmo PAM. Este método implica calcular el PAM con distintos valores de grupos, seguido por la representación gráfica de la silueta promedio en función del número de conglomerados. La silueta media es un indicador de la calidad del agrupamiento, donde un ancho de silueta promedio elevado denota una agrupación efectiva. La identificación del número óptimo de conglomerados, k, se basa en maximizar la silueta promedio dentro de un rango de valores posibles para k, siguiendo la propuesta de Kaufman y Rousseeuw (1990). La función R fviz_nbclust() del paquete factoextra brinda una solución práctica para esta estimación.
library(cluster)
library(factoextra)
fviz_nbclust(df,pam, method="silhouette")+
theme_classic()
5.3.4 Computing PAM clustering
El siguiente código R calcula el algoritmo PAM con k=2.
pam.res <- pam(df,2)
print(pam.res)## Medoids:
## ID Murder Assault UrbanPop Rape
## New Mexico 31 0.8292944 1.3708088 0.3081225 1.1603196
## Nebraska 27 -0.8008247 -0.8250772 -0.2445636 -0.5052109
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 1 1 2 1
## Colorado Connecticut Delaware Florida Georgia
## 1 2 2 1 1
## Hawaii Idaho Illinois Indiana Iowa
## 2 2 1 2 2
## Kansas Kentucky Louisiana Maine Maryland
## 2 2 1 2 1
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 1 2 1 1
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 1 2 2
## New Mexico New York North Carolina North Dakota Ohio
## 1 1 1 2 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 1 2 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 2 2 2
## Objective function:
## build swap
## 1.441358 1.368969
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"dd <- cbind(USArrests, cluster=pam.res$cluster)
head(dd, n=3)## Murder Assault UrbanPop Rape cluster
## Alabama 13.2 236 58 21.2 1
## Alaska 10.0 263 48 44.5 1
## Arizona 8.1 294 80 31.0 15.3.5 Accessing to the results of the pam() function
La función pam() produce un objeto de clase “pam”, cuyos componentes fundamentales son la lista de medoides que representan los grupos y el vector de agrupamiento que asigna cada objeto a su respectivo grupo. Estos resultados brindan una visión detallada de la estructura de los grupos y la asignación de objetos en el conjunto de datos analizado.
#Cluster medoids:New Mexico, Nebraska
pam.res$medoids## Murder Assault UrbanPop Rape
## New Mexico 0.8292944 1.3708088 0.3081225 1.1603196
## Nebraska -0.8008247 -0.8250772 -0.2445636 -0.5052109#Cluster numbers
head(pam.res$clustering)## Alabama Alaska Arizona Arkansas California Colorado
## 1 1 1 2 1 15.3.6 Visualizing PAM clusters
Para representar gráficamente los resultados de la partición en clústeres PAM, emplearemos la función fviz_cluster() del paquete factoextra. Esta función genera un diagrama de dispersión que coloreará los puntos de datos según los números de grupo asignados. En situaciones donde los datos abarcan más de dos variables, se implementa el algoritmo de Análisis de Componentes Principales (PCA) para reducir la dimensionalidad de los datos. En este contexto, las dos primeras dimensiones principales se utilizan para visualizar los datos de manera efectiva.
fviz_cluster(pam.res,
palette=c("#00AFBB","#FC4E07"), #color palette
ellipse.type="t", #concentration ellipse
repel=TRUE, #Avoid label overplotting(slow)
ggtheme=theme_classic()
)
Chapter 6 CLARA-Clustering large application
CLARA (Clustering Large Applications, Kaufman y Rousseeuw, 1990) constituye una extensión de los métodos k-medoids, como se detalla en el Capítulo 5, diseñada para abordar conjuntos de datos que implican una cantidad considerable de objetos, superando varias miles de observaciones. El objetivo principal es mitigar la carga computacional y los desafíos de almacenamiento de la RAM. Para lograrlo, CLARA implementa un enfoque de muestreo.
6.3 Computing CLARA in R
6.3.1 Data format and preparation
Para calcular el algoritmo CLARA en R, los datos deben prepararse como se indica en el Capítulo 2. Aquí, generaremos un conjunto de datos aleatorios. Para que el resultado sea reproducible, comenzamos usando la función set.seed().
set.seed(1234)
#Generate 500 objects,divided into 2 clusters
df <- rbind(cbind(rnorm(200,0,8), rnorm(200,0,8)),
cbind(rnorm(300,50,8),rnorm(300,50,8)))
#Specify column and rows names
colnames(df) <- c("x","y")
rownames(df) <- paste0("S",1:nrow(df))
#Previewing the data
head(df,nrow=6)## x y
## S1 -9.656526 3.881815
## S2 2.219434 5.574150
## S3 8.675529 1.484111
## S4 -18.765582 5.605868
## S5 3.432998 2.493448
## S6 4.048447 6.0836996.3.2 Required R packages and functions
Para estimar el número óptimo de clusters en sus datos, es posible utilizar el método de silueta promedio como se describe en el capítulo de clustering de PAM (Capítulo 5). La función R fuiz_nbclust() [paquete factoextra] proporciona una solución para facilitar este paso.
library(factoextra)
library(cluster)
clara(df,2,metric="euclidean", stand=FALSE,
samples=5, pamLike=FALSE)## Call: clara(x = df, k = 2, metric = "euclidean", stand = FALSE, samples = 5, pamLike = FALSE)
## Medoids:
## x y
## S168 -0.5495492 2.458514
## S397 47.5964047 50.892735
## Objective function: 9.9971
## Clustering vector: Named int [1:500] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "names")= chr [1:500] "S1" "S2" "S3" "S4" "S5" "S6" "S7" ...
## Cluster sizes: 200 300
## Best sample:
## [1] S6 S45 S51 S67 S75 S85 S90 S94 S97 S110 S111 S160 S168 S170 S176
## [16] S181 S201 S219 S249 S260 S264 S275 S296 S304 S317 S319 S337 S361 S362 S369
## [31] S370 S374 S379 S397 S398 S411 S420 S422 S424 S436 S448 S458 S465 S489
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"fviz_nbclust(df, clara, method="silhouette") +
theme_classic()
6.3.4 Computing CLARA
El siguiente código R calcula el algoritmo PAM con k=2.
#Compute CLARA
clara.res <- clara(df,2, samples=50, pamLike=TRUE)
#Print components of clara.res
print(clara.res)## Call: clara(x = df, k = 2, samples = 50, pamLike = TRUE)
## Medoids:
## x y
## S121 -1.531137 1.145057
## S455 48.357304 50.233499
## Objective function: 9.87862
## Clustering vector: Named int [1:500] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "names")= chr [1:500] "S1" "S2" "S3" "S4" "S5" "S6" "S7" ...
## Cluster sizes: 200 300
## Best sample:
## [1] S37 S49 S54 S63 S68 S71 S76 S80 S82 S101 S103 S108 S109 S118 S121
## [16] S128 S132 S138 S144 S162 S203 S210 S216 S231 S234 S249 S260 S261 S286 S299
## [31] S304 S305 S312 S315 S322 S350 S403 S450 S454 S455 S456 S465 S488 S497
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"dd <- cbind(df,cluster=clara.res$cluster)
head(dd,n=4)## x y cluster
## S1 -9.656526 3.881815 1
## S2 2.219434 5.574150 1
## S3 8.675529 1.484111 1
## S4 -18.765582 5.605868 1#Medoids
clara.res$medoids## x y
## S121 -1.531137 1.145057
## S455 48.357304 50.233499#Clustering
head(clara.res$clustering,10)## S1 S2 S3 S4 S5 S6 S7 S8 S9 S10
## 1 1 1 1 1 1 1 1 1 16.3.5 Visualizing CLARA clusters
Para visualizar los resultados de la partición, usaremos la función fviz_cluster() [paquete factoextra]. Dibuja un diagrama de dispersión de puntos de datos coloreados por números de grupo.
fviz_cluster(clara.res,
palette=c("#00AFBB","#FC4E07"), #color palette
ellipse.type="t", #Concentration ellipse
geom="point", pointsize=1,
ggtheme=theme_classic()
)
Chapter 7 Agglomeration Clustering
7.2.1 Data structure and preparation
Los datos deben ser una matriz numérica con:
• Filas que representan observaciones (individuos).
• Columnas que representan variables.
#Load the data
data(USArrests)
#Standardize the data
df <- scale(USArrests)
#Show the first 6 rows
head(df, nrow=6)## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.8649672077.2.2 Similarity measures
En el contexto de la agrupación jerárquica, es crucial evaluar la similitud entre objetos o grupos para determinar las fusiones o divisiones adecuadas. Se dispone de diversos métodos para medir la (des)similitud, como las distancias euclidiana y de Manhattan detalladas en el Capítulo 3. En el entorno R, la función dist() permite calcular estas distancias entre cada par de objetos en un conjunto de datos, generando así una matriz de distancia o disimilitud que sirve como base para las decisiones de agrupamiento.
#Compute the dissimilarity matrix
#df=the standardized data
res.dist <- dist(df, method="euclidean")
as.matrix(res.dist)[1:6, 1:6]## Alabama Alaska Arizona Arkansas California Colorado
## Alabama 0.000000 2.703754 2.293520 1.289810 3.263110 2.651067
## Alaska 2.703754 0.000000 2.700643 2.826039 3.012541 2.326519
## Arizona 2.293520 2.700643 0.000000 2.717758 1.310484 1.365031
## Arkansas 1.289810 2.826039 2.717758 0.000000 3.763641 2.831051
## California 3.263110 3.012541 1.310484 3.763641 0.000000 1.287619
## Colorado 2.651067 2.326519 1.365031 2.831051 1.287619 0.0000007.2.3 Linkage
La información de distancia proporcionada por la función dist() se utiliza en la vinculación para agrupar objetos en pares similares, formando grupos que luego se vinculan para crear estructuras jerárquicas más extensas. Este procedimiento se repite hasta que todos los objetos del conjunto de datos original están interconectados en un árbol jerárquico. hclust() puede ser utilizado de la siguiente manera:
res.hc <- hclust(d=res.dist,method="ward.D2")7.2.4 Dendrogram
Los dendrogramas corresponden a la representación gráfica del árbol jerárquico generado por la función hclust(). El dendrograma se puede producir en R usando la función base plot(res.hc), donde res.hc es la salida de hclust(). Aquí, usaremos la función fuiz_dend()[ en el paquete factoextra R] para producir un hermoso dendrograma. Primero instale factoextra escribiendo esto: install.packages(“factoextra”); A continuación visualice el dendrograma de la siguiente manera:
#cex: label size
library("factoextra")
fviz_dend(res.hc, cex=0.5)## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
7.3 Verify the cluster tree
Las diferentes técnicas de vinculación en la agrupación jerárquica incluyen la vinculación máxima, mínima, media, centroide y el método de varianza mínima de Ward. La vinculación máxima utiliza la distancia máxima entre elementos de dos conglomerados, produciendo grupos más compactos. La vinculación mínima utiliza la distancia mínima, generando grupos más largos y dispersos. La vinculación media utiliza la distancia promedio, mientras que la vinculación centroide utiliza la distancia entre los centroides de los conglomerados. El método de varianza mínima de Ward minimiza la varianza total dentro de cada conglomerado, fusionando en cada paso el par con la distancia mínima.
#Compute cophentic distance
res.coph <- cophenetic(res.hc)
#Correlation between cophenetic distance and
#the original distance
cor(res.dist, res.coph)## [1] 0.6975266res.hc2 <- hclust(res.dist,method="average")
cor(res.dist, cophenetic(res.hc2))## [1] 0.71803827.4 Cut the dendrogram into different groups
Cortar el dendrograma para formar grupos en la agrupación jerárquica plantea el desafío de determinar el número de grupos y la altura adecuada para el corte. La función cutree() de R, aplicada al árbol generado por hclust(), facilita este proceso al permitir especificar el número deseado de grupos o la altura de corte. Al ejecutar cutree(), se obtiene un vector que asigna a cada observación el número correspondiente al grupo al que pertenece.
# Cut tree into 4 groups
grp <- cutree(res.hc,k=4)
head(grp, n=4)## Alabama Alaska Arizona Arkansas
## 1 2 2 3#Number of members in each cluster
table(grp)## grp
## 1 2 3 4
## 7 12 19 12#Get the names for the members of cluster 1
rownames(df)[grp==1]## [1] "Alabama" "Georgia" "Louisiana" "Mississippi"
## [5] "North Carolina" "South Carolina" "Tennessee"#Cut in 4 groups and color by groups
fviz_dend(res.hc, k=4, #Cut in four groups
cex=0.5, #label size
k_colors=c("#2E9FDF","#00AFBB","#E7B800","#FC4E07"),
color_labels_by_k=TRUE, #color labels by groups
rect=TRUE #Add rectangle around groups
)
fviz_cluster(list(data=df, cluster=grp),
palette=c("#2E9FDF","#00AFBB","#E7B800","#FC4E07"),
ellipse.type="convex", #Concentration ellipse
repel=TRUE, #Avoid label overplotting(slow)
show.clust.cent=FALSE, ggtheme=theme_minimal())
7.5 Cluster R package
El paquete cluster de R facilita la realización de análisis de clusters en R. Proporciona las funciones agnes() y diana() para calcular clustering aglomerativo y divisivo, respectivamente. Estas funciones realizan todos los pasos necesarios para usted. No es necesario ejecutar las funciones scale(), dist() y hclust() por separado. Las funciones se pueden ejecutar de la siguiente manera:
library("cluster")
#Agglomerative Nesting(Hierarchical Clustering)
res.agnes <- agnes(x=USArrests, #data matrix
stand=TRUE, #Standardize the data
metric="euclidean", #metric for distance matrix
method="ward" #Linkage method
)
#DIvisive ANAlysis Clustering
res.diana <- diana(x=USArrests, #data matrix
stand=TRUE, #standardize the data
metric="euclidean" #metric for distance matrix
)
fviz_dend(res.agnes, cex=0.6,k=4)
8. Comparing Dendrograms
Tras abordar el cálculo del agrupamiento jerárquico en el Capítulo 7, se presenta ahora la comparación de dos dendrogramas mediante el uso del paquete dendextend en R. Este paquete ofrece diversas funciones para esta tarea, destacando especialmente dos de ellas: tanglegram(), que facilita la comparación visual de dos dendrogramas, y cor.dendlist(), destinada a calcular una matriz de correlación entre dendrogramas.
8.1 Data preparation
Usaremos los conjuntos de datos de USArrests base R y comenzaremos estandarizando las variables usando la función scale() de la siguiente manera
#Subset containing 10 rows
set.seed(123)
ss <- sample(1:50,10)
df<- df[ss,]8.2 Comparing dendrograms
Comenzamos creando una lista de dos dendrogramas calculando la agrupación jerárquica (HC) utilizandodos métodos de vinculación diferentes (“promedio” y “ward.D2”). A continuación, transformamoslos resultados como dendrogramas y creamos una lista para contener los dos dendrogramas.
library(dendextend)##
## ---------------------
## Welcome to dendextend version 1.17.1
## Type citation('dendextend') for how to cite the package.
##
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags:
## https://stackoverflow.com/questions/tagged/dendextend
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------##
## Adjuntando el paquete: 'dendextend'## The following object is masked from 'package:stats':
##
## cutree# Compute distance matrix
res.dist <- dist(df, method="euclidean")
#Compute 2 hierarchical clusterings
hc1 <- hclust(res.dist, method="average")
hc2 <- hclust(res.dist, method="ward.D2")
#Create two dendrograms
dend1 <- as.dendrogram(hc1)
dend2 <- as.dendrogram(hc2)
#Create a list to hold dendrograms
dend_list <- dendlist(dend1, dend2)8.2.1 Visual comparison of two dendrograms
En el contexto de la comparación visual de dos dendrogramas, emplearemos la función tanglegram() del paquete dendextend. Esta función traza ambos dendrogramas uno al lado del otro, conectando sus etiquetas con líneas. La calidad de la alineación de los dos árboles puede evaluarse mediante la función entrelazamiento(), donde un coeficiente más bajo indica una mejor alineación, variando entre 1 (entrelazamiento total) y 0 (sin entrelazamiento).
tanglegram(dend1,dend2)
tanglegram(dend1,dend2,
highlight_distinct_edges = FALSE, #Turn-off dashed lines
common_subtrees_color_lines = FALSE, #Turn-off line colors
common_subtrees_color_branches = TRUE, #Color common branches
main=paste("entanglement=", round(entanglement(dend_list),2))
)
8.2.2 Correlation matrix between a list of dendrograms
La función cor.dendlist() se utiliza para calcular la matriz de correlación “Baker” o “Cophenetic” entre una lista de árboles. El valor puede oscilar entre -1 y 1. Un valor cercano a 0 significa que los dos árboles no son estadísticamente similares.
#Cophenetic correlation matrix
cor.dendlist(dend_list,method = "cophenetic")## [,1] [,2]
## [1,] 1.0000000 0.9925544
## [2,] 0.9925544 1.0000000#Baker correlation matrix
cor.dendlist(dend_list,method="baker")## [,1] [,2]
## [1,] 1.0000000 0.9895528
## [2,] 0.9895528 1.0000000#Cophenetic correlation coefficient
cor_cophenetic(dend1, dend2)## [1] 0.9925544#Baker correlation coefficient
cor_bakers_gamma(dend1,dend2)## [1] 0.9895528#Create multiple dendrograms by chaining
dend1 <- df %>% dist %>% hclust("complete") %>% as.dendrogram
dend2 <- df %>% dist %>% hclust("single") %>% as.dendrogram
dend3 <- df %>% dist %>% hclust("average") %>% as.dendrogram
dend4 <- df %>% dist %>% hclust("centroid") %>% as.dendrogram
#Compute correlation matrix
dend_list <- dendlist("Complete"=dend1,"Single"=dend2,
"Average"=dend3, "Centroid"=dend4)
cors <- cor.dendlist(dend_list)
#Print correlation matrix
round(cors,2)## Complete Single Average Centroid
## Complete 1.00 0.46 0.45 0.30
## Single 0.46 1.00 0.23 0.17
## Average 0.45 0.23 1.00 0.31
## Centroid 0.30 0.17 0.31 1.00#Visualize the correlation matrix using corrplot package
library(corrplot)## corrplot 0.92 loadedcorrplot(cors,"pie","lower")
Chapter 9 Visualizing Dendrograms
En capítulos previos, se explicó que un dendrograma es una representación en forma de árbol de datos generados mediante métodos de agrupamiento jerárquico, como se detalló en el Capítulo 7. Este artículo presenta código en R para la visualización y personalización de dendrogramas, además de ofrecer información sobre cómo guardar y ampliar un dendrograma extenso. El proceso comienza con el cálculo de la agrupación jerárquica utilizando conjuntos de datos de USArrests.
#Load data
data(USArrests)
#Compute distances and hierarchical clustering
dd<- dist(scale(USArrests),method="euclidean")
hc<- hclust(dd,method="ward.D2")9.1 Visualizing dendrograms
Usaremos la función fviz_dend()[en el paquete R factoextra] para crear fácilmente un dendrograma usando el gráfico base de R o ggplot2. También proporciona una opción para dibujar dendrogramas circulares y árboles de tipo filogénico. Para crear dendrogramas básicos, escriba esto:
library(factoextra)
fviz_dend(hc, cex=0.5)
fviz_dend(hc, cex=0.5,
main="Dendrogram-ward.D2",
xlab="Objects",ylab="Distance",sub="")
fviz_dend(hc, cex=0.5, horiz=TRUE)
fviz_dend(hc, k=4,
cex=0.5)
fviz_dend(hc, k=4,
cex=0.5,
k_colors=c("#2E9FDF","#00AFBB","#E7B800","#FC4E07"),
color_labels_by_k = TRUE,
ggtheme=theme_gray()
)
fviz_dend(hc, cex=0.5, k=4,
k_colors="jco")
fviz_dend(hc,k=4,cex=0.4,horiz=TRUE,k_colors="jco",
rect=TRUE,rect_border = "jco",rect_fill = TRUE)
fviz_dend(hc,cex=0.5,k=4,
k_colors="jco",type="circular")
require("igraph")## Cargando paquete requerido: igraph##
## Adjuntando el paquete: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionfviz_dend(hc, k=4,k_colors="jco",
type="phylogenic",repel=TRUE)
require("igraph")
fviz_dend(hc, k=4,
k_colors="jco",
type="phylogenic",
repel=TRUE,
phylo_layout = "layout.gem")
9.2 Case of dendrogram with lage data sets
Cuando realiza el análisis de agrupación jerárquica en un conjunto de datos extenso, puede resultar beneficioso ampliar el dendrograma o representar únicamente un subconjunto de éste. Otra opción viable es generar el dendrograma en una página de gran tamaño en formato PDF, lo cual permite su ampliación sin sufrir pérdida de resolución.
9.2.1 Zooming in the dendrogram
Si desea ampliar los primeros grupos, es posible utilizar la opción xlim e ylim para limitar el área de trazado. Por ejemplo, escriba el siguiente código:
fviz_dend(hc,xlim=c(1,20), ylim=c(1,8))
9.2.2 Plotting a sub-tree of dendrograms
Para trazar un subárbol, seguiremos el procedimiento siguiente:
Cree el dendrograma completo usando fviz_dend() y guarde el resultado en un objeto, llamado dend_plot, por ejemplo.
Utilice la función base de R cut.dendrogram() para cortar el dendrograma, a una altura determinada (h), en múltiples subárboles. Esto devuelve una lista con los componentes Supper y $lower, el primero es una versión truncada del árbol original, también de clase dendrograma, el segundo una lista con las ramas obtenidas al cortar el árbol, cada una un dendrograma.
Visualice subárboles usando fviz_dend().
dend_plot <- fviz_dend(hc,k=4,
cex=0.5,
k_colors="jco"
)
dend_data <- attr(dend_plot,"dendrogram")
dend_cuts <- cut(dend_data, h=10)
fviz_dend(dend_cuts$upper)## Warning in min(-diff(our_dend_heights)): ningún argumento finito para min;
## retornando Inf
print(dend_plot)
fviz_dend(dend_cuts$lower[[1]], main="Subtree 1")
fviz_dend(dend_cuts$lower[[2]],main="Subtree 2")
fviz_dend(dend_cuts$lower[[2]],type="circular")
9.3 Manipulating dendrograms using dendextend
El paquete dendextend facilita la modificación visual de dendrogramas y la comparación entre ellos. En esta sección, utilizamos el operador de encadenamiento (%>%) para simplificar el código, convirtiendo múltiples operaciones en una forma más legible y estructurada. Este operador transforma expresiones del tipo x %>% f(y) en f(x, y), permitiendo la lectura secuencial de las operaciones de izquierda a derecha y de arriba a abajo. Los siguientes códigos R ilustran de manera equivalente esta idea.
data <- scale(USArrests)
dist.res <- dist(data)
hc <- hclust(dist.res, method="ward.D2")
dend <- as.dendrogram(hc)
plot(dend)
library(dendextend)
mycols <- c("#2E9FDF","#00AFBB","#E7B800","#FC4E07")
dend <- as.dendrogram(hc) %>%
set("branches_lwd",1) %>%
set("branches_k_color",mycols,k=4) %>%
set("labels_colors",mycols,k=4) %>%
set("labels_cex",0.5)
fviz_dend(dend)
Chapter 10 Heatmap:Static and Interactive
El mapa de calor, también conocido como imagen de color falso, ofrece una representación visual de la agrupación jerárquica en datos. Este método transforma los valores de los datos en una escala de colores, permitiendo la observación simultánea de grupos de muestras y características. La agrupación jerárquica se aplica a las filas y columnas de la matriz de datos, reorganizando estas filas y columnas según los resultados de la agrupación. Este proceso posiciona las observaciones similares en proximidad. Los bloques de valores “altos” y “bajos” aparecen contiguos en la matriz de datos. Posteriormente, se utiliza un esquema de colores para la visualización final de la matriz de datos. Esta representación facilita la identificación de variables que pueden ser características distintivas de cada grupo de muestras.
10.3 R base heatmap:heatmap()
Se puede utilizar la función R heatmap() incorporada [en el paquete de estadísticas]. Un formato simplificado es:
df <- scale(mtcars)
heatmap(df,scale="row")
heatmap(df,scale="none")
col <- colorRampPalette(c("red","white","blue"))(256)
library("RColorBrewer")
col <- colorRampPalette(brewer.pal(10,"RdYlBu"))(256)
heatmap(df,scale="none",col= col,
RowSideColors = rep(c("blue","pink"), each=16),
ColSideColors = c(rep("purple",5),rep("orange",6)))
10.4 Enhaced heat maps:heatmap.2()
La función heatmap.2() [en el paquete gplots] proporciona muchas extensiones a la función estándar R heatmap() presentada en la sección anterior.
#install.packages("gplots")
library("gplots")##
## Adjuntando el paquete: 'gplots'## The following object is masked from 'package:stats':
##
## lowessheatmap.2(df,scale="none", col=bluered(100),
trace="none", density.info = "none")
10.5 Pretty heat maps:pheatmap()
Primero, instale el paquete pheatmap: install.packages(“pheatmap”); luego escribe esto:
library("pheatmap")
pheatmap(df, cutree_rows=4)
10.6 Interactive heat maps:d3heatmap()
if (!require("devtools")) install.packages("devtools")## Cargando paquete requerido: devtools## Cargando paquete requerido: usethisdevtools::install_github("rstudio/d3heatmap")## Skipping install of 'd3heatmap' from a github remote, the SHA1 (c8a3c64a) has not changed since last install.
## Use `force = TRUE` to force installationlibrary("d3heatmap")##
## ======================
## Welcome to d3heatmap version 0.9.0
##
## Type citation('d3heatmap') for how to cite the package.
## Type ?d3heatmap for the main documentation.
##
## The github page is: https://github.com/talgalili/d3heatmap/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/d3heatmap/issues
## You may ask questions at stackoverflow, use the r and d3heatmap tags:
## https://stackoverflow.com/questions/tagged/d3heatmap
## ======================##
## Adjuntando el paquete: 'd3heatmap'## The following objects are masked from 'package:base':
##
## print, saved3heatmap(scale(mtcars), colors="RdYlBu",
k_row=4, #Number of groups in rows
k_col=2 #Number of groups in columns
)## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette RdYlBu is 11
## Returning the palette you asked for with that many colors
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette RdYlBu is 11
## Returning the palette you asked for with that many colors10.7 Enhancing heatmaps using dendextend
El paquete dendextend se puede utilizar para mejorar funciones de otros paquetes. Los datos de mtcars se utilizan en las siguientes secciones. Comenzaremos definiendo el orden y la apariencia de las filas y columnas usando dendextend. Estos resultados se utilizan en otras funciones de otros paquetes.
library(dendextend)
#order for rows
Rowv <- mtcars %>% scale %>% dist %>% hclust %>% as.dendrogram %>%
set("branches_k_color",k=3) %>% set("branches_lwd",1.2) %>%
ladderize
#Order for columns:We must transpose the data
Colv <- mtcars %>% scale %>% t %>% dist %>% hclust %>% as.dendrogram %>%
set("branches_k_color", k=2, value=c("orange","blue")) %>%
set("branches_lwd",1.2) %>%
ladderizeheatmap(scale(mtcars),Rowv=Rowv,Colv=Colv,
scale="none")
library(gplots)
heatmap.2(scale(mtcars),scale="none",col=bluered(100),
Rowv=Rowv, Colv=Colv,
trace="none",density.info = "none")
10.8 Complex heatmap
ComplexHeatmap es un paquete R/bioconductor, desarrollado por Zuguang Gu, que proporciona una solución flexible para organizar y anotar múltiples mapas de calor. También permite visualizar la asociación entre diferentes datos de diferentes fuentes.
library(ComplexHeatmap)## Cargando paquete requerido: grid## ========================================
## ComplexHeatmap version 2.20.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
## in the original pheatmap() are identically supported in the new function. You
## can still use the original function by explicitly calling pheatmap::pheatmap().##
## Adjuntando el paquete: 'ComplexHeatmap'## The following object is masked from 'package:pheatmap':
##
## pheatmapheatmap(df,
name="mtcars", #title of legend
column_title="Variables",row_title="Samples",
row_names_gp=gpar(fontsize=7) #Text size for row names
)## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.window(...): "column_title" es un parámetro gráfico inválido## Warning in plot.window(...): "row_title" es un parámetro gráfico inválido## Warning in plot.window(...): "row_names_gp" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "column_title" es un parámetro gráfico
## inválido## Warning in plot.xy(xy, type, ...): "row_title" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "row_names_gp" es un parámetro gráfico
## inválido## Warning in title(...): "name" es un parámetro gráfico inválido## Warning in title(...): "column_title" es un parámetro gráfico inválido## Warning in title(...): "row_title" es un parámetro gráfico inválido## Warning in title(...): "row_names_gp" es un parámetro gráfico inválido
library(grid)
library(dendextend)
row_dend=hclust(dist(df)) #row clustering
col_dend=hclust(dist(t(df))) #column clustering
heatmap(df,name="mtcars",
row_names_gp=gpar(fontsize=6.5),
cluster_rows=color_branches(row_dend,k=4),
cluster_columns=color_branches(col_dend,k=2))## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.window(...): "row_names_gp" es un parámetro gráfico inválido## Warning in plot.window(...): "cluster_rows" es un parámetro gráfico inválido## Warning in plot.window(...): "cluster_columns" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "row_names_gp" es un parámetro gráfico
## inválido## Warning in plot.xy(xy, type, ...): "cluster_rows" es un parámetro gráfico
## inválido## Warning in plot.xy(xy, type, ...): "cluster_columns" es un parámetro gráfico
## inválido## Warning in title(...): "name" es un parámetro gráfico inválido## Warning in title(...): "row_names_gp" es un parámetro gráfico inválido## Warning in title(...): "cluster_rows" es un parámetro gráfico inválido## Warning in title(...): "cluster_columns" es un parámetro gráfico inválido
#Divide into 2 groups
set.seed(2)
heatmap(df,name="mtcars",k=2)## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in title(...): "name" es un parámetro gráfico inválido
#split by a vector specifying rowgroups
heatmap(df,name="mtcars",split=mtcars$cyl,
row_names_gp=gpar(fontsize=7))## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.window(...): "split" es un parámetro gráfico inválido## Warning in plot.window(...): "row_names_gp" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "split" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "row_names_gp" es un parámetro gráfico
## inválido## Warning in title(...): "name" es un parámetro gráfico inválido## Warning in title(...): "split" es un parámetro gráfico inválido## Warning in title(...): "row_names_gp" es un parámetro gráfico inválido
10.8.3.1 Simple annotation
Puedes diseñar un heatmap simple de la siguiente manera:
heatmap(df, name="mtcars", split=mtcars$cyl,
row_names_gp=gpar(fontsize=7))## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.window(...): "split" es un parámetro gráfico inválido## Warning in plot.window(...): "row_names_gp" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "split" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "row_names_gp" es un parámetro gráfico
## inválido## Warning in title(...): "name" es un parámetro gráfico inválido## Warning in title(...): "split" es un parámetro gráfico inválido## Warning in title(...): "row_names_gp" es un parámetro gráfico inválido
18.8.3.2 Complex annotation
En esta sección veremos cómo combinar mapas de calor y algunos gráficos básicos para mostrar la distribución de datos. Para gráficos de anotaciones simples, se pueden utilizar las siguientes funciones: anno_points(), anno_barplot(), anno_boxplot(), anno_density() y anno_histogram().
``` r
.hist = anno_histogram(df, gp = gpar(fill = "lightblue"))
.density = anno_density(df, type = "line", gp = gpar(col = "blue"))
ha_mix_top = HeatmapAnnotation(
hist = .hist, density = .density,
height = unit(3.8, "cm")
)
# Definir algunos gráficos
.violin = anno_density(df, type = "violin",
gp = gpar(fill = "lightblue"), which = "row")
.boxplot = anno_boxplot(df, which = "row")
ha_mix_right = HeatmapAnnotation(violin = .violin, bxplt = .boxplot,
which = "row", width = unit(4, "cm"))
# Combinar la anotación con el heatmap
Heatmap(df, name = "mtcars",
column_names_gp = gpar(fontsize = 8),
top_annotation = ha_mix_top) + ha_mix_right
10.8.3.3 Combining multiple heatmaps
Para combinar múltiples heatmaps, sigue los pasos que se muestran a continuación.
#Heatmap 1
ht1=heatmap(df,name="ht1",km=2,
column_names_gp=gpar(fontsize=9))## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.window(...): "km" es un parámetro gráfico inválido## Warning in plot.window(...): "column_names_gp" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "km" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "column_names_gp" es un parámetro gráfico
## inválido## Warning in title(...): "name" es un parámetro gráfico inválido## Warning in title(...): "km" es un parámetro gráfico inválido## Warning in title(...): "column_names_gp" es un parámetro gráfico inválido
#Heatmap 2
ht2=heatmap(df,name="ht2",
col=circlize::colorRamp2(c(-2,0,2),c("green","white","red")),
column_names_gp=gpar(fontsize=9))## Warning in plot.window(...): "name" es un parámetro gráfico inválido## Warning in plot.window(...): "column_names_gp" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "name" es un parámetro gráfico inválido## Warning in plot.xy(xy, type, ...): "column_names_gp" es un parámetro gráfico
## inválido## Warning in title(...): "name" es un parámetro gráfico inválido## Warning in title(...): "column_names_gp" es un parámetro gráfico inválido## Warning in image.default(1L:nc, 1L:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 + :
## supplied color is neither numeric nor character
10.9 Application to gene expression matrix
En los datos de expresión genética, las filas son genes y las columnas son muestras. Se puede adjuntar más información sobre los genes después del mapa de calor de expresión, como la longitud del gen y el tipo de genes.
expr <- readRDS(paste0(system.file(package = "ComplexHeatmap"),
"/extdata/gene_expression.rds"))
mat <- as.matrix(expr[, grep("cell", colnames(expr))])
type <- gsub("s\\d+_", "", colnames(mat))
ha = HeatmapAnnotation(
df = data.frame(type = type),
annotation_height = unit(4, "mm")
)
Heatmap(mat, name = "expression", km = 5, top_annotation = ha,
show_row_names = FALSE, show_column_names = FALSE) +
Heatmap(expr$length, name = "length", width = unit(5, "mm"),
col = circlize::colorRamp2(c(0, 100000), c("white", "orange"))) +
Heatmap(expr$type, name = "type", width = unit(5, "mm")) +
Heatmap(expr$chr, name = "chr", width = unit(5, "mm"),
col = circlize::rand_color(length(unique(expr$chr))))## There are 23 unique colors in the vector `col` and 23 unique values in
## `matrix`. `Heatmap()` will treat it as an exact discrete one-to-one
## mapping. If this is not what you want, slightly change the number of
## colors, e.g. by adding one more color or removing a color.
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.