Casos y Ejemplos de Análisis Multivariante con R

Author

Alex Sánchez y Francesc Carmona

Published

October 23, 2024

Presentación

Este documento contiene ejemplos de análisis multivariante con R.

Los ejemplos están pensados para ilustrar los capítulos del documento “Notas de Análisis Multivariante” pero se mantienen en un documento aparte para poder ser usados autónomamente.

Algunos de los ejemplos se han inspirado en un documento no publicado del Dr. Carles Cuadras nuestro mentor en Análisis Multivariante por lo que deseamos hacer agradecimiento explícito.

Instalación de paquetes utilizados

Muestra el código

if (!require(FactoMineR)) install.packages("FactoMineR", dep=TRUE)
if (!require(factoextra)) install.packages("factoextra", dep=TRUE)
if (!require(ade4))
  install.packages("ade4", dep=TRUE)
if (!require(BiocManager))
  install.packages("BiocManager", dep=TRUE)
if (!require(made4))
  BiocManager::install("made4")

Análisis de componentes principales

Ejemplo PCA-1: Búsqueda de factores latentes en datos ecológicos¹

Un estudio ecológico recogió datos de 49 gorriones hembras que fueron recogidos casi moribundos después de un temporal. Los investigadores midieron 5 magnitudes de los animales (“length”,“wing”,“head”,“humerus”,“sternum”) y anotaron si los animales sobrevivían o no a los pocos días de ser recogidos (“survival”). Entre otros objetivos, el estudio perseguía encontrar factores, es decir algunas variables latentes, no medibles pero presentes en los datos, que explicaran la supervivencia de los mismos ante tempestades como la registrada. Como veremos, las dos primeras componentes principales resultan ser estos factores.

Muestra el código

load("datos/gorriones.RData")
colnames(gorriones) <-c("length","wing","head","humerus","sternum","survival")
gorrionesNum<-as.matrix(gorriones[,1:5])

Empezamos con una exploración de los datos. Aunque existen muchos paquetes que permiten hacerla nos restringiremos a las funciones de R Base (o de tidyverse) para facilitar la reproducibilidad.

Muestra el código

str(gorriones)

'data.frame':   49 obs. of  6 variables:
 $ length  : num  156 154 153 153 155 163 157 155 164 158 ...
 $ wing    : num  245 240 240 236 243 247 238 239 248 238 ...
 $ head    : num  31.6 30.4 31 30.9 31.5 32 30.9 32.8 32.7 31 ...
 $ humerus : num  18.5 17.9 18.4 17.7 18.6 19 18.4 18.6 19.1 18.8 ...
 $ sternum : num  20.5 19.6 20.6 20.2 20.3 20.9 20.2 21.2 21.1 22 ...
 $ survival: Factor w/ 2 levels "N","S": 2 2 2 2 2 2 2 2 2 2 ...

Muestra el código

summary(gorriones)

     length         wing            head          humerus         sternum     
 Min.   :152   Min.   :230.0   Min.   :30.10   Min.   :17.20   Min.   :18.60  
 1st Qu.:155   1st Qu.:238.0   1st Qu.:30.90   1st Qu.:18.10   1st Qu.:20.20  
 Median :158   Median :242.0   Median :31.50   Median :18.50   Median :20.70  
 Mean   :158   Mean   :241.3   Mean   :31.46   Mean   :18.47   Mean   :20.83  
 3rd Qu.:161   3rd Qu.:245.0   3rd Qu.:32.00   3rd Qu.:18.80   3rd Qu.:21.50  
 Max.   :165   Max.   :252.0   Max.   :33.40   Max.   :19.80   Max.   :23.10  
 survival
 N:28    
 S:21

Muestra el código

f<- function(x){
  ifelse (is.numeric(x), 
          hist(x, breaks=5),
          barplot(table(x))
  )
}
par(mfrow=c(3,2))
apply(gorriones,2,f)

  length     wing     head  humerus  sternum survival 
     0.7      0.7      0.7      0.7      0.7      0.7

Muestra el código

par(mfrow=c(1,1))

La escala de las variables numéricas es heterogénea, pero los datos son del mismo tipo (magnitudes biométricas) por lo que nos basaremos en la matriz de covarianzas de los datos centrados.

Muestra el código

gorrionesNum <- scale(gorrionesNum, center = TRUE, scale=FALSE)
apply(gorrionesNum,2, mean)

       length          wing          head       humerus       sternum 
-4.640173e-15 -1.102063e-14  1.232600e-15  4.350306e-16 -2.175331e-16

Calculamos la matriz de varianzas ajustando a dividir por n en vez de por (n-1) para que los resultados sean comparables en los dos casos que haremos.

Muestra el código

n<- dim(gorriones)[1]
S<-cov(gorrionesNum)*(n-1)/n
show(S)

           length      wing      head   humerus   sternum
length  13.081216 13.333195 1.8828405 1.3034569 2.1474802
wing    13.333195 25.158684 2.6582257 2.1528530 2.6035818
head     1.882840  2.658226 0.6187422 0.3352811 0.4061849
humerus  1.303457  2.152853 0.3352811 0.3119200 0.3324448
sternum  2.147480  2.603582 0.4061849 0.3324448 0.9627655

Calculamos la matriz de correlaciones.

Muestra el código

R<-cor(gorrionesNum)
show(R)

           length      wing      head   humerus   sternum
length  1.0000000 0.7349642 0.6618119 0.6452841 0.6051247
wing    0.7349642 1.0000000 0.6737411 0.7685087 0.5290138
head    0.6618119 0.6737411 1.0000000 0.7631899 0.5262701
humerus 0.6452841 0.7685087 0.7631899 1.0000000 0.6066493
sternum 0.6051247 0.5290138 0.5262701 0.6066493 1.0000000

A modo de ilustración empezamos calculando las componentes principales con funciones básicas:

Mediante diagonalización de la matriz de covarianzas

Muestra el código

EIG <- eigen(S)
show(EIG)

eigen() decomposition
$values
[1] 34.60482313  4.52812344  0.61804191  0.30640029  0.07593901

$vectors
            [,1]        [,2]        [,3]        [,4]        [,5]
[1,] -0.53650052  0.82809990  0.15649065  0.04020969 -0.01765243
[2,] -0.82901535 -0.55051223  0.05774395  0.06902156  0.03964203
[3,] -0.09649615  0.03356237 -0.23751487 -0.89762653  0.35695288
[4,] -0.07435219 -0.01459529 -0.20324541 -0.30724056 -0.92658150
[5,] -0.10030441  0.09923405 -0.93512262  0.30575979  0.11021920

Los vectores propios, es decir las columnas de la matriz “eigen$vectors” son las coordenadas de las cinco componentes principales.

La transformación de los datos asociada a las componentes principales se obtiene multiplicando la matriz original por la matriz de vectores propios

Muestra el código

eigenVecs1 <- EIG$vectors
PCAS1 <- gorrionesNum %*% eigenVecs1
head(PCAS1)

           [,1]      [,2]        [,3]        [,4]        [,5]
[1,] -1.9664223 -3.689713  0.16801141 -0.06169590  0.16647852
[2,]  3.5023362 -2.714180  0.81488580  0.49908927  0.03198139
[3,]  3.8434585 -3.430206 -0.52085909  0.07244318 -0.08926600
[4,]  7.2613378 -1.260991 -0.21176257 -0.02111595  0.32098994
[5,]  0.2503842 -3.441451  0.08648433 -0.24206207 -0.04555039
[6,] -7.4958530  1.051783  0.80825616 -0.03255193 -0.15422636

Con esto podemos hacer una primera representación.

Muestra el código

plot(PCAS1[,1], PCAS1[,2], main = "Gorriones. 2 primeras PCs")

Los valores propios representan la varianza de las componentes por lo que cada valor dividido por la suma de éstos es el porcentaje de variabilidad explicado.

Muestra el código

vars1<- EIG$values/sum(EIG$values)
round(vars1,3)

[1] 0.862 0.113 0.015 0.008 0.002

Podemos usar estos porcentajes para etiquetar los ejes indicando la importancia de cada componente.

Muestra el código

xlabel <- paste("PCA1 ", round(vars1[1]*100, 2),"%" )
ylabel <- paste("PCA2 ", round(vars1[2]*100,2),"%" )
plot(PCAS1[,1], PCAS1[,2], main = "Gorriones. 2 primeras PCs",
     xlab=xlabel, ylab=ylabel)

Finalmente podemos visualizar si el gorrión ha sobrevivido o no representándolo en el gráfico con colores distintos.

Muestra el código

bgSurv<- colSurv <- ifelse(gorriones$survival=="N", "red", "blue")
pchSurv <- ifelse(gorriones$survival=="N",1,2)
plot(PCAS1[,1], PCAS1[,2], main = "Gorriones. 2 primeras PCs",
     xlab=xlabel, ylab=ylabel, 
     pch=pchSurv, col=colSurv, bg=bgSurv)
legend( "topright"  , inset = c(0.01,0.01), cex =1, 
        bty = "y", legend = c("N", "S"), 
        text.col = c("red", "blue"),
        col = c("red", "blue"), 
        pt.bg = c("red","blue"), 
        pch = c(1,2)
        )

Calculo de las componentes principales mediante las funcionesn `princomp` y `prcomp`

La función princomp calcula las componentes principales recurriendo básicamente al mismo método, a diferencia de la función prcomp que calcula las componentes principales mediante la descomposición en valores singulares de la matriz de datos.

Muestra el código

PCAS2 <- princomp(gorrionesNum)
names(PCAS2)

[1] "sdev"     "loadings" "center"   "scale"    "n.obs"    "scores"   "call"

Muestra el código

PCAS3 <- prcomp(gorrionesNum)
names(PCAS3)

[1] "sdev"     "rotation" "center"   "scale"    "x"

Observemos como algunos resultados coinciden mientras que otros tendremos que ajustarlos.

EL argumento sdev contiene las desviaciones estándard es decir las raíces cuadradas de los valores propios.

Muestra el código

PCAS2$sdev

   Comp.1    Comp.2    Comp.3    Comp.4    Comp.5 
5.8825864 2.1279388 0.7861564 0.5535344 0.2755703

Muestra el código

PCAS3$sdev

[1] 5.9435475 2.1499905 0.7943033 0.5592706 0.2784261

Muestra el código

sqrt(EIG$values)

[1] 5.8825864 2.1279388 0.7861564 0.5535344 0.2755703

Hay una ligera diferencia en los resultados calculados por prcomp que podemos atribuir a la forma de cálculo.

