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

knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(tidyverse)
library(magrittr)
library(hrbrthemes)

El objetivo de este cuaderno es el de analizar las estructuras experimentales disponibles para FSCN1 y analizarlas en base a la literatura científica publicada.

Revisión bibliográfica

En primer lugar, fruto de la revisión bibliográfica, hemos podido determinar que existe cierto consenso en cuanto a que FSCN1 tiene 3 sitios de unión a actina, los cuales pueden visualizarse en la Figura 1. Además, se sabe que es una estructura altamente dinámica y que ha sido cristalizada tanto en su conformación activa como inactiva (Wu et al., 2021). Se ha propuesto que en la unión a actina sería muy importante tanto la carga neta como la participación de residuos cargados que podrían formar puentes salinos en la propia fascina (estabilizando una conformación inactiva) y con la actina (acercándola y favoreciendo la unión)(Wu et al., 2021). Se ha propuesto, además, un modelo según el cual la fascina interaccionaría con la actina de 2 modos distintos, uniéndose a un único filamento de F-actina con los sitios de unión 1 y 2 y con un segundo haz con el sitio de unión 3 (Aramaki et al., 2016).

En ese sentido, se han estudiado el efecto que podría tener la fosforilización de distintos residuos sobre la actividad de la fascina. La literatura parece apuntar que fosforilaciones en puntos clave, como S39, sitio de unión de PKC, están conservados evolutivamente, disminuyen la capacidad de unión a la actina, probablemente al estabilizar una conformación inactiva, y que tienen asociadas el favorecer el resto de funciones "no canónicas" de la fascina como la interacción con otras proteínas o la localización celular, donde se conoce como modulador epigenético (Lamb & Tootle, 2020).

knitr::include_graphics(path = "../reports/fscn1_binding_sites/fscn1_binding_sites_S39.png",)

Actualmente existen un total de 18 estructuras experimentales publicadas. Después de realizar una revisión bibliográfica de los artículos donde fueron publicadas y, atendiendo a los resultados de Yang et al. (2013), hemos elaborado la siguiente tabla donde se muestra la información relevante al respecto.

info.tabla <- readxl::read_excel(path = "../reports/pbds/pdbs.xlsx") 
info.tabla%>%
  select(`PDB`, `Año de publicación`, `Resolución (Å)`,`Anotación`, `Conformación`) %>%
  knitr::kable(align = 'c',caption = "Estructuras experimentales disponibles de FSCN1 y su conformación según Yang et al. (2013). En los ensayos de actividad 6I0Z no mostro una reduccion de la misma.")

Análisis estructural usando bio3d

Incluyendo estructuras acomplejadas con inhibidores

A continuación, vamos a descargar las estructuras del PDB y a separarlas por cadenas:

library(bio3d)
raw.files <- get.pdb(info.tabla$PDB, path = "raw_pdbs", gzip=TRUE, verbose = FALSE)
files <- pdbsplit(raw.files, ids = info.tabla$PDB, path = "raw_pdbs/split_chain",
  ncore=4, verbose = FALSE)
files

pdbs.all.chains <-  pdbaln(files)

Procedemos a comprobar si existen confórmeros redundantes entre sí, para lo cual filtraremos usando RMSD. Como se puede observar a continuación, con el umbral que habíamos establecido, todas las estructuras son significativamente distintas. Esto es consistente con la literatura científica consultada.

rd <- pdbs.all.chains %$%
  filter.rmsd(xyz, cutoff=0.1, fit=TRUE)
pdbs <- trim(pdbs.all.chains, row.inds=rd$ind)
pdbs$id <- pdbs %$%
  unlist(strsplit(basename(id[rd$ind]), split=".pdb"))
pdbs$id

Por último, vamos a anotar las estructuras, para facilitar luego la interpretación de los resultados, según la Tabla 1.

