1 Objetivo

En este documento se presenta la estrategia y el codigo necesario para identificar subpoblaciones espermáticas cinemáticas, partiendo de datos de espermatozoides individuales, evaluados con un sistema CASA (Computer Assisted Sperm Analysis). Se presenta la estrategia y el código utilizado en el artículo de Rivas et al (2022).

2 Introducción

Para identificación de subpoblaciones espermáticas cinemáticas, generalmente se utilizan los valores de los parámetros de movilidad obtenidos con un sistema CASA. Los valores se calculan para cada espermatozoide. El número de parámetros de movilidad depende del fabricante del sistema CASA; aunque generalemnte se evalúan VCL, VAP, VSL, LIN, STR, BCF, WOB y ALH.

Existen diferentes estrategías estadísticas que utilizan como entrada los valores de los parámetros de movilidad, evaluados en un sistema CASA (Martínez-Pastor et al., 2011; Ramón and Martínez-Pastor, 2018), para identificar subpoblaciones espermáticas cinemáticas. Es común hacer la identificación en dos etapas. En la primera etapa se trata de reducir la dimensionalidad del dataset; esto se realiza mediante un análisis de componentes principales. En el análisis de componentes principales se retienen aquellos componentes que tienen un eigenvalue mayor a 1. En la segunda etapa, los componentes principales sirven como entrada para algoritmos de agrupamiento, como Kmeans o aglomerativo jerárquico.

En el caso del agrupamiento aglomerativo jerárquico, se trata de agrupar los datos mas cercanos entre sí, y mas alejados de otros grupos utilizando diferentes tipos de distancia, aunque la más usada es la distancia euclidiana. El número de grupos es un aspecto importante a considerar, por ello se busca utilizar un criterio objetivo de elección, como puede serlo la ganancia de inercia. En este sentido, el análisis gráfico de la ganancia de inercia es de gran ayuda.

En este documento presentaremos la estrategia y el procedimiento para identificar subpoblaciones espermáticas cinemáticas, a partir de datos de espermatozoides de cerdo expuestos a difefrentes valores de pH (Rivas et al., 2022).

3 Requerimientos

Librerías: Para replicar el procedimiento presentado en este documento se requiere tener instaladas las librerías FactoMineR, y dplyr. Datos: Es necesario descargar el dataset “230819_pluginnuevo3.tab” en su directorio de trabajo. El dataset se localiza en https://doi.org/10.7910/DVN/CRVVSI (Rivas et al., 2022).

Enseguida descargue el archivo:

download.file(url = "https://dataverse.harvard.edu/api/access/datafile/4549661", destfile="motility_results_ph.csv")

Crearemos un objeto llamado mr y meteremos en él el contenido del dataset descargado:

mr<-read.csv("motility_results_ph.csv", header=TRUE, sep=",", stringsAsFactors=TRUE)

El archivo motility_results_ph.csv tiene 11620 lineas de datos, pero para el propósito de este documento, se utilizará una muestra con 500 líneas de datos.

mr_sample<-sample(mr, 500, replace=FALSE, prob=NULL)

4 Procedimiento

4.1 Etapa 1, Adquisición y creación de objetos

El dataset utilizado en este flujo de trabajo puede descargarse desde nuestra cuenta en el sitio de Harvard Dataverse https://dataverse.harvard.edu, como se indicó arriba, o directamente con el siguiente código:

setwd("/Users/andresammx/Documents/RStudio/Markdown/Identificación_de_subpoblaciones")

Descargar el dataset.

download.file(url = "https://dataverse.harvard.edu/api/access/datafile/4549661", destfile="230819_pluginnuevo.tab")

Leemos el archivo:

mr<-read.csv("230819_pluginnuevo3.tab", header=TRUE, sep="\t", stringsAsFactors=TRUE)

Tomamos una muestra de 500 filas del dataset y las guardamos en un nuevo objeto llamado mr_sample:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
mr_sample<- mr %>% group_by(ID2) %>% slice_sample(n=100, replace=FALSE)

Convertimos el objeto mr_sample a dataframe y verificamos su estructura:

