EJERCICIO DE NANCY CATS

A. Librería e importación de paquetería

Primero cargamos la libreria adegenet y la base de datos nancycats que esta incrustada en el paquete adegenet y que podemos llamarlo a nuestro entrono de trabajo con la función data(). Note, que no es necesario asociar a un nombre el resultado de la función data(). El objeto que se llama igual que la base de datos original nancycats

library(adegenet)
## Loading required package: ade4
## 
##    /// adegenet 2.1.10 is loaded ////////////
## 
##    > overview: '?adegenet'
##    > tutorials/doc/questions: 'adegenetWeb()' 
##    > bug reports/feature requests: adegenetIssues()
data("nancycats")

B. Manejo de objetos genInd y genPop

Luego podemos prenguntarnos si es realmente un objeto genInd usando la función is.genind() del paquete adegenet

adegenet::is.genind(nancycats)
## [1] TRUE

Luego que corroboramos que es un objeto genInd, debemos convertirlo a un obeto de poblaciones genPop usando la función adegenet::genind2genpop() definida en el mismo paquete.

nancypop <- adegenet::genind2genpop(nancycats)
## 
##  Converting data from a genind to a genpop object... 
## 
## ...done.
adegenet::is.genpop(nancypop)
## [1] TRUE

C. Analisis de PCA con objeto genPop nancypop

C.1. Extracción de la base de datos

Vamos a extraer la base de datos (tabla) del objeto genPop usando la función tab(). Como tenemos algunos datos que son NA, y para posteriores calculos necesitamos una matriz de datos sin NA, vamos agregar el parámetro NA.method="mean" a la fucnción tab()

x <- tab(nancypop, NA.method = "mean")   # extraemos la tabla con los datos por alelo

C.2. Generación de la matriz de correlación

Ahora vamos a hacer un reduccion de variables usando el estadistico de PCA. Hay varias tecnicas para realizar este analisis en R, tanto con las funciones del paquete adegenet, como funciones basicas de R u otras paqueterías. Aqui vamos a realizar una forma general para el PCA en R. Lo primero que debemos hacer es generar una matriz de correlacion con la función básica de R cor()

mat.cor <- cor(t(x))              # generación de la matriz de correlación

C.3. Gráfica de la matriz de correlación

Para graficar esta matriz, vamos a instalar el paquete ggcorrplot que esta en el CRAN de R, por lo que podemos usar el instalador clásico de R install.packages. (Nota: no olvide definier el parámetro dependencies=T, para descargar todos los paquetes que se necesite). Sin embargo, la versión mas actualizada se encuentra en el repositorio de github. Para ello tenemos que instalar el paquete devtools y luego instalar la libreria ggcorrplot desde github con la función devtools::install_github(). No se olvide cargar la librería a su entorno de trabajo con la función library()

# instalación clasica
# install.packages("ggcorrplot", dependencies = T)

# instalición alternativa desde el repositorio de Github
if(!require(devtools)) install.packages("devtools")
## Loading required package: devtools
## Loading required package: usethis
devtools::install_github("kassambara/ggcorrplot")
## Skipping install of 'ggcorrplot' from a github remote, the SHA1 (0a85456d) has not changed since last install.
##   Use `force = TRUE` to force installation
# Cargamos al libreria al entorno de trabajo
library(ggcorrplot)
## Loading required package: ggplot2

Ahora que ya tenemos esta libraría desplegada en nuestro entorno de trabajo, podemos gráficar el resultado de la matriz de correlación.

ggcorrplot(mat.cor)                # graficando la matriz de correlación

C.4. Reducción de variables, análisis de PCA

En esta matriz de correlacion la relación que hay entre los alelos no queda tan clara y puede ser que el patron que estamos buscando se diluya o se enmascare con el enorme numero de variables que describen una muestra. En este caso, podemos definir al numero de muestras como el numero de poblaciones del objeto gendPop “nancyPop”. Para ello podemos usar la función npop del paquete adegenet, y para definir el numero de variables, que vendría a ser el número de alelos posibles podemos usar la fucnión nAll del mismo paquete