pdbs$ann <- pdbs %$%
  tibble(type = case_when(
  str_detect(id, "1DFC")~ "wild-type",
  str_detect(id, "3LLP")~ "wild-type",
  str_detect(id, "3P53")~ "wild-type",
  str_detect(id, "4GO")~ "mutante",
  str_detect(id, "4GP")~ "mutante",
  str_detect(id, "6B0T")~ "complejo_inhibidor",
  str_detect(id, "6I0Z")~ "complejo_inhibidor",
  str_detect(id, "6I1")~ "complejo_inhibidor",
  TRUE ~ "NA"
  ),
  conf = case_when(
  str_detect(id, "1DFC")~ "activa",
  str_detect(id, "6I0Z")~ "activa",
  id == "3LLP_A" ~ "activa",
  TRUE ~ "inactiva"
  ))
cbind(pdb = pdbs$id, pdbs$ann) %>% head()

Análisis de componente principal

El primer análisis que vamos a realizar es un análisis de componente principal. Para que los resultados que obtengamos sean significativos, el primer paso será superponer todas las estructuras en el "núcleo invariante" de las mismas (Grant et al., 2006).

core <- core.find(pdbs, verbose = FALSE)
# superimpose all structures to core
pdbs$xyz = pdbfit(pdbs, core$c0.5A.xyz)
plot(core)
write.pdb(xyz=pdbs$xyz[3,core$xyz],
          file="3LLP_A_fit-core_all.pdb")

En la Figura 2 se muestra estructura tridimensional de 3LLP A superpuesta con la región invariante de las estructuras experimentales de FSCN1.

knitr::include_graphics(path = "../reports/core_all/core_all.png",)

A continuación, mostramos los resultados del PCA:

# identify gaps, and perform PCA
gaps.pos <- gap.inspect(pdbs$xyz)
gaps.res <- gap.inspect(pdbs$ali)
pc.xray <- pca.xyz(pdbs$xyz[,gaps.pos$f.inds])
plot(pc.xray)

Como podemos observar, parece que las estructuras se agrupan de forma distanta. En la siguiente imagen mostramos la proyección de las 31 estructuras sobre las 2 primeras dimensiones del PCA (varianza total acumulada superior al 90%) anotadas según el tipo de estructura (mutante, wild-type o acomplejada) y su conformación (activa o inactiva).

d1 <- project.pca(pdbs$xyz[, gaps.pos$f.inds], pc.xray)
pdbs %$%
  tibble(pdb.id = id,dim1 = d1[,1], dim2 = d1[,2],
    type = ann$type, conf = ann$conf) %>%
  ggplot(aes(x=dim1, y=dim2, color = type, shape = conf)) +
    geom_point(size=3)+
  scale_colour_ipsum()+
  labs(title = "PCA sobre las estucturas experimentales de FSCN1", x = "PC1 (65.31%)", y = "PC2 (24.97%)")

Los resultados son muy interesantes, se puede observar que:

1. Aquellas estructuras que corresponden a estructuras donde la FSCN1 está acomplejadas con el inhibidor forman un grupo muy separado de las FSCN1 mutantes o wild-type.
2. Las estructuras que corresponden a las dos subunidades de 6I0Z se agrupan con las estructuras mutantes o wild type. Esto es muy interesante porque en los ensayos bioquimicos esta estructura, aunque acomplejada con un analogo que, al menos, pretendia ser un inhibidor, mostro actividad de union a FSCN1.

Para interpretar estos resultados nos será de utilidad visualizar la contribución de cada residuo a cada una de las 2 dimensiones principales (en negro la conttribución a la primera dimension y en azul a la segunda).

plot.bio3d(pc.xray$au[,1], ylab="PC (A)", xlab="Residue Position", typ="l")
points(pc.xray$au[,2], typ="l", col="blue")

Para visualizar las principales variaciones estructurales podemos usar la función mktrj() para generar un archivo PDB de trayectoria interpolando a lo largo del eigenvector 1 y que podría visualizarse usando un software específico.

