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

En este cuaderno se comparan las estructuras de fascina en humanos: 3p53 (una versión "limpiada") de fascina 1 experimental, una estructura de fascina 2 canónica predecida por AlphaFold2 (por Deepmind) y dos estructuras de la isoforma larga de fascina 2 predecidas por AlphaFold2 y RoseTTAfold. El objetivo es explorar las distintas posibilidades de la librería bio3d para el análisis estructural y de la librería r3dmol para realizar representaciones interactivas de proteínas dentro del entorno de los cuadernos de R markdown. Se han excluido aquellas funcionalidades relativas al análisis de trayectorias de simulaciones de dinámica molecular y de Normal Mode Analysis.

knitr::opts_chunk$set(echo = TRUE)
library(bio3d)
library(r3dmol)

En primer lugar, vamos a cargar los 4 archivos PDB.

fascin2b.AlphaFold <- read.pdb("../inst/extdata/structural_modelling/processed/AlphaFold2/O14926-2/O14926_AlphaFold2.pdb")
fascin2b.RoseTTaFold <- read.pdb("../inst/extdata/structural_modelling/processed/RoseTTAFold/O14926-2/O14926-2_RoseTTAFold.pdb")
fascin2a.AlphaFold <- read.pdb("../inst/extdata/structural_modelling/external/O14926-1/AF-O14926-F1-model_v2.pdb")
fascin1 <- read.pdb("../inst/extdata/structural_modelling/external/fascin1/3p53_clean.pdb")

Empezaremos evaluando la calidad de los mismos. Para ello, graficamos el factor B (para estructuras experimentales) y el pLDDT (medida de la confianza por residuo en una escala de 1-100 en estructuras obtenidas in silico). Estas gráficas nos pueden indicar tanto zonas intrínsecamente desestructuradas como zonas donde la predicción / determinación no es buena.

plot.bio3d(fascin1$atom$b[fascin1$calpha], typ="l", ylab="Factor B",
  main= "Factor B por residuo para Fascina 1",
  sub = "Estructura experimental 3P53 limpiada manualmente.")

plot.bio3d(fascin2a.AlphaFold$atom$b[fascin2a.AlphaFold$calpha], typ="l", ylab="pLDDT",
  main= "Estimación de la confianza por residuo de Fascina 2 canónica",
  sub="Estructura predecida por AlphaFold y depositada en base de datos.")

plot.bio3d(fascin2b.AlphaFold$atom$b[fascin2b.AlphaFold$calpha], typ="l", ylab="pLDDT",
  main= "Estimación de la confianza por residuo de Fascina 2 isoforma larga",
  sub = "Estructura predecida por AlphaFold2 usando Google Collab. Elaboración propia.")

plot.bio3d(fascin2b.RoseTTaFold$atom$b[fascin2b.RoseTTaFold$calpha], typ="l", ylab="pLDDT",
  main= "Estimación de la confianza por residuo de Fascina 2 isoforma larga",
  sub = "Estructura predecida por RoseTTAFold usando Google Collab. Elaboración propia.")

Comenzamos explorando las funciones más básicas relativas a comparar estructuras.

Sequence and structural similarity

En primer lugar, vamos a alinear los archivos.

pdbs_list <- list(fascin1, fascin2a.AlphaFold, fascin2b.RoseTTaFold, fascin2b.AlphaFold)
names(pdbs_list) <- c("Fascin1.Exp", "Fascin2a.AlphaFold","Fascin2b.RoseTTAFold", "Fascin2b.AlphaFold")
pdbs <- pdbaln(pdbs_list,fit = FALSE)
pdbs

Calculamos el grado de identidad de las secuencias por pares.

seqidentity(pdbs)

Procedemos a calcular una matriz con los valores de RMSD por pares (después de realizar una superposición usando el algoritmo por defecto, Kabsch).

rmsd(pdbs, fit=TRUE)

A continuación, realizamos un PCA. En él se puede observar cómo la estructura predecida por RoseTTAFold está muy alejada de las otras 3 (valores muy distintos al ser proyectada en el eje del PC1, que explica más del 95% de la varianza). En este gráfico, y en el resto del cuaderno, se emplea la siguiente leyenda de colores:

3P53: negro
O14926-1 (AlphaFold2): rojo
O14926-2 (RoseTTAFold): verde
O14926-2 (AlphaFold2): azul

pc <- pca(pdbfit(pdbs), rm.gaps=TRUE)
plot(pc, col=c("black", "red", "green", "blue"))

A continuación, y para ilustrar la potencia de la librería en cuanto al NMA, mostramos unos resultados rápidos:

modes <- nma(pdbs)
plot(modes, pdbs, spread=TRUE)

