In currocam/FascinRSCA: R and Python code for SCA with fascin protein

knitr::opts_chunk$set(echo = TRUE,
  dpi=200)
library("FactoMineR")
library("factoextra")

library(tidyverse)
library(magrittr)
library(FascinRSCA)
data <- FascinRSCA::freezing %>% distinct(fasta_seq, .keep_all = TRUE) %>%
  mutate(annotation = fct_collapse(annotation, fascin2 = c("fascin2","fascin2a")))

Estandarización de los datos

Seleccionamos datos númericos de interés sobre las secuencias y los escalamos para que tengan desviación estándar uno y media cero. Realizamos PCA sobre las r dim(data)[[1]] secuencias homólogas a fascina no redundantes (individuals) y 39 variables numéricas de interés. Guardamos en resultados las primeras 15 dimensiones (explicación después).

#summary(data)
x<-data %>% select(where(is.numeric))%>% select(-c (ndom, Length, nregions,
  niseqs, pvalue, nincluded, flags, nreported))

#x<- cbind(x, acc =data$acc, ann = data$annotation_rev)
#res.pca<-PCA(x, scale.unit = TRUE, ncp = 13, graph = FALSE, quali.sup = c(40, 41))
res.pca<-PCA(x, scale.unit = TRUE, ncp = 15, graph = FALSE)

Análisis objeto PCA

El objeto generado por PCA contiene toda la información necesaria para nuestro analisis.

Varianza de las dimensiones

Extraemos y visualizamos los eigenvalues correspondientes a cada dimensión. Los eigenvalues corresponden a la cantidad de varianz aexplicada por cada componente principal. Seleccionamos las primera 15 dimensiones al haber establecido como criterio descartar aquellas dimensiones que explican un porcentaje de la varianza inferior al 1 %.

eig.val <- get_eigenvalue(res.pca)
knitr::kable(eig.val[1:16,])

Graficamos el porcentaje de varianza explicada frente al número de dimensiones.

fviz_eig(res.pca,ncp=15, addlabels = TRUE, ylim = c(0, 60))

Análisis de las variables

Calidad de la representación de las variables(cos2)

A continuación, vamos a visualizar el cos2 de las variables de las primeras 15 dimensiones utilizando un corrplot:

library("corrplot")
corrplot(res.pca$var$cos2, is.corr=FALSE)

# Total cos2 of variables on Dim.1 and Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2,xtickslab.rt = 90,top=10)
fviz_cos2(res.pca, choice = "var", axes = c(1, 3),xtickslab.rt = 90,top=10)
fviz_cos2(res.pca, choice = "var", axes = c(2, 3),xtickslab.rt = 90,top=10)
fviz_cos2(res.pca, choice = "var", axes = 1:5,xtickslab.rt = 90,top=10)

Gráfico contribución de las variables

A continuación, vamos a visualizar el la contribución de las variables a las primeras 13 dimensiones utilizando un corrplot:

corrplot(res.pca$var$contrib, is.corr=FALSE)

fviz_pca_var(res.pca, col.var = "contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

A continuación, mostramos las 10 variables que más contribuyen a explicar las primeras 15 dimensiones (porcentaje de la varianza acumulada del 90%)..

fviz_contrib(res.pca, choice = "var", axes = 1:15,xtickslab.rt = 90,top=10)

Análisis de los individuos

Análisis coordenadas, calidad y contribución de los individuos

A continuación, mostramos los valores de coordenadas, cos2 y contribución de los individuos para cada dimensión.

ind <- get_pca_ind(res.pca)  #Extraemos resultados para los individuos
# # Coordinates of individuals
#head(ind$coord)
knitr::kable(head(ind$coord), caption = "Coordenadas individuos")
# # Quality of individuals
#head(ind$cos2)
knitr::kable(head(ind$cos2), caption = "Calidad individuos")
# # Contributions of individuals
# head(ind$contrib)
knitr::kable(head(ind$coord), caption = "Contribucion individuos")

A continuación, graficamos la calidad de la representación para los individuos en las dimensiones 1 y 2, y 1 y 3.

fviz_pca_ind(res.pca, col.ind = "cos2", repel = TRUE)
fviz_pca_ind(res.pca, axes = c(1, 3),col.ind = "cos2", repel = TRUE)
fviz_pca_ind(res.pca, col.ind = data$annotation, pointsize = "cos2",
             repel = TRUE # Avoid text overlapping (slow if many points)
             )

A continuación, graficamos la contribución para los individuos en las dimensiones 1 y 2, y 1 y 3.

fviz_pca_ind(res.pca, col.ind = "contrib")+
  scale_color_gradient2(low="white", mid="blue", high="red", midpoint=0.5)
fviz_pca_ind(res.pca, axes = c(1, 3),col.ind = "contrib")+
  scale_color_gradient2(low="white", mid="blue", high="red", midpoint=0.5)

A continuación, mostramos los individuos en orden decreciente de contribución a las dimensiones 1 y 2, mostrando información sobre las 5 primeras.

# Total contribution on PC1 and PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2, xtickslab.rt = 90,top=50, )

knitr::kable(data %>%
               slice(c(193, 192, 272, 285))%>%
               select(phylum, class, species, extlink))

Biplot

#Biplot
fviz_pca_biplot(res.pca, label = "var", habillage=as.factor(data$annotation),
               addEllipses=TRUE, ellipse.level=0.95,
               ggtheme = theme_minimal())

PCA 3-D

library(plotly)
df_pca <- data.frame(dim1 = res.pca$ind$coord[,1],
                     dim2 = res.pca$ind$coord[,2],
                     dim3 = res.pca$ind$coord[,3],
                     ann = as.factor(data$annotation),
                     text = paste("Phylum: ", data$phylum, 
                        "\nClass: ", data$class,
                       "\nSpecie: ", data$species,
                        #"\nCharge: ", data$charge,
                        #"\nSeq_length: ", data$seq_length,
                        "\nacc: ", data$acc,
                        "\nAnnotation: ", data$annotation,
                        "\nPreviuos Annotation: ", data$annotation_original, sep=""))%>%
  mutate(ann_color = case_when(ann == "fascin" ~ "blue",
          ann == "fascin1" ~ "red",
          ann == "fascin2" ~ "green",
          ann == "singed" ~ "gray",
          ann == "uncharacterized" ~ "black"),)


fig <- plot_ly(df_pca, x = ~dim1, y = ~dim2, z = ~dim3, color = ~ann, text=~text)
fig <- fig %>% add_markers()
fig <- fig %>% layout(scene = list(xaxis = list(title = 'Dim.1'),
                     yaxis = list(title = 'Dim.2'),
                     zaxis = list(title = 'Dim.3')))
fig

library(rgl)
library(glue)
plot3d(df_pca$dim1, df_pca$dim2, df_pca$dim3, col = df_pca$ann_color, size = 4,
       xlab = glue("Dimensión 1 ({axis1}%)", axis1 = eig.val[1,2] %>% round(1)),
       ylab = glue("Dimensión 2 ({axis2}%)", axis2 = eig.val[2,2]%>% round(1)),
       zlab = glue("Dimensión 3 ({axis3}%)", axis3 = eig.val[3,2]%>% round(1)))
rglwidget()