mr_sample<-as.data.frame(mr_sample)
str(mr_sample)
## 'data.frame':    500 obs. of  11 variables:
##  $ ID1      : Factor w/ 75 levels "220219-1","220219-10",..: 58 68 58 68 46 69 35 70 7 1 ...
##  $ ID2      : Factor w/ 5 levels "pH.70","pH.74",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ID3      : Factor w/ 6 levels "","balder","berthrand",..: 2 2 2 2 4 2 5 2 6 5 ...
##  $ VCL      : num  213.2 84.1 163.6 293.8 533.7 ...
##  $ VAP      : num  117.7 37.1 68 121 206.5 ...
##  $ VSL      : num  76.1 14.9 23.1 34.2 71.8 ...
##  $ LIN      : num  0.357 0.177 0.141 0.116 0.134 ...
##  $ STR      : num  0.646 0.402 0.34 0.283 0.348 ...
##  $ WOB      : num  0.552 0.442 0.416 0.412 0.387 ...
##  $ BeatCross: num  27.1 38.3 29 18.6 29.1 ...
##  $ ALH      : num  8.38 3.25 4.74 10.21 17.55 ...

El objeto mr_sample tiene 500 filas de datos y 11 columnas. La columna (ID1) corresponde al identificador del video de la rutina de captura; la columna 2 (ID2) corresponde al tratamiento, los diferetes valores de pH; la columna 3 (ID3) corresponde al identificador del semental donante.

4.2 Etapa 2, Análisis de componentes principales (PCA)

Se cargan las librerías necesarias:

library(FactoMineR)

Ahora debemos crear un nuevo objeto mr_sample_2 que contenga únicamente las variables (columnas) que vamos a utilizar:

mr_sample_2<-mr_sample[,c(2,4:11)]

Se crea un objeto res.pca que contenga los resultados de la función de componentes principales. Los argumentos para la función PCA indican que los datos deben escalarse, que se esperan tres componentes principales, que se está ingresando una variable cualitativa que no tiene peso en el resultado y que pedimos que se grafiquen los resultados.

res.pca = PCA(mr_sample_2, scale.unit=TRUE, ncp=3, quali.sup=1, graph=TRUE)

La funsión PCA roduce dos gráficos, en el PCA graph of individuals se encuentran graficados los valores de los componentes principales 1 y 2; al alado de cada punto se observa el numero de línea correspondiente. La gráfica PCA graph of variables muestra un circulo de correlaciones, en donde cada flecha indica la magnitud y el signo de la correlación de una variable con respecto a los componentes principales 1 y 2. Así, por ejemplo vemos que VCL tiene una alta correlación positiva con el componente principal 1 y una baja correlación con el componente principal 2.

res.pca es un objeto de tipo lista, y se puede acceder a su contenido mediante el operador $. Escribir res.pca en la consola nos mostrará primero a los objetos que componen la lista y enseguida mostrará el significado del contenido de cada objeto:

res.pca
## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 500 individuals, described by 9 variables
## *The results are available in the following objects:
## 
##    name                description                                          
## 1  "$eig"              "eigenvalues"                                        
## 2  "$var"              "results for the variables"                          
## 3  "$var$coord"        "coord. for the variables"                           
## 4  "$var$cor"          "correlations variables - dimensions"                
## 5  "$var$cos2"         "cos2 for the variables"                             
## 6  "$var$contrib"      "contributions of the variables"                     
## 7  "$ind"              "results for the individuals"                        
## 8  "$ind$coord"        "coord. for the individuals"                         
## 9  "$ind$cos2"         "cos2 for the individuals"                           
## 10 "$ind$contrib"      "contributions of the individuals"                   
## 11 "$quali.sup"        "results for the supplementary categorical variables"
## 12 "$quali.sup$coord"  "coord. for the supplementary categories"            
## 13 "$quali.sup$v.test" "v-test of the supplementary categories"             
## 14 "$call"             "summary statistics"                                 
## 15 "$call$centre"      "mean of the variables"                              
## 16 "$call$ecart.type"  "standard error of the variables"                    
## 17 "$call$row.w"       "weights for the individuals"                        
## 18 "$call$col.w"       "weights for the variables"

Los valores eigen y la varianza acumulada presentados en Rivas et al., 2022) se encuentran en el objeto res.pca$eig:

