---
title: "Dopamine aptamer structure determination"
author: 
  - name: Eric Largy
    orcid: 0000-0002-6140-9788
    email: eric.largy@u-bordeaux.fr
    affiliations: ARNA, INSERM U1212, CNRS UMR 5320, Université de Bordeaux
license: CC BY-SA 4.0
abstract: |
  This page presents the determination of an aptamer-dopamine complex structure based on NMR data.
draft: false
date: today
bibliography: biblio/references.bib
number-sections: true
number-depth: 6
echo: false
execute:
  warning: false
  cache: true
# molstar: embed
# embed-resources: true
---

```{r}
#| label: dependencies 
#| cache: false
#| output: false

# load packages (install where necessary)
if (!requireNamespace("librarian", quietly = TRUE)) {
  install.packages("librarian")
}

librarian::shelf(tidyverse, data.table, dtplyr, patchwork,
                 ggsci, ggthemes, ggtext, ggrepel,
                 reticulate, rmarkdown, readxl, writexl, DT, 
                 arrow, purrr, httr2, bio3d, NbClust, magrittr,
                 martinasladek/convertr)

# source files and convert to notebooks for display in article
lapply(
  list('R/theme.R', 'R/rst.error.R', 'R/dssr.R',
       'R/chiR.R', 'R/nuR.R', 'R/Contact_analyzer.R'),
  function(x) {
    source(x)
    r_to_qmd(x, paste0('notebooks/', gsub('R/', '', x), '.qmd'))
  }
)

scaling <- 1

color.d3 <- c(
  "#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
  "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf"
)
```

```{r}
#| label: info

info.short <- read_table(
  'D:/Amber/dop/pmemd/out/4_prod_doprstslt01.out',
  col_names = FALSE,
  skip = 37,
  n_max = 8
) %>% 
  select(1) %>% 
  separate_wider_delim(delim = '=', col = 1, names = c('param', 'value')) %>% 
  mutate(value = as.numeric(gsub(',', '', value))) %>% 
  pivot_wider(names_from = param, values_from = value) %>% 
  summarise(
    time.per.frame = ntpr * dt / 1000, #ns
    total.t = time.per.frame * nstlim/ntpr
  )

info.long <- read_table(
  'D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200.out', 
  col_names = FALSE,
  skip = 37,
  n_max = 8
) %>% 
  select(1) %>% 
  separate_wider_delim(delim = '=', col = 1, names = c('param', 'value')) %>% 
  mutate(value = as.numeric(gsub(',', '', value))) %>% 
  pivot_wider(names_from = param, values_from = value) %>% 
  summarise(
    time.per.frame = ntpr * dt / 1000, #ns
    total.t = time.per.frame * nstlim/ntpr #ns
  )

info.long.unrst <- read_table(
  'D:/Amber/dop/pmemd/out/4_prod_dopunrstslt01_250.out', 
  col_names = FALSE,
  skip = 37,
  n_max = 8
) %>% 
  select(1) %>% 
  separate_wider_delim(delim = '=', col = 1, names = c('param', 'value')) %>% 
  mutate(value = as.numeric(gsub(',', '', value))) %>% 
  pivot_wider(names_from = param, values_from = value) %>% 
  summarise(
    time.per.frame = ntpr * dt / 1000, #ns
    total.t = time.per.frame * nstlim/ntpr #ns
  )

```

```{r}
#| label: minpdbdata
#| output: false

# Minimized structures---

#remove previous file
if (file.exists('images/multi.min.pdb')) {
  file.remove('images/multi.min.pdb')
}

# read all final minimization pdb files
min.pdb.list <- lapply(
  list.files(path = 'data', pattern = "5_min_dop.*\\.pdb$", full.names = TRUE), 
  read.pdb
)

# concatenate xyz coordinates
for (i in 1:length(min.pdb.list)) {
  if (i == 1) {
    xyz <- min.pdb.list[[i]]$xyz[1,]
  } else {
    xyz <- rbind(
      xyz,
      min.pdb.list[[i]]$xyz[1,]
    )
  }
}

# export multimodel pdb file
write.pdb(min.pdb.list[[1]], xyz = xyz, file = "data/input_multi.min.pdb")

# read multimodel pdb file, renumber, remove hydrogens and waters
min.pdb <- read.pdb("data/input_multi.min.pdb", multi = TRUE) %>% 
  convert.pdb(., renumber = TRUE,
              first.resno = 11,
              rm.h = FALSE,
              rm.wat = TRUE) %>% 
  trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))


# clean up of pdb file
min.pdb$atom <- min.pdb$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    )
  )

#align all frames
min.pdb$xyz <- fit.xyz(
  fixed = min.pdb$xyz[1,],
  mobile = min.pdb$xyz,
  fixed.inds = atom.select(min.pdb, elety = unique(min.pdb$atom$elety))$xyz,
  mobile.inds = atom.select(min.pdb, elety = unique(min.pdb$atom$elety))$xyz
)

# write pdb file for deposition
if (!file.exists('images/multi.min.pdb')) {
  write.pdb(min.pdb, file = "images/multi.min.pdb")
}

# re-import without hydrogends for data processing below
pdb <- read.pdb('images/multi.min.pdb', multi = TRUE) %>% 
  convert.pdb(., rm.h = TRUE)
```

```{r}
#| label: rMDdata
#| output: false

for (i in 1:15) {
  
  if (i < 10) {
    j <- paste0('0', i)
  } else {
    j <- i
  }
  
  # Load the minimized structures
  # remove water and hydrogen atoms
  pdb.object <- as.pdb(
    prmtop = read.prmtop(paste0('data/doprstslt', j, '.prmtop')), 
    crd = read.ncdf(paste0('D:/Amber/dop/pmemd/out/4_prod_doprstslt', j, '.nc'))
  ) %>% 
    convert.pdb(., renumber = TRUE,
                first.resno = 11,
                rm.h = TRUE,
                rm.wat = TRUE) %>% 
    trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  
  # clean up pdb
  pdb.object$atom <- pdb.object$atom %>% 
    mutate(
      # clean up residue names
      resid = case_when(
        resid == 'DT3' ~ 'DT',
        resid == 'DT5' ~ 'DT',
        TRUE ~ resid
      ),
      # move dopamine to chain B
      chain = if_else(resid == 'LDP', 'B', 'A'),
      # rename dopamine atoms
      elety = case_when(
        resid == 'LDP' & elety == 'OAH' ~ 'O1',
        resid == 'LDP' & elety == 'HAH' ~ 'HO1',
        resid == 'LDP' & elety == 'OAI' ~ 'O2',
        resid == 'LDP' & elety == 'HAI' ~ 'HO2',
        resid == 'LDP' & elety == 'CAB' ~ 'C1',
        resid == 'LDP' & elety == 'CAD' ~ 'C2',
        resid == 'LDP' & elety == 'HAD' ~ 'H2',
        resid == 'LDP' & elety == 'CAG' ~ 'C3',
        resid == 'LDP' & elety == 'CAC' ~ 'C4',
        resid == 'LDP' & elety == 'CAF' ~ 'C5',
        resid == 'LDP' & elety == 'HAF' ~ 'H5',
        resid == 'LDP' & elety == 'CAE' ~ 'C6',
        resid == 'LDP' & elety == 'HAE' ~ 'H6',
        resid == 'LDP' & elety == 'CAA' ~ 'C7',
        resid == 'LDP' & elety == 'HAA' ~ 'H71',
        resid == 'LDP' & elety == 'HAB' ~ 'H72',
        resid == 'LDP' & elety == 'CAJ' ~ 'C8',
        resid == 'LDP' & elety == 'HAJ' ~ 'H81',
        resid == 'LDP' & elety == 'HAK' ~ 'H82',
        resid == 'LDP' & elety == 'NAK' ~ 'N1',
        resid == 'LDP' & elety == 'HAL' ~ 'HN11',
        resid == 'LDP' & elety == 'HAM' ~ 'HN13',
        resid == 'LDP' & elety == 'HAN' ~ 'HN12',
        TRUE ~ elety
      ),
      elesy = if_else(
        is.na(elesy) | elesy == 'not found',
        case_when(
          elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
          elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
          elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
          elety %in% c(paste0("N", 1:99)) ~ 'N',
          elety %in% c("P", paste0("P", 1:99)) ~ 'P',
          TRUE ~ 'not found'
        ),
        elesy
      )
    )
  
  #align all frames
  pdb.object$xyz <- fit.xyz(
    fixed = pdb.object$xyz[1,],
    mobile = pdb.object$xyz,
    fixed.inds = atom.select(pdb.object, elety = unique(pdb.object$atom$elety))$xyz,
    mobile.inds = atom.select(pdb.object, elety = unique(pdb.object$atom$elety))$xyz
  )
  
  #assign to numbered pdb object
  assign(paste0('pdb.traj.', j), pdb.object)
  
  #export a pdb trajectory with only 1000 frames
  #to avoid very big files
  if (!file.exists(paste0('images/4_prod_doprstslt', j, '.traj.1000frames.pdb'))) {
    write.pdb(
      pdb = pdb.object,
      xyz = pdb.object$xyz[seq(1, 5000, by = 5),],
      file = paste0('images/4_prod_doprstslt', j, '.traj.1000frames.pdb')
    )
  }
  
}

# create a xyz.list that contains all xyz 
xyz.list <- lapply(1:15, function(i) {
  pdb.traj <- get(paste0('pdb.traj.', sprintf('%02d', i)))
  return(pdb.traj$xyz)
})

```

```{r}
#| label: rMDdatalong
#| output: false

pdb.long.files <- lapply(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_doprstslt01_200(\\_restart\\d*)?.nc$", 
    full.names = TRUE), 
  function(prmtop = 'data/doprstslt01.prmtop', crd) {
    as.pdb(
      prmtop = read.prmtop(prmtop),
      crd = read.ncdf(crd)
    ) %>%
      convert.pdb(., renumber = TRUE,
                  first.resno = 11,
                  rm.h = TRUE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }, 
  prmtop = 'data/doprstslt01.prmtop')

pdb.long <- pdb.long.files[[1]]
pdb.long$xyz <- do.call(rbind, lapply(pdb.long.files, function(x) x$xyz))

divider <- 5

pdb.long$xyz <- pdb.long$xyz[1:(5000*divider),] #trim the first 25000 frames (exactly one microsecond)

# clean up pdb
pdb.long$atom <- pdb.long$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

#align all frames
pdb.long$xyz <- fit.xyz(
  fixed = pdb.long$xyz[1,],
  mobile = pdb.long$xyz,
  fixed.inds = atom.select(pdb.long, elety = unique(pdb.long$atom$elety))$xyz,
  mobile.inds = atom.select(pdb.long, elety = unique(pdb.long$atom$elety))$xyz
)

###

#for reference, these are the imported files
# [1] "D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200.nc"         
# [2] "D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200_restart.nc" 
# [3] "D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200_restart2.nc"
# [4] "D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200_restart3.nc"
# [5] "D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200_restart4.nc"
# [6] "D:/Amber/dop/pmemd/out/4_prod_doprstslt01_200_restart5.nc"
```

```{r}
#| label: rMDdatalongh
#| output: false

pdb.long.h.files <- lapply(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_doprstslt01_200(\\_restart\\d*)?.nc$", 
    full.names = TRUE), 
  function(prmtop = 'data/doprstslt01.prmtop', crd) {
    as.pdb(
      prmtop = read.prmtop(prmtop),
      crd = read.ncdf(crd)
    ) %>%
      convert.pdb(., renumber = TRUE,
                  first.resno = 11,
                  rm.h = FALSE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }, 
  prmtop = 'data/doprstslt01.prmtop')

pdb.long.h <- pdb.long.h.files[[1]]
pdb.long.h$xyz <- do.call(rbind, lapply(pdb.long.h.files, function(x) x$xyz))

divider <- 5

pdb.long.h$xyz <- pdb.long.h$xyz[1:(5000*divider),] #trim the first 25000 frames (exactly one microsecond)

# clean up pdb
pdb.long.h$atom <- pdb.long.h$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

#align all frames
pdb.long.h$xyz <- fit.xyz(
  fixed = pdb.long.h$xyz[1,],
  mobile = pdb.long.h$xyz,
  fixed.inds = atom.select(pdb.long.h, elety = unique(pdb.long.h$atom$elety))$xyz,
  mobile.inds = atom.select(pdb.long.h, elety = unique(pdb.long.h$atom$elety))$xyz
)

```

```{r}
#| label: uMDdatalong
#| output: false


pdb.unrst.long.files <- lapply(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_dopunrstslt01_250(\\_restart\\d*)?.nc$", 
    full.names = TRUE), 
  function(prmtop = 'data/doprstslt01.prmtop', crd) {
    as.pdb(
      prmtop = read.prmtop(prmtop),
      crd = read.ncdf(crd)
    ) %>%
      convert.pdb(., renumber = TRUE,
                  first.resno = 11,
                  rm.h = TRUE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }, 
  prmtop = 'data/doprstslt01.prmtop')

pdb.unrst.long <- pdb.unrst.long.files[[1]]
pdb.unrst.long$xyz <- do.call(rbind, lapply(pdb.unrst.long.files, function(x) x$xyz))

unrst.divider <- length(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_dopunrstslt01_250(\\_restart\\d*)?.nc$")
)

pdb.unrst.long$xyz <- pdb.unrst.long$xyz[1:(5000*25/20*unrst.divider),] #trim the first 25000 frames (exactly one microsecond)

# clean up pdb
pdb.unrst.long$atom <- pdb.unrst.long$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

#align all frames
pdb.unrst.long$xyz <- fit.xyz(
  fixed = pdb.unrst.long$xyz[1,],
  mobile = pdb.unrst.long$xyz,
  fixed.inds = atom.select(pdb.unrst.long, elety = unique(pdb.unrst.long$atom$elety))$xyz,
  mobile.inds = atom.select(pdb.unrst.long, elety = unique(pdb.unrst.long$atom$elety))$xyz
)

###
###
```

```{r}
#| label: uMDdatalongh
#| output: false

pdb.unrst.h.long.files <- lapply(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_dopunrstslt01_250(\\_restart\\d*)?.nc$", 
    full.names = TRUE), 
  function(prmtop = 'data/doprstslt01.prmtop', crd) {
    as.pdb(
      prmtop = read.prmtop(prmtop),
      crd = read.ncdf(crd)
    ) %>%
      convert.pdb(., renumber = TRUE,
                  first.resno = 11,
                  rm.h = FALSE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }, 
  prmtop = 'data/doprstslt01.prmtop')

pdb.unrst.h.long <- pdb.unrst.h.long.files[[1]]
pdb.unrst.h.long$xyz <- do.call(rbind, lapply(pdb.unrst.h.long.files, function(x) x$xyz))

unrst.divider <- length(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_dopunrstslt01_250(\\_restart\\d*)?.nc$")
)

pdb.unrst.h.long$xyz <- pdb.unrst.h.long$xyz[1:(5000*25/20*unrst.divider),] #trim the first 25000 frames (exactly one microsecond)

# clean up pdb
pdb.unrst.h.long$atom <- pdb.unrst.h.long$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

#align all frames
pdb.unrst.h.long$xyz <- fit.xyz(
  fixed = pdb.unrst.h.long$xyz[1,],
  mobile = pdb.unrst.h.long$xyz,
  fixed.inds = atom.select(pdb.unrst.h.long, elety = unique(pdb.unrst.h.long$atom$elety))$xyz,
  mobile.inds = atom.select(pdb.unrst.h.long, elety = unique(pdb.unrst.h.long$atom$elety))$xyz
)

```

```{r}
#| label: noldpMDdata
#| output: false

pdb.noldp.long.files <- lapply(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_noligand_250(\\_restart\\d*)?.nc$", 
    full.names = TRUE), 
  function(prmtop = 'D:/Amber/dop/leap/noligand.prmtop', crd) {
    as.pdb(
      prmtop = read.prmtop(prmtop),
      crd = read.ncdf(crd)
    ) %>%
      convert.pdb(., renumber = TRUE,
                  first.resno = 11,
                  rm.h = TRUE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }, 
  prmtop = 'D:/Amber/dop/leap/noligand.prmtop')

pdb.noldp.long <- pdb.noldp.long.files[[1]]
pdb.noldp.long$xyz <- do.call(rbind, lapply(pdb.noldp.long.files, function(x) x$xyz))

noldp.divider <- length(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_noligand_250(\\_restart\\d*)?.nc$")
)

pdb.noldp.long$xyz <- pdb.noldp.long$xyz[1:(5000*25/20*noldp.divider),] #trim the first 25000 frames (exactly one microsecond)

# clean up pdb
pdb.noldp.long$atom <- pdb.noldp.long$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

#align all frames
pdb.noldp.long$xyz <- fit.xyz(
  fixed = pdb.noldp.long$xyz[1,],
  mobile = pdb.noldp.long$xyz,
  fixed.inds = atom.select(pdb.noldp.long, elety = unique(pdb.noldp.long$atom$elety))$xyz,
  mobile.inds = atom.select(pdb.noldp.long, elety = unique(pdb.noldp.long$atom$elety))$xyz
)

```

```{r}
#| label: noldpMDh
#| output: false

pdb.noldp.h.long.files <- lapply(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_noligand_250(\\_restart\\d*)?.nc$", 
    full.names = TRUE), 
  function(prmtop = 'D:/Amber/dop/leap/noligand.prmtop', crd) {
    as.pdb(
      prmtop = read.prmtop(prmtop),
      crd = read.ncdf(crd)
    ) %>%
      convert.pdb(., renumber = TRUE,
                  first.resno = 11,
                  rm.h = FALSE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }, 
  prmtop = 'D:/Amber/dop/leap/noligand.prmtop')

pdb.noldp.h.long <- pdb.noldp.h.long.files[[1]]
pdb.noldp.h.long$xyz <- do.call(rbind, lapply(pdb.noldp.h.long.files, function(x) x$xyz))

noldp.divider <- length(
  list.files(
    path = 'D:/Amber/dop/pmemd/out/', 
    pattern = "4_prod_noligand_250(\\_restart\\d*)?.nc$")
)

pdb.noldp.long$xyz <- pdb.noldp.long$xyz[1:(5000*25/20*noldp.divider),] #trim the first 25000 frames (exactly one microsecond)

# clean up pdb
pdb.noldp.h.long$atom <- pdb.noldp.h.long$atom %>% 
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DT3' ~ 'DT',
      resid == 'DT5' ~ 'DT',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

#align all frames
pdb.noldp.h.long$xyz <- fit.xyz(
  fixed = pdb.noldp.h.long$xyz[1,],
  mobile = pdb.noldp.h.long$xyz,
  fixed.inds = atom.select(pdb.noldp.h.long, elety = unique(pdb.noldp.h.long$atom$elety))$xyz,
  mobile.inds = atom.select(pdb.noldp.h.long, elety = unique(pdb.noldp.h.long$atom$elety))$xyz
)

```

```{r}
#| label: doplong
#| output: false

source('R/doplong.R')
```

```{r}
#| label: doplongminstructures
#| output: false

doplong7.finalmin <- read.pdb('D:/Amber/dop/pmemd/out/5_min_doplong7.pdb') %>%
  convert.pdb(
    .,
    rm.h = FALSE,
    rm.wat = TRUE
  ) %>%
  trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))

doplong7.finalmin$atom <- doplong7.finalmin$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    )
  )

doplong8.finalmin <- read.pdb('D:/Amber/dop/pmemd/out/5_min_doplong8.pdb') %>%
  convert.pdb(
    .,
    rm.h = FALSE,
    rm.wat = TRUE
  ) %>%
  trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))

doplong8.finalmin$atom <- doplong8.finalmin$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    )
  )

### Export----
if (!file.exists('data/doplong/5_min_doplong7.pdb')) {
  write.pdb(doplong7.finalmin, file = 'data/doplong/5_min_doplong7.pdb')
  write.pdb(doplong8.finalmin, file = 'data/doplong/5_min_doplong8.pdb')
} else {print('file exists')}


```