adegenet::nAll(nancypop)            # numero de alelos por locus
##  fca8 fca23 fca43 fca45 fca77 fca78 fca90 fca96 fca37 
##    16    11    10     9    12     8    12    12    18
sum(adegenet::nAll(nancypop))       # numero de alelos por todos los loci
## [1] 108
adegenet::nPop(nancypop)            # numero de poblaciones del objeto genPop
## [1] 17

Con esto vemos que el número total de alelos es 108 y que el número de poblaciones es 17. Entonces, lo que queremos hacer con este analisis es reducir las 108 variables a solo dos (los dos primeros compnentes principales) para describir/caracterizar nuestras 17 poblaciones. Para ello usaremos la función princomp del

pca1 <- princomp(mat.cor)

Para explorar este resultado, podemos usar la función summary()

summary(pca1)
## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     0.2464258 0.1998988 0.1885052 0.15464922 0.14507554
## Proportion of Variance 0.2505844 0.1648929 0.1466319 0.09869093 0.08685006
## Cumulative Proportion  0.2505844 0.4154773 0.5621092 0.66080009 0.74765016
##                            Comp.6     Comp.7     Comp.8     Comp.9    Comp.10
## Standard deviation     0.13150591 0.09907121 0.09636279 0.08934614 0.07183266
## Proportion of Variance 0.07136287 0.04050202 0.03831779 0.03294074 0.02129245
## Cumulative Proportion  0.81901302 0.85951504 0.89783283 0.93077357 0.95206603
##                           Comp.11    Comp.12     Comp.13     Comp.14
## Standard deviation     0.06300533 0.05086877 0.043834983 0.038462077
## Proportion of Variance 0.01638085 0.01067786 0.007929092 0.006104459
## Cumulative Proportion  0.96844687 0.97912473 0.987053824 0.993158284
##                            Comp.15     Comp.16      Comp.17
## Standard deviation     0.033700484 0.022853285 1.477755e-09
## Proportion of Variance 0.004686558 0.002155158 9.011278e-18
## Cumulative Proportion  0.997844842 1.000000000 1.000000e+00

El resultado de summary() sobre el output de princomp es una tabla de doble entrada donde la tercera fila corresponde a la proporción acumulada y las columnas a los componentes. Vemos que con los dos primeros componentes podemos recoger el 46.09% de toda la variación de nuestros datos.

C.5. Graficando PCA

Para tener una mejor apreciación de como los datos se “re-distribuyen” considerando solo los dos primeros componentes, podemos proyectarlos en un gráfica de dos dimensiones. Para ello vamos a instalar un paquete llamado “factoexgtra” que se encuentra en el respositorio de github.

library(devtools)          # cargamos nuevamente la libraria de devtools (just in case)
devtools::install_github("kassambara/factoextra")     # instalamos paqueteria
## Skipping install of 'factoextra' from a github remote, the SHA1 (1689fc74) has not changed since last install.
##   Use `force = TRUE` to force installation
library(factoextra)       # desplegamos libria en entorno de trabajo
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Ahora podemos hacer una grafica de la proporción de la variación explicada por cada componente

fviz_eig(pca1, addlabels = TRUE)

El siguiente paso sería graficar los los resultados basados en los dos primeros componentes

fviz_pca_var(pca1, col.var = "black")

PREGUNTAS

  1. Basado en la última gráfica de PCA, que puede concluir?, que poblaciones son las mas cercanas, es decir, que poblaciones se agrupan?.

  2. Si asume que estas poblaciones tienen una distribución geográfica propia, y que las que se encuentran más alejadas tienen menor oportunidad de intercambiar material genético, Como imaginaría la distribución de estas poblaciones?

  3. Porque debemos convertir el objeto genInd a un objeto genPop?. Haga el mismo analisis de PCA, pero con el objeto genInd. El resultado es diferente?. Que ventajas nos da una analisis sobre el otro.