In epongpipat/ppi: Psychophysiological Interaction (PPI)

I believe it is helpful to think of PPI as a essentially a lot of data manipulation/wrangling followed by multi-level modeling.

This example will focus on how to perform the data manipulation portion of PPI.

This example assumes that the reader:

1. is familiar with standard task fMRI analyses and jargon
2. has a conceptual understanding of psychophysiological interaction (PPI)
3. has pre-processed the fMRI data

R Packages

First, let's load some packages.

packages <- c("tidyverse", "eepR", "furrr", "ggpubr")
xfun::pkg_attach(packages, install = T, message = F)

script_dir <- "~/Documents/GitHub/ppi/R/"
script_files <- list.files(script_dir)
script_path <- paste0(script_dir, script_files)
script_ppi <- lapply(script_path, source)

Define Parameters

Let's define some MRI parameters and an upsample factor for convolution and deconvolution.

In this example, we are going to use SPM's first gamma variate (spmg1) for the HRF. This is SPM's default HRF other current options include gam and block.

tr <- 1.5
n_volumes <- 260
upsample_factor <- 16
hrf_name <- "spmg1"

events_file <- "~/Box/my_mri/behavioral/nback_3tb4188_2017_Sep_13_1341_soa.csv"

psy_contrast_table <- cbind(nback_vs_baseline = c(1, 1, 1, 1, -4)/5,
                                              task_vs_control = c(-3, 1, 1, 1, 0)/4,
                                              linear_task = c(0, -1, 0, 1, 0),
                                              quadratic_task = c(0, -1, 2, -1, 0)/3)
psy_contrast_table

out_file <- NULL

Ideal HRF

Let's first obtain the hemodynamic response function (HRF) that we are going to use during convolution and deconvolution steps.

hrf <- get_hrf_afni(hrf_name, tr, upsample_factor)

ggplot(hrf %>% mutate(volume = row_number()), aes(volume, hrf)) +
  geom_line() +
  theme_minimal() +
  labs(title = paste0("HRF: ", toupper(hrf_name)),
       subtitle = "Upsample: ", upsample_factor)

Psychological Variables

Label Volume by Trial Type

df_psy <- read.csv(events_file) %>%
  rename(onset = start, trial_type = block) %>%
  mutate(duration = ifelse(trial_type == "0-back", 10 * 2.5, 20 * 2.5)) %>%
  group_by(run) %>%
  nest() %>%
  mutate(data = future_map(data, function(x) x %>% select(onset, duration, trial_type)),
         data_trial_type = future_pmap(list(data, tr, n_volumes),  create_trial_type_by_volume_list)) %>%
  select(-data) %>%
  unnest(data_trial_type) %>%
  ungroup()
head(df_psy)

Contrast Code

func_fig_contrast <- function(data) {
  data %>%
    as.data.frame() %>%
    mutate(volume = row_number()) %>%
    gather(., psy_contrast, value, -volume) %>%
    ggplot(., aes(volume, value)) +
    geom_line() +
    facet_wrap(~ psy_contrast) +
    theme_minimal()
}

func_add_fig_title <- function(fig, text) {
  fig + 
    labs(title = text)
}