El argumento loadings contiene los vectores propios en el objeto calculado por `princomp.

En el caso de prcomp estan en el objeto rotation.

Muestra el código

PCAS2$loadings


Loadings:
        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
length   0.537  0.828  0.156              
wing     0.829 -0.551                     
head                  -0.238 -0.898 -0.357
humerus               -0.203 -0.307  0.927
sternum  0.100        -0.935  0.306 -0.110

               Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
SS loadings       1.0    1.0    1.0    1.0    1.0
Proportion Var    0.2    0.2    0.2    0.2    0.2
Cumulative Var    0.2    0.4    0.6    0.8    1.0

Muestra el código

PCAS3$rotation[,1]

    length       wing       head    humerus    sternum 
0.53650052 0.82901535 0.09649615 0.07435219 0.10030441

Muestra el código

EIG$vectors

            [,1]        [,2]        [,3]        [,4]        [,5]
[1,] -0.53650052  0.82809990  0.15649065  0.04020969 -0.01765243
[2,] -0.82901535 -0.55051223  0.05774395  0.06902156  0.03964203
[3,] -0.09649615  0.03356237 -0.23751487 -0.89762653  0.35695288
[4,] -0.07435219 -0.01459529 -0.20324541 -0.30724056 -0.92658150
[5,] -0.10030441  0.09923405 -0.93512262  0.30575979  0.11021920

Observamos como en este caso coinciden todos los valores, salvo por el signo del primer valor propio del cálculo basado en la diagonalización. Esto no es un error sino que se debe a que el vector propio no es único y está indeterminado por un producto por -1.

Finalmente el argumento scores contiene las componentes principales calculadas por princomp. en el caso de prcomp se encuentran en el campo x. Para hacerlas comparables centraremos las componentes que hemos calculado nosotros.

Muestra el código

head(PCAS2$scores)

         Comp.1    Comp.2      Comp.3      Comp.4      Comp.5
[1,]  1.9664223 -3.689713  0.16801141 -0.06169590 -0.16647852
[2,] -3.5023362 -2.714180  0.81488580  0.49908927 -0.03198139
[3,] -3.8434585 -3.430206 -0.52085909  0.07244318  0.08926600
[4,] -7.2613378 -1.260991 -0.21176257 -0.02111595 -0.32098994
[5,] -0.2503842 -3.441451  0.08648433 -0.24206207  0.04555039
[6,]  7.4958530  1.051783  0.80825616 -0.03255193  0.15422636

Muestra el código

head(PCAS3$x)

            PC1       PC2         PC3         PC4         PC5
[1,]  1.9664223  3.689713 -0.16801141  0.06169590  0.16647852
[2,] -3.5023362  2.714180 -0.81488580 -0.49908927  0.03198139
[3,] -3.8434585  3.430206  0.52085909 -0.07244318 -0.08926600
[4,] -7.2613378  1.260991  0.21176257  0.02111595  0.32098994
[5,] -0.2503842  3.441451 -0.08648433  0.24206207 -0.04555039
[6,]  7.4958530 -1.051783 -0.80825616  0.03255193 -0.15422636

Muestra el código

head(PCAS1)

           [,1]      [,2]        [,3]        [,4]        [,5]
[1,] -1.9664223 -3.689713  0.16801141 -0.06169590  0.16647852
[2,]  3.5023362 -2.714180  0.81488580  0.49908927  0.03198139
[3,]  3.8434585 -3.430206 -0.52085909  0.07244318 -0.08926600
[4,]  7.2613378 -1.260991 -0.21176257 -0.02111595  0.32098994
[5,]  0.2503842 -3.441451  0.08648433 -0.24206207 -0.04555039
[6,] -7.4958530  1.051783  0.80825616 -0.03255193 -0.15422636

De nuevo los valores coinciden aunque los dos métodos prcomp y princomp han centrado las componentes principales antes de almacenarlas.

Interpretación de las componentes

A partir de los resultados obtenidos podemos escribir así las dos primeras componentes:

$Y_1 = 0.5365\times length + 0.8290\times wing + 0.0965\times head + 0.0743\times humerus + 0.1003\times sternum$

$Y_2 = -0.8281\times length + 0.5505\times wing -0.0336\times head + 0.0146\times humerus -0.0992\times sternum$

La primera componente tiene todos los coeficientes positivos del mismo signo y se puede interpretar como el tamaño del gorrión.

El peso principal de las variables originales se da en las dos primeras. El mayor peso se da a la extensión de las alas, con un coeficiente de 0.83. Le sigue la longitud del pájaro con un coeficiente de 0.54. Las demás variables presentan una contribución notablemente menor a la PC1. Esta primera componente explica un 82% de la variabilidad.

La segunda componente tiene coeficientes positivos y negativos y se puede interpretar como un factor forma. En este caso por contraposición de las dos primeras variables. En el eje PC2 las variables con mayor peso en valor absoluto son las mismas pero aquí una aparece con signo positivo y la otra con signo negativo: 0.55 para las alas, −0.83 para la longitud. Esta segunda componente explica un 11.2% de la variabilidad.

Valores positivos de PC2 corresponderán a pájaros cortos con gran envergadura de alas, y valores negativos de PC2 corresponderán a pájaros más proporcionados o incluso de mayor longitud que envergadura.

Tamaño y forma son pues dos componentes interrelacionadas que, entre ambas explican el 97.51% de la variabilidad. Si atendemos a la distribución de los “S” y “N” en el gráfico podríamos sugerir que los gorriones “extremos” ya sea en tamaño o en forma tienen menos posibilidades de sobrevivir.

¹ El ejemplo de los gorriones se ha tomado del documento Exemples d’Anàlisi Multivariant del Dr. Carles Cuadras disponible en su página web.

Ejemplo PCA-2: Detección del efecto batch en datos de microarrays.

Muchos datos de alto rendimiento como los microarrays o los de proteómica presentan un tipo particular de complicación técnica que se conoce como efecto . Dicho efecto consiste en que muestras producidas en el mismo ``lote’’ (mismo día, misma tanda de hibridación, mismo operario, …) se parecen más entre ellas que muestras producidas en lotes distintos pudiendo llegar a causar confusión en estudios comparativos en los que el efecto batch enmascare el efecto de los tratamientos en estudio.

Si la presencia del efecto batch se conoce o espera puede, en principio, ser controlado o cancelado mediante un adecuado diseño experimental, en el que el lote se tratará usualmente como un bloque. En este caso un adecuado balanceo entre bloques y tratamientos puede cancelar el efecto batch. Alternativamente, dicho efecto puede incluirse como un factor en el modelo de análisis de la varianza, lo que eliminara del análisis la variación atribuible al lote.

En muchas (demasiadas) ocasiones el efecto batch no ha sido considerado de antemano, o incluso, cuando sí lo ha sido, puede que no se hayan tenido en cuenta todos los posible efectos (a lo mejor se ha tenido en cuenta el día pero no que los reactivos utilizados provenían de dos lotes o que las piezas se han procesado en grupos, o …).

En prevención de los problemas que esto puede generar es conveniente realizar algún tipo de análisis exploratorio que permita detectar la posible presencia de estos efectos. En caso de detectarse efectos indeseados puede plantearse su eliminación mediante alguna de las técnicas disponibles para ello.

Los datos para el análisis

Los datos para este ejemplo consisten en datos de microarrays de expresión génica, tipo hgu95-a utilizados para analizar un estudio de cáncer de mama en el que se analiza el efecto de distintos tratamientos y distintos tiempos de exposición a éstos, sobre la expresión de los genes.

Los datos y la información del estudio se hallan disponibles en la base de datos Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo con número de acceso ``GSE848’’.

Para este ejemplo hemos utilizado un subconjunto de estos datos basado en 18 muestras, que hemos guardado en un objeto binario datosMicroarrays.Rda con el fin de simplificar el proceso y centrar la atención en el análisis.

Muestra el código

load(file.path("datos", "datosMicroarrays.Rda"))
head(x)

            E2.8.A   E2.8.B  E2.48.A  E2.48.B E2.ICI.8.A E2.ICI.8.B E2.ICI.48.A
X100_g_at 7.933691 8.098032 8.256445 7.619120   7.800253   7.886916    8.015694
X101_at   6.877744 6.811214 6.803485 6.989820   6.816344   7.043301    6.716991
X102_at   4.491853 3.121015 4.655352 3.877744   3.536053   2.765535    4.343408
X103_at   5.244126 3.350497 3.017922 4.277985   3.419539   2.536053    3.916477
X104_at   5.173927 4.842979 5.350497 4.596935   4.209453   4.990955    5.185867
X105_at   2.104337 1.807355 2.263034 2.536053   3.485427   1.584963    4.035624
          E2.ICI.48.B E2.Ral.8.A E2.Ral.8.B E2.Ral.48.A E2.Ral.48.B E2.TOT.8.A
X100_g_at    7.878971   7.660353   7.613237    7.798958    7.675957   8.126188
X101_at      6.772150   6.515700   6.774787    6.205549    6.841722   6.623516
X102_at      3.017922   3.350497   2.963474    4.153805    3.169925   4.321928
X103_at      3.954196   3.201634   2.459432    4.566815    2.765535   2.722466
X104_at      4.781360   5.314697   4.921246    5.228819    4.807355   4.672425
X105_at      4.321928   2.807355   3.722466    4.153805    3.459432   4.232661
          E2.TOT.8.B E2.TOT.48.A E2.TOT.48.B
X100_g_at   8.037821    8.249351    7.720415
X101_at     6.888743    6.675251    7.130313
X102_at     1.536053    4.035624    3.857981
X103_at     2.510962    3.104337    3.104337
X104_at     5.255501    5.057450    5.185867
X105_at     2.321928    3.765535    2.201634

Muestra el código

rownames(targets) <- targets$ShortName

La tabla targetscontiene información sobre el tratamiento y el lote de cada muestra analizada. los individuos analizados.

Muestra el código

kableExtra::kable(targets)

	Filename	Treatment	Time	Batch	ShortName
E2.8.A	GSM13099.txt	E2	8	A	E2.8.A
E2.8.B	GSM13138.txt	E2	8	B	E2.8.B
E2.48.A	GSM13139.txt	E2	48	A	E2.48.A
E2.48.B	GSM13140.txt	E2	48	B	E2.48.B
E2.ICI.8.A	GSM15900.txt	E2.ICI	8	A	E2.ICI.8.A
E2.ICI.8.B	GSM15901.txt	E2.ICI	8	B	E2.ICI.8.B
E2.ICI.48.A	GSM15902.txt	E2.ICI	48	A	E2.ICI.48.A
E2.ICI.48.B	GSM15903.txt	E2.ICI	48	B	E2.ICI.48.B
E2.Ral.8.A	GSM15904.txt	E2.Ral	8	A	E2.Ral.8.A
E2.Ral.8.B	GSM15905.txt	E2.Ral	8	B	E2.Ral.8.B
E2.Ral.48.A	GSM15906.txt	E2.Ral	48	A	E2.Ral.48.A
E2.Ral.48.B	GSM15907.txt	E2.Ral	48	B	E2.Ral.48.B
E2.TOT.8.A	GSM15908.txt	E2.TOT	8	A	E2.TOT.8.A
E2.TOT.8.B	GSM15909.txt	E2.TOT	8	B	E2.TOT.8.B
E2.TOT.48.A	GSM15910.txt	E2.TOT	48	A	E2.TOT.48.A
E2.TOT.48.B	GSM15911.txt	E2.TOT	48	B	E2.TOT.48.B

Obsérvese que los nombres de las columnas de los datos coinciden con los de las filas de targets

Muestra el código

rownames(targets)==colnames(x)

 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[16] TRUE

Análisis del efecto batch

La detección de efectos batch puede hacerse de diversas formas pero esencialmente, de lo que se trata es de ver si las muestras se agrupan por causas distintas a las que se podría esperar, por ejemplo si en vez de parecerse más entre si muestras que han recibido un mismo tratamiento, lo hacen las que han recibido algún tratamiento el mismo día.

La detección de estos efectos puede hacerse mediante distintas técnicas de las que, la más popular es el análisis de componentes principales. Una vez detectado si existe efecto batch es posible mirar de eliminarlo. Si cada tratamiento en estudio está presente en cada lote suele ser posible separar ambos efectos. Si, no es así no suele ser posible evitar cierto grado de tratamiento–lote.

En este ejercicio nos centraremos únicamente en la detección del efecto, no en su eliminación.

Cabe destacar una diferencia entre esta aplicación del PCA y la anterior.

En el caso de los gorriones hemos utilizado el PCA para buscar y explicar variables latentes que nos ayuden a interpretar los datos.
En este caso no pretendemos detectar variables asociadas al efecto batch, es decir no buscaremos interpretarlos sino tan sólo ponerlo en evidencia.

Esto no es porque no sea importante saber que causa, sino porque muchas veces puede tener orígenes desconocidos y sobretodo es plausible que no hayamos recogido información de las variables que lo explican. Esto ha llevado al desarrollo de métodos multivariantes, como el análisis de variables subrogadas del paquete SVA https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3307112/ orientadas a modelizar y quitar el efecto batch globalmente sin buscar su interpretación.

Detectando el efecto batch

La detección de efectos batch puede hacerse de diversas formas pero esencialmente, de lo que se trata es de ver si las muestras se agrupan por causas distintas a las que se podría esperar, por ejemplo si en vez de parecerse más entre si muestras que han recibido un mismo tratamiento, lo hacen las que han recibido cualquier tratamiento el mismo día.

La detección de estos efectos puede hacerse mediante distintas técnicas de las que, la más popular es el análisis de componentes principales.

Una vez detectado si existe efecto batch es posible mirar de eliminarlo. Si cada tratamiento en estudio está presente en cada lote suele ser posible separar ambos efectos usando el análisis de la varianza. Si, no es así no suele ser posible evitar cierto grado de tratamiento-lote.

Muestra el código

colores <- c(rep("black", 4), rep("blue", 4), rep("red", 4), rep("green", 4))

Un diagrama de cajas de todos los datos no muestra ninguna tendencia clara que separe los lotes A y B

Muestra el código

sampleNames <- targets$ShortName
boxplot(as.data.frame(x), cex.axis=0.6, col=colores, las=2, names=sampleNames, 
        main="Signal distribution for selected chips")

Para simplificar la visualización, la función siguiente combina la realización de un PCA, la extracción de las dos primeras componentes y su visualización, debidamente etiquetada, por tratamientos y por lotes.

Muestra el código

plotPCA <- function ( X, labels=NULL, colors=NULL, dataDesc="", scale=FALSE, cex4text=0.8)
{
  pcX<-prcomp(X, scale=scale) # o prcomp(t(X))
  loads<- round(pcX$sdev^2/sum(pcX$sdev^2)*100,1)
  xlab<-c(paste("PC1",loads[1],"%"))
  ylab<-c(paste("PC2",loads[2],"%"))
  if (is.null(colors)) colors=1
  plot(pcX$x[,1:2],xlab=xlab,ylab=ylab, col=colors,  
       xlim=c(min(pcX$x[,1])-10, max(pcX$x[,1])+10),
       ylim=c(min(pcX$x[,2])-10, max(pcX$x[,2])+10))
  text(pcX$x[,1],pcX$x[,2], labels, pos=3, cex=cex4text)
  title(paste("Plot of first 2 PCs for expressions in", dataDesc, sep=" "), cex=0.8)
}

Con esta función podemos realizar y visualizar de forma sencilla el PCA de los datos. Obsérvese que realizamos el PCA de los datos traspuestos puesto que los datos de microarrays se presentan traspuestos a la forma habitual, es decir las filas son las variables y las columnas los individuos.

Muestra el código

sampleNames <- targets$ShortName
plotPCA(t(x), labels=sampleNames, colors=colores, dataDesc="selected samples", cex4text=0.6)

Las dos primeras componentes explican muy poca variabilidad, pero es imediato ver que

La primera componente se asocia bien al factor “time” (valor 8 en general a la derechay valor 48 a la izquierda)
La segunda componente parece asociada al factor “Batch”, las muestras del Batch A, arriba y las del batch “B” abajo.

Esto sugiere que el efcto “Batch” puede enmascarar los efectos de los tratamientos por lo qyue será recomendable eliminarlo para poder estudiar bien el efecto Tratamiento.

En este caso esto será posible, al menos en parte, porque el diseño se encuentra balanceado entre Tratamiento, Tiempo y Batch, por lo que un modelo de Análisis de la Varianza nos permitirá descomponer la variabilidad asociada a cada causa y eliminar la debida al efecto batch.

Análisis basado en distancias

Un enfoque alternativo aunque relacionado es realizar un análisis basado en distancias. Podemos hacerlo calculando la matriz de distancias y visualizándola mediante un mapa de colores o usando escalamiento multidimensional que se discute más adelante.

Muestra el código

manDist <- dist(t(x))
heatmap (as.matrix(manDist), col=heat.colors(16))

Muestra el código

require(MASS)
sam1<-sammon (manDist, trace=FALSE)
plot(sam1$points)
text(sam1$points, targets$Batch, pos=4)

Todas las visualizaciones coinciden en mostrar una separación asociada al factor batch A o B.

Análisis de proximidades

Breve repaso sobre medidas de similaridad

Supongamos que tenemos $n$ variables dicotómicas, basadas en ausencia (-) o presencia (+) de caracteres cualitativos. Un individuo queda entonces caracterizado por la presencia o la ausencias de los caracteres, y dos individuos seran tano más similares cuanto más parecido sea su perfil de presencias-ausencias.

	$x_1$	$x_2$	$x_3$	$\cdots$	$x_n$
$i$	$+$	$-$	$+$	$\cdots$	$+$

Podemos determinar a asociación entre dos individuos $i$, $j$ a partir de la tabla de frecuencias

			$i$
		$+$		$-$
	$+$	$a$		$b$
$j$
	$-$	$c$		$d$

donde: $n = a + b + c + d$, que contiene:

el número de caracteres comunes a $i$ y a $j$, $a$,
el número de caracteres presentes en $i$ pero no en $j$, $b$,
el número de caracteres presentes en $j$ pero no en $i$, $c$,
el número de caracteres no comunes a $i$ y a $j$, $d$,

La asociación se mide mediante un coeficiente de similaridad $s_{ij}$ que verifica, en general:

$0 \leq s_{ij} \leq 1$,
$s_{ij} = 0$ si $c+b = n$ (discrepancia total)
$s_{ij} = 1$ si $a+d = n$ (concordancia total)

La similaridad $s_{ij}$ da el grado de semejanza entre $i$, $j$ en relación a los $n$ caracteres. Ejemplos de coeficientes de similaridad son

Ejemplo 1: Distancias geográficas²

En muchas ocasiones se dispone de datos en forma de una matriz de distancias, de similaridades o de disimilariddes. En algún sentido, que puede variar de un caso a otro, estas matrices cuantifican el parecido, o la diferencia, o incluso la distancia geográfica por parejas de un grupo de variables. Nuestro objetivo, en estos casos, suele ser intentar recuperar, a partrir de las maytrices anteriores, las coordenadas de los puntos, de forma que podamos visualizarlos mostrando las relaciones entre ellos.

Este ejemplo utiliza las distancias aéreas, en millas, entre 10 ciudades de Estados Unidos. A partir de éstas deseamos recuperar las coordenadas de cada dciudad de forma que podamos representarlas en un plano que refleje lo mejor posible las distancias entre ellas. Hemos de tener en cuenta que puesto que se trata de duistancias aéreas, no son distancias euclídeas (los aviones no se desplazan en línea recta) por lo que la recuperación tan ssólo será aproximada.

Muestra el código

library(readxl)
DistSim <- as.data.frame(read_excel("datos/Dist_Sim_Disim.xlsx", sheet = "GeoDistances"))
rownames(DistSim) <- DistSim[,1]
DistSim <-as.matrix(DistSim [,-1])

Muestra el código

library(magrittr)
DistSim %>%
  kableExtra::kbl() %>%
  kableExtra::kable_paper("hover", full_width = F, font_size = 9)

	Atl	Chi	Den	Hou	LA	Mia	NY	SF	Sea	Was
Atlanta	0	NA	NA	NA	NA	NA	NA	NA	NA	NA
Chicago	587	0	NA	NA	NA	NA	NA	NA	NA	NA
Denver	1212	920	0	NA	NA	NA	NA	NA	NA	NA
Houston	701	940	879	0	NA	NA	NA	NA	NA	NA
Los Angeles	1936	1745	831	1734	0	NA	NA	NA	NA	NA
Miami	604	1188	1726	968	2339	0	NA	NA	NA	NA
New York	748	713	1631	1420	2451	1092	0	NA	NA	NA
San Francisco	2139	1858	949	1645	347	2594	2571	0	NA	NA
Seattle	2182	1737	1021	1891	959	2734	2408	678	0	NA
Washington	543	597	1494	1220	2300	923	205	2442	2329	0

Nuestro objetivo es lograr una visualización de los datos que represente lo mejor posible las distancias entre éstas, aunque no conozcamos, como en este caso, las coordenadas de cada ciudad. Además, puesto que se trata de distancias aéreas, no será posible recuperar directamente las distancias euclídeas entre las ciudades, por lo que procederemos a realizar un escalamiento multidimensional clásico

Siguiendo a @Everitt2011, capítulo 4, procederemos a:

Diagonalizar la matriz de distancias.
Obtener las coordenadas de cada individuo (ciudad).

\[ X = V_1\cdot \Lambda_1^{\frac 12} , \] donde $V_1$ contiene los vectores propios no nulos de $B$ y \[ \Lambda_1^{\frac 12}= \mbox{diag}(\lambda_1^{\frac 12}, \lambda_2^{\frac 12},..., \lambda_q^{\frac 12}), \] contiene la raíz cuadrada de los primeros $q$ valores propios no nulos de $B$.

Para realizar un multidimensonal scaling clásico nos basaremos en la instrucción cmdscale del paquete stats:

Muestra el código

airdist <- as.dist(as.matrix(DistSim))
airline_mds <- cmdscale(airdist, k = 9, eig = TRUE)
airline_mds$points

                    [,1]       [,2]        [,3]        [,4]         [,5]
Atlanta        -712.4178  157.34246   67.519636  33.9648977  -8.45853360
Chicago        -374.1124 -332.99012    2.928936  28.7749927 -15.84330898
Denver          486.6607   14.78013  -41.304984  45.9929928  14.75126177
Houston        -228.9426  515.95616 -465.504161  -0.5668829   0.04794259
Los Angeles    1219.0246  277.99731  458.396262  -0.1879388   0.02297080
Miami         -1128.4262  598.95103  142.601024 -32.8128340   1.93168561
New York      -1062.3316 -534.19580   39.122338 -35.6808373   0.22573428
San Francisco  1423.4269  181.43505 -110.535264 -24.0833935  -4.86588885
Seattle        1347.8203 -536.47617 -140.929698 -18.5126346   0.87459912
Washington     -970.7018 -342.80004   47.705911   3.1116379  11.31353725

Si inspeccionamos los valores propios veremos que algunos son negativos:

Muestra el código

lam <- airline_mds$eig
show(lam)

 [1]  9.609694e+06  1.561695e+06  4.893155e+05  7.379541e+03  6.963801e+02
 [6] -8.149073e-10 -4.056578e+02 -9.605809e+02 -2.371259e+04 -2.945694e+05

Esto se explica porque, tal como se ha comentado, las distancias no son euclídeas. Para determinar si podemos prescindir de algunas componentes calculamos, como hacíamos en el análisis de componentes principales, el porcentaje de la suma de los valores propios (en valor absoluto o al cuadrado) que representa cada valor propio individual.

Muestra el código

cumsum(abs(lam)) / sum(abs(lam))

 [1] 0.8015807 0.9318476 0.9726633 0.9732788 0.9733369 0.9733369 0.9733708
 [8] 0.9734509 0.9754289 1.0000000

Muestra el código

cumsum(lam^2) / sum(lam^2)

 [1] 0.9709215 0.9965638 0.9990812 0.9990818 0.9990818 0.9990818 0.9990818
 [8] 0.9990818 0.9990877 1.0000000

El elevado valor del primer término en ambos casos, muestra como, bastará con las dos primeras coordenadas para obtener una buena representación.

Muestra el código

lim <- range(airline_mds$points[,1] * (-1)) * 1.2
plot(airline_mds$points[,1] * (-1), airline_mds$points[,2] * (-1),
     type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2",
     xlim = lim, ylim = lim)
text(airline_mds$points[,1] *(-1), airline_mds$points[,2] * (-1), 
     labels(airdist), cex = 0.7)

Como puede verse en el gráfico la primera coordenada representa la dirección Este Oeste, y la seguna la Norte-Sur. Las posiciones de los puntos reproducen relativamente bien las posiciones relativas de las ciudades en el mapa.

² El ejemplo de la visualización de distancias entre ciudades se ha tomado del libro de @Everitt2011 An introduction to aplied multivariate analysis using R y del repositorio de github en donde se encuentra el código R de dicho libro: https://rdrr.io/cran/MVA/

³ El ejemplo de la relación entre partidos políticos se ha tomado del documento @Cuadras2007, Exemples d’Anàlisi Multivariant del Dr. Carles Cuadras disponible en su página web.

Ejemplo 2: Similaridades entre partidos políticos³

La tabla siguiente contiene una matriz de distancias entre partidos políticos españoles en la primera década del siglo XXI obtenida a partir de 11 características sociológicas y económicas (sobre constitución del estado, poder judicial y ejecutivo, financiación autonómica, etc.) aplicando la similaridad de Jaccard. Las distancias estan elevadas al cuadrado.

Muestra el código

require(readxl)
simPartidos <- as.data.frame(read_excel("datos/Dist_Sim_Disim.xlsx", sheet = "SimilarJaccard-Partidos"))
rownames(simPartidos) <- simPartidos[,1]
simPartidos <-as.matrix(simPartidos [,-1])

Muestra el código

library(magrittr)
simPartidos %>%
  kableExtra::kbl() %>%
  kableExtra::kable_paper("hover", full_width = F, font_size = 9)

	PP	PSOE	CU	IU	ERC	ICV
PP	0.000	NA	NA	NA	NA	NA
PSOE	1.400	0.000	NA	NA	NA	NA
CU	1.714	2.000	0.000	NA	NA	NA
IU	1.273	0.444	1.800	0.000	NA	NA
ERC	2.000	1.250	1.600	1.111	0.000	NA
ICV	1.818	1.111	1.714	0.667	0.667	0

Para trabajar con los datos los convertimos en una matriz simétrica:

Muestra el código

makeSymm <- function(m) {
   m[upper.tri(m)] <- t(m)[upper.tri(m)]
   return(m)
}
distPartidos <- makeSymm(as.matrix(simPartidos))
show(distPartidos)

        PP  PSOE    CU    IU   ERC   ICV
PP   0.000 1.400 1.714 1.273 2.000 1.818
PSOE 1.400 0.000 2.000 0.444 1.250 1.111
CU   1.714 2.000 0.000 1.800 1.600 1.714
IU   1.273 0.444 1.800 0.000 1.111 0.667
ERC  2.000 1.250 1.600 1.111 0.000 0.667
ICV  1.818 1.111 1.714 0.667 0.667 0.000

En este caso tenemos una matriz de distancias $\Delta=\delta_{ij}$. Siguiendo a @Cuadras2019, cap 8, podemos transformarla en una matriz de similaridad mediante la transformación: $S_{ij} = -\frac 12 \delta_{ij}^2$.

En este caso se nos indica que las distancias ya estan al cuadrado por lo que tan sólo multiplicaremos por $-1/2$:

Muestra el código

simPartidos <-  S<- -0.5 *distPartidos

Una vez tenemos la matriz de similaridades $S$ la transformaremos en \[ S^*=s_{ij}^* = s_{ij} - \overline{s}_{i}- \overline{s}_{j} + \overline{s}, \]

Esto puede hacerse multiplicando por la izquierda y poir la derecha por la matriz $H$, es decir: \[ S*= HSH, \mbox{ donde $H$ es la matriz de centrado de datos: } H = I_n- 1_n 1_n'. \]

Aquí:

Muestra el código

n = nrow(S)
ones = rep(1, n)
H = diag(ones) - (1/n) * (ones %*% t(ones))
Sprime <- H %*% S %*% H

A continuación diagonalizamos esta matriz:

Muestra el código

diagS <- eigen(Sprime)
diagS

eigen() decomposition
$values
[1] 1.272211e+00 1.122880e+00 5.198830e-01 3.642334e-01 1.489596e-01
[6] 1.554312e-15

$vectors
           [,1]       [,2]       [,3]        [,4]         [,5]       [,6]
[1,]  0.4716114 -0.5709237  0.5305741 -0.03758486  0.045174070 -0.4082483
[2,] -0.2713463 -0.3554746 -0.5219773 -0.38458580  0.461493025 -0.4082483
[3,]  0.6693854  0.4654397 -0.4027922  0.07985754  0.001916061 -0.4082483
[4,] -0.2682520 -0.2465235 -0.2559949  0.32458983 -0.727810701 -0.4082483
[5,] -0.2614706  0.4537676  0.3837563 -0.59891639 -0.230416406 -0.4082483
[6,] -0.3399280  0.2537145  0.2664340  0.61663969  0.449643951 -0.4082483

Las coordenadas principales seran:

Muestra el código

U <- diagS$vectors
Dhalf <- sqrt(diagS$values)
coords<- U %*% diag(Dhalf)[,1:2]

Y la representación gráfica de los partidos según las dos primeras CP es:

Muestra el código

xlim <- range(coords[,1] * (-1)) * 1.2
ylim <- range(coords[,2] * (-1)) * 1.2
plot(coords[,1] * (-1), coords[,2] * (-1),
     type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2",
     xlim = xlim, ylim = ylim)
text(coords[,1] * (-1), coords[,2] * (-1), 
     rownames(simPartidos), cex = 0.7)

Como puede verse surge una interpretación bastante natural, en tanto que el primer eje representa si el partido es (pretendidamente) de izquierdas o de derecha, ientrtas que el segundo eje alinea partidos nacionalistas en la parte superior y partidos de ámbito más estatal en la inferior. Como era de esperar, en el centro no hay nadie.

Análisis de correspondencias

El análisis de correspondencias es una técnica que resulta especialmente adecuada para analizar y visualizar datos de tablas de contingencia con datos de frecuencias numéricas.
Como en el caso del PCA y el MDS este análisis proporciona una representación de datos en dimensión reducida que facilita la interpretación y la comprensión de los datos.
En este caso, además, es posible una representación simultánea muy intuitiva de los atributos de las filas y las columnas, lo que, en las circunstancias adecuadas, permite visualizar el grado de asociación entre unas y otras.

Ejemplo 1: Representación de mutaciones cromosómicos entre poblaciones

Este primer ejemplo es un análisis clásico basado en estudios de genética de poblaciones de los años 70. En primer lugar se presenta la realización del análisis “empaquetado” para después ilustrarl su realización paso a paso.

@Alonso1975 hace un amplio estudio de la distribución geográfica del polimorfismo cromosómico de Drosophila subobscura utilizando análisis factorial de correspondencias.

Sobre una tabla de frecuencias de 66 poblaciones y $8 + 3 + 7 + 23 + 11$ ordenaciones de los $5$ cromosomas $A$, $J$, $E$, $U$ y $O$ respectivamente, hace un análisis de correspondencias global, y varios análisis parciales. Utilizaremos, como ilustración uno de los análisis parciales, concretamente la que coge $13$ poblaciones y $3$ inversiones del cromosoma $A$. Los datos, se dan en forma de porcentaje en la Tabla 1.

TABLA.- Porcentajes de frecuencias de tres inversiones del cromosoma A para 13 poblaciones de Drosophila subobscura:

Muestra el código

kableExtra::kable(datos)

	A-ST	A-1	A-2
HELSINKI	96.0	4.0	0.0
DROBACK	78.4	16.2	5.4
HARIOT	100.0	0.0	0.0
DALKEITH	100.0	0.0	0.0
GRONINGEN	80.0	16.0	4.0
FONTAINEBLEAU	88.5	7.7	3.8
VIENA	56.9	35.8	7.4
ZURICH	67.8	24.4	7.8
FRUSKA-GORA	36.1	55.6	8.3
LAGRASSE	72.5	17.5	10.0
MONTPELLIER	60.2	24.3	15.5
CARASCO	50.0	31.8	18.2

Obsérvese como un gráfico sencillo puede dar buena idea de qué relaciones entre filas y columnas esperamos encontrar:

Muestra el código

library("gplots")
# 1. convert the data as a table
dt <- as.table(datos)
# 2. Graph
balloonplot(t(dt), main ="Polimorfismos en Europa", xlab ="Cromosoma A", ylab="Ciudades",
            label = FALSE, show.margins = FALSE)

Si realizamos el análisis usando la función CA del paquete FactoMineR obtenemos los siguientes resultados:

Muestra el código

library(FactoMineR)
datos.ca <- CA(datos, graph = FALSE)
plot.CA(datos.ca)

Observamos que Heriot y Dalkeith quedan representadas en un mismo punto, dada su distribución idéntica. También queda reflejada la similaridad entre Drobak y Fontainebleau. La población Fruska-Gora queda apartada del resto, debido a la influencia de la ordenación $A-1$, menos frecuente en las demás poblaciones. Las proporciones de las 3 ordenaciones quedan muy bien reflejadas en la gráfica; se ve que las poblaciones quedan más cercanas de las ordenaciones que se presentan más.

Realización del análisis paso por paso

Aunque, en la práctica, solemos utilizar librerías de R que “empaquetan” en funciones propias las transformaciones y análisis de los datos, resulta instructivo realizar algún análisis paso a paso, ya sea para tener una mejor comprensión del análisis realizado o, incluso, para saber cuando puede ser preciso buscar alternativas porque se dan condiciones particulares.

En esta sección realizaremos paso a paso el análisis de los datos anteriores comparando los resultados que obtenemos con los proporcionados por el paquete utilizado.

Empezamos transformando los datos en frecuencias relativas y calculando las frecuencias relativas marginales.

Muestra el código

tabla.N <- datos
N <- sum(tabla.N)

El análisis de esta tabla de frecuencias mediante análisis de correspondencias puede plantearse de dos formas:

1- Como un escalado multidimensional de filas y de columnas con la distancia ji-cuadrado. 2- Como un análisis de componentes principales sobre la matriz Z estandarizada.

Como un escalado multidimensional

Calculamos las tablas de frecuencias relativas y marginales:

Muestra el código

tabla.F <-   tabla.N/N
margin.f <- margin.table(tabla.F,1)
margin.c <- margin.table(tabla.F,2)
show(round(addmargins(tabla.F),4))

                A-ST    A-1    A-2    Sum
HELSINKI      0.0800 0.0033 0.0000 0.0833
DROBACK       0.0653 0.0135 0.0045 0.0833
HARIOT        0.0833 0.0000 0.0000 0.0833
DALKEITH      0.0833 0.0000 0.0000 0.0833
GRONINGEN     0.0667 0.0133 0.0033 0.0833
FONTAINEBLEAU 0.0737 0.0064 0.0032 0.0833
VIENA         0.0474 0.0298 0.0062 0.0834
ZURICH        0.0565 0.0203 0.0065 0.0833
FRUSKA-GORA   0.0301 0.0463 0.0069 0.0833
LAGRASSE      0.0604 0.0146 0.0083 0.0833
MONTPELLIER   0.0502 0.0202 0.0129 0.0833
CARASCO       0.0417 0.0265 0.0152 0.0833
Sum           0.7386 0.1944 0.0670 1.0000

La solución mediante escalado multidimensional pasa por calcular la matriz de distancias jui-cuadrado por filas y/o por columnas y, sobre esta matriz, realizamos un escalado multidimensional. Es decir realizamos dos cálculos separados que podremos superponer después

Análisis de la tabla basado en los perfiles “fila”

A partir de la tabla de frecuencias relativas podemos calcular la tabla de frecuencias relativas condicionada por filas (haciendo que éstas sumen 1).

Muestra el código

Y<- tabla.P <- diag(1/margin.f) %*% tabla.F
show(round(addmargins(Y),4))

      A-ST    A-1    A-2 Sum
    0.9600 0.0400 0.0000   1
    0.7840 0.1620 0.0540   1
    1.0000 0.0000 0.0000   1
    1.0000 0.0000 0.0000   1
    0.8000 0.1600 0.0400   1
    0.8850 0.0770 0.0380   1
    0.5684 0.3576 0.0739   1
    0.6780 0.2440 0.0780   1
    0.3610 0.5560 0.0830   1
    0.7250 0.1750 0.1000   1
    0.6020 0.2430 0.1550   1
    0.5000 0.3180 0.1820   1
Sum 8.8634 2.3326 0.8039  12

Calculamos la matriz de distancias ji cuadrado entre las filas. Obsérvese que para el cálculo de la distancia entre dos perfiles filas aplicamos la fórmula:

\[\begin{equation} D^{2}\left(\mathbf{p}_{u}, \mathbf{p}_{v}\right)=\left(\mathbf{p}_{u}-\mathbf{p}_{v}\right)^{\prime} \mathbf{D}_{c}^{-1}\left(\mathbf{p}_{u}-\mathbf{p}_{v}\right) \end{equation}\]

Muestra el código

nf <- nrow(tabla.F)
D2.chisq.f <- matrix(0,nf,nf)
for(i in 1:(nf-1))
  for(j in i:nf)
    D2.chisq.f[i,j] <-  t(Y[i,]-Y[j,]) %*% diag(1/margin.c) %*% (Y[i,]-Y[j,])
   D2.chisq.f <- D2.chisq.f + t(D2.chisq.f)
rownames(D2.chisq.f) <- colnames(D2.chisq.f) <- rownames(tabla.F)
show(round(D2.chisq.f,6))

              HELSINKI  DROBACK   HARIOT DALKEITH GRONINGEN FONTAINEBLEAU
HELSINKI      0.000000 0.162028 0.010397 0.010397  0.132616      0.036212
DROBACK       0.162028 0.000000 0.241693 0.241693  0.003293      0.054798
HARIOT        0.010397 0.241693 0.000000 0.000000  0.209726      0.069958
DALKEITH      0.010397 0.241693 0.000000 0.000000  0.209726      0.069958
GRONINGEN     0.132616 0.003293 0.209726 0.209726  0.000000      0.045279
FONTAINEBLEAU 0.036212 0.054798 0.069958 0.069958  0.045279      0.000000
VIENA         0.808178 0.265734 0.991703 0.991703  0.290720      0.560091
ZURICH        0.412555 0.058399 0.537446 0.537446  0.078002      0.225358
FRUSKA-GORA   1.958238 1.053343 2.245859 2.245859  1.095190      1.582225
LAGRASSE      0.317785 0.037167 0.409191 0.409191  0.062509      0.141441
MONTPELLIER   0.744114 0.230863 0.876825 0.876825  0.285920      0.454512
CARASCO       1.178466 0.478943 1.353089 1.353089  0.551247      0.808971
                 VIENA   ZURICH FRUSKA-GORA LAGRASSE MONTPELLIER  CARASCO
HELSINKI      0.808178 0.412555    1.958238 0.317785    0.744114 1.178466
DROBACK       0.265734 0.058399    1.053343 0.037167    0.230863 0.478943
HARIOT        0.991703 0.537446    2.245859 0.409191    0.876825 1.353089
DALKEITH      0.991703 0.537446    2.245859 0.409191    0.876825 1.353089
GRONINGEN     0.290720 0.078002    1.095190 0.062509    0.285920 0.551247
FONTAINEBLEAU 0.560091 0.225358    1.582225 0.141441    0.454512 0.808971
VIENA         0.000000 0.082935    0.261880 0.214932    0.167245 0.188767
ZURICH        0.082935 0.000000    0.637165 0.034706    0.096325 0.232512
FRUSKA-GORA   0.261880 0.637165    0.000000 0.930412    0.659970 0.463832
LAGRASSE      0.214932 0.034706    0.930412 0.000000    0.089422 0.274098
MONTPELLIER   0.167245 0.096325    0.659970 0.089422    0.000000 0.053903
CARASCO       0.188767 0.232512    0.463832 0.274098    0.053903 0.000000

Finalmente realizamos escalado multidimensional sobre la matriz de distancias entre filas.

Muestra el código

mds.f <- cmdscale(sqrt(D2.chisq.f),eig=TRUE)
mds.f

$points
                     [,1]        [,2]
HELSINKI      -0.49930839 -0.08153341
DROBACK       -0.10258763 -0.01341023
HARIOT        -0.59336881 -0.04217206
DALKEITH      -0.59336881 -0.04217206
GRONINGEN     -0.13576678 -0.06022835
FONTAINEBLEAU -0.33501017  0.01447700
VIENA          0.39799723 -0.13649197
ZURICH         0.13905244 -0.01046753
FRUSKA-GORA    0.88289376 -0.30006242
LAGRASSE       0.02154641  0.13409481
MONTPELLIER    0.29331962  0.25884055
CARASCO        0.52460114  0.27912567

$eig
 [1]  2.413616e+00  2.858896e-01  3.877601e-16  3.737291e-17 -4.164152e-18
 [6] -1.369120e-17 -2.166838e-17 -3.447503e-17 -4.854955e-17 -5.148866e-17
[11] -1.013651e-16 -1.593154e-16

$x
NULL

$ac
[1] 0

$GOF
[1] 1 1

El resultado del MDS se puede representar para visualizar la relación entre las filas.

Muestra el código

plot(mds.f$points,ty="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(mds.f$points[,1],mds.f$points[,2],labels=substr(rownames(tabla.N),1,6),cex=0.6)

Análisis de la tabla basado en los perfiles “columna”

Muestra el código

Y<- tabla.C <-   tabla.F %*% diag(1/margin.c)
show(round(addmargins(Y),4))

                                      Sum
HELSINKI      0.1083 0.0171 0.0000 0.1254
DROBACK       0.0884 0.0694 0.0672 0.2251
HARIOT        0.1128 0.0000 0.0000 0.1128
DALKEITH      0.1128 0.0000 0.0000 0.1128
GRONINGEN     0.0903 0.0686 0.0498 0.2086
FONTAINEBLEAU 0.0998 0.0330 0.0473 0.1801
VIENA         0.0642 0.1535 0.0920 0.3097
ZURICH        0.0765 0.1046 0.0970 0.2781
FRUSKA-GORA   0.0407 0.2383 0.1032 0.3823
LAGRASSE      0.0818 0.0750 0.1244 0.2812
MONTPELLIER   0.0679 0.1042 0.1928 0.3649
CARASCO       0.0564 0.1363 0.2264 0.4191
Sum           1.0000 1.0000 1.0000 3.0000

A continuación calculamos la matriz de distancias ji cuadrado entre las columnas.

Muestra el código

nc <- ncol(tabla.F)
D2.chisq.c <- matrix(0,nc,nc)
for(i in 1:(nc-1))
  for(j in i:nc)
    D2.chisq.c[i,j] <-  t(Y[,i]-Y[,j]) %*% diag(1/margin.f) %*% (Y[,i]-Y[,j])
   D2.chisq.c <- D2.chisq.c + t(D2.chisq.c)
rownames(D2.chisq.c) <- colnames(D2.chisq.c) <- colnames(tabla.F)
show(round(D2.chisq.c,6))

         A-ST      A-1      A-2
A-ST 0.000000 1.135265 1.121353
A-1  1.135265 0.000000 0.496044
A-2  1.121353 0.496044 0.000000

Finalmente realizamos escalado multidimensional sobre la matriz de distancias entre columnas.

Muestra el código

mds.c <- cmdscale(sqrt(D2.chisq.c),eig=TRUE)
mds.c

$points
           [,1]         [,2]
A-ST  0.6680866 -0.003873443
A-1  -0.3395554 -0.350171946
A-2  -0.3285312  0.354045388

$eig
[1]  6.695703e-01  2.479835e-01 -3.330669e-16

$x
NULL

$ac
[1] 0

$GOF
[1] 1 1

El resultado del MDS se puede representar para visualizar la relación entre las filas.

Muestra el código

plot(mds.c$points,ty="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(mds.c$points[,1],mds.c$points[,2],labels=substr(colnames(tabla.F),1,6),cex=0.6)

Finalmente observemos que podemos superponer los dos gráficos:

Muestra el código

xinf<- min(min(mds.c$points[,1]), min(mds.f$points[,1]))
xsup <- max(max(mds.c$points[,1]), max(mds.f$points[,1]))
yinf<- min(min(mds.c$points[,2]), min(mds.f$points[,2]))
ysup <- max(max(mds.c$points[,2]), max(mds.f$points[,2]))
xdelta<- (xsup-xinf)/10
ydelta<- (ysup-yinf)/10
xLim <- c(xinf-xdelta, xsup+xdelta)
yLim <- c(yinf-ydelta, ysup+ydelta)

plot(mds.c$points,ty="n",xlab="PC1",ylab="PC2", xlim=xLim, ylim=yLim)
abline(v=0,h=0, col="gray",lty=4)
text(mds.f$points[,1],mds.f$points[,2],labels=substr(rownames(tabla.F),1,6),cex=0.6, col="blue")
text(mds.c$points[,1],mds.c$points[,2],labels=substr(colnames(tabla.F),1,6),cex=0.6, col="red")

AC como un ACP sobre la matriz estandarizada

El análisis se realiza sobre la matriz $\mathbf{Z}$:

\[\begin{equation} \mathbf{Z}=\mathbf{D}_{f}^{1 / 2} \mathbf{Y}=\mathbf{D}_{f}^{-1 / 2} \mathbf{F} \mathbf{D}_{c}^{-1 / 2}=\left(\frac{f_{i j}}{\sqrt{f_{i} f_{\cdot j}}}\right) \end{equation}\]

Sin embargo y para evitar que tengamos un valor propio 1, inútil ya que es una solución trivial, es mejor estandarizar la matriz de correspondencias y utilizar la matriz

\[\begin{equation} \mathbf{Z}=\mathbf{D}_{f}^{-1 / 2}\left(\mathbf{F}-\mathbf{f} \mathbf{c}^{\prime}\right) \mathbf{D}_{c}^{-1 / 2} \end{equation}\]

Muestra el código

Dfmh <- diag(1/sqrt(margin.f))
Dcmh <- diag(1/sqrt(margin.c))
Z <-  Dfmh %*% (tabla.F -margin.f%*% t(margin.c)) %*% Dcmh
rownames(Z) <- rownames(tabla.F)
colnames(Z) <- colnames(tabla.F)
round(Z,4)

                 A-ST     A-1     A-2
HELSINKI       0.0744 -0.1011 -0.0747
DROBACK        0.0152 -0.0212 -0.0145
HARIOT         0.0878 -0.1273 -0.0747
DALKEITH       0.0878 -0.1273 -0.0747
GRONINGEN      0.0206 -0.0225 -0.0301
FONTAINEBLEAU  0.0492 -0.0769 -0.0323
VIENA         -0.0572  0.1069  0.0077
ZURICH        -0.0204  0.0325  0.0123
FRUSKA-GORA   -0.1268  0.2367  0.0179
LAGRASSE      -0.0046 -0.0127  0.0368
MONTPELLIER   -0.0459  0.0318  0.0981
CARASCO       -0.0801  0.0809  0.1283

Realizaremos el ACP mediante la descomposición en valores singulares de Z:

Muestra el código

Z.svd <- svd(Z)
Z.svd

$d
[1] 4.484765e-01 1.543493e-01 9.676243e-17

$u
             [,1]        [,2]        [,3]
 [1,] -0.32140123  0.15248587  0.05827654
 [2,] -0.06605194  0.02506336 -0.60314261
 [3,] -0.38194412  0.07887700 -0.37662848
 [4,] -0.38194412  0.07887700 -0.37662848
 [5,] -0.08740702  0.11262391 -0.11173040
 [6,] -0.21565180 -0.02708007  0.12161619
 [7,]  0.25628095  0.25535378 -0.09044851
 [8,]  0.08948031  0.01954840 -0.33124219
 [9,]  0.56826058  0.56111153 -0.21770114
[10,]  0.01384486 -0.25080562 -0.22405817
[11,]  0.18877038 -0.48411733 -0.01438563
[12,]  0.33763503 -0.52206548 -0.32579981

$v
           [,1]        [,2]      [,3]
[1,] -0.5099814  0.03624682 0.8594214
[2,]  0.7751769  0.45243825 0.4409087
[3,]  0.3728536 -0.89105881 0.2588328

De aquí se obtienen las coordenadas principales y estándar:

Muestra el código

f.sc <- diag(1/sqrt(margin.f)) %*% Z.svd$u
c.sc <- diag(1/sqrt(margin.c)) %*% Z.svd$v
f.pc <- f.sc %*% diag(Z.svd$d)
c.pc <- c.sc %*% diag(Z.svd$d)

Podemos finalmente representar las distancias con filas y columnas de Z

Muestra el código

library(MASS)
# solución simétrica
eqscplot(f.pc[,1:2],type="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(f.pc[,1],f.pc[,2],labels=rownames(tabla.F),cex=0.8,font=2,col="blue")
text(c.pc[,1],c.pc[,2],labels=colnames(tabla.F),cex=0.8,font=2,col= "red")
title(main="Solución simétrica",line=1)

Muestra el código

#
# solución asimétrica
eqscplot(c.sc[,1:2],type="n",xlab="PC1",ylab="PC2")
abline(v=0,h=0, col="gray",lty=4)
text(f.pc[,1],f.pc[,2],labels=rownames(tabla.F),cex=0.8,font=2,col="blue")
text(c.sc[,1],c.sc[,2],labels=colnames(tabla.F),cex=0.8,font=2,col= "red")
title(main="Solución asimétrica",line=1)

Ejemplo 2: Análisis de correspondencias de datos de microarrays

El AC se presenta tradicionalmente como una técnica para el análisis de tablas de frecuencias, pero estríctamente hablando, si podemos expresar unos datos cuantitativos como perfiles fila y perfiles columna, tambien se puede aplicar a otros tipos de datos.

Un ejemplo de esto fue la aplicación a principios de este siglo a datos de microarrays por @Fellenberg2001. El interes de esta aplicación reside sobretodo en la posibilidad de visualizar simultaneamente expresión de genes e individuos/condiciones experimentales.

El ejemplo se lleva a cabo usando el paquete de R/Bioconductor made4 que es una extensión, para ser usada en datos ómicos, del paquete ade4, un programa muy popular en francia para el análisis exploratorio de datos.

El paquete made4 (@Culhane2005) contiene algunos datasets de ejemplo. Usaremos el subconjunto de datos data.train del dataset khan. La función overview permite una rápida visualización de un dataset dado.

Muestra el código

library(made4)
data("khan")
names(khan)

[1] "train"                "test"                 "train.classes"       
[4] "test.classes"         "annotation"           "gene.labels.imagesID"
[7] "cellType"

Muestra el código

k.data<-khan$train
k.class<-khan$train.classes
overview(k.data)

EL paquete contiene la función ord que permite realizar distintos métodos de reducción de la dimensión cambiando un argumento. Para análisis de correspondencias es coa. Para saber más puede hacerse help (ord).

Muestra el código

k.coa<- ord(k.data, type="coa")
show(k.coa)

$ord
Duality diagramm
class: coa dudi
$call: dudi.coa(df = data.tr, scannf = FALSE, nf = ord.nf)

$nf: 63 axis-components saved
$rank: 63
eigen values: 0.1713 0.1383 0.1032 0.05995 0.04965 ...
  vector length mode    content       
1 $cw    64     numeric column weights
2 $lw    306    numeric row weights   
3 $eig   63     numeric eigen values  

  data.frame nrow ncol content             
1 $tab       306  64   modified array      
2 $li        306  63   row coordinates     
3 $l1        306  63   row normed scores   
4 $co        64   63   column coordinates  
5 $c1        64   63   column normed scores
other elements: N 

$fac
NULL

attr(,"class")
[1] "coa" "ord"

La forma más sencilla de ver los resultados producidos por ord es utilizar plot. plot (k.ord) dibujará un gráfico de los valores propios, junto con los gráficos de las variables (genes) y un gráfico de los casos (muestras de microarrays). En este ejemplo, las muestras de microarrays están codificadas por colores utilizando la variable k.classes que se guarda como k.class.

Muestra el código

plot(k.coa, classvec=k.class, genecol="grey3")

La imagen muestra arriba a la izquierda (A) el gráfico de los valores propios, lo que permite hacerse una idea de cuantas dimensiones retener.
Arriba a la derecha (B) se muestra la proyección de las muestras de microarrays de pacientes coloreadas por tipos de tumores EWS (rojo), BL (azul), NB (verde) o RMS (marrón)
Abajo a la izquierda (C) se muestra la proyección de los genes (círculos grises rellenos).
Finalmente, abajo a la derecha (D) se muestra el biplot con los genes superpuestos con la muestras.

Las muestras y genes con una fuerte asociado se proyectan en la misma dirección desde el origen. Cuanto mayor sea la distancia desde el origen, más fuerte será la asociación.

Análisis de conglomerados (“Cluster Analysis”)

El análisis de conglomerados pertenece al conjunto de técnicas conocidas com “análisis no supervisado” porque su objetivo es descubrir grupos en los datos, es decir grupos de individuos (o eventualmente de variables) que pueden considerarse más similares entre ellos que a los individuos de los otros grupos, según el criterio de similaridad que se decida utilizar.

Presentaremos dos ejemplos distintos:

Un estudio sobre como se agrupan los países europeos por su consumo de carne
Un análisis sobre la clasificación de los estados del ciclo celular y los genes que influyen en cada uno de ellos basado en un estudio de expresión génica.

Agrupación de genes y de muestras en el análisis del ciclo celular

Clustering y visualización

En este caso se presentan distintos métodos de agrupación ("clustering") disponibles en R y Bioconductor. Se utilizará el conjunto de datos yeast generado por el proyecto "Yeast Cell Cycle Project" (Proyecto de estudio del ciclo celular de la levadura).

Estos datos ya no se encuentran disponibles en la web, pero se han descargado para su uso, al subidrectorio datos de este proyecto en github https://github.com/ASPteaching/AMVCasos

Pre-procesamiento de los datos

La agrupación de datos se suele aplicar después del pre-procesamiento y filtrado de los mismos. El conjunto de datos de ejemplo ya ha sido normalizado y filtrado por lo que no nos ocuparemos de este aspecto.

El archivo de datos contiene información de varios experimentos (factor alfa, cdc15 y cinéticas de decantación). Para simplificar en este ejercicio nos ocuparemos únicamente de los datos denominados "cdc15" ("cell division control protein 15").

EL primer paso consiste en extraer los valores del grupo cdc15 del conjunto de datos pre-normalizados y eliminar los genes con valores faltantes:

Muestra el código

workingDir <- getwd()
dataDir <- file.path(workingDir, "datos")
d <- read.table(file.path(dataDir, "combined.txt"), sep="\t", header=T)
names(d)

 [1] "X"         "cln3.1"    "cln3.2"    "clb"       "clb2.2"    "clb2.1"   
 [7] "alpha"     "alpha0"    "alpha7"    "alpha14"   "alpha21"   "alpha28"  
[13] "alpha35"   "alpha42"   "alpha49"   "alpha56"   "alpha63"   "alpha70"  
[19] "alpha77"   "alpha84"   "alpha91"   "alpha98"   "alpha105"  "alpha112" 
[25] "alpha119"  "cdc15"     "cdc15_10"  "cdc15_30"  "cdc15_50"  "cdc15_70" 
[31] "cdc15_80"  "cdc15_90"  "cdc15_100" "cdc15_110" "cdc15_120" "cdc15_130"
[37] "cdc15_140" "cdc15_150" "cdc15_160" "cdc15_170" "cdc15_180" "cdc15_190"
[43] "cdc15_200" "cdc15_210" "cdc15_220" "cdc15_230" "cdc15_240" "cdc15_250"
[49] "cdc15_270" "cdc15_290" "cdc28"     "cdc28_0"   "cdc28_10"  "cdc28_20" 
[55] "cdc28_30"  "cdc28_40"  "cdc28_50"  "cdc28_60"  "cdc28_70"  "cdc28_80" 
[61] "cdc28_90"  "cdc28_100" "cdc28_110" "cdc28_120" "cdc28_130" "cdc28_140"
[67] "cdc28_150" "cdc28_160" "elu"       "elu0"      "elu30"     "elu60"    
[73] "elu90"     "elu120"    "elu150"    "elu180"    "elu210"    "elu240"   
[79] "elu270"    "elu300"    "elu330"    "elu360"    "elu390"

Muestra el código

cdc15 <- which(substr(names(d),1,6)=="cdc15_")
dat <- d[,cdc15]
names(dat)

 [1] "cdc15_10"  "cdc15_30"  "cdc15_50"  "cdc15_70"  "cdc15_80"  "cdc15_90" 
 [7] "cdc15_100" "cdc15_110" "cdc15_120" "cdc15_130" "cdc15_140" "cdc15_150"
[13] "cdc15_160" "cdc15_170" "cdc15_180" "cdc15_190" "cdc15_200" "cdc15_210"
[19] "cdc15_220" "cdc15_230" "cdc15_240" "cdc15_250" "cdc15_270" "cdc15_290"

Muestra el código

rownames(dat) <- d$X
dat <- na.omit(dat)
colnames(dat)

 [1] "cdc15_10"  "cdc15_30"  "cdc15_50"  "cdc15_70"  "cdc15_80"  "cdc15_90" 
 [7] "cdc15_100" "cdc15_110" "cdc15_120" "cdc15_130" "cdc15_140" "cdc15_150"
[13] "cdc15_160" "cdc15_170" "cdc15_180" "cdc15_190" "cdc15_200" "cdc15_210"
[19] "cdc15_220" "cdc15_230" "cdc15_240" "cdc15_250" "cdc15_270" "cdc15_290"

Los datos descargados de la web habían sido previamente normalizados. Para comprobarlo podemos realizar un boxplot y comprobar que los datos presentan una distribución similar y centrada en el cero como sería de esperar con valores de expresión relativa normalizados.

Muestra el código

boxplot(dat, las=2, cex.axis=0.7)

A continuación seleccionaremos sólo aquellos genes que se encuentran entre los que poseen el 1% de las desviaciones estándar más altas. Esto es una práctica habitual en estudios de clasificación puesto que se considera que solo aquellos individuos con algo variabilidad resultan relevantes para la construcción o la distinción de grupos.

Muestra el código

percentage <- c(0.975)
sds <- apply(dat, MARGIN=1, FUN="sd")
sel <- (sds>quantile(sds,percentage))
dat.sel <- dat[sel, ]
dim(dat.sel)

[1] 110  24

Agrupación jerárquica de las muestras

La primera aproximación que suele hacerse para establecer la relación entre entre muestras es la aplicación de algun método de agrupamiento jerárquico. En este caso se realizará un cluster jerárquico basado en distancias euclídeas y enlaces promedio (“average linkage’’).

Muestra el código

distmeth <- c("euclidian")
Distan <- dist(t(dat.sel), method=distmeth)
treemeth <- c("average")
hc <- hclust(Distan, method=treemeth)

Para representar la estructura jerárquica de un agrupamiento se utiliza el dendrograma, un gráfico que tiene forma de árbol invertido, donde los nombres de los genes equivaldrían a las hojas.

Esta no es la única forma de realizar este tipo de agrupación y , en general se recomienda probar con distintos métodos "razonables" y ver hasta que punto cambian los resultados que se obtiene.

Ejercicio

Dibujar un dendrograma utilizando el método de distancias mínimas ("single linkage") y el de distancias máximas ("complete linkage").

A pesar que en los ejemplos anteriores se han utilizado distancias euclídeas, los perfiles de expresión génica acostumbran a agruparse en función de los coeficientes de correlación.

<<<<<<< HEAD:index.Rmd Para ello es necesario calcular la correlación entre los genes y cambiar el formato de la matriz de correlación para pueda contener las distancias.

Muestra el código

cor.pe <- cor(as.matrix(dat.sel), method=c("pearson"))
cor.pe.rows <- cor(t(as.matrix(dat.sel)), method=c("pearson"))
cor.sp <- cor(as.matrix(dat.sel), method=c("spearman"))
dist.pe <- as.dist(1-cor.pe)
dist.pe.rows <- as.dist(1-cor.pe.rows)
dist.sp <- as.dist(1-cor.sp)
hc.cor <- hclust(dist.pe, method=treemeth)
hc.cor.rows <- hclust(dist.pe.rows, method=treemeth)

Para ello es necesario calcular la correlación entre los genes y cambiar el formato de la matriz de correlación para pueda contener las distancias. A continuación, el árbol puede dibujarse normalmente.

Repetir el ejercicio anterior cambiando el coeficiente de Pearson por el de Spearman. ?Los resultados son diferentes?

Agrupación de genes por el método de las `k-means`

Si en lugar de muestras deseamos agrupar por genes es mejor usar métodos no jerárquicos partitivos como k-means.

Un inconveniente de l’agrupación por k-means es que hace falta escoger un número el número de clusters ($k$) antes de realizar la agrupación.

Por ejemplo para producir un agrupamiento k-means con 5 grupos haremos:

Muestra el código

k <- c(4)
km <- kmeans(dat.sel, k, iter.max=1000)

El método de las k-means permite estudiar la consistencia de la agrupación a partir del cálculo de las sumas de cuadrados dentro de cada grupo (Within SS). El promedio de estas sumas de cuadrados para todos los grupos creados (variable "withinss") es una medida de la variabilidad dentro de los grupos, es decir qué tan parecidos son los genes dentro de los grupos.

Muestra el código

mean(km$withinss)

Una manera de evaluar las agrupaciones generadas mediante "k-means" es ejecutar el mismo análisis varias veces – guardando cada vez el resultado como un nuevo objeto – y seleccionar aquella agrupación que ofrece un menor valor de "withinss".

Dibujo de una agrupación `k-means`

Los resultados de los métodos divisivos como k-meansk-means no pueden visualizarse mediante un dendrograma por lo que es preciso recurrir a otras representaciones gráficas.

Podemos, por ejemplo, pintar los miembros de cada grupo con un color y un símbolo distintos.

Una alternativa es utilizar de nuevo el paquete factoextra

Muestra el código

fviz_cluster(km, data = dat.sel,
palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
ellipse.type = "euclid", # Concentration ellipse
star.plot = TRUE, # Add segments from centroids to items
repel = TRUE, # Avoid label overplotting (slow)
ggtheme = theme_minimal()
)

Para visualizar un sólo grupo, y sus valores de expresión, seleccionaremos los genes que lo constituyen y a continuación, dibujaremos los perfiles de expresión de estos genes utilizando diferentes colores.

Si hay más de un grupo podemos visualizarlos todos en el mismo gráfico utilizando un peque~no bucle for

Muestra el código

km <- kmeans(dat.sel, 4, iter.max=1000)

Anotación de resultados

Una vez que se han identificado algunos grupos de interés, se puede comprobar que genes constituyen de cada uno de los "clusters". Esto puede hacerse con el paquete de anotaciones para levadura de donde podremos extraer información sobre nombres, descripciones y otras anotaciones en múltiples bases de datos:

Muestra el código

if(!require(org.Sc.sgd.db)){
 BiocManager::install("org.Sc.sgd.db")
}
require("org.Sc.sgd.db")
anotData <- capture.output(org.Sc.sgd())
print(anotData[1:15])

 [1] "Quality control information for org.Sc.sgd:"               
 [2] ""                                                          
 [3] ""                                                          
 [4] "This package has the following mappings:"                  
 [5] ""                                                          
 [6] "org.Sc.sgdALIAS has 4468 mapped keys (of 8061 keys)"       
 [7] "org.Sc.sgdCHR has 8061 mapped keys (of 8061 keys)"         
 [8] "org.Sc.sgdCHRLENGTHS has 17 mapped keys (of 17 keys)"      
 [9] "org.Sc.sgdCHRLOC has 7993 mapped keys (of 8061 keys)"      
[10] "org.Sc.sgdCHRLOCEND has 7993 mapped keys (of 8061 keys)"   
[11] "org.Sc.sgdCOMMON2ORF has 12361 mapped keys (of 12361 keys)"
[12] "org.Sc.sgdDESCRIPTION has 8052 mapped keys (of 8061 keys)" 
[13] "org.Sc.sgdENSEMBL has 5796 mapped keys (of 8061 keys)"     
[14] "org.Sc.sgdENSEMBL2ORF has 5802 mapped keys (of 5802 keys)" 
[15] "org.Sc.sgdENSEMBLPROT has 5796 mapped keys (of 8061 keys)"

Muestra el código

cat ("... output continues until ", length(anotData), " lines.\n")

... output continues until  46  lines.

Podemos utilizar las funciones de anotacion para convertir los identificadores propios a identificadores más generals (“Entrezs”) enriqueciéndolos con información sobre dichos genes.

Muestra el código

genes<-rownames(dat.sel) 
columns(org.Sc.sgd.db)

 [1] "ALIAS"        "COMMON"       "DESCRIPTION"  "ENSEMBL"      "ENSEMBLPROT" 
 [6] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL" 
[11] "GENENAME"     "GO"           "GOALL"        "INTERPRO"     "ONTOLOGY"    
[16] "ONTOLOGYALL"  "ORF"          "PATH"         "PFAM"         "PMID"        
[21] "REFSEQ"       "SGD"          "SMART"        "UNIPROT"

Muestra el código

geneAnots <-select(org.Sc.sgd.db, genes,c("SGD", "ENTREZID", "GENENAME"))
head(geneAnots)

      ORF        SGD ENTREZID GENENAME
1 YAR007C S000000065   851266     RFA1
2 YBL003C S000000099   852283     HTA2
3 YBL051C S000000147   852229     PIN4
4 YBR038W S000000242   852326     CHS2
5 YBR054W S000000258   852343     YRO2
6 YBR092C S000000296   852389     PHO3

Muestra el código

###################################################
### code chunk number 15: genesInClusters
###################################################
genesInClusters <- cbind(geneAnots, cl)
kableExtra::kable(genesInClusters[order(genesInClusters$cl),])

	ORF	SGD	ENTREZID	GENENAME	cl
YCL027W	YCL027W	S000000532	850330	FUS1	1
YCL040W	YCL040W	S000000545	850317	GLK1	1
YCL055W	YCL055W	S000000560	850303	KAR4	1
YCR098C	YCR098C	S000000695	850462	GIT1	1
YER103W	YER103W	S000000905	856840	SSA4	1
YGL089C	YGL089C	S000003057	852791	MF(ALPHA)2	1
YJR004C	YJR004C	S000003764	853460	SAG1	1
YKL177W	YKL177W	S000001660	NA	NA	1
YKL178C	YKL178C	S000001661	853677	STE3	1
YLL024C	YLL024C	S000003947	850636	SSA2	1
YLL026W	YLL026W	S000003949	850633	HSP104	1
YLR247C	YLR247C	S000004237	850949	IRC20	1
YMR072W	YMR072W	S000004676	855094	ABF2	1
YMR186W	YMR186W	S000004798	855224	HSC82	1
YMR189W	YMR189W	S000004801	855227	GCV2	1
YNL007C	YNL007C	S000004952	855725	SIS1	1
YNL057W	YNL057W	S000005002	NA	NA	1
YNR044W	YNR044W	S000005327	855780	AGA1	1
YOL058W	YOL058W	S000005419	854096	ARG1	1
YPR158W	YPR158W	S000006362	856281	CUR1	1
YBL051C	YBL051C	S000000147	852229	PIN4	2
YBR110W	YBR110W	S000000314	852407	ALG1	2
YBR158W	YBR158W	S000000362	852455	AMN1	2
YCL024W	YCL024W	S000000529	850334	KCC4	2
YDR097C	YDR097C	S000002504	851671	MSH6	2
YGL055W	YGL055W	S000003023	852825	OLE1	2
YGR044C	YGR044C	S000003276	852935	RME1	2
YGR109C	YGR109C	S000003341	853003	CLB6	2
YJL159W	YJL159W	S000003695	853281	HSP150	2
YKL163W	YKL163W	S000001646	853693	PIR3	2
YKL164C	YKL164C	S000001647	853692	PIR1	2
YKL172W	YKL172W	S000001655	853682	EBP2	2
YKL185W	YKL185W	S000001668	853650	ASH1	2
YLR049C	YLR049C	S000004039	850738	NA	2
YLR079W	YLR079W	S000004069	850768	SIC1	2
YLR194C	YLR194C	S000004184	850891	NCW2	2
YMR193W	YMR193W	S000004806	855231	MRPL24	2
YNL082W	YNL082W	S000005026	855642	PMS1	2
YNR067C	YNR067C	S000005350	855804	DSE4	2
YOL101C	YOL101C	S000005461	854052	IZH4	2
YOR308C	YOR308C	S000005835	854485	SNU66	2
YAR007C	YAR007C	S000000065	851266	RFA1	3
YBL003C	YBL003C	S000000099	852283	HTA2	3
YDL003W	YDL003W	S000002161	851561	MCD1	3
YDL055C	YDL055C	S000002213	851504	PSA1	3
YDR224C	YDR224C	S000002632	851810	HTB1	3
YDR225W	YDR225W	S000002633	851811	HTA1	3
YER124C	YER124C	S000000926	856861	DSE1	3
YGR189C	YGR189C	S000003421	853102	CRH1	3
YHR143W	YHR143W	S000001186	856546	DSE2	3
YIL066C	YIL066C	S000001328	854744	RNR3	3
YIL140W	YIL140W	S000001402	854666	AXL2	3
YJL078C	YJL078C	S000003614	853367	PRY3	3
YKR013W	YKR013W	S000001721	853882	PRY2	3
YLR121C	YLR121C	S000004111	850812	YPS3	3
YLR183C	YLR183C	S000004173	850880	TOS4	3
YLR286C	YLR286C	S000004276	850992	CTS1	3
YML027W	YML027W	S000004489	854981	YOX1	3
YMR199W	YMR199W	S000004812	855239	CLN1	3
YNL030W	YNL030W	S000004975	855701	HHF2	3
YNL031C	YNL031C	S000004976	855700	HHT2	3
YOL007C	YOL007C	S000005367	854155	CSI2	3
YOL090W	YOL090W	S000005450	854063	MSH2	3
YPL127C	YPL127C	S000006048	855976	HHO1	3
YPL256C	YPL256C	S000006177	855819	CLN2	3
YBR038W	YBR038W	S000000242	852326	CHS2	4
YBR054W	YBR054W	S000000258	852343	YRO2	4
YBR092C	YBR092C	S000000296	852389	PHO3	4
YBR138C	YBR138C	S000000342	852435	NA	4
YBR202W	YBR202W	S000000406	852501	MCM7	4
YCL064C	YCL064C	S000000569	850295	CHA1	4
YCR005C	YCR005C	S000000598	850361	CIT2	4
YCR024C-A	YCR024C-A	S000000619	850389	PMP1	4
YDL037C	YDL037C	S000002195	851524	BSC1	4
YDL039C	YDL039C	S000002197	851522	PRM7	4
YDL077C	YDL077C	S000002235	851482	VAM6	4
YDL234C	YDL234C	S000002393	851364	GYP7	4
YDR033W	YDR033W	S000002440	851597	MRH1	4
YDR039C	YDR039C	S000002446	851609	ENA2	4
YDR273W	YDR273W	S000002681	851866	DON1	4
YFL011W	YFL011W	S000001883	850536	HXT10	4
YGL008C	YGL008C	S000002976	852876	PMA1	4
YGR108W	YGR108W	S000003340	853002	CLB1	4
YGR279C	YGR279C	S000003511	853196	SCW4	4
YHL028W	YHL028W	S000001020	856357	WSC4	4
YHR137W	YHR137W	S000001179	856539	ARO9	4
YJL079C	YJL079C	S000003615	853366	PRY1	4
YKL096W	YKL096W	S000001579	853766	CWP1	4
YKR039W	YKR039W	S000001747	853912	GAP1	4
YLR126C	YLR126C	S000004116	850817	NA	4
YLR190W	YLR190W	S000004180	850887	MMR1	4
YML052W	YML052W	S000004516	854953	SUR7	4
YML058W	YML058W	S000004523	854945	SML1	4
YMR001C	YMR001C	S000004603	855013	CDC5	4
YMR031C	YMR031C	S000004633	855047	EIS1	4
YMR032W	YMR032W	S000004635	855048	HOF1	4
YMR058W	YMR058W	S000004662	855080	FET3	4
YMR203W	YMR203W	S000004816	855243	TOM40	4
YMR215W	YMR215W	S000004828	855255	GAS3	4
YNL160W	YNL160W	S000005104	855562	YGP1	4
YNL270C	YNL270C	S000005214	855451	ALP1	4
YNR009W	YNR009W	S000005292	855743	NRM1	4
YOR256C	YOR256C	S000005782	854430	TRE2	4
YOR298W	YOR298W	S000005824	854473	MUM3	4
YPL021W	YPL021W	S000005942	856086	ECM23	4
YPL061W	YPL061W	S000005982	856044	ALD6	4
YPL132W	YPL132W	S000006053	855971	COX11	4
YPL265W	YPL265W	S000006186	855863	DIP5	4
YPR119W	YPR119W	S000006323	856236	CLB2	4
YPR149W	YPR149W	S000006353	856272	NCE102	4

Algunas fuentes con ejemplos adicionales sobre análisis de conglomerados:

Ejemplos de cluster analysis con R del profesor Gabriel Martos El paquete cluster.datasets, con los datasets de un clásico sobre análisisd e conglomerados (Hartigan) Materiales de Thomas Girke, bioinformatico de la UCR, sobre Clustering

	\(x_1\)	\(x_2\)	\(x_3\)	\(\cdots\)	\(x_n\)
\(i\)	\(+\)	\(-\)	\(+\)	\(\cdots\)	\(+\)

			\(i\)
		\(+\)		\(-\)
	\(+\)	\(a\)		\(b\)
\(j\)
	\(-\)	\(c\)		\(d\)

Presentación

Instalación de paquetes utilizados

Análisis de componentes principales

Ejemplo PCA-1: Búsqueda de factores latentes en datos ecológicos1

Mediante diagonalización de la matriz de covarianzas

Calculo de las componentes principales mediante las funcionesn princomp y prcomp

Interpretación de las componentes

Ejemplo PCA-2: Detección del efecto batch en datos de microarrays.

Los datos para el análisis

Análisis del efecto batch

Detectando el efecto batch

Análisis basado en distancias

Análisis de proximidades

Breve repaso sobre medidas de similaridad

Ejemplo 1: Distancias geográficas2

Ejemplo 2: Similaridades entre partidos políticos3

Análisis de correspondencias

Ejemplo 1: Representación de mutaciones cromosómicos entre poblaciones

Realización del análisis paso por paso

Como un escalado multidimensional

AC como un ACP sobre la matriz estandarizada

Ejemplo 2: Análisis de correspondencias de datos de microarrays

Análisis de conglomerados (“Cluster Analysis”)

Agrupación de genes y de muestras en el análisis del ciclo celular

Clustering y visualización

Pre-procesamiento de los datos

Agrupación jerárquica de las muestras

Ejercicio

Agrupación de genes por el método de las k-means

Dibujo de una agrupación k-means

Anotación de resultados

Algunas fuentes con ejemplos adicionales sobre análisis de conglomerados:

Referencias

Ejemplo PCA-1: Búsqueda de factores latentes en datos ecológicos¹

Calculo de las componentes principales mediante las funcionesn `princomp` y `prcomp`

Ejemplo 1: Distancias geográficas²

Ejemplo 2: Similaridades entre partidos políticos³

Agrupación de genes por el método de las `k-means`

Dibujo de una agrupación `k-means`