Superposición de coordenadas y alineamiento estructural

A continuación, vamos a explorar distintos métodos de superposición propuestos en la documentación.

Superposicion usando Kabsch algorithm

A continuación, vamos a emplear la función pdbfit para superponer las estructuras. Esta funcion emplea el Kabsch algorithm, es una superposición sequence based y es la empleada por defecto cuando fit = TRUE.

xyz.Kabsch <- pdbfit(pdbs)
xyz.Kabsch

Superposición con ciclos de refinamiento y empleando MUSCLE

A continuación, vamos a emplear la función struct.aln() que realiza un alineamiento de secuencias empleando MUSCLE y refina la superposición eliminando átomos con una desviación estructural muy alta. Este método solo se puede realizar con 2 archivos pdb y, hasta donde he podido entender, no hay una forma sencilla de emplear otro programa para los alineamientos (que MUSCLE sería, por lo que me comentó Coral, el más adecuado para estas secuencias).

RoseTTAFold Vs AlphaFold2.

Vamos a emplearlo para superponer la estructura de RoseTTAFold a la de AlphaFold2.

aln.fascin2b <- struct.aln(fixed = fascin2b.AlphaFold,
  mobile = fascin2b.RoseTTaFold)
fascin2b.RoseTTaFold$xyz <- aln.fascin2b$xyz

Ahora vamos a escribir la estructura resultante en un archivo que podríamos visualizar con precisión usando un programa específico.

w <- write.pdb(cat.pdb(fascin2b.RoseTTaFold, fascin2b.AlphaFold, rechain= TRUE,
  renumber=TRUE), file="O14926-2_RoseTTAFold_at_O14926_AlphaFold2.pdb", )

No obstante, y para explorar las posibilidades de la librería r3dmol para la visualización dentro del entorno de los cuadernos de R Markdown, vamos a construir un imagen interactiva.

r3dmol(                         # Set up the initial viewer
  viewer_spec = m_viewer_spec(
    cartoonQuality = 10,
    lowerZoomLimit = 50,
    upperZoomLimit = 350
  )
) %>%
  m_add_model(                  # Add model to scene
    data = "O14926-2_RoseTTAFold_at_O14926_AlphaFold2.pdb",
    format = "pdb"
  )%>%
  m_zoom_to() %>%               # Zoom to encompass the whole scene
  m_add_style(
    style = c(
      m_style_cartoon("green")),
    sel = m_sel(chain = "A", )
  ) %>%
  m_add_label(
    text = "RoseTTaFold",
    sel = m_sel(resi = 1:10, chain = "A"), 
    style = m_style_label(fontSize = 8,fontColor = "green",
                          inFront = TRUE, backgroundOpacity = 0.9)
  )%>%
  m_add_style(
    style = c(
      m_style_cartoon(color = "blue")),
    sel = m_sel(chain = "B", )
  )%>%
  m_add_label(
    text = "AlphaFold2",
    sel = m_sel(resi = 1:10,chain = "B"), 
    style = m_style_label(fontSize = 8,fontColor = "blue",
                          inFront = TRUE, backgroundOpacity = 0.9)
  )%>%
  m_rotate(                     # Rotate the scene by given angle on given axis
    angle = 90,
    axis = "y"
  )

FSCN2-1 Vs FSCN2-2.

A continuación, vamos a superponer la fascina 2 canónica y la isoforma larga.

aln.fascin2 <- struct.aln(fixed = fascin2a.AlphaFold,
  mobile = fascin2b.AlphaFold)
fascin2b.AlphaFold$xyz <- aln.fascin2$xyz
w <- write.pdb(cat.pdb(fascin2a.AlphaFold, fascin2b.AlphaFold, rechain= TRUE),
          file="O14926-2_at_O14926-1_AlphaFold2.pdb")

Y procedemos a visualizar la estructura resultante:

r3dmol(                         # Set up the initial viewer
  viewer_spec = m_viewer_spec(
    cartoonQuality = 10,
    lowerZoomLimit = 50,
    upperZoomLimit = 350
  )
) %>%
  m_add_model(                  # Add model to scene
    data = "O14926-2_at_O14926-1_AlphaFold2.pdb",
    format = "pdb"
  )%>%
  m_zoom_to() %>%               # Zoom to encompass the whole scene
  m_add_style(
    style = c(
      m_style_cartoon("red")),
    sel = m_sel(chain = "A", )
  ) %>%
  m_add_label(
    text = "Fascin2a",
    sel = m_sel(resi = 100:110, chain = "A"), 
    style = m_style_label(fontSize = 8,fontColor = "red",
                          inFront = TRUE, backgroundOpacity = 0.9)
  )%>%
  m_add_style(
    style = c(
      m_style_cartoon(color = "blue")),
    sel = m_sel(chain = "B", )
  )%>%
  m_add_label(
    text = "Fascin2b",
    sel = m_sel(resi = 1:10,chain = "B"), 
    style = m_style_label(fontSize = 8,fontColor = "blue",
                          inFront = TRUE, backgroundOpacity = 0.9)
  )%>%
  m_rotate(                     # Rotate the scene by given angle on given axis
    angle = 90,
    axis = "y"
  )