df_psy_contrast <- contrast_code_categorical_variable(df_psy, psy_contrast_table) %>%
  group_by(run) %>%
  nest() %>%
  mutate(data = future_map(data, function(x) x %>% select(contains("psy"))),
         fig = future_map(data, func_fig_contrast),
         fig = future_map2(fig, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_psy_contrast$fig, ncol = 1)

Upsample

df_psy_upsample <- df_psy_contrast %>%
  mutate(data_upsample = future_map(data, function(x) apply(x, 2, function(x) upsample(x, upsample_factor))),
         fig_upsample = future_map(data_upsample, func_fig_contrast),
         fig_upsample = future_map2(fig_upsample, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_psy_upsample$fig_upsample, ncol = 1)

Convolve

df_psy_convolve <- df_psy_upsample %>%
  mutate(data_convolve = future_map(data_upsample, function(x) apply(x, 2, function(x) convolve_afni(x, hrf, tr, n_volumes, upsample_factor))),
         fig_convolve = future_map(data_convolve, func_fig_contrast),
         fig_convolve = future_map2(fig_convolve, paste0("Run:", run), func_add_fig_title))


ggarrange(plotlist = df_psy_convolve$fig_convolve, ncol = 1)

Downsample

df_psy_downsample <- df_psy_convolve %>%
  mutate(data_downsample = future_map(data_convolve, function(x) apply(x, 2, function(x) downsample(x, upsample_factor))),
         fig_downsample = future_map(data_downsample, func_fig_contrast),
         fig_downsample = future_map2(fig_downsample, paste0("Run:", run), func_add_fig_title))


ggarrange(plotlist = df_psy_downsample$fig_downsample, ncol = 1)

Physiological Variable

Create Sphere(s)

This step can be skipped if have you have your own mask.

Extract ROI

func_fig_phys <- function(data) {
  colnames(data) <- "data"

  data <- as.data.frame(data)

  fig <- data %>%
    mutate(volume = row_number()) %>%
    ggplot(., aes(volume, data)) +
    geom_line() +
    theme_minimal()
  return(fig)
}

dir_phys <- "/Volumes/shared/KK_KR_JLBS/Wave1/MRI/FMRI/PPI/Nback_individual_covariates/PHYS_MNI_42_-42_42/3tb1780/"
files_phys <- list.files(dir_phys, "_mean.1D")
path_phys <- paste0(dir_phys, files_phys)
df_phys <- tibble(file = files_phys,
                  path = path_phys) %>%
  mutate(run = str_remove(file, "_mean.1D"),
         run = str_remove(run, "r"),
         data = future_map(path, function(x) read.csv(x, header = F, col.names = "data")),
         fig = future_map(data, func_fig_phys),
         fig = future_map2(fig, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_phys$fig, ncol = 1)

Detrend

df_phys_detrend <- df_phys %>%
  mutate(data_detrend = future_map(data, function(x) detrend(x, 2)),
         fig_detrend = future_map(data_detrend, func_fig_phys),
         fig_detrend = future_map2(fig_detrend, paste0("Run: ", run), func_add_fig_title))


ggarrange(plotlist = df_phys_detrend$fig_detrend, ncol = 1)

Upsample

df_phys_upsample <- df_phys_detrend %>%
  mutate(data_upsample = future_map(data_detrend, function(x) apply(x, 2, function(x) upsample(x, upsample_factor))),
         fig_upsample = future_map(data_upsample, func_fig_phys),
         fig_upsample = future_map2(fig_upsample, paste0("Run: ", run), func_add_fig_title))

ggarrange(plotlist = df_phys_upsample$fig_upsample, ncol = 1)

Deconvolve

df_phys_deconvolve <- df_phys_upsample %>%
  mutate(data_deconvolve = future_map(data_upsample, function(x) apply(x, 2, function(x) deconvolve_afni(x, hrf))),
         fig_deconvolve = future_map(data_deconvolve, func_fig_phys),
         fig_deconvolve = future_map2(fig_deconvolve, paste0("Run:", run), func_add_fig_title))

ggarrange(plotlist = df_phys_deconvolve$fig_deconvolve, ncol = 1)

Psychophysiological Interaction Variable

Multiply/Create Interaction

df_ppi <- tibble(run = c(1:3),
                 data = NA)

for (i in 1:nrow(df_ppi)) {
  df_ppi$data[i] <- list(apply(df_psy_upsample$data_upsample[[i]], 2, function(x) x * df_phys_deconvolve$data_deconvolve[[i]]))
}

df_ppi <- df_ppi %>%
  mutate(fig = future_map(data, func_fig_contrast),
         fig = future_map2(fig, paste0("Run:", run), func_add_fig_title))

ggarrange(plotlist = df_ppi$fig, ncol = 1)

Convolve

df_ppi_convolve <- df_ppi %>%
  mutate(data_convolve = future_map(data, function(x) apply(x, 2, function(x) convolve_afni(x, hrf, tr, n_volumes, upsample_factor))),
         fig_convolve = future_map(data_convolve, func_fig_contrast),
         fig_convolve = future_map2(fig_convolve, paste0("Run:", run), func_add_fig_title))

ggarrange(plotlist = df_ppi_convolve$fig_convolve, ncol = 1)

Downsample

df_ppi_downsample <- df_ppi_convolve %>%
  mutate(data_downsample = future_map(data_convolve, function(x) apply(x, 2, function(x) downsample(x, upsample_factor))),
         fig_downsample = future_map(data_downsample, func_fig_contrast),
         fig_downsample = future_map2(fig_downsample, paste0("Run:", run), func_add_fig_title))

ggarrange(plotlist = df_ppi_downsample$fig_downsample, ncol = 1)

Design Matrix

Initialize

df_design_mat <- tibble(run = c(1:3),
                        data = NA)

for (i in 1:nrow(df_design_mat)) {

  temp_psy <- df_psy_downsample$data_downsample[[i]]
  temp_phys <- df_phys_detrend$data_detrend[[i]] %>% rename(phys = residual)
  temp_ppi <- df_ppi_downsample$data_downsample[[i]]
  colnames(temp_ppi) <- str_replace(colnames(temp_ppi), "psy_", "ppi_")

  df_design_mat$data[i] <- list(cbind(temp_psy, temp_phys, temp_ppi))
}

df_design_mat <- df_design_mat %>%
  unnest()

func_fig_data_heat_map <- function(data) {
  colnames(data) <- paste0(str_pad(1:ncol(data), 2, "left", "0"), "_", colnames(data))
  apply(data, 2, scale_min_max) %>%
    as.data.frame() %>%
    mutate(volume = row_number()) %>%
    gather(., variable, value, -volume) %>%
    ggplot(., aes(variable, volume, fill = value)) +
    geom_raster() + 
    scale_fill_distiller(palette = "Greys", direction = 1) + 
    theme_minimal() +
    theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
    labs(x = NULL)
}

func_fig_data_ts_long <- function(data) {
  colnames(data) <- paste0(str_pad(1:ncol(data), 2, "left", "0"), "_", colnames(data))
  data %>%
  mutate(volume = row_number()) %>%
  gather(., variable, value, -volume) %>%
  ggplot(., aes(volume, value)) +
  geom_line() +
  facet_wrap(~ variable, scales = "free_y", ncol = 1) +
  theme_minimal()
}

func_fig_data_heat_map(df_design_mat)

func_fig_data_ts_long(df_design_mat)

Add nuisance variables

Final

df_design_mat_final <- df_design_mat %>%
  group_by(run) %>%
  nest() %>%
  mutate(data = future_map(data, function(x) x %>% mutate(time_linear = row_number(),
                                                          time_linear = scale(time_linear),
                                                          time_quadratic = time_linear^2))) %>%
  unnest() %>%
  ungroup()

func_fig_data_heat_map(df_design_mat_final)

func_fig_data_ts_long(df_design_mat_final)

Save

if (!is.null(out_file)) {
  write.csv(df_design_mat_final, out_file)
}

epongpipat/ppi documentation built on Jan. 31, 2024, 1:02 p.m.