```{r}
#| label: doplongshorttraj
#| output: false

doplong7.traj <- as.pdb(
  prm = read.prmtop('D:/Amber/dop/leap/doplong7.prmtop'),
  crd = read.ncdf('D:/Amber/dop/pmemd/out/4_prod_doplong7.nc'),
) %>%
  convert.pdb(
    .,
    rm.h = FALSE,
    rm.wat = TRUE
  ) %>%
  trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))

doplong7.traj$xyz <- fit.xyz(
  fixed = doplong7.traj$xyz[1,],
  mobile = doplong7.traj$xyz,
  fixed.inds = atom.select(doplong7.traj, elety = unique(doplong7.traj$atom$elety))$xyz,
  mobile.inds = atom.select(doplong7.traj, elety = unique(doplong7.traj$atom$elety))$xyz
)

doplong7.traj$atom <- doplong7.traj$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    )
  )

doplong8.traj <- as.pdb(
  prm = read.prmtop('D:/Amber/dop/leap/doplong8.prmtop'),
  crd = read.ncdf('D:/Amber/dop/pmemd/out/4_prod_doplong8.nc'),
) %>%
  convert.pdb(
    .,
    rm.h = FALSE,
    rm.wat = TRUE
  ) %>%
  trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))

doplong8.traj$xyz <- fit.xyz(
  fixed = doplong8.traj$xyz[1,],
  mobile = doplong8.traj$xyz,
  fixed.inds = atom.select(doplong8.traj, elety = unique(doplong8.traj$atom$elety))$xyz,
  mobile.inds = atom.select(doplong8.traj, elety = unique(doplong8.traj$atom$elety))$xyz
)

doplong8.traj$atom <- doplong8.traj$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    )
  )

## Export----

doplong7.traj.1000f <- doplong7.traj
doplong8.traj.1000f <- doplong8.traj

doplong.traj.1000f <- lapply(
  list(doplong7.traj.1000f, doplong8.traj.1000f),
  function(x) {
    x$xyz <- x$xyz[seq(1, 5000, by = 5),]
    return(x)
  }
)

names(doplong.traj.1000f) <- c('doplong7', 'doplong8')

if (!file.exists('data/doplong/4_prod_doplong7_1000f.pdb')) {
  lapply(
    seq_along(doplong.traj.1000f),
    function(x) {
      write.pdb(
        pdb = doplong.traj.1000f[[x]],
        file = paste0('data/doplong/4_prod_', names(doplong.traj.1000f[x]), '_1000f.pdb')
      )
    }
  )
} else {print('file exists')}


```

```{r}
#| label: doplonglongtraj
#| output: false

doplong7_longtrajfiles <- lapply(
  paste0('D:/Amber/dop/pmemd/out/4_prod_doplong7_250_', 1:4, '.nc'),
  function(x) {
    as.pdb(
      prm = read.prmtop('D:/Amber/dop/leap/doplong7.prmtop'),
      crd = read.ncdf(x)
    ) %>%
      convert.pdb(.,
                  rm.h = TRUE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }
)

doplong7_longtraj <- doplong7_longtrajfiles[[1]]

# rbind all $xyz items
doplong7_longtraj$xyz <- do.call(rbind, lapply(doplong7_longtrajfiles, function(x) x$xyz))

doplong7_longtraj$xyz <- fit.xyz(
  fixed = doplong7_longtraj$xyz[1,],
  mobile = doplong7_longtraj$xyz,
  fixed.inds = atom.select(doplong7_longtraj, elety = unique(doplong7_longtraj$atom$elety))$xyz,
  mobile.inds = atom.select(doplong7_longtraj, elety = unique(doplong7_longtraj$atom$elety))$xyz
)

doplong7_longtraj$atom <- doplong7_longtraj$atom %>%
    mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
      elesy = if_else(
        is.na(elesy) | elesy == 'not found',
        case_when(
          elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
          elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
          elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
          elety %in% c(paste0("N", 1:99)) ~ 'N',
          elety %in% c("P", paste0("P", 1:99)) ~ 'P',
          TRUE ~ 'not found'
        ),
        elesy
    )
  )

doplong8_longtrajfiles <- lapply(
  paste0('D:/Amber/dop/pmemd/out/4_prod_doplong8_250_', 1:4, '.nc'),
  function(x) {
    as.pdb(
      prm = read.prmtop('D:/Amber/dop/leap/doplong8.prmtop'),
      crd = read.ncdf(x)
    ) %>%
      convert.pdb(.,
                  rm.h = TRUE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }
)

doplong8_longtraj <- doplong8_longtrajfiles[[1]]

# rbind all $xyz items
doplong8_longtraj$xyz <- do.call(rbind, lapply(doplong8_longtrajfiles, function(x) x$xyz))

doplong8_longtraj$xyz <- fit.xyz(
  fixed = doplong8_longtraj$xyz[1,],
  mobile = doplong8_longtraj$xyz,
  fixed.inds = atom.select(doplong8_longtraj, elety = unique(doplong8_longtraj$atom$elety))$xyz,
  mobile.inds = atom.select(doplong8_longtraj, elety = unique(doplong8_longtraj$atom$elety))$xyz
)

doplong8_longtraj$atom <- doplong8_longtraj$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
      elesy = if_else(
        is.na(elesy) | elesy == 'not found',
        case_when(
          elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
          elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
          elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
          elety %in% c(paste0("N", 1:99)) ~ 'N',
          elety %in% c("P", paste0("P", 1:99)) ~ 'P',
          TRUE ~ 'not found'
        ),
        elesy
    )
  )

gc()
```

```{r}
#| label: doplonglongtrajh
#| output: false

doplong7h_longtrajfiles <- lapply(
  paste0('D:/Amber/dop/pmemd/out/4_prod_doplong7_250_', 1:4, '.nc'),
  function(x) {
    as.pdb(
      prm = read.prmtop('D:/Amber/dop/leap/doplong7.prmtop'),
      crd = read.ncdf(x)
    ) %>%
      convert.pdb(.,
                  rm.h = FALSE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }
)

doplong7h_longtraj <- doplong7h_longtrajfiles[[1]]

# rbind all $xyz items
doplong7h_longtraj$xyz <- do.call(rbind, lapply(doplong7h_longtrajfiles, function(x) x$xyz))

doplong7h_longtraj$xyz <- fit.xyz(
  fixed = doplong7h_longtraj$xyz[1,],
  mobile = doplong7h_longtraj$xyz,
  fixed.inds = atom.select(doplong7h_longtraj, elety = unique(doplong7h_longtraj$atom$elety))$xyz,
  mobile.inds = atom.select(doplong7h_longtraj, elety = unique(doplong7h_longtraj$atom$elety))$xyz
)

doplong7h_longtraj$atom <- doplong7h_longtraj$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

doplong8h_longtrajfiles <- lapply(
  paste0('D:/Amber/dop/pmemd/out/4_prod_doplong8_250_', 1:4, '.nc'),
  function(x) {
    as.pdb(
      prm = read.prmtop('D:/Amber/dop/leap/doplong8.prmtop'),
      crd = read.ncdf(x)
    ) %>%
      convert.pdb(.,
                  rm.h = FALSE,
                  rm.wat = TRUE) %>%
      trim.pdb(., atom.select(., elety = c('Na+', 'Cl-'), inverse = TRUE))
  }
)

doplong8h_longtraj <- doplong8h_longtrajfiles[[1]]

# rbind all $xyz items
doplong8h_longtraj$xyz <- do.call(rbind, lapply(doplong8h_longtrajfiles, function(x) x$xyz))

doplong8h_longtraj$xyz <- fit.xyz(
  fixed = doplong8h_longtraj$xyz[1,],
  mobile = doplong8h_longtraj$xyz,
  fixed.inds = atom.select(doplong8h_longtraj, elety = unique(doplong8h_longtraj$atom$elety))$xyz,
  mobile.inds = atom.select(doplong8h_longtraj, elety = unique(doplong8h_longtraj$atom$elety))$xyz
)

doplong8h_longtraj$atom <- doplong8h_longtraj$atom %>%
  mutate(
    # clean up residue names
    resid = case_when(
      resid == 'DG3' ~ 'DG',
      resid == 'DC5' ~ 'DC',
      TRUE ~ resid
    ),
    # move dopamine to chain B
    chain = if_else(resid == 'LDP', 'B', 'A'),
    # rename dopamine atoms
    elety = case_when(
      resid == 'LDP' & elety == 'OAH' ~ 'O1',
      resid == 'LDP' & elety == 'HAH' ~ 'HO1',
      resid == 'LDP' & elety == 'OAI' ~ 'O2',
      resid == 'LDP' & elety == 'HAI' ~ 'HO2',
      resid == 'LDP' & elety == 'CAB' ~ 'C1',
      resid == 'LDP' & elety == 'CAD' ~ 'C2',
      resid == 'LDP' & elety == 'HAD' ~ 'H2',
      resid == 'LDP' & elety == 'CAG' ~ 'C3',
      resid == 'LDP' & elety == 'CAC' ~ 'C4',
      resid == 'LDP' & elety == 'CAF' ~ 'C5',
      resid == 'LDP' & elety == 'HAF' ~ 'H5',
      resid == 'LDP' & elety == 'CAE' ~ 'C6',
      resid == 'LDP' & elety == 'HAE' ~ 'H6',
      resid == 'LDP' & elety == 'CAA' ~ 'C7',
      resid == 'LDP' & elety == 'HAA' ~ 'H71',
      resid == 'LDP' & elety == 'HAB' ~ 'H72',
      resid == 'LDP' & elety == 'CAJ' ~ 'C8',
      resid == 'LDP' & elety == 'HAJ' ~ 'H81',
      resid == 'LDP' & elety == 'HAK' ~ 'H82',
      resid == 'LDP' & elety == 'NAK' ~ 'N1',
      resid == 'LDP' & elety == 'HAL' ~ 'HN11',
      resid == 'LDP' & elety == 'HAM' ~ 'HN13',
      resid == 'LDP' & elety == 'HAN' ~ 'HN12',
      TRUE ~ elety
    ),
    elesy = if_else(
      is.na(elesy) | elesy == 'not found',
      case_when(
        elety %in% c(paste0("C", 1:99), paste0("C", 1:99, "'")) ~ 'C',
        elety %in% c(paste0("H", 1:99), paste0("H", 1:99, "'"), 'HN11', 'HN12', 'HN13', 'HO1', 'HO2') ~ 'H',
        elety %in% c(paste0("O", 1:99), paste0("O", 1:99, "'"), paste0("OP", 1:99)) ~ 'O',
        elety %in% c(paste0("N", 1:99)) ~ 'N',
        elety %in% c("P", paste0("P", 1:99)) ~ 'P',
        TRUE ~ 'not found'
      ),
      elesy
    )
  )

gc()
```

```{r}
#| label: trajlongexport

if (!file.exists('images/longrMD.1000frames.pdb')) {
  write.pdb(
    pdb = pdb.long,
    xyz = pdb.long$xyz[seq(1, nrow(pdb.long$xyz), by = nrow(pdb.long$xyz)/1000),],
    file = 'images/longrMD.1000frames.pdb'
  )
} 

if (!file.exists('images/longuMD.1000frames.pdb')) {
  write.pdb(
    pdb = pdb.unrst.long,
    xyz = pdb.unrst.long$xyz[seq(1, nrow(pdb.unrst.long$xyz), by = nrow(pdb.unrst.long$xyz)/1000),],
    file = 'images/longuMD.1000frames.pdb'
  )
} 


if (!file.exists('images/longnoLDPMD.1000frames.pdb')) {
  write.pdb(
    pdb = pdb.noldp.long,
    xyz = pdb.noldp.long$xyz[seq(1, nrow(pdb.noldp.long$xyz), by = nrow(pdb.noldp.long$xyz)/1000),],
    file = 'images/longnoLDPMD.1000frames.pdb'
  )
} 

if (!file.exists('images/longdoplong7.1000frames.pdb')) {
  write.pdb(
    pdb = doplong7_longtraj,
    xyz = doplong7_longtraj$xyz[seq(1, nrow(doplong7_longtraj$xyz), by = nrow(doplong7_longtraj$xyz)/1000),],
    file = 'images/longdoplong7.1000frames.pdb'
  )
} 

if (!file.exists('images/longdoplong8.1000frames.pdb')) {
  write.pdb(
    pdb = doplong8_longtraj,
    xyz = doplong8_longtraj$xyz[seq(1, nrow(doplong8_longtraj$xyz), by = nrow(doplong8_longtraj$xyz)/1000),],
    file = 'images/longdoplong8.1000frames.pdb'
  )
} 
```

```{css}
.justify {
text-align: justify !important
}
```

## Methods

### System preparation

::: justify
The system was prepared using the `leap` program from the Amber24 suite [@case2023]. The OL21 and AMBER General Force Field for organic molecules (Version 1.81) were used to describe the complex [@Wang2004; @Zgarbová2021; @Love2023]. The structure was explicitly solvated in a truncated octahedral box of water molecules, using the OPC model with the ad hoc Li/Merz ion parameters of atomic ions (12-6 set)[@Izadi2014; @Li2020], with a minimum of 14 Å between the solute and the box edge. For simulations in absence of salts (case of unrestrained simulations), the system was neutralized by adding 25 Na^+^ counter-cations. In cases where 140 mM NaCl were added to better reflect the SELEX conditions (restrained simulations), the number of required ions ($N_{\pm}$; here around 35 Na<sup>+</sup> and 10 Cl<sup>-</sup>) was determined following the SLTCAP method [@Schmit2018], using @eq-reqions simplified by Machado and Pantano into @eq-requionssimp, where $\nu_w$ is the water volume of the simulation box (around $3\times 10^5$ Å<sup>3</sup> here) in reduced units, $c_0$ the salt concentration, $Q$ the total charge of the complex (here -25: 26 phosphates on the aptamer and 1 ammonium group on the dopamine), and $N_0 = \frac{N_w \times c_0}{55.5}$, with $N_w$. the number of water molecules in the simulation box (here, around $9\times10^4$) [@Machado2020].

