Puckers

Function to calculate the pseudorotation angle nu for a single residue

In [1]:

res.nuR <- function(pdb.in, res) {

  # trim pdb to residue
  x <- trim.pdb(pdb.in, resno = res)

  # number of frames
  frames <- nrow(x$xyz)

  # residue name
  residue <- unique(x$atom$resid)

  nu0.sel <- atom.select(x, elety = c("C4'", "O4'", "C1'", "C2'"))
  nu1.sel <- atom.select(x, elety = c("O4'", "C1'", "C2'", "C3'"))
  nu2.sel <- atom.select(x, elety = c("C1'", "C2'", "C3'", "C4'"))
  nu3.sel <- atom.select(x, elety = c("C2'", "C3'", "C4'", "O4'"))
  nu4.sel <- atom.select(x, elety = c("C3'", "C4'", "O4'", "C1'"))

  pdb.nu0 <- trim.pdb(x, trim = nu0.sel)
  pdb.nu1 <- trim.pdb(x, trim = nu1.sel)
  pdb.nu2 <- trim.pdb(x, trim = nu2.sel)
  pdb.nu3 <- trim.pdb(x, trim = nu3.sel)
  pdb.nu4 <- trim.pdb(x, trim = nu4.sel)


  atom.match.0 <- match(pdb.nu0$atom$elety, c("C4'", "O4'", "C1'", "C2'"))
  atom.match.1 <- match(pdb.nu1$atom$elety, c("O4'", "C1'", "C2'", "C3'"))
  atom.match.2 <- match(pdb.nu2$atom$elety, c("C1'", "C2'", "C3'", "C4'"))
  atom.match.3 <- match(pdb.nu3$atom$elety, c("C2'", "C3'", "C4'", "O4'"))
  atom.match.4 <- match(pdb.nu4$atom$elety, c("C3'", "C4'", "O4'", "C1'"))

  # function to calculate the torsion angle nu of frame frame
  nux <- function(frame, nu){

    xyz <- get(paste0('pdb.nu', nu))$xyz[frame,]

    xyz.reord <- numeric()
    for (frame in get(paste0('atom.match.', nu))) { #add the coordinates in the correct order
      start_index <- 3 * (frame - 1) + 1
      xyz.reord <- c(xyz.reord, xyz[start_index:(start_index + 2)])
    }

    return(torsion.xyz(xyz.reord))

  }

  expand.grid(frame = 1:frames, nu = 0:4) %>%
    mutate(angle = unlist(map2(frame, nu, ~ nux(.x, .y)))) %>%
    pivot_wider(names_from = nu, values_from = angle) %>%
    mutate(residue = paste0(residue, res), .before = 1) %>%
    rename(
      nu0 = `0`,
      nu1 = `1`,
      nu2 = `2`,
      nu3 = `3`,
      nu4 = `4`
    )

}

Function to calculate the pseudorotation angle nu for a range of residues

In [2]:

nuR <- function(pdb.in, chain = NULL, res.start = NULL, res.end = NULL){

  # pull chain names if none supplied
  if (is.null(chain)) {
    chain <- unique(pdb.in$atom$chain)
  }

  # trim pdb to chain
  x <- trim.pdb(pdb.in, chain = chain)

  # pull residue min from selected chain if none supplied
  if (is.null(res.start)) {
    res.start <- min(x$atom$resno)
  }

  # pull residue max from selected chain if none supplied
  if (is.null(res.end)) {
    res.end <- max(x$atom$resno)
  }

  map(res.start:res.end, ~ res.nuR(x, .)) %>%
    bind_rows()

}

Function to format the output of nuR

In [3]:

puck.formatR <- function(nuR.output) {
  nuR.output %>%
    group_by(residue, frame) %>%
    mutate(
      P0 = atan2(nu4 + nu1 - nu3 - nu0, 2.0 * nu2 * (sin(pi/5) + sin(pi/2.5))),
      amplitude = nu2/cos(P0),
      angle = 180/pi * P0,
      angle = (angle + 360) %% 360
    ) %>%
    group_by(residue) %>%
    mutate(
      residue = str_replace_all(residue, 'D|3_|5_|_', ''),
      mean.angle = mean(angle),
      mod.angle = angle,
      # mod.angle = if_else(
      #   angle > 324 & mean.angle > 0, angle - 360, angle
      # ),
      mean.angle = median(mod.angle),
      mean.amplitude = median(amplitude),
      puck = case_when(
        mean.angle > 0 & mean.angle < 36 ~ "C3'-endo",
        mean.angle > 36 & mean.angle < 72 ~ "C4'-exo",
        mean.angle > 72 & mean.angle < 108 ~ "O4'-endo",
        mean.angle > 108 & mean.angle < 144 ~ "C1'-exo",
        mean.angle > 144 & mean.angle < 180 ~ "C2'-endo",
        mean.angle > 180 & mean.angle < 216 ~ "C3'-exo",
        mean.angle > 216 & mean.angle < 252 ~ "C4'-endo",
        mean.angle > 252 & mean.angle < 288 ~ "O4'-exo",
        mean.angle > 288 & mean.angle < 324 ~ "C1'-endo",
        mean.angle > 324 & mean.angle < 360 ~ "C2'-exo"
      )
    )
}

Function to plot the output of puck.formatR as a polar plot

In [4]:

puck.plotR <- function(puck.formatR.output, alpha = 0.01) {
  puck.formatR.output %>%
    ggplot(., aes(
      angle, amplitude,
      color = factor(
        puck,
        levels = c(
          "C3'-endo", "C4'-exo", "O4'-endo", "C1'-exo", "C2'-endo",
          "C3'-exo", "C4'-endo", "O4'-exo", "C1'-endo", "C2'-exo"
        )
      )
    )
    ) +
    facet_wrap(
      ~factor(
        residue, levels = puck.formatR.output %>%
          arrange(mean.angle) %>%
          pull(residue) %>%
          unique())
    ) +
    geom_point(alpha = alpha) +
    geom_point(
      data = puck.formatR.output %>%
        select(residue, mean.amplitude, mean.angle, puck) %>%
        unique(),
      aes(mean.angle, mean.amplitude, fill = factor(
        puck,
        levels = c(
          "C3'-endo", "C4'-exo", "O4'-endo", "C1'-exo", "C2'-endo",
          "C3'-exo", "C4'-endo", "O4'-exo", "C1'-endo", "C2'-exo"
        )
      )),
      size = 2, shape = 21, color = 'black'
    ) +
    scale_x_continuous(limits = c(0,360), breaks = 36*(0:10)) +
    scale_y_continuous(expand = c(0,0)) +
    scale_color_manual(name = 'pucker', values = color.d3, drop = FALSE) +
    scale_fill_manual(name = 'pucker', values = color.d3, drop = FALSE) +
    ggthemes::theme_pander() +
    theme(
      axis.text.x = element_text(size = 12 * scaling, face = 'bold'),
      axis.text.y = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      strip.text.x = element_markdown(size = 14 * scaling, face = 'bold'),
      strip.text.y = element_markdown(size = 14 * scaling, face = 'bold'),
      strip.background = element_blank(),
      panel.grid.major.x = element_line(
        linewidth = 0.5 * scaling, color = 'grey',
        linetype = 'dotted'
      ),
      panel.grid.minor.x = element_blank(),
      panel.grid.major.y = element_line(
        linewidth = 0.5 * scaling, color = 'grey',
        linetype = 'dotted'
      ),
      panel.grid.minor.y = element_blank(),
      panel.border = element_blank(),
      panel.background = element_blank(),
      plot.title = element_markdown(size = 16 * scaling, face = 'bold'),
      plot.subtitle = element_markdown(size = 14 * scaling),
      plot.caption = element_markdown(hjust = 1, size = 10 * scaling)
    ) +
    coord_polar(clip = 'off')
}