res.pca$eig
##         eigenvalue percentage of variance cumulative percentage of variance
## comp 1 4.076974926             50.9621866                          50.96219
## comp 2 2.355910038             29.4488755                          80.41106
## comp 3 0.926365554             11.5795694                          91.99063
## comp 4 0.455068897              5.6883612                          97.67899
## comp 5 0.107369499              1.3421187                          99.02111
## comp 6 0.048517430              0.6064679                          99.62758
## comp 7 0.021014711              0.2626839                          99.89026
## comp 8 0.008778945              0.1097368                         100.00000

De acuerdo al criterio de Kaiser, se eligen los componentes principales con valores eigen mayores a 1.

Se pueden guardar los valores de las coordenadas de cada componente principal, con el propósito de generar graficos limpios y personalizados:

dind<-as.data.frame(res.pca$ind$coord)
names(dind)
## [1] "Dim.1" "Dim.2" "Dim.3"

Y se pueden pegar esas columnas al objeto mr_sample_2, por ejemplo:

tab2<-cbind(mr_sample_2$ID2,dind)
head(tab2)
##   mr_sample_2$ID2       Dim.1     Dim.2       Dim.3
## 1           pH.70  0.99670454  1.160061 -0.04373894
## 2           pH.70 -1.68269163 -2.027713 -0.36329757
## 3           pH.70 -0.06122691 -2.176256 -0.78562889
## 4           pH.70  2.49214245 -2.014954 -1.13065128
## 5           pH.70  4.79649558 -1.056998  1.05761245
## 6           pH.70  0.65297406 -2.020568 -1.49055675

Preparamos tab2 para producir la gráfica personalizada:

tab2<-as.data.frame(tab2)
names(tab2)
## [1] "mr_sample_2$ID2" "Dim.1"           "Dim.2"           "Dim.3"
head(tab2)
##   mr_sample_2$ID2       Dim.1     Dim.2       Dim.3
## 1           pH.70  0.99670454  1.160061 -0.04373894
## 2           pH.70 -1.68269163 -2.027713 -0.36329757
## 3           pH.70 -0.06122691 -2.176256 -0.78562889
## 4           pH.70  2.49214245 -2.014954 -1.13065128
## 5           pH.70  4.79649558 -1.056998  1.05761245
## 6           pH.70  0.65297406 -2.020568 -1.49055675
colnames(tab2)[1]<-"pH"
names(tab2)
## [1] "pH"    "Dim.1" "Dim.2" "Dim.3"

Producimos la gráfica:

plot(tab2$Dim.1, tab2$Dim.2, col=tab2$pH, xlab = "PC1", ylab = "PC2", xlim = c(-4,8), ylim = c(-4,8))
legend("topright", legend = levels(tab2$pH), col = 1:length(levels(tab2$pH)), pch = 1, cex = 1.2)

4.3 Etapa 3, Agrupamiento jerárquico sobre los componentes principales

Una vez que se ha realizado el análisis de componentes principales, el objeto res.pca servirá como entrada para la función HCPC; se crea un nuevo objeto llamado res.hcpc en donde se guardarán los resultados de HCPC.

res.hcpc <- HCPC(res.pca, nb.clust=-1, conso=0, min=3, max=10, metric="euclidean", method="ward")

La función HCPCgenera tres gráficos, el primero es “Hierarchical clustering” que tiene una línea horizontal sugiriendo el punto de corte para la formación de grupos en el dendograma. Esta gráfico tiene un inserto con una gráfica de barras que indica la ganancia de inercia. Una vez que se hace clic sobre la barra o en algún otro nivel de la grafica, se generan otras dos gráficas (Para los propoósitos de este documento el argumento nb.clust se ajustó a -1 para hacer el corte automático, pero el ajuste igual a 0 permite al usuario elegir la altura de corte). La gráfica “Hierarchical clustering on the factor map” es una representación en tres dimensiones de los componentes principales 1 y 2 en los ejex x y y, respectivamente; mientras que z se dibuja el dendograma presentado bidimensionalemnte en “Hierarchical clustering”. El grafico “Factor map” presenta los valores de los componentes principales 1 y 2 ahora coloreados de acuerdo a los grupos formados o clusters. Estos grupos son las subpoblaciones cinemáticas.

El objeto res.hcpc es un objeto de tipo lista, y se puede acceder a su contenido mediante el operador $. Escribir res.hcpc en la consola nos mostrará primero a los objetos que componen la lista y enseguida mostrará el significado del contenido de cada objeto:

res.hcpc
## **Results for the Hierarchical Clustering on Principal Components**
##    name                   
## 1  "$data.clust"          
## 2  "$desc.var"            
## 3  "$desc.var$quanti.var" 
## 4  "$desc.var$quanti"     
## 5  "$desc.axes"           
## 6  "$desc.axes$quanti.var"
## 7  "$desc.axes$quanti"    
## 8  "$desc.ind"            
## 9  "$desc.ind$para"       
## 10 "$desc.ind$dist"       
## 11 "$call"                
## 12 "$call$t"              
##    description                                             
## 1  "dataset with the cluster of the individuals"           
## 2  "description of the clusters by the variables"          
## 3  "description of the cluster var. by the continuous var."
## 4  "description of the clusters by the continuous var."    
## 5  "description of the clusters by the dimensions"         
## 6  "description of the cluster var. by the axes"           
## 7  "description of the clusters by the axes"               
## 8  "description of the clusters by the individuals"        
## 9  "parangons of each clusters"                            
## 10 "specific individuals"                                  
## 11 "summary statistics"                                    
## 12 "description of the tree"

El objeto res.hcpc$data.clust contiene los datos de rm_sample_2 y al final tiene una columna llamada clust, que indica el cluster al cual pertenece cada línea de datos:

head(res.hcpc$data.clust)
##     ID2       VCL       VAP      VSL      LIN      STR      WOB BeatCross
## 1 pH.70 213.24402 117.73412 76.07623 0.356757 0.646170 0.552110  27.07317
## 2 pH.70  84.07959  37.14471 14.92351 0.177493 0.401767 0.441780  38.27586
## 3 pH.70 163.56094  67.98609 23.09495 0.141201 0.339701 0.415662  28.96552
## 4 pH.70 293.77853 121.00965 34.22053 0.116484 0.282792 0.411908  18.62069
## 5 pH.70 533.67860 206.52632 71.76946 0.134481 0.347508 0.386986  29.09091
## 6 pH.70 176.40210  82.96847 23.80340 0.134938 0.286897 0.470337  24.31035
##         ALH clust
## 1  8.377957     2
## 2  3.246029     1
## 3  4.740679     2
## 4 10.212711     2
## 5 17.550550     3
## 6  6.082295     2

Este objeto es muy importante ya que nos permitirá analizar el dataset en función de los tratamientos y de los clusters. El contenido de res.hcpc$data.clust se puede guardar en un nuevo objeto para analizarlo posteriormente:

mr_sub<-res.hcpc$data.clust

El objeto res.hcpc$desc.axes contiene los valores de cada cluster en función de las variables cuantitativas (los parámetros de movilidad).

4.4 Etapa 4, Subpoblaciones en el dataset

Uno de los objetivos de la identificación de subpobaciones espermáticas cinemáticas, es describir el efecto de los tratamientos experimentales en la estrucura de las subpoblaciones, iniciando por conocer la proporción de espermatozoides de cada tratamiento en cada subpoblación. Así mismo es relevante describir el efecto de cada tratamiento en cada cluster. Este trabajo será abordado en otro documento de Rpubs.

5 Conclusiones

Con el flujo de trabajo propuesto en este documento se puede manejar un archivo de parámetros de movilidad, obtenidos con un sistema CASA, para identificar subpoblaciones cinemáticas en dos pasos. La resultante es un objeto que contiene los mismos datos de entrada mas una columna que indica a que clusyter pertenece cada línea de datos. El manejo del objeto resultante permitirá identificar la proporción de espermatozoides en cada subpoblación y el efecto de cada tratamiento en cada subpoblación.

Bibliografía

Martínez-Pastor F, Tizado EJ, Garde JJ, Anel L and Paz P de (2011) Statistical Series: Opportunities and challenges of sperm motility subpopulation analysis. Theriogenology 75 783–795.
Ramón M and Martínez-Pastor F (2018) Implementation of novel statistical procedures and other advanced approaches to improve analysis of CASA data. Reproduction, Fertility, and Development 30 860–866.
Rivas AC, Ayala EME and Aragon MA (2022) Effect of various pH levels on the sperm kinematic parameters of boars. South African Journal of Animal Science 52 693–704.