FSCN1 Vs FSCN2-1.

A continuación, realizamos la superposición para la fascina 2a y la fascina 1.

aln.fascin <- struct.aln(fixed = fascin1,
  mobile = fascin2a.AlphaFold)
fascin2a.AlphaFold$xyz <- aln.fascin$xyz
w <- write.pdb(cat.pdb(fascin1, fascin2a.AlphaFold, rechain= TRUE),
            file="O14926-1_at_3P53.pdb")

y procedemos a visualizar el archivo resultante:

r3dmol(                         # Set up the initial viewer
  viewer_spec = m_viewer_spec(
    cartoonQuality = 10,
    lowerZoomLimit = 50,
    upperZoomLimit = 350
  )
) %>%
  m_add_model(                  # Add model to scene
    data = "O14926-1_at_3P53.pdb",
    format = "pdb"
  )%>%
  m_zoom_to() %>%               # Zoom to encompass the whole scene
  m_add_style(
    style = c(
      m_style_cartoon("gray")),
    sel = m_sel(chain = "A", )
  ) %>%
  m_add_label(
    text = "Fascin1",
    sel = m_sel(resi = 100:110, chain = "A"), 
    style = m_style_label(fontSize = 8,fontColor = "gray",
                          inFront = TRUE, backgroundOpacity = 0.9)
  )%>%
  m_add_style(
    style = c(
      m_style_cartoon(color = "red")),
    sel = m_sel(chain = "B", )
  )%>%
  m_add_label(
    text = "Fascin2a",
    sel = m_sel(resi = 1:10,chain = "B"), 
    style = m_style_label(fontSize = 8,fontColor = "red",
                          inFront = TRUE, backgroundOpacity = 0.9)
  )%>%
  m_rotate(                     # Rotate the scene by given angle on given axis
    angle = 90,
    axis = "y"
  )

Constructing Experimental Structure Ensembles for a Protein Family

Por último, vamos a explorar las posibilidades relativas a la construcción y al análisis de ensembles de estructuras de una misma familia de proteínas.

Comenzamos expandiendo el número de estructuras de proteínas disponibles empleando blast y buscando proteínas proteínas similares a fascina 1 (3P53).

pdb.3p53 <- read.pdb("3p53")
seq.3p53 <- pdbseq(pdb.3p53)
blast.3p53 <- blast.pdb(seq.3p53)

Examinamos los evalue de los hits y establecemos un valor de umbral.

hits.3p53 <- plot.blast(blast.3p53, cutoff=250)

Las estructuras seleccionadas son las siguientes.

unique(hits.3p53$pdb.id)

Procedemos a descargarlos y a alinear las secuencias:

files.3p53 <- get.pdb(hits.3p53, path="raw_pdbs", split = TRUE)
pdbs.3p53 <- pdbaln(files.3p53)

A continuación, vamos a emplear la función core.find() para identificar las posiciones que conforman la región rígida invariante de las estructuras. Para ello se realiza un proceso iterativo de superposición donde se excluyen progresivamente los residuos con mayores diferencias en cada iteración. La lista de residuos que conforman este core podemos emplearla, por ejemplo, para realizar la superposición únicamente en este core.

core <- core.find(pdbs.3p53)
col=rep("black", length(core$volume))
col[core$volume<2]="pink"; col[core$volume<1]="red"
plot(core, col=col)
core.inds <- print(core, vol=1.0)
xyz.3p53 <- pdbfit(pdbs.3p53, core.inds )

A continuación, mostramos el core "invariante" para fascina 1:

core A continuación, podemos repetir el análisis de RMSD pero únicamente para este core.

rd <- rmsd(xyz.3p53)
hist(rd, breaks=40, xlab="RMSD (Å)", main="Histograma de RMSD para el core de proteínas homólogas a fascina 1")

Realizamos un clustering a partir de la matriz de distancias RMSD y lo representamos como un dendograma.

hc.rd <- hclust(as.dist(rd))
pdbs.3p53$id <- substr(basename(pdbs.3p53$id), 1, 6)
hclustplot(hc.rd, k = 1, labels=pdbs.3p53$id, cex=0.5,
          ylab="RMSD (Å)", main="RMSD Cluster Dendrogram", fillbox=FALSE)