$$N_{\pm}=\nu_w c_0 e^{\mp ArcSinh(\frac{Q}{2 \nu_w c_0})}$$ {#eq-reqions} $$N_{\pm}=N_0\sqrt{1+(\frac{Q}{2N_0})^2 \mp \frac{Q}{2}}$$ {#eq-requionssimp}

Note that the SPLIT method describe by Machado and Pantano cannot be applied as our system does not satisfy the $N_0 \gg Q$ condition; however it yields identical values in most cases or deviate by a single ion.

### Restrained minimization, heating and molecular dynamics

Using an in-house R script [@R2023; @tidyverse], the distance restraints obtained in ARIA were converted to an 8-column format suitable for its processing by the `makeDIST_RST` function from the `sander` module of Amber. A custom map file was prepared to define common names for groups of protons sharing a given restraints (e.g. 3 H from a same methyl group). The resulting DISANG restraints file was then applied to all steps below.

All simulation steps were performed with `pmemd.cuda` (v. 18.0.0) from the CUDA version of AMBER [@Götz2012; @Salomon-Ferrer2013; @legrand2013], on an NVIDIA H100 PCIe Tensor core GPU (CUDA version: 12.4) from the DOREMI CALI v3 cluster of the *Mésocentre de Calcul Intensif Aquitain* (Université de Bordeaux). The system was minimized for 20000 cycles using the steepest descent algorithm for the first 4000 steps and the conjugate gradient for the next 16000 steps. The weights of the restraints were kept constant at 100 kcal.mol^-1^.Å^-2^. The system was then heated at constant volume from 0 to 298 K over 18 ps then kept for 2 ps at the final temperature, using a time step of 2 fs, the Langevin thermostat with a 2.0 ps^-1^ collision frequency and a different seed for the pseudo-random number generation for every run to avoid synchronization artifacts [@Sindhikara2009], an 8 Å non-bonded cutoff, and the bonds involving hydrogen were constrained with the SHAKE algorithm. The system was further equilibrated five times at 298 K with the parameters above and the restraint weights ramping down from 100 to 5 kcal.mol^-1^.Å^-2^. The pressure was kept at 1 bar with the Berendsen barostat [@berendsen1984].

Restrained molecular dynamics for subsequent minimization were ran for 10 ns with the parameters above, and restraints weights set at 20 kcal.mol^-1^.Å^-2^. The final coordinates were further minimized, as described above except for the restraint weights set at 20 kcal.mol^-1^.Å^-2^. Remaining NMR violations were summarized with the `sviol` function from the Amber package.

Production MD were run over a microsecond with and without restraints, as well as in absence of dopamine, with the parameters above.

### Data analysis

The `bio3d` package was used for minimized structure and trajectory files cleanup, alignment, filtering, averaging and analysis [@bio3d]. The determination of RMSD and RMSF was performed with the `rmsd` and `rmsf` functions. Dihedral angles, sugar pucker angles and amplitudes $\theta_M$ were obtained with in-house R functions leveraging the `torsion.xyz` function, following @eq-amp, where the pucker P is determined by @eq-tanp and the sugar torsion angles $\nu_i$ are defined by four atoms as shown below.

$$
\theta_M = \frac{\nu_2}{\cos P}
$$ {#eq-amp}

$$
\tan P = \frac{(\nu_4 + \nu_1)-(\nu_3 + \nu_0)}{2\nu_2(\sin(\frac{\pi}{5}) + \sin(\frac{2\pi}{5}))}
$$ {#eq-tanp}

$$\nu_0=C4'-O4'-C1'-C2'$$ $$\nu_1=O4'-C1'-C2'-C3'$$ $$\nu_2=C1'-C2'-C3'-C4'$$ $$\nu_3=C2'-C3'-C4'-O4'$$ $$\nu_4=C3'-C4'-O4'-C1'$$

All atom-atom and ring-ring distances and angles were measured with the `dist.xyz` and `angle.xyz` functions. In-house R scripts were used to infer the formation of H-bonds from these values on full trajectories, by first selecting realistic donor/acceptor atoms for each residue/ligand, then verifying that the donor/acceptor distances and angles were compatible with H-bond formation.

Base pairs were further characterized with DSSR 2.3.2 [@Lu2015; @Lu2020], using the web API through the `httr2` R package [@httr2].

Principal component analysis on the final structures and trajectories (all atoms except hydrogen) was performed with the `pca` function from `bio3d`. The results were clustered with the k-means method using Euclidean distance and the best number of clusters determined by `Nbclust` [@NbClust]. Cross-correlation analysis was carried out with the `dccm` function of `bio3d`.

Molecular structures images were created in PyMOL 2.5 [@PyMOL], interactive molecular structures with mol\*, and ligand binding site diagram with LigPlot<sup>+</sup> [@Laskowski2011]. All further data processing and all plotting was performed in R 4.3. Apache Arrow was used to write/read processed data files in feather format [@arrow].
:::

## Results

::: justify
Hereafter, are described the results for the restrained simulations over 10 ns for the fifteen starting structures, the first of which was also simulated for a microsecond. As a control, the latter was also performed in absernce of restraints.
:::

### Trajectory

#### Visualization

::: justify
##### Microsecond simulation

Below can be seen a comparison between the restrained and unrestrained microsecond simulations. One frame every 10 ns was extracted.

::: {layout-ncol="2"}
{{< mol-url images/longrMD.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/longuMD.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}
:::

##### Production simulation for minimization

All twenty restrained molecular dynamic trajectories (10 ns) can be visualized below. Only one every 5 frames were extracted. Click on the top-left play button to animate the trajectory.
:::

::: {#traj20 layout-nrow="5" layout-ncol="3"}
{{< mol-url images/4_prod_doprstslt01.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt02.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt03.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt04.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt05.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt06.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt07.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt08.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt09.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt10.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt11.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt12.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt13.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt14.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

{{< mol-url images/4_prod_doprstslt15.traj.1000frames.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}
:::

#### Average structures

```{r}
#| label: AvStructR
#| output: false

# split pdb.long$xyz into 20 matrices of n.frames/20 rows
mean.xyz.list <- split(
  as.data.table(pdb.long$xyz), 
  rep(1:20, each = nrow(pdb.long$xyz) %/% 20)
)
# average each coordinate (per atom/dimension) for each matrix
mean.xyz <- lapply(mean.xyz.list, function(x) as.numeric(colMeans(x)))
# Collapse all numerics into a single matrix of 20 sets of coordinates
mean.xyz <- matrix(unlist(mean.xyz), nrow = 20, byrow = TRUE)

if (!file.exists('images/mean_long_rMD.pdb')) {
  write.pdb(pdb = pdb.long, xyz = mean.xyz, file = 'images/mean_long_rMD.pdb')
}

###
```

::: justify
Twenty structures of the microsecond restrained simulation were averaged over `r nrow(pdb.long$xyz) %/% 20` frames of the simulation, without further minimization. It is not incredibly useful as it was not further minimized, but it can be done. The minimized structures can be found in @sec-minimized-structures.
:::

::: {layout-nrow="1"}
{{< mol-url images/mean_long_rMD.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}
:::

#### RMSD

##### Restrained simulations

::: justify
The trajectory of almost all fifteen simulations seem to converge after the first hundreds of ps (@fig-rmsd20). Only structure 2 requires around 5 ns to converge. Minimization of the coordinates after ten nanoseconds of simulation is appropriate.
:::

```{r}
#| label: fig-rmsd20
#| fig-width: 10
#| fig-height: 6
#| fig-cap: "RMSD on all residues (including 5' and 3'-ends) without hydrogens across the restrained simulations (10 ns), where the first frame is used as reference."


#apply rmsd on all objects of xyz.list
rmsd.list <- lapply(1:15, function(i) {
  rmsd(
    a = xyz.list[[i]][1,],
    b = xyz.list[[i]][seq(1, nrow(xyz.list[[i]]), by = 5),],
    fit = FALSE #already aligned at the import stage
  )
})


# bind rows of rmsd.list to make a data.frame, and add a column for the frame number, and a column for the list number
rmsd.df <- do.call(rbind, lapply(1:15, function(i) {
  data.frame(
    frame = seq(1, nrow(xyz.list[[i]]), by = 5),
    rmsd = rmsd.list[[i]],
    structure = i
  ) 
})) %>% 
  mutate(t = frame * info.short$time.per.frame)


#plot rmsds
rmsd.df %>% 
  ggplot(aes(x = t, y = rmsd)) +
  geom_line() +
  custom.theme(0.8) +
  facet_wrap(~structure,
             ncol = 5,
             scales = 'free_x') +
  # scale_x_continuous(n.breaks = 3)  +
  labs(
    x = 't (ns)',
    y = 'RMSD (&angst;)',
  )

```

::: justify
The microsecond simulation has a very stable RMSD vs. the first frame.

```{r}
#| label: fig-rmsdlong
#| fig-width: 6
#| fig-height: 4
#| fig-cap: RMSD on all residues (including 5' and 3'-ends) without hydrogens across the restrained simulation (1 microsecond), where the first frame is used as reference.


data.frame(
  rmsd = rmsd(
    a = pdb.long$xyz[1,],
    b = pdb.long$xyz[seq(1, nrow(pdb.long$xyz), by = divider), ],
    fit = FALSE
  ),
  frame = (seq(1, nrow(pdb.long$xyz), by = divider))
) %>% 
  mutate(t = frame * info.long$time.per.frame) %>% 
  ggplot(aes(x = t, y = rmsd)) +
  geom_line() +
  custom.theme(1) +
  labs(
    x = 't (ns)',
    y = 'RMSD (&angst;)',
  )

```

Two frames with similar one-dimensional RMSD are not necessarily similar. Pairwise RMSD better captures the extent of structural diversity across the simulation and avoids the bias from comparing all frames to a single (first) conformer. A few clusters of structures are visible in the RMSD matrices, suggesting that the simulations have converged to a few distinct states (@fig-pairwise20). However, there are no very significant variations along the trajectories (RMSD never exceeds \~ 2 Å). The nature of the visited states is investigated for the structure 1 in the sections below.
:::

```{r}
#| label: pairwisermsdr
#| output: false
#| cache-lazy: false

#calculate only for 1000 frames
rmsd.pairwise <- lapply(1:15, function(i) {
  rmsd(
    a = xyz.list[[i]][seq(1, nrow(xyz.list[[i]]), by = 5),], #keep one in every 5 frames
    fit = FALSE
  ) %>%
    lazy_dt() %>%
    mutate(var1 = 1:nrow(.), .before = 1) %>%
    pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
    mutate(var2 = as.numeric(gsub("V", "", var2)),
           structure = i) %>%
    as.data.table()
}) %>%
  rbindlist()

arrow::write_feather(rmsd.pairwise, "data/rmsd.pairwise.feather")

```

```{r}
#| label: pairwisermsdrlong
#| output: false
#| cache-lazy: false
#calculate only for 1000 frames

rmsd.pairwise.long <- rmsd(
  a = pdb.long$xyz[seq(1, nrow(pdb.long$xyz), by = nrow(pdb.long$xyz)/1000),], #keep 1000 frames only
  fit = FALSE
) %>%
  lazy_dt() %>%
  mutate(var1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
  mutate(var2 = as.numeric(gsub("V", "", var2))) %>%
  as.data.table()

arrow::write_feather(rmsd.pairwise.long, "data/rmsd.pairwise.long.feather")

###
```

```{r}
#| label: fig-pairwise20
#| fig-cap: All-to-all RMSD (all residues and ligand atoms except H) for the fifteen simulations
#| fig-width: 12
#| fig-height: 9

scaling <- 0.75

arrow::read_feather('data/rmsd.pairwise.feather') %>%
  mutate(t1 = var1 * info.short$time.per.frame * 5,
         t2 = var2 * info.short$time.per.frame * 5) %>%
  ggplot(aes(x = t1, y = t2, fill = value)) +
  geom_raster() +
  facet_wrap(~structure,
             ncol = 5,
             scales = 'free') +
  scale_fill_viridis_c(name = 'RMSD (&angst;)') +
  scale_x_continuous(expand = c(0, 0), limits = c(0, NA), 
                     breaks = c(0, 5, 10),
                     name = 't (ns)') +
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA), 
                     breaks = c(0, 5, 10),
                     name = 't (ns)') +
  ggthemes::theme_pander() +
  theme(
    axis.text.x = element_markdown(size = 12 * scaling),
    axis.text.y = element_markdown(size = 12 * scaling),
    axis.title.x = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.title.y = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.line = element_line(linewidth = 0.75 * scaling),
    axis.ticks = element_line(linewidth = 0.75 * scaling),
    legend.text = element_markdown(size = 12 * scaling),
    legend.title = element_markdown(size = 14 * scaling, face = 'bold'),
    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 = element_blank(),
    panel.grid.minor = 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)
  ) +
  theme(legend.position = 'bottom') + 
  theme(plot.margin = margin(1, 1, 0, 0, 'lines'))

```

```{r}
#| label: fig-pairwiselong
#| fig-cap: All-to-all RMSD (all residues and ligand atoms except H) for the microsecond simulation
#| fig-width: 6
#| fig-height: 7

scaling <- 1.0

arrow::read_feather('data/rmsd.pairwise.long.feather') %>%
  mutate(t1 = var1 * info.long$time.per.frame * nrow(pdb.long$xyz)/1000,
         t2 = var2 * info.long$time.per.frame * nrow(pdb.long$xyz)/1000) %>%
  ggplot(aes(x = t1, y = t2, fill = value)) +
  geom_raster() +
  scale_fill_viridis_c(name = 'RMSD (&angst;)') +
  scale_x_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  ggthemes::theme_pander() +
  theme(
    axis.text.x = element_markdown(size = 12 * scaling),
    axis.text.y = element_markdown(size = 12 * scaling),
    axis.title.x = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.title.y = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.line = element_line(linewidth = 0.75 * scaling),
    axis.ticks = element_line(linewidth = 0.75 * scaling, color = 'black'),
    legend.text = element_markdown(size = 12 * scaling),
    legend.title = element_markdown(size = 14 * scaling, face = 'bold'),
    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 = element_blank(),
    panel.grid.minor = 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)
  ) +
  theme(legend.position = 'bottom') +
  theme(plot.margin = margin(1, 1, 0, 0, 'lines'))

###
```

##### Unrestrained MD

The trajectory is relatively stable across the unrestrained microsecond simulation, although with a larger RMSD (@fig-rmsd, @fig-pairwisermsd).

```{r}
#| label: pairwiseumd
#| output: false

rmsd.pairwise.unrst.long <- rmsd(
  a = pdb.unrst.long$xyz[seq(1, nrow(pdb.unrst.long$xyz), by = nrow(pdb.unrst.long$xyz)/1000),], #keep 1000 frames only
  fit = FALSE
) %>%
  lazy_dt() %>%
  mutate(var1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
  mutate(var2 = as.numeric(gsub("V", "", var2))) %>%
  as.data.table()

arrow::write_feather(rmsd.pairwise.unrst.long, "data/uMD.rmsd.pairwise.feather")

```

```{r}
#| label: fig-rmsd 
#| fig-width: 6
#| fig-height: 4 
#| fig-cap: "RMSD on all residues (including 5' and 3'-ends) without hydrogens across the unrestrained simulation, where the first frame is used as reference." 
#| cache: false   

data.frame(
  rmsd = rmsd(
    a = pdb.unrst.long$xyz[1,],
    b = pdb.unrst.long$xyz[seq(1, nrow(pdb.unrst.long$xyz), by = unrst.divider), ],
    fit = FALSE
  ),
  frame = (seq(1, nrow(pdb.unrst.long$xyz), by = unrst.divider))
) %>% 
  mutate(t = frame * info.long$time.per.frame) %>% 
  ggplot(aes(x = t, y = rmsd)) +
  geom_line() +
  custom.theme(1) +
  labs(
    x = 't (ns)',
    y = 'RMSD (&angst;)',
  )


```

```{r}
#| label: fig-pairwisermsd
#| fig-width: 6
#| fig-height: 6.5
#| fig-cap: All-to-all RMSD (all residues and ligand atoms except H) of the unrestrained simulation
#| 
scaling <- 1.0

arrow::read_feather('data/uMD.rmsd.pairwise.feather') %>%
  mutate(t1 = var1 * info.long.unrst$time.per.frame * nrow(pdb.unrst.long$xyz)/1000,
         t2 = var2 * info.long.unrst$time.per.frame * nrow(pdb.unrst.long$xyz)/1000) %>%
  ggplot(aes(x = t1, y = t2, fill = value)) +
  geom_raster() +
  scale_fill_viridis_c(name = 'RMSD (&angst;)') +
  scale_x_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  ggthemes::theme_pander() +
  theme(
    axis.text.x = element_markdown(size = 12 * scaling),
    axis.text.y = element_markdown(size = 12 * scaling),
    axis.title.x = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.title.y = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.line = element_line(linewidth = 0.75 * scaling),
    axis.ticks = element_line(linewidth = 0.75 * scaling, color = 'black'),
    #legend.text = element_markdown(size = 12 * scaling),
    legend.title = element_markdown(size = 14 * scaling, face = 'bold'),
    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 = element_blank(),
    panel.grid.minor = 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)
  ) +
  theme(legend.position = 'bottom') +
  theme(plot.margin = margin(1, 1, 0, 0, 'lines'))

##
```

#### PCA

::: justify
Given the RMSD heatmaps shown above, it is likely that the aptamer visits different conformations. To isolate discrete states from the microsecond simulation, principal component analysis was performed.
:::

```{r}
#| label: pcatraj
#| output: false
#| cache-lazy: false

#PCA
pca_pdb_traj <- pca.pdbs(
  pdb.long, 
  use.svd = FALSE, 
  rm.gaps = TRUE, 
  fit = FALSE #pdb models are already aligned at the import stage
)

scree <- data.frame(
  pc = 1:length(pca_pdb_traj$L),
  L = pca_pdb_traj$L
) %>% 
  mutate(
    var = L/sum(L) * 100,
    cum.var = cumsum(var)
  ) %>% 
  filter(pc <= 15) %>%
  select(-L) %>%
  mutate(
    label = ifelse(
      pc %in% 1:3 | 
        pc %in% seq(4, nrow(.), 3) | 
        pc == nrow(.), 
      signif(cum.var, 3), NA)
  ) %>% 
  ggplot(., aes(x = pc, y = var)) +
  geom_text_repel(aes(label = label), size = 5, fontface = 'bold', force = 100) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 2) +
  custom.theme(scaling) +
  labs(
    x = 'Number of Principal Components',
    y = 'Proportion of variance (%)'
  ) 

###
```

```{r}
#| label: nbclusttraj
#| output: false
#| cache-lazy: false

pca.nb.traj <- NbClust(
  data.frame(pca_pdb_traj$z) %>% select(all_of(1:4)),
  distance = "euclidean",
  min.nc = 2, max.nc = 5,
  method = "kmeans"
)

###
```

```{r}
#| label: pcaplotrtrajfunc
#| output: false

pca.plotr.traj <- function(pca_object, dim.1, dim.2, 
                           nb = pca.nb.traj,
                           scale = scaling){
  
  scores <- pca_object$z %>% 
    as_tibble() %>% 
    set_names(paste0('PC', 1:ncol(.))) %>%
    mutate(State = 1:n(), .before = 1) %>% 
    mutate(cluster = nb$Best.partition) %>% 
    select(State, cluster, !!sym(dim.1), !!sym(dim.2))
  
  # print(scores)
  
  # Function to calculate convex hull coordinates for each cluster
  get_hull_coordinates <- function(cluster_data) {
    hull_indices <- chull(cluster_data$PC1, cluster_data$PC2)
    hull_coordinates <- cluster_data[hull_indices, ]
    return(hull_coordinates)
  }
  
  # Calculate hull coordinates for each cluster
  renamed_scores <- scores %>% magrittr::set_colnames(c('State', 'cluster', 'PC1', 'PC2'))
  hull_coordinates <- do.call(rbind, lapply(split(renamed_scores, renamed_scores$cluster), get_hull_coordinates)) %>% 
    magrittr::set_colnames(c('State', 'cluster', dim.1, dim.2))
  
  
  # Calculate the variance explained by each PC
  var <- signif(pca_object$L/sum(pca_object$L) * 100, 2)[c(as.numeric(gsub('PC', '', dim.1)), as.numeric(gsub('PC', '', dim.2)))]
  
  representative_states <- renamed_scores %>%
    left_join(
      renamed_scores %>%
        group_by(cluster) %>%
        summarize(
          centroid_PC1 = mean(PC1),
          centroid_PC2 = mean(PC2)
        ), 
      by = "cluster") %>%
    mutate(
      distance_to_centroid = sqrt((PC1 - centroid_PC1)^2 + (PC2 - centroid_PC2)^2)
    ) %>%
    group_by(cluster) %>%
    slice_min(distance_to_centroid) %>% 
    pull(State)
  
  # Plot the scores and color by cluster assignment
  pca_plot <- hull_coordinates %>% 
    ggplot(., aes_string(x = dim.1, y = dim.2)) +
    geom_polygon(
      aes(group = cluster, fill = factor(cluster), color = factor(cluster)),
      alpha = 0.05,
      linewidth = 0.75
    ) +
    geom_point(
      mapping = aes(color = factor(cluster)),
      data = scores,
      alpha = 0.2
    ) +
    geom_text_repel(
      data = scores %>% 
        filter(State %in% representative_states),
      aes(label = State),
      size = 5,
      fontface = 'bold',
      force = 100
    ) +
    geom_point(
      data = scores %>% 
        filter(State %in% representative_states),
      aes(fill = factor(cluster)),
      size = 3, shape = 21, color = 'black', linewidth = 0.75
    ) +
    custom.theme(scale) +
    labs(
      x = glue::glue(dim.1, ' (', var[1], '%)'),
      y =  glue::glue(dim.2, ' (', var[2], '%)')
    ) +
    scale_color_d3(name = 'Cluster') +
    scale_fill_d3(name = 'Cluster')
  
  return(pca_plot)
  
}
```

```{r}
#| label: pcaloadings
#| output: false
#| cache-lazy: false

p.loads.traj <- data.frame(
  u = c(pca_pdb_traj$U[,1], pca_pdb_traj$U[,2], pca_pdb_traj$U[,3]),
  eleno = rep(pdb.long$atom$eleno, 3),
  PC = c(rep('PC1', length(pca_pdb_traj$U[,1])),
         rep('PC2', length(pca_pdb_traj$U[,2])),
         rep('PC3', length(pca_pdb_traj$U[,3])))
) %>% 
  left_join(
    pdb.long$atom,
    by = c("eleno" = "eleno")
  ) %>% 
  group_by(resno, resid, PC, chain) %>%
  summarise(u = sum(abs(u))) %>% 
  mutate(resno = if_else(resid == 'LDP', 38, resno),
         resid = gsub('^D', '', resid),
         label = if_else(
           resid == 'LDP', 
           #make LDP bold and pink
           paste0('<span style="color:pink"><b>', resid, '</b></span>'),
           paste0('<b>', resid, resno, '</b>')
         )
  ) %>%
  filter(PC %in% paste0('PC', 1:2)) %>% 
  ggplot(aes(
    x = factor(
      resno,
      levels = unique(resno),
      labels = unique(label)
    ),
    y = u,
    fill = chain
  )
  ) +
  facet_wrap(~PC, ncol = 1) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c('grey', 'pink')) +
  scale_y_continuous(expand = c(0, 0)) +
  custom.theme(1) +
  theme(axis.title.x = element_blank(),
        legend.position = 'none') +
  labs(y = 'sum of absolute loadings')

###
```

```{r}
#| label: pcafigprep
#| output: false
#| cache-lazy: false

p.pca <- scree +
  pca.plotr.traj(pca_pdb_traj, 'PC1', 'PC2', nb = pca.nb.traj) +
  p.loads.traj +
  plot_layout(
    design = '
              AB
              CC
              ') &
  plot_annotation(tag_levels = c('A', 'B', 'C', 'D'))

###
```

```{r}
#| label: fig-pcatraj
#| fig-cap: 'Principal component analysis on the coordinates of the restrained microsecond trajectory. Scree plot (A) showing the contribution of each principal component on the total variance and the cumulative variance labelled on selected data points. (B) Score plots along the first two principal components, colored by kmeans clusters. (C) Sum of absolute loadings of residues for the first two first principal components'
#| fig-width: 12
#| fig-height: 8

p.pca

###
```

```{r}
#| label: representativestate
#| output: false

#if (!file.exists('data/representative_states_pca_microsec.pdb')) {
write.pdb(
  pdb = pdb.long,
  xyz = pdb.long$xyz[c(10062, 8099, 12207),],
  file = paste0('data/representative_states_pca_microsec.pdb')
)
# }

###

```

::: justify
From the first four principal components, the data is best clustered in three groups. The states closest to the centroids of each cluster are shown below. The largest differences are at the 3'-end region (G36:T37 and G13:A14, different in cluster 3), the T23 chain reversal loop (different in cluster 1), and - to a lower extent - the G25. This correlates well with the root mean squared flucturations measured during the simulation (@fig-rmsfclust)

![Representative structures of the two PCA clusters. Differences are highlighted.](images/cluster%20centroids.png){fig-align="center" width="800"}
:::

```{r}
#| label: fig-rmsfclust
#| fig-cap: Root mean squared fluctuations (RMSF), highlighting the residues with the largest mobility during the restrained simulation.
#| fig-width: 12
#| fig-height: 6

rmsf <- rmsf(pdb.long$xyz)

labels <- pdb.long$atom %>% 
  select(resid, resno, eleno) %>%
  group_by(resno) %>%
  filter(eleno == round(mean(eleno), 0)) %>% 
  mutate(residue = if_else(
    resno < 38,
    paste0(gsub('D', '', resid), resno),
    resid
  ))

rmsf.df <- tibble(
  rmsf = rmsf
) %>% 
  mutate(eleno = 1:nrow(.)) %>%
  left_join(pdb.long$atom, by = c('eleno' = 'eleno'))

rmsf.df %>% 
  filter(resno < 38) %>% 
  ggplot(aes(x = eleno, y = rmsf, fill = factor(resno), color = factor(resno))) +
  geom_vline(
    data = pdb.long$atom %>%
      select(resno, eleno) %>%
      filter(resno < 38) %>% 
      group_by(resno) %>%
      filter(eleno == max(eleno)),
    aes(xintercept = eleno),
    color = 'grey', linetype = 'dotted'
  ) +
  geom_text(
    data = labels %>% filter(resno < 38),
    aes(label = residue, y = 1.05*max(rmsf)),
    fontface = 'bold',
    show.legend = FALSE
  ) +
  geom_text(
    data = labels %>% filter(resno == 38),
    aes(label = residue, y = 1.05*max(rmsf)),
    fontface = 'bold',
    show.legend = FALSE,
    color = 'pink'
  ) +
  geom_col(
    show.legend = FALSE,
    width = 1.2,
    alpha = 0.9
  ) +
  geom_col(
    data = rmsf.df %>% filter(resno == 38),
    color = 'pink', fill = 'pink',
    show.legend = FALSE,
    width = 1.2,
    alpha = 0.9
  ) +
  scale_x_continuous(expand = c(0, 0), limits = c(0,NA)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0,NA)) +
  scale_color_viridis_d(option = 'D') +
  scale_fill_viridis_d(option = 'D') +
  labs(
    x = 'Atom number',
    y = 'RMSF (Å)'
  ) +
  custom.theme(1) +
  coord_cartesian(clip = 'off')

###
```

#### Cross-correlation analysis

The correlation of atomic movements was assessed by calculation of the atom-wise cross-correlation matrix, shown in @fig-dccm as a heatmap. Negative values indicate that the atoms/residues move in opposite directions.

Here, it is possible to observe a positive correlation of the dopamine binding site residues T18 and T19 H-bonded to the cathecol, A27 and G28 to the ammonium, and the stacked C26, T32 and G33 (see @sec-dopamine-binding below).

```{r}
#| label: fig-dccm
#| fig-cap: Dynamical cross-correlation map (DCCM) of the microsecond simulation across atoms (A) or residues (B)
#| fig-width: 12
#| fig-height: 7.5

cij <- dccm(pdb.long$xyz)

cij.atoms <- cij %>% 
  lazy_dt() %>% 
  mutate(atom1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -atom1, names_to = "atom2", values_to = "value") %>% 
  mutate(atom2 = as.numeric(gsub("V", "", atom2))) %>% 
  left_join(pdb.long$atom %>% 
              select(eleno, resno, resid), 
            by = c("atom1" = "eleno")) %>%
  rename(resno1 = resno, resid1 = resid) %>%
  left_join(pdb.long$atom %>% 
              select(eleno, resno, resid), 
            by = c("atom2" = "eleno")) %>%
  rename(resno2 = resno, resid2 = resid) %>% 
  filter(resid1 != 'LDP', resid2 != 'LDP') %>%
  as.data.table()



p.cij.atoms <- ggplot(cij.atoms, aes(x = atom1, y = atom2, fill = value)) +
  geom_raster() +
  scale_fill_gradientn(
    colors = hcl.colors(n = 1E5, palette = "Tropic"),
    limits = c(-1, 1),
    name = ''
  ) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  custom.theme(1) +
  coord_cartesian(clip = "off") +
  labs(x = "Atom", y = "Atom")


cij.res <- cij.atoms %>%
  group_by(resno1, resno2) %>%
  summarise(value = mean(value)) 

p.cij.res <- ggplot(cij.res, aes(x = resno1, y = resno2, fill = value)) +
  geom_raster(show.legend = FALSE) +
  scale_fill_gradientn(
    colors = hcl.colors(n = 1E5, palette = "Tropic"),
    limits = c(-1, 1)
  ) +
  scale_x_continuous(expand = c(0,0), breaks = 11:37) +
  scale_y_continuous(expand = c(0,0), breaks = 11:37) +
  custom.theme(1) +
  coord_cartesian(clip = "off") +
  labs(x = "Residue", y = "Residue")



p.cij.atoms + p.cij.res +
  plot_layout(guides = "collect", ncol = 2) +
  plot_annotation(
    tag_levels = c('A', 'B')
  ) &
  theme(legend.position = 'bottom')

###

```

#### Domain analysis

The identification of geometrically rigid regions of the aptamer in the microsecond simulation was performed with GeoStaS algorithm [@Romanowska2012]. In @fig-gs, we see that most of the aptamer is rigid, with the exception of the T23 loop and, to a lower extent, the nearby T30 loop.

```{r}
#| label: gs
#| output: false
#| cache-lazy: false

gs <- geostas(
  pdb.long$xyz[seq(1, nrow(pdb.long$xyz), by = nrow(pdb.long$xyz)/1000),], 
  fit = FALSE
)

###

```

```{r}
#| label: fig-gs
#| fig-cap: Atomic movement similarity matrix obtained with the GeoStaS analysis for the microsecond simulation across atoms (A) or residues (B). The calculation was carried out on 1000 evenly-spaced frames.
#| fig-width: 12
#| fig-height: 7.5
#| cache-lazy: false

gs.amsm <- gs$amsm %>%
  lazy_dt() %>% 
  mutate(atom1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -atom1, names_to = "atom2", values_to = "value") %>% 
  mutate(atom2 = as.numeric(gsub("V", "", atom2))) %>% 
  left_join(pdb.long$atom %>% 
              select(eleno, resno, resid), 
            by = c("atom1" = "eleno")) %>%
  rename(resno1 = resno, resid1 = resid) %>%
  left_join(pdb.long$atom %>% 
              select(eleno, resno, resid), 
            by = c("atom2" = "eleno")) %>%
  rename(resno2 = resno, resid2 = resid) %>% 
  filter(resid1 != 'LDP', resid2 != 'LDP') %>%
  as.data.table()

p.gs.amsm <- ggplot(gs.amsm, aes(x = atom1, y = atom2, fill = value)) +
  geom_raster() +
  scale_fill_viridis_c(direction = -1, name = '') +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  custom.theme(1) +
  coord_cartesian(clip = "off") +
  labs(x = "Atom", y = "Atom")


gs.amsm.res <- gs.amsm %>%
  group_by(resno1, resno2) %>%
  summarise(value = mean(value)) 

p.gs.amsm.res <- ggplot(gs.amsm.res, aes(x = resno1, y = resno2, fill = value)) +
  geom_raster(show.legend = FALSE) +
  scale_fill_viridis_c(direction = -1, name = '') +
  scale_x_continuous(expand = c(0,0), breaks = 11:37) +
  scale_y_continuous(expand = c(0,0), breaks = 11:37) +
  custom.theme(1) +
  coord_cartesian(clip = "off") +
  labs(x = "Residue", y = "Residue") 

p.gs.amsm + p.gs.amsm.res +
  plot_layout(guides = "collect", ncol = 2) +
  plot_annotation(
    tag_levels = c('A', 'B')
  ) &
  theme(legend.position = 'bottom')

###
```

### Minimized structures {#sec-minimized-structures}

```{r}
#| label: rmsd15
#| cache: false

rmsd.interstates <- rmsd(pdb, 
                         a.inds = atom.select(pdb),
                         b.inds = atom.select(pdb),
                         fit = TRUE) %>% 
  as.data.frame() %>%
  magrittr::set_colnames(as.character(1:15)) %>%
  mutate(State = 1:15, .before = 1) %>% 
  #replace all 0 with NA in columns 2:16
  mutate(across(2:16, ~ifelse(. == 0, NA_real_, .))) %>%
  rbind(c(NA_real_, sapply(.[, -1], mean, na.rm = TRUE))) %>% 
  mutate(State = ifelse(is.na(State), 'Mean', as.character(State)))

# Mean RMSD
mean.rmsd.doprstslt <- rmsd.interstates[rmsd.interstates$State == "Mean", 2:16] %>% 
  #collapse into a single vector
  unlist() %>% 
  mean(.) %>% 
  round(1)

# Std. dev. of RMSD
sd.rmsd.doprstslt <- rmsd.interstates[rmsd.interstates$State == "Mean", 2:16] %>% 
  #collapse into a single vector
  unlist() %>% 
  sd(.) %>% 
  round(2)
```

#### Structures and energies

::: justify
The fifteen minimized structures of the dopamine aptamer obtained by restrained molecular dynamics are displayed in @fig-20static (mean RMSD across structures: `r mean.rmsd.doprstslt` $\pm$ `r sd.rmsd.doprstslt` Å). The longer sequence (with duplex stem at the bottom) are shown alongside (also minimized after 10 ns rMD ; see @sec-long-aptamer).

![Fifteen minimized structures obtained by restrained molecular dynamics and minimization (left) and corresponding structures for the longer sequence calculated with and without A38•T37 base pairing. The dopamine ligand is shown in yellow, guanines in tan, thymines in green, adenines in blue and cytosine in purple.](images/doprstslt_allstates_vs_doplong.png){#fig-20static fig-align="center" width="1200"}

The summary of the simulation energies is given in @tbl-energies.

The inter-states RMSD are given in @tbl-rmsd
:::

{{< mol-url images/multi.min.pdb layoutShowSequence=true layoutShowControls=false viewportShowSelectionMode=true >}}

```{r}
#| label: tbl-energies
#| tbl-cap: Energies
#| cache: false

# list all files in sander/out with the format 5_min_doprstslt**.mdinfo where ** are numbers
files <- list.files('data/', pattern = "5_min_doprstslt[0-9]+.mdinfo", full.names = TRUE)

# function to read the mdinfo file
mdinfo.reader <- function(input.file) {
  
  blob <- fread(input.file, skip = 5, fill = TRUE) %>% 
    as.data.frame()
  
  mdinfo <- data.frame(
    file = basename(input.file),
    bond = blob[1,3],
    angle = blob[1,6],
    dihedral = blob[1,9],
    vdw = blob[2, 3],
    eel = blob[2, 6],
    hbond = blob[2, 9],
    vdw14 = blob[3, 4],
    eel14 = blob[3, 8],
    restraint = blob[3, 11],
    eamber = blob[4, 3],
    bond.rst = blob[5, 5],
    angle.rst = blob[5, 8],
    torsion.rst = blob[5, 11]
  ) %>% 
    mutate(
      structure = gsub('.mdinfo', '', gsub("5_min_doprstslt", "", file)),
      .before = 1
    ) %>% 
    select(structure, eamber, restraint, bond, angle, dihedral, vdw, eel,
           vdw14, eel14) %>% 
    mutate(across(2:ncol(.), as.double))
  
  return(mdinfo)
  
}

# read and bind all files
energies <- do.call(rbind, lapply(files, mdinfo.reader))

digits <- energies %>% 
  pivot_longer(2:ncol(.), names_to = "energy", values_to = "value") %>%
  group_by(energy) %>%
  summarise(rst = sd(value)/mean(value)) %>% 
  mutate(
    digits = case_when(
      abs(rst) < 0.001 ~ 5,
      abs(rst) < 0.01 ~ 4,
      abs(rst) < 0.1 ~ 3,
      abs(rst) < 1 ~ 2,
      TRUE ~ 1
    )
  )

DT::datatable(
  energies %>% 
    filter(structure <= 15),
  filter = 'top', 
  extensions = c('Buttons', 'FixedHeader', 'Scroller'),   
  options = list(     
    dom = 'Bfrtip',     
    buttons = c('copy', 'csv', 'excel', 'pdf'),     
    scrollX = TRUE,     
    scrollY = 450,     
    pageLength = 20   
  ),
  #change column names
  colnames = c(
    'Structure', 'Amber energy', 'NMR restraints', 'Bond', 'Angle', 
    'Dihedral', 'VdW', 'EEL', '1-4 VdW', '1-4 EEL'
  ),
  style = "bootstrap"
) %>%  
  formatSignif(columns = 'eamber', digits = digits %>% filter(energy == 'eamber') %>% pull(digits)) %>%
  formatSignif(columns = 'restraint', digits = digits %>% filter(energy == 'restraint') %>% pull(digits)) %>%
  formatSignif(columns = 'vdw', digits = digits %>% filter(energy == 'vdw') %>% pull(digits)) %>%
  formatSignif(columns = 'eel', digits = digits %>% filter(energy == 'eel') %>% pull(digits)) %>%
  formatSignif(columns = 'bond', digits = digits %>% filter(energy == 'bond') %>% pull(digits)) %>%
  formatSignif(columns = 'angle', digits = digits %>% filter(energy == 'angle') %>% pull(digits)) %>%
  formatSignif(columns = 'dihedral', digits = digits %>% filter(energy == 'dihedral') %>% pull(digits)) %>%
  formatSignif(columns = 'vdw14', digits = digits %>% filter(energy == 'vdw14') %>% pull(digits)) %>%
  formatSignif(columns = 'eel14', digits = digits %>% filter(energy == 'eel14') %>% pull(digits))

####
```

```{r}
#| label: tbl-rmsd
#| tbl-cap: Pairwise RMSD between the top 15 structures
#| cache: false

brks <- quantile(rmsd.interstates[-1], probs = seq(0.005, 0.995, .0025), na.rm = TRUE)

rmsd.interstates %>% 
  DT::datatable(
    extensions = c('Buttons', 'FixedHeader', 'Scroller'),  
    rownames = FALSE,
    options = list(     
      dom = 'Brt',
      buttons = c('copy', 'csv', 'excel', 'prmsd.interstates'),     
      scrollX = TRUE,   
      pageLength = 16   
    ),   
    style = "bootstrap"
  ) %>% 
  formatRound(columns = 2:16, digits = 2) %>%
  formatStyle(
    columns = 1,
    fontWeight = 'bold'
  ) %>% 
  formatStyle(
    columns = 1:16,
    target = 'row',
    fontWeight = styleEqual(c('Mean'), c('bold'))
  ) %>% 
  formatStyle(
    columns = 2:16,
    background = styleInterval(
      quantile(rmsd.interstates[-1], probs = seq(0.005, 0.995, .0025), na.rm = TRUE), 
      c('steelblue', 
        colorRampPalette(c('steelblue', 'white', 'tomato'))(length(brks) - 1),
        'tomato')
    )
  ) %>% 
  formatStyle(
    columns = 1:16,
    textAlign = 'center'
  ) 
###
```

::: justify
#### PCA

It is apparent that some structures deviate more from others, e.g. 1, 4, and 8 and, to a lower extent, 9,12 and 14. We aimed at identifying clusters of structures to summarize the conformational space explored by the aptamer. Below are shown the results of a principal component analysis based on all atoms except hydrogen atoms.
:::

```{r}
#| label: pca
#| echo: false
#| output: false
#| cache: false


pca_pdb <- pca.pdbs(
  pdb, 
  use.svd = FALSE, 
  rm.gaps = TRUE, 
  fit = FALSE #pdb models are already aligned at the import stage
)

scree <- data.frame(
  pc = 1:length(pca_pdb$L),
  L = pca_pdb$L
) %>% 
  mutate(
    var = L/sum(L) * 100,
    cum.var = cumsum(var)
  ) %>% 
  filter(pc <= 15) %>%
  select(-L) %>%
  mutate(
    label = ifelse(
      pc %in% 1:3 | 
        pc %in% seq(4, nrow(.), 3) | 
        pc == nrow(.), 
      signif(cum.var, 3), NA)
  ) %>% 
  ggplot(., aes(x = pc, y = var)) +
  geom_text_repel(aes(label = label), size = 5, fontface = 'bold', force = 100) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 2) +
  custom.theme(scaling) +
  labs(
    x = 'Number of Principal Components',
    y = 'Proportion of variance (%)'
  ) 


pca.plotr <- function(pca_object, dim.1, dim.2, 
                      max.pc = 4, meth = "kmeans", dist = "euclidean", 
                      min_nc = 3, max_nc = 7,
                      scale = scaling){
  
  # Run NbClust to get the best number of clusters and cluster assignment
  nb <- NbClust(
    data.frame(pca_object$z) %>% select(all_of(1:max.pc)), 
    distance = dist, 
    min.nc = min_nc, max.nc = max_nc, 
    method = meth
  )
  
  # Extract and format scores
  scores <- pca_object$z %>% 
    as_tibble() %>% 
    set_names(paste0('PC', 1:ncol(.))) %>%
    mutate(State = 1:n(), .before = 1) %>% 
    mutate(cluster = nb$Best.partition) %>% 
    select(State, cluster, !!sym(dim.1), !!sym(dim.2))
  
  # print(scores)
  
  # Function to calculate convex hull coordinates for each cluster
  get_hull_coordinates <- function(cluster_data) {
    hull_indices <- chull(cluster_data$PC1, cluster_data$PC2)
    hull_coordinates <- cluster_data[hull_indices, ]
    return(hull_coordinates)
  }
  
  # Calculate hull coordinates for each cluster
  renamed_scores <- scores %>% magrittr::set_colnames(c('State', 'cluster', 'PC1', 'PC2'))
  hull_coordinates <- do.call(rbind, lapply(split(renamed_scores, renamed_scores$cluster), get_hull_coordinates)) %>% 
    magrittr::set_colnames(c('State', 'cluster', dim.1, dim.2))
  
  
  # Calculate the variance explained by each PC
  var <- signif(pca_object$L/sum(pca_object$L) * 100, 2)[c(as.numeric(gsub('PC', '', dim.1)), as.numeric(gsub('PC', '', dim.2)))]
  
  # Plot the scores and color by cluster assignment
  pca_plot <- hull_coordinates %>% 
    ggplot(., aes_string(x = dim.1, y = dim.2)) +
    geom_polygon(
      aes(group = cluster, fill = factor(cluster), color = factor(cluster)),
      alpha = 0.2
    ) +
    geom_point(
      mapping = aes(color = factor(cluster)),
      data = scores
    ) +
    geom_text_repel(
      mapping = aes(label = State, color = factor(cluster)),
      data = scores,
      size = 5,
      fontface = 'bold',
      show.legend = FALSE
    ) +
    custom.theme(scale) +
    labs(
      x = glue::glue(dim.1, ' (', var[1], '%)'),
      y =  glue::glue(dim.2, ' (', var[2], '%)')
    ) +
    scale_color_d3(name = 'Cluster') +
    scale_fill_d3(name = 'Cluster')
  
  return(pca_plot)
  
}

p.pca <- scree +
  pca.plotr(pca_pdb, 'PC1', 'PC2') +
  pca.plotr(pca_pdb, 'PC1', 'PC3') +
  pca.plotr(pca_pdb, 'PC2', 'PC3') +
  plot_layout(guides = 'collect') &
  plot_annotation(tag_levels = c('A', 'B', 'C', 'D'))
```

```{r}
#| label: fig-pca
#| cache: false
#| fig-width: 10
#| fig-heigth: 18
#| fig-cap: 'Principal component analysis on the coordinates of the 15 minimized structures. Scree plot (A) showing the contribution of each principal component on the total variance and the cumulative variance labelled on selected data points, and (B) to (D): score plots along two dimensions for the first three principal components. Clusters determined by kmeans are colored.'

#plot
p.pca

```

```{r}
#| label: fig-loadings
#| fig-cap: Sum of absolute loadings of residues for the first three first principal components
#| fig-width: 14

data.frame(
  u = c(pca_pdb$U[,1], pca_pdb$U[,2], pca_pdb$U[,3]),
  eleno = rep(pdb$atom$eleno, 3),
  PC = c(rep('PC1', length(pca_pdb$U[,1])),
         rep('PC2', length(pca_pdb$U[,2])),
         rep('PC3', length(pca_pdb$U[,3])))
) %>% 
  left_join(
    pdb$atom,
    by = c("eleno" = "eleno")
  ) %>% 
  group_by(resno, resid, PC, chain) %>%
  summarise(u = sum(abs(u))) %>% 
  mutate(resno = if_else(resid == 'LDP', 38, resno),
         resid = gsub('^D', '', resid),
         label = if_else(
           resid == 'LDP', 
           #make LDP bold and pink
           paste0('<span style="color:tomato"><b>', resid, '</b></span>'),
           paste0('<b>', resid, resno, '</b>')
         )
  ) %>%
  ggplot(aes(
    x = factor(
      resno,
      levels = unique(resno),
      labels = unique(label)
    ),
    y = u,
    fill = chain
  )
  ) +
  facet_wrap(~PC, ncol = 1) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c('grey', 'pink')) +
  custom.theme(1) +
  theme(axis.title.x = element_blank(),
        legend.position = 'none') +
  labs(y = 'sum of absolute loadings')

###
```

As expected from the analysis of trajectories, the key differences between the clusters are the loop thymine T23, the 3'end, and G25.

![Minimized structures colored by clusters](images/clusters_15minstruct.png){width="707"}

#### NMR violations

```{r}
#| label: violations

rstslt.violations <- rst.error.batchr(input.file = 'data/nmr_viol_doprstslt.txt', error.thresh = 0.15, occurence.thresh = 0, states = 15)

writexl::write_xlsx(rstslt.violations$violations, 'data/rstslt_all_violations.xlsx')
writexl::write_xlsx(rstslt.violations$summary, 'data/rstslt_summary_violations.xlsx')

```

::: justify
The average violations observed across the 20 minimized structures, with a cut-off of 0.15 Å are presented in @tbl-violations. No violation exceeds a mean error of `r rstslt.violations$summary %>% pull(mean.error) %>%  max(.)` Å. Out of the `r rstslt.violations$summary %>% nrow(.)` violations above the cutoff, `r rstslt.violations$summary %>% filter(occurence > 7) %>% nrow(.)` are present in more than half of the structures.

Dopamine is involved in `r rstslt.violations$summary %>% filter(str_detect(atom.1, 'LDP28') | str_detect(atom.2, 'LDP28')) %>% nrow(.)` violation(s).
:::

```{r}
#| label: tbl-violations
#| tbl-cap: "NMR restraint violations (mean error > 0.15 &angst;)"  
#| cache: false

DT::datatable(   
  rstslt.violations$summary %>%     
    filter(mean.error > 0.15) %>% 
    mutate(
      atom.1 = gsub('HAD', 'H2', atom.1)
    ),   
  filter = 'top',
  extensions = c('Buttons', 'FixedHeader', 'Scroller'),   
  options = list(     
    dom = 'Bfrtip',     
    buttons = c('copy', 'csv', 'excel', 'pdf'),     
    scrollX = TRUE,     
    scrollY = 400,     
    pageLength = 10   
  ),   
  #change column names   
  colnames = c(
    'Atom 1', 'Atom 2', 'Target', 'Mean error (&angst;)',
    'SD error (&angst;)', 'Occurence'
  ),
  style = "bootstrap"
) 
```

#### Base pairs

::: justify
The base pairs identified in the 20 minimized structures are named, where possible, in @tbl-dssrbp, and further described using three nomenclatures.

The <u>Saenger</u> nomeclature is a list of 28 (plus one added later) base pairs identified by Roman numerals [@saenger1984c; @gesteland2006rna]. A canonical Watson-Crick base pair is XIX for GC, and XX for AT.

The <u>Leontis-Westhof</u> classification [@Leontis2001], latter expanded by Lemieux and Major [@Lemieux2002], describes base-pairing by considering whether H-bonds are established on the Watson-Crick, Hoogsteen or sugar edges, and what is the glycosidic orientation (*cis*/*trans*: c/t) of the nucleotides. A canonical Watson-Crick base pair would then be cWW. Note that the sugar edge has been defined with RNA in mind, and therefore with 2'-OH that are absent here. This sugar edge is absent for *syn* glycosidic dihedral angles (see section @sec-tetrad below), which is only the case of G31 here.

The <u>DSSR</u> classification is similar to the Leontis-Westhof, as it indicates the *cis*/*trans* (c/t) orientation followed by the interacting edges, which can be Watson-Crick (W), major (M; \~Hoogsteen) and minor grooves (m; \~sugar). It additionally incorporate the +/- base 'face indication as in the Saenger nomenclature. Here, a canonical Watson-Crick base pair is called cW-W.
:::

```{r}
#| label: dssr_request
#| output: false

c(1:15) %>% # initialize structure numbers
  as.character() %>%
  # if the number is less than 10, add a leading zero
  sapply(function(x) {
    if (nchar(x) == 1) {
      paste0('0', x)
    } else {
      x
    }
  }) %>%
  # paste the url
  sapply(function(x) {
    paste0('https://raw.githubusercontent.com/EricLarG4/EricLarG4.github.io/master/indiv_pdb/5_min_doprstslt', x, '.pdb')
  }) %>% 
  # feed to dssr function
  lapply(., dssr)
```

```{r}
#| label: tbl-dssrbp
#| tbl-cap: "Base pairs identified by DSSR in the 15 minimized structures. Abbreviations. +/-: the paired nucleotides have same/opposite faces (as in canonical antiparallel Watson-Crick pairing); WC and W: Watson-Crick pairing/edge; H: Hoogsteen edge; S: sugar edge; M/m: Major/minor groove edges; .: undefined edge; r: reverse; c/t: *cis*/*trans* orientation along the glycosidic bond." 
#| cache: false

dssrbp <- lapply(list.files('data/dssr', full.names = TRUE), dssr.bp) %>% 
  bind_rows() %>% 
  # filter(structure <= 15) %>% 
  #offset the numbers from nt1 and nt2 strings by 10
  mutate(
    nt1 = paste0(str_extract(nt1, '[A-Z]+'), as.numeric(str_extract(nt1, '\\d+')) + 10),
    nt2 = paste0(str_extract(nt2, '[A-Z]+'), as.numeric(str_extract(nt2, '\\d+')) + 10)
  )

DT::datatable(
  dssrbp,
  filter = 'top',
  extensions = c('Buttons', 'FixedHeader', 'Scroller'),
  options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf'),
    scrollX = TRUE,
    scrollY = 450,
    pageLength = 30
  ),
  colnames = c(
    'Structure', 'ID', 'nucleotide 1', 'nucleotide 2', 'base pair', 'Name', 'Saenger', 'Leontis-Westhof', 'DSSR'
  ),
  style = "bootstrap"
)

###
```

::: justify
Remarkably, the whole structure only contains two canonical Watson-Crick base pairs (C24-G29 and C26-G33; @fig-WC). One GC (C20-G28) and two AT ((T12-A35 and T19-A27) base pairs have rWC configurations (@fig-rWC). Note that T19 is bound by the dopamine ligand as well (see @sec-dopamine-binding).
:::

::: {#fig-bp layout-ncol="2"}
![Watson-Crick GC base pairs](images/WC.png){#fig-WC}

![Reverse Watson-Crick AT base pairs](images/rWC.png){#fig-rWC}

Watson-Crick and reverse Watson-Crick base pairs
:::

### Dihedral angles {#sec-tetrad}

::: justify
All residues are in the *anti* configuration except for G31 in [*syn*]{style="color:tomato"} (@tbl-miniangles). G31 forms a trans Watson-Crick base pair with G21 (tW+W; @fig-G2111).
:::

```{r}
#| label: tbl-miniangles
#| tbl-cap: Dihedral &chi; angles of the 15 minimized structures
#| cache: false

chi20 <- chiR(pdb, chain = 'A') %>% 
  rename(structure = frame) 


chi20.sum <- chi20 %>%
  group_by(residue) %>% 
  summarise(
    mean.chi = mean(chi),
    rsd.chi = sd(chi)/mean.chi*100
  )

chi20 %>%
  pivot_wider(
    names_from = structure,
    values_from = chi
  ) %>%
  left_join(chi20.sum, by = 'residue') %>% 
  select(residue, mean.chi, rsd.chi, everything()) %>%
  mutate(gba = ifelse(mean.chi > 60 & mean.chi < 90, 'syn', 'anti'),
         .before = mean.chi) %>%
  DT::datatable(
    extensions = c('Buttons', 'FixedHeader', 'Scroller'),   
    options = list(     
      dom = 'Bfrtip',     
      buttons = c('copy', 'csv', 'excel', 'pdf'),     
      scrollX = TRUE,     
      scrollY = 450,     
      pageLength = 30   
    ),
    #change column names
    colnames = c(
      'Residue', 'GBA', 'Mean', 'RSD (%)', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15'
    ),
    style = "bootstrap"
  ) %>%
  formatSignif(columns = c(3, 5:19), digits = 3) %>%
  formatSignif(columns = 4, digits = 2) %>% 
  formatStyle(
    columns = 2,
    color = styleEqual(c('syn', 'anti'), c('lightcoral', 'black'))
  ) 


```

![Non-canonical base pair between G31 and G21 (tW+W), within a triplet with G29 (tm+m)](images/G21G11.png){#fig-G2111 fig-align="center" width="500"}

.

```{r}
#| label: chirstslt
#| output: false
# cache-lazy: false

chi.rstslt <- chiR(pdb.long, chain = 'A', res.start = 11, res.end = 37)
chi.unrstslt <- chiR(pdb.unrst.long, chain = 'A', res.start = 11, res.end = 37)

###
```

::: justify
Dihedral angles are stable along the simulation (@fig-chirstslt). The same is observed during the unrestrained simulation, albeit with more variance (@fig-chi).
:::

```{r}
#| label: fig-chirstslt
#| fig-width: 9
#| fig-height: 9
#| fig-cap: "Dihedral angles of all residues for the rMD, wrapped in [-180;180], along the trajectory (from center to outside of circle)"

scaling <- 0.7

dihedrals_rstslt <- chi.rstslt %>%
  group_by(residue) %>% 
  mutate(
    chi = ifelse(chi > 180, chi - 360, chi),
    mean.chi = mean(chi),
    mod.chi = if_else(
      mean.chi < 20 & mean.chi > 0 & chi < 0, chi + 360,
      if_else(
        mean.chi < -160 & chi > 0, chi - 360, chi
      )
    ),
    mean.chi = mean(mod.chi),
    gba = case_when(
      mean.chi > 60 & mean.chi < 90 ~ 'syn',
      mean.chi < -60 | mean.chi > -180 ~ 'anti',
      TRUE ~ 'other'
    ),
    residue = str_replace_all(residue, 'D|3_|5_|_', ''),
    #add html colors to the residue as function of their gba
    residue.label = case_when(
      gba == 'syn' ~ paste0('<span style="color:tomato">', residue, '</span>'),
      gba == 'anti' ~ paste0('<span style="color:grey10">', residue, '</span>'),
      TRUE ~ '<span style="color:steelblue4">other</span>'
    ),
    residue.number = as.numeric(str_extract(residue, '\\d+'))
  ) 

ggplot(dihedrals_rstslt, aes(chi, frame, color = gba)) +
  facet_wrap(~factor(
    residue.label,
    levels = dihedrals_rstslt %>%
      arrange(residue.number) %>%
      pull(residue.label) %>%
      unique()
  )) +
  geom_point(alpha = 0.005) +
  scale_y_continuous(expand = c(0,0), breaks = c(99999)) +
  scale_x_continuous(limits = c(-180, 180), breaks = c(-90, 0, 90, 180)) +
  scale_color_manual(values = rev(c('tomato', 'grey10'))) +
  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(),
    axis.ticks = element_blank(),
    strip.text.x = element_markdown(size = 16 * scaling, face = 'bold'),
    strip.text.y = element_markdown(size = 16 * scaling, face = 'bold'),
    strip.background = element_rect(),
    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(),
    panel.spacing.y = unit(1, 'lines'),
    panel.spacing.x = unit(2, 'lines'),
    plot.title = element_markdown(size = 16 * scaling, face = 'bold', hjust = 0.5),
    plot.subtitle = element_markdown(size = 14 * scaling),
    plot.caption = element_markdown(hjust = 1, size = 14 * scaling),
    legend.position = 'none'
  ) +
  coord_polar(clip = 'off')

##

```

```{r}
#| label: fig-chi
#| fig-width: 9
#| fig-height: 9
#| fig-cap: "Dihedral angles of all residues for the unrestrained MD, wrapped in [-180;180], along the trajectory (from center to outside of circle)"


scaling <- 0.7

dihedrals_unrstslt <- chi.unrstslt %>%
  group_by(residue) %>% 
  mutate(
    chi = ifelse(chi > 180, chi - 360, chi),
    mean.chi = mean(chi),
    mod.chi = if_else(
      mean.chi < 20 & mean.chi > 0 & chi < 0, chi + 360,
      if_else(
        mean.chi < -160 & chi > 0, chi - 360, chi
      )
    ),
    mean.chi = mean(mod.chi),
    gba = case_when(
      mean.chi > 60 & mean.chi < 90 ~ 'syn',
      mean.chi < -60 | mean.chi > -180 ~ 'anti',
      TRUE ~ 'other'
    ),
    residue = str_replace_all(residue, 'D|3_|5_|_', ''),
    #add html colors to the residue as function of their gba
    residue.label = case_when(
      gba == 'syn' ~ paste0('<span style="color:tomato">', residue, '</span>'),
      gba == 'anti' ~ paste0('<span style="color:grey10">', residue, '</span>'),
      TRUE ~ '<span style="color:steelblue4">other</span>'
    ),
    residue.number = as.numeric(str_extract(residue, '\\d+'))
  ) 

ggplot(dihedrals_unrstslt, aes(chi, frame, color = gba)) +
  facet_wrap(~factor(
    residue.label,
    levels = dihedrals_unrstslt %>%
      arrange(residue.number) %>%
      pull(residue.label) %>%
      unique()
  )) +
  geom_point(alpha = 0.005) +
  scale_y_continuous(expand = c(0,0), breaks = c(99999)) +
  scale_x_continuous(limits = c(-180, 180), breaks = c(-90, 0, 90, 180)) +
  scale_color_manual(values = rev(c('tomato', 'grey10'))) +
  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(),
    axis.ticks = element_blank(),
    strip.text.x = element_markdown(size = 16 * scaling, face = 'bold'),
    strip.text.y = element_markdown(size = 16 * scaling, face = 'bold'),
    strip.background = element_rect(),
    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(),
    panel.spacing.y = unit(1, 'lines'),
    panel.spacing.x = unit(2, 'lines'),
    plot.title = element_markdown(size = 16 * scaling, face = 'bold', hjust = 0.5),
    plot.subtitle = element_markdown(size = 14 * scaling),
    plot.caption = element_markdown(hjust = 1, size = 14 * scaling),
    legend.position = 'none'
  ) +
  coord_polar(clip = 'off')

```

### Puckers {#sec-puckers}

::: justify
Most nucleotides have their sugar adopting a C2'-endo (typical in B DNA) or close to it (C1'-exo, C3'-exo; @tbl-puckers, @fig-pucker20). Among those, T23 also visits the C1'-endo pucker during a large portion of the simulation (@fig-puckstab), and can be found in several of the minimized structures; however this is of little significance given that this residue does not form any intra- or inter-molecular contacts.

G25, whose amino group binds to the N7 of G28 (tm-M; @fig-G15), is found in half the structures in C1'-exo or C2'-endo but adopts a variety of puckers in others. During the microsecond simulation, it mostly exist in C1'-exo, and only visits the other configurations for a short period of time.

T18, which is stacked on G17, is tightly bound by the dopamine ligand (@fig-T8) and forms an off-plane tW-M base pair with G33, has a C1'-endo configuration in most conformers, which has been observed for dsDNA around intercalator drugs or when gaining increasing flexibility [@wang1987; @shen2011]. It is therefore likely that the backbone of T18 adopts a non-canonical conformation to allow dopamine binding. It remains stable during the simulation.

G16, C20, and G34 mostly exist in the C4'-exo/O4'-endo configurations, which is not very common for DNA [@Anosova2016]. G16 pairs with the Hoogsteen face of G34 (@fig-G6; tW-M) and their puckering is quite unstable during the simulation (@fig-puckstab). C20 is involved in a triplet made of a reverse Watson-Crick base pair with G28 (tW+W) and a single H-bond to T32 (tW-.; @fig-C10).

G17 has an extremely varied pucker depending on the minimized structure (9 unique puckers observed!), and has a remarkably diffuse pucker angle and amplitude during the course of the microsecond simulation. G17 is involved in a non coplanar triplet with the Hoogsteen face of G33 (tW-M; stacked under the dopamine) and T11 (cW-m; @fig-G7). It also binds to G34 (tW-M) and is stacked on G16, whose puckering is also fairly unstable.

A15 is the only nucleotide that mostly adopts a C4'-endo pucker. It forms an asymetric homopurine base pair with A35 (@fig-A5; tW-M) , and the sugar pucker is stable during the simulation (@fig-puckstab).
:::

```{r}
#| label: fig-pucker20
#| fig-cap: Sugar puckering of the 15 minimized structures. Puckers are presented from most to least frequent.
#| fig-width: 10
#| fig-height: 10

scaling <- 0.7 

puckers.rstslt20 <- nuR(pdb, chain = 'A', res.start = 11, res.end = 37) %>% 
  puck.formatR(.) 


puckers.rstslt20 %>% 
  ungroup() %>% 
  mutate(
    puck = case_when(
      angle > 0 & angle < 36 ~ "C3'-endo",
      angle > 36 & angle < 72 ~ "C4'-exo",
      angle > 72 & angle < 108 ~ "O4'-endo",
      angle > 108 & angle < 144 ~ "C1'-exo",
      angle > 144 & angle < 180 ~ "C2'-endo",
      angle > 180 & angle < 216 ~ "C3'-exo",
      angle > 216 & angle < 252 ~ "C4'-endo",
      angle > 252 & angle < 288 ~ "O4'-exo",
      angle > 288 & angle < 324 ~ "C1'-endo",
      angle > 324 & angle < 360 ~ "C2'-exo"
    )
  ) %>% 
  ggplot(aes(x = angle, y = amplitude, fill = puck)) +
  geom_point(size = 3, alpha = 0.5, 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') +
  facet_wrap(~residue)


```

```{r}
#| label: tbl-puckers
#| tbl-cap: Sugar puckering of the 15 minimized structures. Puckers are presented from most to least frequent.
#| cache: false


puckers.rstslt20 <- nuR(pdb, chain = 'A', res.start = 11, res.end = 37) %>% 
  puck.formatR(.) %>% 
  group_by(residue) %>% 
  mutate(
    mean.angle = mean(angle),
    rsd.angle = sd(angle)/mean.angle*100
  ) %>% 
  ungroup() %>% 
  mutate(
    puck = case_when(
      angle > 0 & angle < 36 ~ "C3'-endo",
      angle > 36 & angle < 72 ~ "C4'-exo",
      angle > 72 & angle < 108 ~ "O4'-endo",
      angle > 108 & angle < 144 ~ "C1'-exo",
      angle > 144 & angle < 180 ~ "C2'-endo",
      angle > 180 & angle < 216 ~ "C3'-exo",
      angle > 216 & angle < 252 ~ "C4'-endo",
      angle > 252 & angle < 288 ~ "O4'-exo",
      angle > 288 & angle < 324 ~ "C1'-endo",
      angle > 324 & angle < 360 ~ "C2'-exo"
    )
  )


puckers.rstslt20 <- puckers.rstslt20 %>% 
  select(-c(puck, nu0, nu1, nu2, nu3, nu4, P0, amplitude, mean.amplitude, mod.angle)) %>% 
  left_join(
    puckers.rstslt20 %>% 
      group_by(residue) %>%
      count(puck) %>%
      arrange(-n) %>% 
      mutate(if_else(n == max(n), paste0('**', puck, '**'), puck)) %>%
      summarise(puck = paste(unique(puck), collapse = ", ")), 
    by = 'residue'
  ) 


puckers.rstslt20 %>%
  pivot_wider(
    names_from = frame,
    values_from = angle
  ) %>%  
  select(residue, puck, mean.angle, rsd.angle, everything()) %>%
  DT::datatable(
    extensions = c('Buttons', 'FixedHeader', 'Scroller'),   
    options = list(     
      dom = 'Bfrtip',     
      buttons = c('copy', 'csv', 'excel', 'pdf'),     
      scrollX = TRUE,     
      scrollY = 450,     
      pageLength = 30   
    ),
    #change column names
    colnames = c(
      'Residue', 'Pucker', 'Mean', 'RSD (%)', 
      '1', '2', '3', '4', '5', '6', '7', '8', 
      '9', '10', '11', '12', '13', '14', '15'
    ),
    style = "bootstrap"
  ) %>%
  formatSignif(columns = c(3, 5:19), digits = 3) %>%
  formatSignif(columns = 4, digits = 2)



```

```{r}
#| label: fig-puckrstslt
#| fig-width: 10
#| fig-height: 10
#| fig-cap: "Sugar puckering calculated with the Altona & Sundaralingam method, wrapped from 0 to 360 for the restrained simulation. The position of the points relative to the center/border of the circles scales with the pucker amplitude (the larger the closer to the border). The circled points represent the median pucker/amplitude. Calculation performed on 5000 evenly spaced frames."

scaling <- 0.8

pdb.long.puck <- pdb.long
pdb.long.puck$xyz <- pdb.long$xyz[seq(1, nrow(pdb.long$xyz), by = nrow(pdb.long$xyz)/5000),]

puckers.rstslt <- nuR(
  pdb.long.puck, chain = 'A', res.start = 11, res.end = 37
) %>% 
  puck.formatR(.)

puck.plotR(puckers.rstslt, alpha = 0.01)




```

```{r}
#| label: fig-puckstab
#| fig-width: 10
#| fig-height: 10
#| fig-cap: "Pucker angles of all residues for the rMD, along the trajectory (from center to outside of circle). Residues with relative standard deviation above 20% are labeled as flexible"

scaling <- 0.8

puckers.rstslt %>% 
  group_by(residue) %>%
  mutate(
    pucker.color = if_else(
      sd(angle)/mean(angle) > 0.20,
      'flexible', 'stable'
    ),
    pucker.color = if_else(residue == 'C10', 'stable', pucker.color),
    pucker.color = factor(pucker.color, levels = c('stable', 'flexible'))
  ) %>% 
  ggplot(aes(y = frame, x = angle, color = pucker.color)) +
  geom_point(alpha = 0.025) +
  geom_point(
    data = . %>% 
      filter(frame == 1)
  ) +
  facet_wrap(~factor(
    residue, levels = puckers.rstslt %>%
      arrange(mean.angle) %>%
      pull(residue) %>%
      unique())) +
  scale_x_continuous(limits = c(0,360), breaks = 36*(0:10)) +
  scale_y_continuous(expand = c(0,0)) +
  scale_color_manual(values = c('grey', 'pink'), name = 'pucker stability') +
  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') 


```

::: {#fig-binding layout-ncol="2" layout-nrow="3"}
![Sheared base pair between the sugar edge of G25 (left) and the Hoogsteen edge of G28 (tm-M)](images/G15.png){#fig-G15}

![Binding of T18 to dopamine, G33 (tW-M), and stacking on G17](images/T8.png){#fig-T8 fig-align="center"}

![Watson-Crick to Hoogsteen base pairing of G16 and G34 (tW-M)](images/G6G24.png){#fig-G6}

![Reverse Watson-Crick base pair between C20 and G28 (tW+W) within a triplet with T32](images/C10G18.png){#fig-C10}

![Triplet of G17 (foreground, top right), T11 (cW-m; background) and G33 (tW-M, foreground, top left) under the dopamine binding site, and above G34 (bottom left; tW-M with G16 and G17, tW+W with T11).](images/G7.png){#fig-G7}

![tW-M base pairing of A15 (right) with A35 (left), and to a lower extent G36 (bottom)](images/A5.png){#fig-A5}

Non-canonical base pairing leading to 'weird' puckers
:::

::: justify
In the unrestrained simulations, almost all nucleotides have their sugar adopting a C2'-endo (typical in B DNA) or close to it (C1'-exo; @fig-puck). Only T18 adopts a C4'-exo configuration. The non-canonical puckers observed above are therefore imposed by the restraints. For the specific case of T18, the more RNA-like pucker observed in the unrestrained simulation still differs from the C1'-endo pucker of its restrained counterpart.
:::

```{r}
#| label: fig-puck
#| fig-width: 10
#| fig-height: 10
#| fig-cap: "Sugar puckering calculated with the Altona & Sundaralingam method, wrapped from 0 to 360 for the unrestrained simulation. The position of the points relative to the center/border of the circles scales with the pucker amplitude (the larger the closer to the border). The circled points represent the median pucker/amplitude."

scaling <- 0.8

pdb.unrst.long.puck <- pdb.unrst.long
pdb.unrst.long.puck$xyz <- pdb.unrst.long$xyz[seq(1, nrow(pdb.unrst.long$xyz), by = nrow(pdb.unrst.long$xyz)/5000),]

puckers.rstslt <- nuR(
  pdb.unrst.long.puck, chain = 'A', res.start = 11, res.end = 37
) %>% 
  puck.formatR(.)

puck.plotR(puckers.rstslt, alpha = 0.01)


###
```

### Dopamine binding {#sec-dopamine-binding}

::: justify
The dopamine is intercalated between T32 and G33/C16, and forms H-bond contacts with the O4s of T18 and T19 with its two hydroxy groups (<i>O1</i>, <i>O2</i>; @fig-dopbindrst). The ammonium group (<i>N1</i> forms two additional H-bonds, with the N3 of of A27 (and O4'; although not always canonically since the bond angle is sometimes below 120°; @fig-anglecanon) and the deoxyribose O4' of G28.

The H-bonds between O1/O2 of the dopamine and T18/T19 are globally stable during the restrained simulation (@fig-dophbonds), in particular the short LDP:O2/T18:O4 (@fig-hbondfreqallexp). The ammonium group (N1) also establishes stable H-bonds with G28 and A27.

In absence of restraints, the H-bond with A27 and G28 are longer and less stable. Conversely, the bond to T19 is shorter and stable. H-bond angles are also significantly much less stable.

```{r}
#| label: contacts
#| output: false

pair.list.rst <- contactr(
 min.pdb,
 dist.thresh = 3.1, freq.thresh = 0.2,
 lgd.resno = 38,
 lgd.donors = c('O2', 'O1', 'N1'),
 lgd.acceptors = NA,
 donor.atoms = c('N', 'O', 'S', 'F'),
 acceptor.atoms = c('O', 'N', 'S', 'F', 'Cl', 'Br', 'I')
  )

pair.dist.rst <- lapply(
  pair.list.rst$pair,
  atom.distr,
  pdb.input = pdb.long,
  span = 0.01
) %>%
  do.call(rbind, .)

pair.angle.rst <- lapply(
  pair.list.rst$pair,
  bond.angler,
  pdb.input = pdb.long.h, 
  info = info.long, 
  span = 0.01, 
  clean = TRUE, 
  rectify = 50
) %>%
  do.call(rbind, .)

###
```

```{r}
#| label: contactdiscovery
#| output: false

# pdb.noldp.5000 <- pdb.noldp.long
# pdb.noldp.5000$xyz <- pdb.noldp.long$xyz[seq(1, nrow(pdb.noldp.long$xyz), by = nrow(pdb.noldp.long$xyz)/5000),]
# 
# pair.list.noldp <- contactr(pdb.noldp.5000, 3.1, freq.thresh = 0.05)

```

```{r}
#| label: ucontacts
#| output: false


pair.dist.unrst <- lapply(
  pair.list.rst$pair, #recycling of same contacts than rst!
  atom.distr,
  pdb.input = pdb.unrst.long,
  span = 0.01
) %>%
  do.call(rbind, .)

pair.angle.unrst <- lapply(
  pair.list.rst$pair, #recycling of same contacts than rst!
  bond.angler,
  pdb.input = pdb.unrst.h.long, 
  info = info.long, 
  span = 0.01, 
  clean = TRUE, 
  rectify = 50
) %>%
  do.call(rbind, .)

###
```

```{r}
#| label: noldpcontacts
#| output: false

pair.dist.noldp <- lapply(
  #recycling of same contacts than rst!
  #but without LDP
  pair.list.rst$pair[!str_detect(pair.list.rst$pair, "LDP")], 
  atom.distr,
  pdb.input = pdb.noldp.long,
  span = 0.01
) %>%
  do.call(rbind, .)

pair.angle.noldp <- lapply(
  #recycling of same contacts than rst!
  #but without LDP
  pair.list.rst$pair[!str_detect(pair.list.rst$pair, "LDP")], 
  bond.angler,
  pdb.input = pdb.noldp.h.long, 
  info = info.long, 
  span = 0.01, 
  clean = TRUE, 
  rectify = 50
) %>%
  do.call(rbind, .)

###
```

```{r}
#| label: fig-dophbonds
#| fig-cap: 'Lengths of H-bonds (donor-acceptor distance) involved in dopamine binding during the restrained (top) and unrestrained (bottom) simulations. The lines are the result of local polynomial regression fitting with a span of 0.01. The canonical C26•G33 bp was used as a reference.'
#| fig-width: 12
#| fig-height: 10

(
  pair.dist.rst %>%
    filter(
      str_detect(pair, 'LDP') | 
        (str_detect(pair, 'G33') & str_detect(pair, 'C26'))
    ) %>% 
    ggplot(., aes(x = t, y = dist.loess, color = pair)) +
    geom_hline(yintercept = 3.0, linetype = 'dashed') +
    geom_text_repel(
      data = . %>% filter(t == max(t)),
      aes(label = pair),
      xlim = c(1020, 1100),
      fontface = 'bold',
      show.legend = FALSE,
      size = 5*scaling
    ) +
    geom_line(show.legend = FALSE,
              linewidth = 0.75 * scaling) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(limits = c(2.5, 3.5), expand = c(0,0)) +
    scale_color_d3() +
    labs(x = 'Time (ns)', y = 'Distance (Å)') +
    custom.theme(scaling) +
    theme(plot.margin = margin(1, 150, 1, 1)) +
    coord_cartesian(clip = 'off')
) +
  (
    pair.dist.unrst %>%
      filter(
        str_detect(pair, 'LDP') | 
          (str_detect(pair, 'G33') & str_detect(pair, 'C26'))
      ) %>% 
      ggplot(., aes(x = t, y = dist.loess, color = pair)) +
      geom_hline(yintercept = 3.0, linetype = 'dashed') +
      geom_text_repel(
        data = . %>% filter(t == max(t)),
        aes(label = pair),
        xlim = c(1020, 1100),
        fontface = 'bold',
        show.legend = FALSE,
        size = 5*scaling
      ) +
      geom_line(show.legend = FALSE,
                linewidth = 0.75 * scaling) +
      scale_x_continuous(expand = c(0, 0)) +
      scale_y_continuous(limits = c(2.5, 3.5), expand = c(0,0)) +
      scale_color_d3() +
      labs(x = 'Time (ns)', y = 'Distance (Å)') +
      custom.theme(scaling) +
      theme(plot.margin = margin(1, 150, 1, 1)) +
      coord_cartesian(clip = 'off')
  ) +
  plot_layout(ncol = 2)


###
```

```{r}
#| label: fig-anglecanon
#| fig-cap: Angles of H-bonds (donor-H-acceptor atoms) established by the dopamine ligand during the microsecond simulation, with the canonical C26•G33 bp used as a reference. When multiple interconverting hydrogen atoms are involved, all are considered and each data point reflect the angle formed by the hydrogen most likely to establish an H-bond. The lines are the result of local polynomial regression fitting with a span of 0.01. The angle stays below 140 &deg; during the simulation in the case of A27:O4', which is suboptimal to establish an H-bond. Note that the hydrogen atoms interconvert frequently during the simulation; all three atoms were considered and .
#| fig-width: 12
#| fig-height: 7


(
  pair.angle.rst %>% 
    filter(
      str_detect(pair, 'LDP') | 
        (str_detect(pair, 'G33') & str_detect(pair, 'C26'))
    ) %>% 
    ggplot(., aes(x = t, y = angle.loess, color = pair)) +
    geom_hline(yintercept = c(120, 140), linetype = 'dashed', linewidth = 0.75 * scaling) +
    geom_text_repel(
      data = . %>% filter(t == max(t)),
      aes(label = pair),
      xlim = c(1020, 1100),
      fontface = 'bold',
      show.legend = FALSE,
      size = 5*scaling
    ) +
    geom_line(linewidth = 0.75, show.legend = FALSE) +
    labs(x = 't (ns)', y = 'Angle (°)') +
    custom.theme(1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), limits = c(100, 180),
                       breaks = c(100, 120, 140, 160, 180)) +
    scale_color_d3() +
    theme(plot.margin = margin(10, 150, 1, 1)) +
    coord_cartesian(clip = 'off')
) + (
  pair.angle.unrst %>% 
    filter(
      str_detect(pair, 'LDP') | 
        (str_detect(pair, 'G33') & str_detect(pair, 'C26'))
    ) %>% 
    ggplot(., aes(x = t, y = angle.loess, color = pair)) +
    geom_hline(yintercept = 120, linetype = 'dashed', linewidth = 0.75 * scaling) +
    geom_text_repel(
      data = . %>% filter(t == max(t)),
      aes(label = pair),
      xlim = c(1020, 1100),
      fontface = 'bold',
      show.legend = FALSE,
      size = 5*scaling
    ) +
    geom_line(linewidth = 0.75, show.legend = FALSE) +
    labs(x = 't (ns)', y = 'Angle (°)') +
    custom.theme(1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), limits = c(100, 180),
                       breaks = c(100, 120, 140, 160, 180)) +
    scale_color_d3() +
    theme(plot.margin = margin(10, 150, 1, 1)) +
    coord_cartesian(clip = 'off')
)

###
```

T18 is off-plane as it also binds to the Hoogsteen face of G33 stacked 'under' the dopamine. G33 is itself base pairing canonically to C26 (see also @fig-WC).

Conversely, T19 is coplanar to the dopamine, and the two form a triplet with A27. T19 and A27 are arranged in a reverse WC configuration (see also @fig-rWC).

::: {#fig-dopbindrst layout-ncol="2"}
![Binding site of dopamine](images/bindingsiterstslt.png){#fig-dopbindrst3d}

![2D interaction diagram. The ligand is shown in yellow, residues in orange, water in light blue, carbon in black, nitrogen in blue and oxygen in red. Hydrogen bonds are shown as dashed green lines. Hydrophobic interactions are highlighted with red lines.](images/ligplotplus_dopamine_aptamer.png){#fig-ligplot fig-align="center" width="500"}

Dopamine binding in the minimized structures obtained by rMD.
:::

In absence of NMR restraints, we noticed that the ammonium side-chain flips between configurations 'above' or 'below' the catechol plane. We performed PCA on the binding pocket (T18:T19-C26:G28-T32:G33 and dopamine) to identify frames representative of the different configurations.

```{r}
#| label: pcaunrst
#| output: false

# Create a binding pocket trimmed pdb object
# from the long unrestrained simulation

core.pdb <- pdb.unrst.long %>% 
  trim.pdb(
    .,
    atom.select(
      .,
      resno = c(18, 19, 26, 27, 28, 32, 33, 38)
    )
  )

# Re-align frames to the binding pocket only
core.pdb$xyz <- fit.xyz(
  fixed = core.pdb$xyz[1,],
  mobile = core.pdb$xyz,
  fixed.inds = atom.select(core.pdb, elety = unique(core.pdb$atom$elety))$xyz,
  mobile.inds = atom.select(core.pdb, elety = unique(core.pdb$atom$elety))$xyz
)


# Perform PCA
pca_unrst_core <- pca.pdbs(
  core.pdb, 
  use.svd = FALSE, 
  rm.gaps = TRUE, 
  fit = FALSE #pdb models are already aligned at the import stage
)

# Scree plot
scree_unrst_core <- data.frame(
  pc = 1:length(pca_unrst_core$L),
  L = pca_unrst_core$L
) %>% 
  mutate(
    var = L/sum(L) * 100,
    cum.var = cumsum(var)
  ) %>% 
  filter(pc <= 15) %>%
  select(-L) %>%
  mutate(
    label = ifelse(
      pc %in% 1:3 | 
        pc %in% seq(4, nrow(.), 3) | 
        pc == nrow(.), 
      signif(cum.var, 3), NA)
  ) %>% 
  ggplot(., aes(x = pc, y = var)) +
  geom_text_repel(aes(label = label), size = 5, fontface = 'bold', force = 100) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 2) +
  custom.theme(scaling) +
  labs(
    x = 'Number of Principal Components',
    y = 'Proportion of variance (%)'
  ) 

# NbClust determination [LONG!]
pca.nb.unrst <- NbClust(
  data.frame(pca_unrst_core$z) %>% select(all_of(1:4)),
  distance = "euclidean",
  min.nc = 3, max.nc = 5,
  method = "kmeans"
)


# count the number of times each number of clusters was selected as the best
tbl.clusters <- table(pca.nb.unrst$Best.partition)
frac.clust.1 <- round(tbl.clusters[1]/sum(tbl.clusters)*100, 1)
frac.clust.4 <- round(tbl.clusters[4]/sum(tbl.clusters)*100, 1)

# plot loadings
p.loads.unrst <- data.frame(
  u = c(pca_unrst_core$U[,1], pca_unrst_core$U[,2], pca_unrst_core$U[,3]),
  eleno = rep(core.pdb$atom$eleno, 3),
  PC = c(rep('PC1', length(pca_unrst_core$U[,1])),
         rep('PC2', length(pca_unrst_core$U[,2])),
         rep('PC3', length(pca_unrst_core$U[,3])))
) %>% 
  left_join(
    core.pdb$atom,
    by = c("eleno" = "eleno")
  ) %>% 
  group_by(resno, resid, PC, chain) %>%
  summarise(u = sum(abs(u))) %>% 
  mutate(resno = if_else(resid == 'LDP', 38, resno),
         resid = gsub('^D', '', resid),
         label = if_else(
           resid == 'LDP', 
           #make LDP bold and pink
           paste0('<span style="color:pink"><b>', resid, '</b></span>'),
           paste0('<b>', resid, resno, '</b>')
         )
  ) %>%
  filter(PC %in% paste0('PC', 1:2)) %>% 
  ggplot(aes(
    x = factor(
      resno,
      levels = unique(resno),
      labels = unique(label)
    ),
    y = u,
    fill = chain
  )
  ) +
  facet_wrap(~PC, ncol = 1) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c('grey', 'pink')) +
  scale_y_continuous(expand = c(0, 0)) +
  custom.theme(1) +
  theme(axis.title.x = element_blank(),
        legend.position = 'none') +
  labs(y = 'sum of absolute loadings')

#full plot
p.pca.unrst <- scree_unrst_core +
  pca.plotr.traj(pca_unrst_core, 'PC1', 'PC2', nb = pca.nb.unrst) +
  p.loads.unrst +
  plot_layout(
    design = '
              AB
              CC
              ') &
  plot_annotation(tag_levels = c('A', 'B', 'C', 'D'))

#remake a core but with hydrogens for figure below
core.pdb.h <- pdb.unrst.h.long %>% 
  trim.pdb(
    .,
    atom.select(
      .,
      resno = c(18, 19, 26, 27, 28, 32, 33, 38)
    )
  )

core.pdb.h$xyz <- core.pdb.h$xyz[c(22927, 13823, 11460, 18656),]

# Re-align frames to the binding pocket only
core.pdb.h$xyz <- fit.xyz(
  fixed = core.pdb.h$xyz[1,],
  mobile = core.pdb.h$xyz,
  fixed.inds = atom.select(core.pdb.h, elety = unique(core.pdb.h$atom$elety))$xyz,
  mobile.inds = atom.select(core.pdb.h, elety = unique(core.pdb.h$atom$elety))$xyz
)

if (!file.exists('data/core.h.pdb')) {
  write.pdb(core.pdb.h, file = 'data/core.h.pdb')
}
```

```{r}
#| label: fig-pcaunrst
#| fig-cap: 'Principal component analysis on the coordinates of the binding pocket in the unrestrained microsecond trajectory. Scree plot (A) showing the contribution of each principal component on the total variance and the cumulative variance labelled on selected data points. (B) Score plots along the first two principal components, colored by kmeans clusters. (C) Sum of absolute loadings of residues for the first two first principal components'
#| fig-width: 12
#| fig-height: 8


p.pca.unrst
```

Four clusters were identified and the structures closest to their centroids are shown in @fig-bindingpocketunrst. The main difference if the configuration of the ammonium chain. The ammonium binds to either:

i.  The O4' and N3 of G33, which is stacked under the dopamine (top-left structure)

ii. the O2 of T32, which is stacked above the dopamine (top-right structure)

iii. the O4' and N3 of G28 (bottom-left structure)

iv. the O2 of T32 and N3 of G28, in a configuration midway between the two above (bottom-right structure)

The first configuration only accounts for `{r} frac.clust.1`% of the simulation. In the three latter case, T32 overlaps more with the catechol moiety, which increases $\pi-\pi$ interactions. The fourth configuration is the more frequent (`{r} frac.clust.4`% of the frames are in this cluster).

Neither of the four configurations follows exactly the restrained one. The 'above' configurations are closer to it, but the ammonium binds differently in each three case. It is folded back above the chain (akin to a scorpio) and does not bind to A27. However, beyond these four discrete examples, we evidence below that the ammonium is still able to bind to A27:N3 in around of third of the frames (@fig-hbondfreqallexp). For reference, it is almost always bound when restrained (@fig-hbondfreqallexp).

![Binding modes in absence of NMR restraints. Four main binding modes have been evidenced by PCA, mostly varying by the configuration of the amino chain.](images/binding_unrst.png){#fig-bindingpocketunrst fig-align="center"}
:::

### H-bond formation and stability

#### Complex

::: justify
Besides those from the binding site, the formation of all H-bonds was monitored over the course of the simulation. Below, the list of monitored bonds is based on H-bonds detected in the minimized structures obtained by rMD; therefore @fig-hbondfreqallexp ignores by design bonds formed exclusively in the uMD simulation and/or in absence of dopamine. The purpose of @fig-hbondfreqallexp is therefore to establish the stability of H-bonds during the rMD simulation and assess whether they 'survived' in absence of restraints and dopamine.

The tight binding of dopamine to T18:O4 by the catechol moiety is confirmed. It is still present in absence of restraints but less frequently. The ammonium group almost always bind the sugar of G28 and the N3 of A27. Similarly, the frequency of H-bond formation is decreased in absence of restraints.

The less canonical binding to the sugar of the A27 is also observed, with low frequency but across the whole simulation, in the rMD. It is almost absent in absence of restraints, as seen above.

Conversely to the observations above, the H-bonding of dopamine to T19 is less frequent and tight in the rMD than in the uMD, as expected from bond length analysis above (@fig-dophbonds).

```{r}
#| label: hbondfreqr
#| output: false

pair.rst.freq <- bond.freqr(
  dist = pair.dist.rst,
  angl = pair.angle.rst,
  dist.cutoff = 3.0,
  angl.cutoff = 140,
  dist.lax.cutoff = 3.5,
  angl.lax.cutoff = 120
)

pair.unrst.freq <- bond.freqr(
  dist = pair.dist.unrst,
  angl = pair.angle.unrst,
  dist.cutoff = 3.0,
  angl.cutoff = 140,
  dist.lax.cutoff = 3.5,
  angl.lax.cutoff = 120
)

pair.noldp.freq <- bond.freqr(
  dist = pair.dist.noldp,
  angl = pair.angle.noldp,
  dist.cutoff = 3.0,
  angl.cutoff = 140,
  dist.lax.cutoff = 3.5,
  angl.lax.cutoff = 120
)

###
```

```{r}
#| label: fig-hbondfreqallexp
#| fig-height: 20
#| fig-width: 12
#| fig-cap: "Formation of H-bonds over the course of the simulations with restraints (rMD), without restraints (uMD) and unrestrained with no dopamine (no LDP). The list of H-bonds was established on the minimized structures from the rMD simulations, with a maximum donor-acceptor distance of maximum 3.1 &angst;, keeping those formed in at least 20% of structures. H-bonds in the 5'- and 3'-end regions (residues 11:17 and 33:37() are colored in blue, and those involving the dopamine in pink. Top: each line is a frame (one frame every 0.04 ns); the darker a frame, the shorter and more linear the H-bond (the transparency follows linearly the bond angle/length ratio). Bottom: Cumulative frequencies. The frequency of bond formation is shown for relaxed (3.5 &angst; and 120°; light colors) and more stringent (3.0 &angst; and 140°; dark colors) cut-offs. Only inter-residue H-bonds with a 'stringent' frequency above 0.10 are shown."

all.exp <- list(
  bond.index = rbind(
    pair.rst.freq$bond.index %>% mutate(exp = 'rMD'),
    pair.unrst.freq$bond.index %>% mutate(exp = 'uMD'),
    pair.noldp.freq$bond.index %>% mutate(exp = 'no LDP')
  ),
  bond.freq = rbind(
    pair.rst.freq$bond.freq %>% mutate(exp = 'rMD'),
    pair.unrst.freq$bond.freq %>% mutate(exp = 'uMD'),
    pair.noldp.freq$bond.freq %>% mutate(exp = 'no LDP')
  )
)

hbond.cumul.plotr(
  all.exp, freq.cutoff = 0.10,
                  facet = TRUE, 
  main.exp = 'rMD', 
  exp.order = c('rMD', 'uMD', 'no LDP'),
  highlight = 'LDP', 
  highlight.region = c(11:16, 34:37)
  )

###
###
###
```

#### (Pre-folded) free DNA

Generally speaking, the residues in the 5' and 3'-end regions (in blue in @fig-hbondfreqallexp) are remarkably well conserved in absence of dopamine, when starting the simulation from the folded structure. This may indicate that the oligonucleotide can prefold into a short stem intermediate before the binding occurs between this stem and the very dynamic loops, after which the latter could finish folding.

PCA on the trajectory yielded three clusters (@fig-pcanoldp), whose centroid coordinates are shown in @fig-noldppdb. The terminus region can be nicely aligned across clusters, but also against the complex, further demonstrating that the stem is very stable during the simulation, and retains the complex structure.

The loop region also maintains the same local folding in all clusters (it can be locally aligned), however its orientation relative to the terminus stem changes dramatically. It can be orthogonal (cluster 1; @fig-noldppdb; left), yielding a linear configuration ('open'), or parallel, forming a more globular structure ('closed'), closer to that of the complex (cluster 3; @fig-noldppdb; right). The second cluster contain intermediate conformations between these two extrema (cluster 2; @fig-noldppdb; center).

```{r}
#| label: pcanoldp
#| output: false

#PCA
pca_noldp <- pca.pdbs(
  pdb.noldp.long, 
  use.svd = FALSE, 
  rm.gaps = TRUE, 
  fit = FALSE #pdb models are already aligned at the import stage
)

scree <- data.frame(
  pc = 1:length(pca_noldp$L),
  L = pca_noldp$L
) %>% 
  mutate(
    var = L/sum(L) * 100,
    cum.var = cumsum(var)
  ) %>% 
  filter(pc <= 15) %>%
  select(-L) %>%
  mutate(
    label = ifelse(
      pc %in% 1:3 | 
        pc %in% seq(4, nrow(.), 3) | 
        pc == nrow(.), 
      signif(cum.var, 3), NA)
  ) %>% 
  ggplot(., aes(x = pc, y = var)) +
  geom_text_repel(aes(label = label), size = 5, fontface = 'bold', force = 100) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 2) +
  custom.theme(scaling) +
  labs(
    x = 'Number of Principal Components',
    y = 'Proportion of variance (%)'
  ) 

```

```{r}
#| label: pcanoldpnbclust
#| output: false

pca.nb.noldp <- NbClust(
  data.frame(pca_noldp$z) %>% select(all_of(1:4)),
  distance = "euclidean",
  min.nc = 3, max.nc = 5,
  method = "kmeans"
)

```

```{r}
#| label: pcanoldpfigprep
#| output: false
#| cache-lazy: false

p.loads.noldp <- data.frame(
  u = c(pca_noldp$U[,1], pca_noldp$U[,2], pca_noldp$U[,3]),
  eleno = rep(pdb.noldp.long$atom$eleno, 3),
  PC = c(rep('PC1', length(pca_noldp$U[,1])),
         rep('PC2', length(pca_noldp$U[,2])),
         rep('PC3', length(pca_noldp$U[,3])))
) %>% 
  left_join(
    pdb.noldp.long$atom,
    by = c("eleno" = "eleno")
  ) %>% 
  group_by(resno, resid, PC, chain) %>%
  summarise(u = sum(abs(u))) %>% 
  mutate(
    resno = if_else(resid == 'LDP', 38, resno),
    resid = gsub('^D', '', resid),
    label = if_else(
      resid == 'LDP', 
      #make LDP bold and pink
      paste0('<span style="color:pink"><b>', resid, '</b></span>'),
      paste0('<b>', resid, resno, '</b>')
    )
  ) %>%
  filter(PC %in% paste0('PC', 1:2)) %>% 
  ggplot(aes(
    x = factor(
      resno,
      levels = unique(resno),
      labels = unique(label)
    ),
    y = u,
    fill = chain
  )
  ) +
  facet_wrap(~PC, ncol = 1) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c('grey', 'pink')) +
  scale_y_continuous(expand = c(0, 0)) +
  custom.theme(1) +
  theme(axis.title.x = element_blank(),
        legend.position = 'none') +
  labs(y = 'sum of absolute loadings')


p.pca.noldp <- scree +
  pca.plotr.traj(pca_noldp, 'PC1', 'PC2', nb = pca.nb.noldp) +
  p.loads.noldp +
  plot_layout(
    design = '
              AB
              CC
              ') &
  plot_annotation(tag_levels = c('A', 'B', 'C', 'D'))

```

```{r}
#| label: fig-pcanoldp
#| fig-width: 12
#| fig-height: 8
#| fig-cap: 'Principal component analysis on the coordinates of the unbound oligonucleotide in the unrestrained microsecond trajectory. Scree plot (A) showing the contribution of each principal component on the total variance and the cumulative variance labelled on selected data points. (B) Score plots along the first two principal components, colored by kmeans clusters. (C) Sum of absolute loadings of residues for the first two first principal components'

p.pca.noldp

```

```{r}
#| label: pcanoldpcentroids
#| output: false

#create pdb object
pdb.noldp.pca <- pdb.noldp.h.long

#keep the coordinates of the centroids of the clusters
pdb.noldp.pca$xyz <- pdb.noldp.pca$xyz[c(4453, 14840, 20352),]

#align on the stem
pdb.noldp.pca$xyz <- fit.xyz(
  fixed = pdb.noldp.pca$xyz[1,],
  mobile = pdb.noldp.pca$xyz,
  fixed.inds = atom.select(
    pdb.noldp.pca, 
    resno = c(11:16, 34:36)
  )$xyz,
  mobile.inds = atom.select(
    pdb.noldp.pca,
    resno = c(11:16, 34:36)
  )$xyz
)

#write pdb file
if (!file.exists('data/pdb.noldp.pca.pdb')) {
  write.pdb(pdb.noldp.pca, file = 'data/pdb.noldp.pca.pdb')
}

```

![Structures of the centroids of each cluster (1 to 3 from left to right) superimposed with the first minimized structure of the complex (grey).](images/noldppdb.png){#fig-noldppdb}

```{r}
#| label: tableopenings
#| output: false

tbl.clusters.noldp <- table(pca.nb.noldp$Best.partition)
frac.clust.no.ldp.1 <- round(tbl.clusters.noldp[1]/sum(tbl.clusters.noldp)*100, 1)
frac.clust.no.ldp.2 <- round(tbl.clusters.noldp[2]/sum(tbl.clusters.noldp)*100, 1)
frac.clust.no.ldp.3 <- round(tbl.clusters.noldp[3]/sum(tbl.clusters.noldp)*100, 1)

open.times <- tibble(
  frame = 1:length(pca.nb.noldp$Best.partition),
  cluster = pca.nb.noldp$Best.partition
) %>% 
  #create a group column that increments every time the cluster changes
  mutate(
    group = cumsum(c(1, diff(cluster) != 0)),
    t = frame * info.long$time.per.frame
  ) %>% 
  group_by(cluster, group) %>%
  #calculate the time spent in each cluster
  summarise(
    t = max(t) - min(t)
  ) %>%
  group_by(cluster) %>%
  summarise(
    max.t = signif(max(t), 2),
    mean.t = signif(mean(t), 2)
  )


```

The unbound aptamer spends most of its time in the 'closed' (`{r} frac.clust.no.ldp.3`%) and 'intermediate' (`{r} frac.clust.no.ldp.2`%) configurations, and only occasionally visits the 'open', linear conformation (`{r} frac.clust.no.ldp.1`%) (@fig-openingevents). The 'closed' conformation has a mean lifetime of `{r} open.times$mean.t[open.times$cluster == 3]` ns and a max residency time of `{r} open.times$max.t[open.times$cluster == 3]` ns. The 'open' configuration has a lower mean lifetime `{r} open.times$mean.t[open.times$cluster == 1]` ns, and is never visited more than `{r} open.times$max.t[open.times$cluster == 1]` ns. The so-called 'intermediate' conformation is aptly named as it has an even lower mean lifetime `{r} open.times$mean.t[open.times$cluster == 2]` ns, and visits never exceed `{r} open.times$max.t[open.times$cluster == 2]` ns.

```{r}
#| label: openingevents
#| output: false
#| cache-lazy: false

p.open <- (
  tibble(
    frame = 1:length(pca.nb.noldp$Best.partition),
    cluster = pca.nb.noldp$Best.partition
  ) %>% 
    mutate(
      t = frame * info.long$time.per.frame,
      cluster = case_when(
        cluster == 1 ~ 'open',
        cluster == 2 ~ 'intermediate',
        cluster == 3 ~ 'closed',
        TRUE ~ 'who knows'
      ),
      cluster = factor(
        cluster,
        levels = c('open', 'intermediate', 'closed')
      )
    ) %>% 
    ggplot(aes(t, cluster)) +
    geom_point(aes(color = cluster),
               show.legend = FALSE,
               alpha = 0.25) +
    geom_line(alpha = 0.1) +
    custom.theme(1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_y_discrete(limits = rev) +
    scale_color_d3(name = 'Cluster') +
    labs(
      x = '',
      y = 'Conformation'
    ) 
) + (
  atom.distr(
    pair = "C24:C1'-A35:C1'", 
    pdb.noldp.long, span = 0.001, 
    info = info.long) %>% 
    ggplot(aes(x = t, y = dist.loess)) +
    geom_hline(
      aes(color = factor(frame), yintercept = dist),
      data = atom.distr(
        pair = "C24:C1'-A35:C1'", 
        pdb.noldp.pca, span = 1, 
        info = info.long),
      linetype = 'dashed', linewidth = 0.75
    ) +
    geom_line(alpha = 0.75) +
    custom.theme(1) +
    scale_x_continuous(expand = c(0, 0)) +
    scale_color_d3(name = 'Cluster') +
    labs(y = "Distance [C24:C1'-A35:C1'] (Å)",
         x = 't (ns)') 
) +
  plot_layout(guides = 'collect', nrow = 2)

###

```

```{r}
#| label: fig-openingevents
#| fig-width: 12
#| fig-height: 10
#| fig-cap: "Opening/closing of the unbound oligonucleotide during the microsecond simulation. Top: Cluster assignment for each frame. Bottom: Distance between the C1' atoms of C24 (loop region) and A35 (stem region). The distance measured for each cluster centroid is shown as a dashed line."

p.open


```
:::

### Long aptamer {#sec-long-aptamer}

The long aptamer was studied from starting structure including (`doplong8`) or devoid of (`doplong7`) base pairing between A8 and T37. The minimized structures shown in @fig-20static were obtained by minimizing the last frame of 10-ns \*restrained\* simulations. The dynamics of these systems was explored with 1 microsecond *unrestrained* simulations.

#### RMSD

##### 10 ns restrained simulations (for minimization)

::: justify
The RMSD (vs. the first frame) is quickly plateauing when the A38•T37 is pre-formed (and restrained), but never does in absence of this base pair (@fig-shortsim, @fig-pairwisedoplong).
:::

```{r}
#| label: fig-shortsim
#| fig-cap: "RMSD on all residues (including 5' and 3'-ends) without hydrogens across the 10 ns restrained simulation of the long aptamer, where the first frame is used as reference." 
#| fig-width: 8
#| fig-height: 5

short.longdop.traj <- list('doplong7' = doplong7.traj, 'doplong8' = doplong8.traj)

lapply(
  seq_along(short.longdop.traj),
  function(x) {
    data.frame(
      rmsd = rmsd(
        a = short.longdop.traj[[x]]$xyz[1,],
        b = short.longdop.traj[[x]]$xyz[seq(1, nrow(short.longdop.traj[[x]]$xyz), by = 1), ],
        fit = FALSE
      ),
      frame = (seq(1, nrow(short.longdop.traj[[x]]$xyz), by = 1)) #divider set to 1 to keep all frames
    ) %>% 
      mutate(t = frame * info.short$time.per.frame,
             traj = x)
  }
) %>% 
  do.call(rbind, .) %>% 
  mutate(traj = case_when(
    traj == 1 ~ 'doplong7 (no 8:37)',
    traj == 2 ~ 'doplong8 (with 8:37)'
  )) %>% 
  ggplot(aes(x = t, y = rmsd, color = traj)) +
  geom_line() +
  custom.theme(1) +
  labs(
    x = 't (ns)',
    y = 'RMSD (&angst;)',
  ) +
  scale_color_manual(values = c('#e62e73', '#003e83')) 
###
```

```{r}
#| label: pairwisedoplongshort
#| output: false

rmsd.pairwise.doplong7 <- rmsd(
  a = doplong7.traj$xyz[seq(1, nrow(doplong7.traj$xyz), by = nrow(doplong7.traj$xyz)/1000),], #keep one out of 5 frames
  fit = FALSE
) %>%
  lazy_dt() %>%
  mutate(var1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
  mutate(var2 = as.numeric(gsub("V", "", var2))) %>%
  mutate(traj = 'doplong7') %>% 
  as.data.table()

rmsd.pairwise.doplong8 <- rmsd(
  a = doplong8.traj$xyz[seq(1, nrow(doplong8.traj$xyz), by = nrow(doplong8.traj$xyz)/1000),], #keep one out of 5 frames
  fit = FALSE
) %>%
  lazy_dt() %>%
  mutate(var1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
  mutate(var2 = as.numeric(gsub("V", "", var2))) %>%
  mutate(traj = 'doplong8') %>% 
  as.data.table()

rbind(rmsd.pairwise.doplong7, rmsd.pairwise.doplong8) %>% 
arrow::write_feather(., "data/doplong.rmsd.pairwise.feather")
###
```

```{r}
#| label: fig-pairwisedoplong
#| fig-cap: All-to-all RMSD (all residues and ligand atoms except H) for the 10 ns restrained simulation of the long aptamer.
#| fig-width: 12
#| fig-height: 7.5

arrow::read_feather('data/doplong.rmsd.pairwise.feather') %>%
  mutate(t1 = var1 * info.short$time.per.frame * nrow(doplong7.traj$xyz)/1000,
         t2 = var2 * info.short$time.per.frame * nrow(doplong7.traj$xyz)/1000) %>%
  ggplot(aes(x = t1, y = t2, fill = value)) +
  facet_wrap(~traj) +
  geom_raster() +
  scale_fill_viridis_c(name = 'RMSD (&angst;)') +
  scale_x_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  ggthemes::theme_pander() +
  theme(
    axis.text.x = element_markdown(size = 12 * scaling),
    axis.text.y = element_markdown(size = 12 * scaling),
    axis.title.x = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.title.y = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.line = element_line(linewidth = 0.75 * scaling),
    axis.ticks = element_line(linewidth = 0.75 * scaling, color = 'black'),
    #legend.text = element_markdown(size = 12 * scaling),
    legend.title = element_markdown(size = 14 * scaling, face = 'bold'),
    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 = element_blank(),
    panel.grid.minor = 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)
  ) +
  theme(legend.position = 'bottom') +
  theme(plot.margin = margin(1, 1, 0, 0, 'lines'))
###
```

##### Unrestrained microsecond simulations

::: justify
In absence of restraints, doplong7 settles at RMSD values close to that of the final frame from the restrained 10 ns simulation (@fig-longsim). The pairwise data also reveals that an event clearly occurs after 56 ns (@fig-pairwisedoplonglonger), that is close to the start of the microsecond simulation but past the endpoint of the restrained 10 ns simulation. After this conformational change, the pairwise RMSD remains apparently stable, but other changes still occur (see below).

In presence of the pre-formed A38•T37 (doplong8), the simulation plateaus much faster (@fig-longsim)., as seen above for the retrained production run, but at about twice the RMSD. A very significant event is visible at 576 ns, resulting in a doubled RMSD vs. the first frame. The RMSD then varies much more significantly than before, with another significant change at 663 ns.
:::

```{r}
#| label: fig-longsim
#| fig-cap: "RMSD on all residues (including 5' and 3'-ends) without hydrogens across the microsecond unrestrained simulation of the long aptamer, where the first frame is used as reference. The arrows point to conformational changes that are not all evident from this metric alone." 
#| fig-width: 8
#| fig-height: 5

doplong7_longtraj.traj <- list('doplong7' = doplong7_longtraj, 'doplong8' = doplong8_longtraj)

lapply(
  seq_along(doplong7_longtraj.traj),
  function(x) {
    data.frame(
      rmsd = rmsd(
        a = doplong7_longtraj.traj[[x]]$xyz[1,],
        b = doplong7_longtraj.traj[[x]]$xyz[seq(1, nrow(doplong7_longtraj.traj[[x]]$xyz), by = 1), ],
        fit = FALSE
      ),
      frame = (seq(1, nrow(doplong7_longtraj.traj[[x]]$xyz), by = 1)) #divider set to 1 to keep all frames
    ) %>% 
      mutate(t = frame * info.long$time.per.frame,
             traj = x)
  }
) %>% 
  do.call(rbind, .) %>% 
  mutate(traj = case_when(
    traj == 1 ~ 'doplong7 (no 8:37)',
    traj == 2 ~ 'doplong8 (with 8:37)'
  )) %>% 
  ggplot(aes(x = t, y = rmsd, color = traj)) +
  geom_line() +
  custom.theme(1) +
  labs(
    x = 't (ns)',
    y = 'RMSD (&angst;)',
  ) +
  scale_color_manual(values = c('doplong7 (no 8:37)' = '#e62e73', 'doplong8 (with 8:37)' = '#003e83')) +
  geom_segment(
    inherit.aes = FALSE,
    data = data.frame(
      traj = c(rep('doplong7 (no 8:37)', 3), rep('doplong8 (with 8:37)', 2)),
      x = c(56.5, 375, 553, 576, 663),
      y = c(13, 13.5, 14, 14, 15)
    ),
    aes(
    x = x, xend = x,
    y = y, yend = (y - 1),
    color = traj
    ),
    arrow = arrow(type = "closed", length = unit(0.02, "npc")),
    linewidth = 1
  )

```

```{r}
#| label: pairwisedoplongtrajlonb
#| output: false

rmsd.pairwise.doplong7_longtraj <- rmsd(
  a = doplong7_longtraj$xyz[seq(1, nrow(doplong7_longtraj$xyz), by = nrow(doplong7_longtraj$xyz)/1000),], #keep one out of 5 frames
  fit = FALSE
) %>%
  lazy_dt() %>%
  mutate(var1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
  mutate(var2 = as.numeric(gsub("V", "", var2))) %>%
  mutate(traj = 'doplong7') %>% 
  as.data.table()

rmsd.pairwise.doplong8_longtraj <- rmsd(
  a = doplong8_longtraj$xyz[seq(1, nrow(doplong8_longtraj$xyz), by = nrow(doplong8_longtraj$xyz)/1000),], #keep one out of 5 frames
  fit = FALSE
) %>%
  lazy_dt() %>%
  mutate(var1 = 1:nrow(.), .before = 1) %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "value") %>%
  mutate(var2 = as.numeric(gsub("V", "", var2))) %>%
  mutate(traj = 'doplong8') %>% 
  as.data.table()

rbind(rmsd.pairwise.doplong7_longtraj, rmsd.pairwise.doplong8_longtraj) %>% 
arrow::write_feather(., "data/doplong.longtraj.rmsd.pairwise.feather")
###
```

```{r}
#| label: fig-pairwisedoplonglonger
#| fig-cap: All-to-all RMSD (all residues and ligand atoms except H) for the unrestrained microsecond simulation. The white arrows point to significant events. The arrows point to conformational changes that are not all evident from this metric alone.
#| fig-width: 12
#| fig-height: 7.5

arrow::read_feather('data/doplong.longtraj.rmsd.pairwise.feather') %>%
  mutate(t1 = var1 * info.long$time.per.frame * nrow(doplong7_longtraj$xyz)/1000,
         t2 = var2 * info.long$time.per.frame * nrow(doplong7_longtraj$xyz)/1000) %>%
  ggplot(aes(x = t1, y = t2, fill = value)) +
  facet_wrap(~traj) +
  geom_raster() +
  scale_fill_viridis_c(name = 'RMSD (&angst;)') +
  scale_x_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA), 
                     name = 't (ns)') +
  ggthemes::theme_pander() +
  theme(
    axis.text.x = element_markdown(size = 12 * scaling),
    axis.text.y = element_markdown(size = 12 * scaling),
    axis.title.x = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.title.y = element_markdown(size = 14 * scaling, face = 'bold'),
    axis.line = element_line(linewidth = 0.75 * scaling),
    axis.ticks = element_line(linewidth = 0.75 * scaling, color = 'black'),
    #legend.text = element_markdown(size = 12 * scaling),
    legend.title = element_markdown(size = 14 * scaling, face = 'bold'),
    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 = element_blank(),
    panel.grid.minor = 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)
  ) +
  theme(legend.position = 'bottom') +
  theme(plot.margin = margin(1, 1, 0, 0, 'lines')) +
  geom_segment(
    inherit.aes = FALSE,
    data = data.frame(
      t1 = c(56.5, 375, 553, 576, 663), 
      t2 = c(56.5, 375, 553, 576, 663), 
      traj = c(rep('doplong7', 3), rep('doplong8', 2))),
    aes(x = (t1 + 50), y = (t2 + 50), xend = t1, yend = t2),
    arrow = arrow(type = "closed", length = unit(0.02, "npc")),
    linewidth = 1, color = 'white'
  )
###
```

#### Base pairing

##### A38•T37 formation

::: justify
Analysis of base pairing along the unrestrained simulation provides a straightforward explanation for the large conformational change of doplong7: the A8 and T37, initially far apart, rapidly base pair (canonical cW-W). This is accompanied by the formation of a base pair between G9 and the Hoogsteen face of G13. The latter is also rapidly formed in doplong8, wherein the A38•T37 bp is initially formed; however the base pairing is 'inverted':

-   doplong7: tW-M (VII in Saenger nomenclature), with G9•N2-G13:O6 and G9•N1-G13-N7 H-bonds

-   doplong8: a cW+M (VI in Saenger nomenclature), with G9•N2-G13-N7 and G9•N1-G13:O6 H-bonds

Note that this base pair is fairly unstable in doplong8, but A38•T37 remains formed nonetheless (reminder: there is no restraint here).
:::

```{r}
#| label: bpdoplong
#| output: false

doplong7h_100frames <- doplong7h_longtraj
doplong7h_100frames$xyz <- doplong7h_100frames$xyz[seq(1, nrow(doplong7h_100frames$xyz), by = nrow(doplong7h_100frames$xyz)/100),]

pair.list.doplong7 <- contactr(
  doplong7h_100frames,
  dist.thresh = 3.1, freq.thresh = 0.2,
  lgd.resno = 45,
  lgd.donors = c('O2', 'O1', 'N1'),
  lgd.acceptors = NA,
  donor.atoms = c('N', 'O', 'S', 'F'),
  acceptor.atoms = c('O', 'N', 'S', 'F', 'Cl', 'Br', 'I')
)

pair.dist.doplong7 <- lapply(
  pair.list.doplong7$pair,
  atom.distr,
  pdb.input = doplong7h_longtraj,
  span = 0.01
) %>%
  do.call(rbind, .)

doplong8h_100frames <- doplong8h_longtraj
doplong8h_100frames$xyz <- doplong8h_100frames$xyz[seq(1, nrow(doplong8h_100frames$xyz), by = nrow(doplong8h_100frames$xyz)/100),]

pair.list.doplong8 <- contactr(
  doplong8h_100frames,
  dist.thresh = 3.1, freq.thresh = 0.2,
  lgd.resno = 45,
  lgd.donors = c('O2', 'O1', 'N1'),
  lgd.acceptors = NA,
  donor.atoms = c('N', 'O', 'S', 'F'),
  acceptor.atoms = c('O', 'N', 'S', 'F', 'Cl', 'Br', 'I')
)

pair.dist.doplong8 <- lapply(
  pair.list.doplong8$pair,
  atom.distr,
  pdb.input = doplong8h_longtraj,
  span = 0.01
) %>%
  do.call(rbind, .)

###
```

```{r}
#| label: fig-bpdoplong
#| fig-cap: Length of potential A38•T37, and G9•G13 H-bonds during the simulation. Note that the G9•G13 pairing differs. The pink arrow shows the onset of basepairing for doplong7.
#| fig-width: 8
#| fig-height: 10

scaling <- 1

pair.dist.doplong <- rbind(
  pair.dist.doplong7 %>% mutate(traj = 'doplong7'),
  pair.dist.doplong8 %>% mutate(traj = 'doplong8')
) %>% 
  filter(
    (str_detect(pair, 'T37') & str_detect(pair, 'A8')) |
      (str_detect(pair, 'G9') & str_detect(pair, 'G13'))
    ) %>% 
  arrange(pair)

pair.levels <- unique(pair.dist.doplong$pair)

pair.dist.doplong %>% 
  ggplot(aes(t, dist.loess, color = factor(pair, levels = pair.levels))) +
  facet_wrap(~traj, ncol = 1, nrow = 2) +
  geom_hline(yintercept = 3.0, linetype = 'dashed') +
  geom_text_repel(
    data = . %>% filter(t == max(t)),
    aes(label = pair),
    xlim = c(1020, 1100),
    fontface = 'bold',
    show.legend = FALSE,
    size = 5*scaling
  ) +
  geom_line(show.legend = FALSE,
            linewidth = 0.75 * scaling) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, NA), expand = c(0,0)) +
  scale_color_d3() +
  labs(x = 'Time (ns)', y = 'Distance (Å)') +
  custom.theme(scaling) +
  theme(plot.margin = margin(1, 150, 1, 1)) +
  coord_cartesian(clip = 'off') +
  geom_segment(
    data = data.frame(traj = 'doplong7'),
    inherit.aes = FALSE,
    x = 56, xend = 56,
    y = 19, yend = 17,
    arrow = arrow(type = "closed", length = unit(0.02, "npc")),
    linewidth = 1, color = '#e62e73'
  )


```

##### Other events

::: justify
The above does not explain what happens in the second part of doplong8's simulation, nor is the only event occuring for doplong7.

-   doplong7: Two discrete events occur at around 375 and 553 ns

    -   First, two of the three H-bonds from C24•G29 bonds get disrupted

    -   Then the third gives up as well. The neighboring G25•G28 base pair is also disrupted at this last onset (@fig-clustpcadoplong).

-   doplong8: Two discrete events as well, at around 576 and 663 ns.

    -   First occurs the formation of the consecutive A15•A35, G16•G34 and G17•G33 base pairs in lieu of G34•T11, which is the expected base pair from NMR data of the short aptamer.

    -   Second, the C24•G29 and G25•G28 base pairs are lost, similar to doplong7 (@fig-clustpcadoplong).

    -   Note also the transient flipping out of G9 after these events. The G36 base pair remains in contact, although it now adopts a co-planar configuration.

All of these conformational changes occur within the region of the aptamer present in the smaller version, and are very likely to be the consequence of the absence of retraints in this microsecond simulation; since all lost base pairs are present in the restrained simulations.
:::

```{r}
#| label: fig-bpdoplongchanges
#| fig-cap: H-bond lengths for selected nucleotides undergoing base pair formation/disruption during the unrestrained microsecond simulation. The vertical lines indicate the onset of formation/disruption and correspond to the arrows in the figures above.
#| fig-width: 12
#| fig-height: 12

rbind(
  pair.dist.doplong7 %>% mutate(traj = 'doplong7'),
  pair.dist.doplong8 %>% mutate(traj = 'doplong8')
) %>% 
  arrange(pair) %>% 
  filter(
    (str_detect(pair, 'A15') & str_detect(pair, 'A35')) |
      (str_detect(pair, 'A35') & str_detect(pair, 'A15') | str_detect(pair, 'T11')) |
      (str_detect(pair, 'C24') & str_detect(pair, 'G29')) |
      (str_detect(pair, 'G16') & str_detect(pair, 'G34')) |
      (str_detect(pair, 'G17') & str_detect(pair, 'G33') | str_detect(pair, 'G34')) |
      (str_detect(pair, 'G25') & str_detect(pair, 'G28')) |
      (str_detect(pair, 'G29') & str_detect(pair, 'C24') | str_detect(pair, 'C24')) |
      (str_detect(pair, 'G34') & str_detect(pair, 'T11')) |
      (str_detect(pair, 'G11')) |
      (str_detect(pair, 'A8')) |
      str_detect(pair, 'G9') |
      str_detect(pair, 'T18')  |
      str_detect(pair, 'G36')
  ) %>% 
  filter(!str_detect(pair, 'LDP')) %>%  
  ggplot(aes(t, dist.loess, color = traj)) +
  facet_wrap(~pair, ncol = 3) +
  geom_line(show.legend = TRUE,
            linewidth = 0.75 * scaling) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, NA), expand = c(0,0)) +
  scale_color_manual(
    name = 'trajectory',
    values = c('doplong7' = '#e62e73', 'doplong8' = '#003e83')
    ) +
  scale_linetype_manual(values = c(2:4), name = 'onset') +
  labs(x = 'Time (ns)', y = 'Distance (Å)') +
  custom.theme(scaling) +
  theme(plot.margin = margin(1, 150, 1, 1)) +
  coord_cartesian(clip = 'off') +
  geom_vline(
    data = data.frame(
      traj = c(rep('doplong7', 3), rep('doplong8', 2)),
      x = c(56.5, 375, 553, 576, 663),
      onset = c('1', '2', '3', '1', '2')
    ),
    aes(xintercept = x, color = traj, linetype = onset),
    show.legend = FALSE,
    linewidth = 0.75 * scaling
  ) +
  theme(legend.position = 'bottom') +
  guides(linetype = 'none')

```

#### PCA

Principal component analysis allows to extract coordinates characteristic of the conformers prior and after the conformational changes.

```{r}
#| label: pcadoplong7
#| output: false

pca_doplong7_longtraj <- pca.pdbs(
  doplong7_longtraj,
  use.svd = FALSE,
  rm.gaps = TRUE,
  fit = FALSE #pdb models are already aligned at the import stage
)

```

```{r}
#| label: pcadoplong8
#| output: false

pca_doplong8_longtraj <- pca.pdbs(
  doplong8_longtraj,
  use.svd = FALSE,
  rm.gaps = TRUE,
  fit = FALSE #pdb models are already aligned at the import stage
)

```

```{r}
#| label: screedoplong
#| output: false

scree7 <- data.frame(
  pc = 1:length(pca_doplong7_longtraj$L),
  L = pca_doplong7_longtraj$L
) %>%
  mutate(
    var = L/sum(L) * 100,
    cum.var = cumsum(var)
  ) %>%
  filter(pc <= 15) %>%
  select(-L) %>%
  mutate(
    label = ifelse(
      pc %in% 1:3 |
        pc %in% seq(4, nrow(.), 3) |
        pc == nrow(.),
      signif(cum.var, 3), NA)
  ) %>%
  ggplot(., aes(x = pc, y = var)) +
  geom_text_repel(aes(label = label), size = 5, fontface = 'bold', force = 100) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 2) +
  custom.theme(scaling) +
  labs(
    x = 'Number of Principal Components',
    y = 'Proportion of variance (%)'
  )

scree8 <- data.frame(
  pc = 1:length(pca_doplong8_longtraj$L),
  L = pca_doplong8_longtraj$L
) %>%
  mutate(
    var = L/sum(L) * 100,
    cum.var = cumsum(var)
  ) %>%
  filter(pc <= 15) %>%
  select(-L) %>%
  mutate(
    label = ifelse(
      pc %in% 1:3 |
        pc %in% seq(4, nrow(.), 3) |
        pc == nrow(.),
      signif(cum.var, 3), NA)
  ) %>%
  ggplot(aes(x = pc, y = var)) +
  geom_text_repel(aes(label = label), size = 5, fontface = 'bold', force = 100) +
  geom_line(linewidth = 0.75) +
  geom_point(size = 2) +
  custom.theme(scaling) +
  labs(
    x = 'Number of Principal Components',
    y = 'Proportion of variance (%)'
  )
```

```{r}
#| label: nbclust7
#| output: false

pca.nb.doplong7_longtraj <- NbClust(
  data.frame(pca_doplong7_longtraj$z) %>% select(all_of(1:4)),
  distance = "euclidean",
  min.nc = 2, max.nc = 5,
  method = "kmeans"
)

```

```{r}
#| label: nblust8
#| output: false

pca.nb.doplong8_longtraj <- NbClust(
  data.frame(pca_doplong8_longtraj$z) %>% select(all_of(1:4)),
  distance = "euclidean",
  min.nc = 2, max.nc = 5,
  method = "kmeans"
)

```

```{r}
#| label: prepfigpcadoplong
#| output: false
#| cache-lazy: false

p.loads.doplong7_longtraj <- data.frame(
  u = c(pca_doplong7_longtraj$U[,1], pca_doplong7_longtraj$U[,2], pca_doplong7_longtraj$U[,3]),
  eleno = rep(doplong7_longtraj$atom$eleno, 3),
  PC = c(rep('PC1', length(pca_doplong7_longtraj$U[,1])),
         rep('PC2', length(pca_doplong7_longtraj$U[,2])),
         rep('PC3', length(pca_doplong7_longtraj$U[,3])))
) %>%
  left_join(
    doplong7_longtraj$atom,
    by = c("eleno" = "eleno")
  ) %>%
  group_by(resno, resid, PC, chain) %>%
  summarise(u = sum(abs(u))) %>%
  mutate(
    resno = if_else(resid == 'LDP', 45, resno),
    resid = gsub('^D', '', resid),
    label = if_else(
      resid == 'LDP',
      #make LDP bold and pink
      paste0('<span style="color:pink"><b>', resid, '</b></span>'),
      paste0('<b>', resid, resno, '</b>')
    )
  ) %>%
  filter(PC %in% paste0('PC', 1:2)) %>%
  ggplot(aes(
    x = factor(
      resno,
      levels = unique(resno),
      labels = unique(label)
    ),
    y = u,
    fill = chain
  )
  ) +
  facet_wrap(~PC, ncol = 1) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c('grey', 'pink')) +
  scale_y_continuous(expand = c(0, 0)) +
  custom.theme(0.75) +
  theme(axis.title.x = element_blank(),
        legend.position = 'none') +
  labs(y = 'sum of absolute loadings')


p.pca.doplong7_longtraj <- scree7 +
  pca.plotr.traj(pca_doplong7_longtraj, 'PC1', 'PC2', nb = pca.nb.doplong7_longtraj) +
  p.loads.doplong7_longtraj +
  plot_layout(
    design = '
              AB
              CC
              ') &
  plot_annotation(tag_levels = c('A', 'B', 'C', 'D'))

p.loads.doplong8_longtraj <- data.frame(
  u = c(pca_doplong8_longtraj$U[,1], pca_doplong8_longtraj$U[,2], pca_doplong8_longtraj$U[,3]),
  eleno = rep(doplong8_longtraj$atom$eleno, 3),
  PC = c(rep('PC1', length(pca_doplong8_longtraj$U[,1])),
         rep('PC2', length(pca_doplong8_longtraj$U[,2])),
         rep('PC3', length(pca_doplong8_longtraj$U[,3])))
) %>%
  left_join(
    doplong8_longtraj$atom,
    by = c("eleno" = "eleno")
  ) %>%
  group_by(resno, resid, PC, chain) %>%
  summarise(u = sum(abs(u))) %>%
  mutate(
    resno = if_else(resid == 'LDP', 45, resno),
    resid = gsub('^D', '', resid),
    label = if_else(
      resid == 'LDP',
      #make LDP bold and pink
      paste0('<span style="color:pink"><b>', resid, '</b></span>'),
      paste0('<b>', resid, resno, '</b>')
    )
  ) %>%
  filter(PC %in% paste0('PC', 1:2)) %>%
  ggplot(aes(
    x = factor(
      resno,
      levels = unique(resno),
      labels = unique(label)
    ),
    y = u,
    fill = chain
  )
  ) +
  facet_wrap(~PC, ncol = 1) +
  geom_bar(stat = "identity") +
  scale_fill_manual(values = c('grey', 'pink')) +
  scale_y_continuous(expand = c(0, 0)) +
  custom.theme(0.75) +
  theme(axis.title.x = element_blank(),
        legend.position = 'none') +
  labs(y = 'sum of absolute loadings')


p.pca.doplong8_longtraj <- scree8 +
  pca.plotr.traj(pca_doplong8_longtraj, 'PC1', 'PC2', nb = pca.nb.doplong8_longtraj) +
  p.loads.doplong8_longtraj +
  plot_layout(
    design = '
              AB
              CC
              ') &
  plot_annotation(tag_levels = c('A', 'B', 'C', 'D'))

```

```{r}
#| label: fig-pcadoplong7
#| fig-width: 12
#| fig-height: 8
#| fig-cap: 'Principal component analysis on the coordinates of the long aptamer doplong7 (starting without A38•T37 bp) in the unrestrained microsecond trajectory. Scree plot (A) showing the contribution of each principal component on the total variance and the cumulative variance labelled on selected data points. (B) Score plots along the first two principal components, colored by kmeans clusters. (C) Sum of absolute loadings of residues for the first two first principal components'

p.pca.doplong7_longtraj

```

```{r}
#| label: fig-pcadoplong8
#| fig-width: 12
#| fig-height: 8
#| fig-cap: 'Principal component analysis on the coordinates of the long aptamer doplong8 (starting with A38•T37 bp) in the unrestrained microsecond trajectory. Scree plot (A) showing the contribution of each principal component on the total variance and the cumulative variance labelled on selected data points. (B) Score plots along the first two principal components, colored by kmeans clusters. (C) Sum of absolute loadings of residues for the first two first principal components'

p.pca.doplong8_longtraj

```

```{r}
#| label: pcadoplongpdbexport
#| output: false

#create pdb object
pca.doplong7_longtraj.pca <- doplong7_longtraj
pca.doplong8_longtraj.pca <- doplong8_longtraj

#keep the coordinates of the centroids of the clusters
pca.doplong7_longtraj.pca$xyz <- pca.doplong7_longtraj.pca$xyz[c(750, 24576),]
pca.doplong8_longtraj.pca$xyz <- pca.doplong8_longtraj.pca$xyz[c(2035, 21687),]

#align
pca.doplong7_longtraj.pca$xyz <- fit.xyz(
  fixed = pca.doplong7_longtraj.pca$xyz[1,],
  mobile = pca.doplong7_longtraj.pca$xyz,
  fixed.inds = atom.select(pca.doplong7_longtraj.pca, elety = unique(pca.doplong7_longtraj.pca$atom$elety))$xyz,
  mobile.inds = atom.select(pca.doplong7_longtraj.pca, elety = unique(pca.doplong7_longtraj.pca$atom$elety))$xyz
)
pca.doplong8_longtraj.pca$xyz <- fit.xyz(
  fixed = pca.doplong8_longtraj.pca$xyz[1,],
  mobile = pca.doplong8_longtraj.pca$xyz,
  fixed.inds = atom.select(pca.doplong8_longtraj.pca, elety = unique(pca.doplong8_longtraj.pca$atom$elety))$xyz,
  mobile.inds = atom.select(pca.doplong8_longtraj.pca, elety = unique(pca.doplong8_longtraj.pca$atom$elety))$xyz
)

#write pdb file
if (!file.exists('data/pca.doplong7_longtraj.pca.pdb')) {
  write.pdb(pca.doplong7_longtraj.pca, file = 'data/pca.doplong7_longtraj.pca.pdb')
}
if (!file.exists('data/pca.doplong8_longtraj.pca.pdb')) {
  write.pdb(pca.doplong8_longtraj.pca, file = 'data/pca.doplong8_longtraj.pca.pdb')
}

```

::: {#fig-clustpcadoplong layout-ncol="2"}
![<b>pcadoplong7:</b> The A38•T37 base pair is initially absent but canonically formed in cluster 2. However, the latter loses the C24•G29 G25•G28 base pairs. The two events are not directly linked as they occur around half a microsecond from one another.](images/pcadoplong7.png)

![<b>pcadoplong8</b>: The second cluster is characterized by the zipping of the 15:17/33:35 region (between the binding site, above, and the duplex stem, below), while the 24:25 nucleotides flip out of the loop region. In each panel, the starting cluster is characterized by the interaction between T34 (green; foreground) with G34, which is the expected base pair with restraints. Note the difference in backbone geometry.](images/pcadoplong8.png)

Superimposed coordinates of the cluster centroids, aligned by their duplex stem. Clusters 1 and 2 and on the left (with blue dopamine) and right (with orange dopamine), respectively.

![Initially, G9 is paired to G13, with G36 forming an orthogonal contact to the latter (left). In the second cluster' centroid coordinates, G9 is flipped out and G36 is coplanar to G13, hence forming a more canonical base pair (right).](images/G36.png){fig-align="center"}
:::