p1 <- mktrj(pc.xray, pc=1,
      resno=pdbs$resno[1, gaps.res$f.inds],
      resid=pdbs$resid[1, gaps.res$f.inds])
write.ncdf(p1, "trj_pc1_all.nc")

Excluyendo estructuras acomplejadas con inhibidores

En este apartado vamos a repetir y ampliar el análisis realizado en el apartado anterior, pero excluyendo a las estructuras acomplejadas con inhibidores:

index <- which(!pdbs$ann$type == "complejo_inhibidor")
pdbs.noInh <- trim(pdbs, row.inds= index)
pdbs.noInh$ann <- pdbs %$%
  tibble(type = ann$type[index], conf = ann$conf[index])

Procedemos a calcular el core de estas secuencias y a visualizarlo. Como era de esperar en este caso el número de residuos que conforman este core es muy superior al caso anterior:

core <- core.find(pdbs.noInh, verbose = FALSE)
# superimpose all structures to core
pdbs.noInh$xyz = pdbfit(pdbs.noInh, core$c0.5A.xyz)
plot(core)
write.pdb(xyz=pdbs.noInh$xyz[3,core$xyz],
          file="3LLP_A_fit-core_sinInh.pdb")

En la Figura 3 mostramos la región invariable sobre la estructura 3LLP A. A diferencia que en el caso anterior, en este caso la región invariante comprende los dominios 3 y 4, mientras que en el core incluyendo la presencia de inhibidores, la región incluía únicamente el dominio 2. Esto podría ser consistente con el hecho de que en literatura las modificaciones que parecen estabilizar una conformación inactivada afectan a los dominios 1 y 2.

knitr::include_graphics(path = "../reports/core_all/core_sinInh.png",)

Ahora, después de haber superpuesto todas las estructuras a la región invariante, vamos a mostrar el RMSF para cada residuo tomando de referencia la estructura 1DFC_A para las estructuras secundarias (alfa-hélices se muestran coloreadas en rojo y hebras beta coloreadas en azul).

# identify gaps, and perform PCA
gaps.pos <- gap.inspect(pdbs.noInh$xyz)
gaps.res <- gap.inspect(pdbs.noInh$ali)

# Tailor the PDB structure to exclude gap positions for SSE annotation
pdb.1DFC_A<- read.pdb("raw_pdbs/split_chain/1DFC_A.pdb")
id <- grep("1DFC_A", pdbs.noInh$id)
inds <- atom.select(pdb.1DFC_A, resno = pdbs.noInh$resno[id, gaps.res$f.inds])
ref.pdb <- trim.pdb(pdb.1DFC_A, inds = inds)

# Plot RMSF with SSE annotation and labeled with residue numbers (Figure 8.)
rf <- rmsf(pdbs.noInh$xyz[, gaps.pos$f.inds])
plot.bio3d(rf, sse=ref.pdb, ylab="RMSF (Å)", 
           xlab="Residue No.", typ="l",
           helix.col = "blue", sheet.col = "red")

Procedemos a realizar el PCA con los resultados:

pc.xray.sinInh <- pca.xyz(pdbs.noInh$xyz[,gaps.pos$f.inds])
plot(pc.xray.sinInh)

De nuevo, parece que las estructuras forman distintos grupos una vez proyectadas sobre las primeras dos dimensiones. Procedemos a representarlas de forma más precisa:

d2 <- project.pca(pdbs.noInh$xyz[, gaps.pos$f.inds], pc.xray.sinInh)
pdbs.noInh %$%
  tibble(pdb.id = id,dim1 = d2[,1], dim2 = d2[,2],
    type = ann$type, conf = ann$conf) %>%
  ggplot(aes(x=dim1, y=dim2, color = type, shape = conf, label = pdb.id)) +
    geom_point(size=3)+
  geom_text(hjust=0, vjust=1)+
  scale_colour_ipsum()+
  labs(title = "PCA sobre las estucturas experimentales de FSCN1 sin inhibidores",
    x = "PC1 (80.45%)", y = "PC2 (11.91%)")


pdbs.noInh %$%
  tibble(pdb.id = id,dim1 = d2[,1], dim3 = d2[,3],
    type = ann$type, conf = ann$conf) %>%
  ggplot(aes(x=dim1, y=dim3, color = type, shape = conf, label = pdb.id)) +
    geom_point(size=3)+
  geom_text(hjust=0, vjust=1)+
  scale_colour_ipsum()+
  labs(title = "PCA sobre las estucturas experimentales de FSCN1 sin inhibidores",
    x = "PC1 (80.45%)", y = "PC3 (3.09%)")

Analizando el gráfico anterior, podemos concluir que:

1. Las subunidades 4GOYB (mutante en sitio de union 1 e inactivada) y 4GOV_B (mutante fosfomimetica de fosforilacion en S39) son muy parecidas entre sí y con 3P53_B (conformacion inactiva wild-type); y parecen formar un grupo cuando se proyectan en las primeras dos dimensiones junto con 4GPO_B (mutante en sitio de union 2 e inactivada) y 3LLP_B (conformacion inactiva wild-type).
2. Las subunidades que corresponden a 1DFC (wild-type conformación activa) parece agruparse entre sí en el cuadrante superior derecho.
3. Sin embargo, la tercera subunidad activa (y que habíamos tomado como referencia, al ser la subunidad activa con mayor resolución), se agrupa cerca de 4GPO_A y 4GP3_A (estructuras inactivadas por una mutación en el sitio de union 2 y 3 respectivamente), y parece estar agrupada con 4GOYA y 4GOV_A,
4. De nuevo, las subunidades 4GOYA y 4GOV_A, respectivamente mutantes en sitio 1 y fosfomimética, vuelven a ser muy similares entre sí.
5. Las estructuras mutantes y wild-type parecen separarse a lo largo de de la 3ª dimensioón del PCA, pero esta varianza explica un porcentaje muy pequeño de la varianza total.

En la siguiente imagen se muestra la contribución de cada residuo a los 3 componentes principales:

par(mfrow = c(3, 1), cex = 0.75, mar = c(3, 4, 1, 1))
plot.bio3d(pc.xray.sinInh$au[,1], sse=ref.pdb, ylab="PC1",
  helix.col = "blue", sheet.col = "red")
plot.bio3d(pc.xray.sinInh$au[,2], sse=ref.pdb, ylab="PC2",
  helix.col = "blue", sheet.col = "red")
plot.bio3d(pc.xray.sinInh$au[,3], sse=ref.pdb, ylab="PC3",
  helix.col = "blue", sheet.col = "red")

Por último, vamos a realizar un clustering jerárquico, usando distancias euclídeas y el método ward.D2, para observar cómo se agrupan las estructuras según las 3 primeras componentes del PCA.

hc <- dist(pc.xray.sinInh$z[,1:3], method = "euclidean") %>%
  hclust(method = "ward.D2")
hc$labels <- pdbs.noInh$id
hcd <- as.dendrogram(hc)

nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), cex = 0.7)
plot(hcd,  xlab = "", horiz = TRUE, nodePar = nodePar)

Conclusión

Hasta este punto, el análisis realizado no parece corroborar lo propuesto por Yang et al. (2013) sobre qué estructuras corresponden a una conformación activa e inactiva. Esto me hace pensar que quizás no estoy utilizando una metodología correcta. Funciona muy bien, sin embargo, con los inhibidores, no solo separando aquellas estructuras acomplejadas con un inhibidor del resto sino, también, separando la estructura de 6I0Z que estaba acomplejada con un análogo de un inhibidor pero sí presentaba actividad del resto de estructuras acomplejadas.

currocam/FascinRSCA documentation built on March 21, 2022, 6:29 a.m.