7 Run RSA Searchlight
7.1 Prepare functional data
We want to run RSA on the multi-voxel patterns corresponding to the object presentations during the picture viewing tasks. For this, we use the data as it was preprocessed by Lorena Deuker. We will do some further cleaning and extract the relevant volumes from the time series to eventually run RSA searchlights. What will be done in the following sections has been done previously for the ROI data only. Now we do it for the whole field of view in preparation of the RSA searchlights.
Data cleaning: Extract residuals from motion parameter regression
In this first section of the script we will create a dataframe that includes the relevant file names for the files that are needed to calculate the clean time series. These files include the motion parameters from FSL FEAT, which was run on the functional data from the picture viewing task. As described above, these data were split into 10 blocks. We will use the preprocessed FEAT output images which were already co-registered to the analysis space (‘samespace’). Further, we will use a mask to only include voxels with data in all blocks. This mask has already been created and is available in the analysis space.
#PREPARE DATA FRAME FOR CLEANING
# create a tibble with filenames etc for data cleaning, start with subject IDs
func_dat_df <- tibble(subject = rep(as.numeric(subjects), each = n_runs * n_blocks))
# add run info (run 1 and 2 equal the pre and post PVT)
func_dat_df <- add_column(func_dat_df, run = rep(rep(1:n_runs, each = n_blocks), n_subs))
# add block information (blocks 1 to 10 are the PVT blocks preprocessed separately)
func_dat_df <- add_column(func_dat_df, block = rep(1:n_blocks, n_subs*n_runs))
# add the motion parameter file
func_dat_df <- add_column(func_dat_df, mc_par_fn =
file.path(dirs$feat_dir,
paste0("VIRTEM_P0", func_dat_df$subject),
sprintf("RSA_%02d_Block%02d.feat",
func_dat_df$run, func_dat_df$block),
"mc", "prefiltered_func_data_mcf.par"))
# add the functional data file
func_dat_df <- add_column(func_dat_df, filtered_func =
file.path(dirs$samespace_dir,
paste0("VIRTEM_P0", func_dat_df$subject),
sprintf("VIRTEM_P%03d_RSA_%02d_Block%02d_func.nii.gz",
func_dat_df$subject, func_dat_df$run, func_dat_df$block)))
# add the brain mask data file
func_dat_df <- add_column(func_dat_df, brain_mask =
file.path(dirs$samespace_dir,
paste0("VIRTEM_P0", func_dat_df$subject),
sprintf("VIRTEM_P%03d_common_mask_perBlock.nii.gz",
func_dat_df$subject)))
# add output directory
func_dat_df <- add_column(func_dat_df, out_dir = file.path(dirs$searchlight, "clean_timeseries"))
if(!dir.exists(file.path(dirs$searchlight, "clean_timeseries"))){
dir.create(file.path(dirs$searchlight, "clean_timeseries"))
}Define function for data cleaning
In this next section, we will first define a function and then run it on all datasets. In this function, the preprocessed data will be loaded and transformed to matrix format. For every voxel within the brain mask, movement correction parameters are used as predictors in a GLM with the voxel time series as dependent variable. Finally, the residuals from this GLM (i.e. what could not be explained by motion) will be written to a text file.
# DEFINE THE FUNCTION TO GET THE RESIDUAL TIME SERIES
run_motion_glm <- function(df = func_dat_df[1,]){
# load the motion params and convert to tibble
mc_pars <- read.table(df$mc_par_fn, header = FALSE)
mp_df <- as_tibble(mc_pars)
colnames(mp_df) <- paste0("mp", 1:6)
# load the brain mask, check it's only zeros and ones, then convert to logical
brain_mask_nii <- readNIfTI2(df$brain_mask)
assertthat::assert_that(all.equal(unique(c(img_data(brain_mask_nii))), c(0,1)))
brain_mask <- array(NA, dim(brain_mask_nii))
brain_mask[] <- as.logical(img_data(brain_mask_nii))
# load the functional data
func_nii <- readNIfTI2(df$filtered_func)
# initialize output
clean_dat <- array(0, dim(func_nii))
counter <- 0
for (i_dim1 in 1:dim(func_nii)[1]){
for (i_dim2 in 1:dim(func_nii)[2]){
for (i_dim3 in 1:dim(func_nii)[3]){
if (brain_mask[i_dim1,i_dim2,i_dim3]){
# extract this voxel's data and merge with motion params
mp_df$vox_dat <- func_nii[i_dim1,i_dim2,i_dim3,]
# run the glm and store the residuals
clean_dat[i_dim1,i_dim2,i_dim3,] <- resid(glm('vox_dat~mp1+mp2+mp3+mp4+mp5+mp6', data = mp_df))
}
# print a message to show progress
counter <-counter+1
if (counter%%(prod(dim(func_nii)[1:3])/100) == 0) {
print(paste0((counter/prod(dim(func_nii)[1:3])*100), "% done"))}
}
}
}
# create nifti with same header info as original data from clean time series
clean_dat_nii = copyNIfTIHeader(img = func_nii, arr = clean_dat)
# create output folder
out_dir_sub <- file.path(df$out_dir, sprintf("%03d",df$subject))
if(!dir.exists(out_dir_sub)){dir.create(out_dir_sub, recursive = TRUE)}
# write cleaned nifti to file
fn <- sprintf("%03d_run%02d_block%02d_clean_func_data.nii.gz",
df$subject, df$run, df$block)
writenii(nim = clean_dat_nii, filename = file.path(out_dir_sub,fn))
}Run data cleaning
Now actually run the function either in parallel or serially.
# next step depends on whether we are in parallel or serial mode
if (!run_parallel){ # run serially
print("Attempting to run motion parameter cleaning serially. This will take a very long time if runnning for all subjects")
# run the function for each row of the data frame,
# i.e. for each block in each run for each subject
for(i in 1:nrow(func_dat_df)){
tic(sprintf("Subject %s, run %d, block %d",
func_dat_df$subject[i], func_dat_df$run[i], func_dat_df$block[i]))
run_motion_glm(df = func_dat_df[i,])
toc()
}
} else if (run_parallel){ # run in parallel, assumes CBS HTCondor is available
# write the data frame to file
fn <- file.path(here("data","mri", "rsa", "searchlight", "clean_timeseries"),
"htc_config_clean_time_series.txt")
fn_def <- cat(sprintf('"fn <- "%s"',fn))
write.table(func_dat_df, fn,
append = FALSE, sep = ",", dec = ".", row.names = FALSE, col.names = TRUE)
# store the function definition as text
func_def <- capture.output(run_motion_glm)
func_def[1] <- paste0("run_motion_glm <- ",func_def[1])
#func_def <- func_def[-length(func_def)]
#write the Rscript that we want to run
rscript_fn <- here("data","mri", "rsa", "searchlight", "clean_timeseries", "run_clean_timeseries.R")
con <- file(rscript_fn)
open(con, "w")
writeLines(c(
"\n# handle input",
"args = commandArgs()",
"i <- as.numeric(args[length(args)])",
"\n#load required packages",
noquote(sprintf('lib_dir <- "%s"',"/data/pt_02261/virtem/virtem_code/R3.6.1/library/Linux")),
'.libPaths(c(lib_dir,.libPaths()))',
'lapply(c("oro.nifti", "assertthat", "dplyr", "neurobase"), library, character.only = TRUE)',
"\n# read the data and transform ROI column back to list",
noquote(sprintf('fn <- "%s"',fn)),
'func_dat_df <- read.table(fn, sep = ",", header = TRUE, stringsAsFactors = FALSE)',
"\n#define the function to run motion GLM",
func_def,
"\n# run the function on one line of the data frame",
"run_motion_glm(df = func_dat_df[i,])"),con)
close(con)
# folder for condor output
htc_dir <- here("htc_logs", "clean_timeseries")
if(!exists(htc_dir)){dir.create(htc_dir, recursive = TRUE)}
# write the submit script
fn <- here("data","mri", "rsa", "searchlight", "clean_timeseries", "run_clean_timeseries.submit")
con <- file(fn)
open(con, "w")
writeLines(c(
"universe = vanilla",
"executable = /afs/cbs.mpg.de/software/scripts/envwrap",
"request_memory = 9000",
"notification = error"
),con)
for (i in 1:nrow(func_dat_df)){
writeLines(c(
sprintf("\narguments = R+ --version 3.6.1 Rscript %s %d", rscript_fn, i),
sprintf("log = %s/%d.log", htc_dir, i),
sprintf("output = %s/%d.out", htc_dir, i),
sprintf("error = %s/%d.err", htc_dir, i),
sprintf("Queue\n")),con)
}
close(con)
# submit to condor
#system("reset-memory-max")
batch_id <- system(paste("condor_submit", fn), intern = TRUE)
batch_id <- regmatches(batch_id[2], gregexpr("[[:digit:]]+", batch_id[2]))[[1]][2]
cat(sprintf("submitted jobs (ID = %s) to clean time series for searchlights. Time to wait...", batch_id))
pause_until_batch_done(batch_id = batch_id, wait_interval = 600)
#system("reset-memory-max")
#system("condor_q")
}Extract relevant volumes
As the presentation of images in the PVT pre and post blocks was locked to the onset of a new volume (see above), the second volume after image onset was selected for every trial (effectively covering the time between 2270-4540 ms after stimulus onset).
In this step, we split the presentation logfile from the picture viewing task into the 10 blocks. We reference the volume count to the first volume of each block and extract the relevant volumes from the cleaned 4D timeseries data, accounting for the temporal offset. This is done for each block separately in both the pre and post runs. The data are written to 3D-niftis.
offset_tr = 2
for (i_sub in subjects){
for (i_run in 1:n_runs){
# load the logfile for this run (pre or post)
log_fn <- file.path(dirs$pvt_log_dir, sprintf('P%s_%svirtem.txt', i_sub, runs[i_run]))
log <- read.table(log_fn)
colnames(log) <- c("pic", "fix_start", "pic_start", "volume", "response", "RT", "trial_end")
# split the log file into the 10 blocks
log_split <- split(log, rep(1:10, each = 21))
# reference volume numbers to the first volume of that block
vol_block <- lapply(log_split, function(x){x$volume - x$volume[1]})
log_split <- mapply(cbind, log_split, "vol_block"=vol_block, SIMPLIFY=FALSE)
for (i_block in 1:n_blocks){
# define the full 4D file
clean_func_fn <- file.path(dirs$searchlight, "clean_timeseries", i_sub,
sprintf("%s_run%02d_block%02d_clean_func_data.nii.gz",
i_sub, i_run, i_block))
# index of relevant volumes when accounting for offset
rel_vols <- log_split[[i_block]]$vol_block + offset_tr
# extract the relevant volumes from the ROI timeseries data
out_dir <- file.path(dirs$searchlight, "rel_vols_3D", i_sub)
if(!dir.exists(out_dir)){dir.create(out_dir, recursive = TRUE)}
for (i_vol in 1:length(rel_vols)){
# subtract 1 from volume number because of 0-based FSL indexing
fslroi(file = clean_func_fn, tmin = rel_vols[i_vol]-1, tsize = 1,
outfile = file.path(out_dir,
sprintf("%s_run%02d_block%02d_pic%02d.nii.gz",
i_sub, i_run, i_block, log_split[[i_block]]$pic[i_vol])),
verbose = FALSE)
}
}
}
}Concatenate relevant volumes
Next, we create 4D-files for both the pre and the post run. Each file will have consist of 200 (20 relevant pictures x 10 blocks) volumes. We take care to order the volumes according to picture identity as this is most convenient for later RSA.
for (i_sub in subjects){
for (i_run in 1:n_runs){
# files to concatenate
in_dir <- file.path(dirs$searchlight, "rel_vols_3D", i_sub)
files <- c(file.path(in_dir, sprintf("%s_run%02d_block%02d_pic%02d.nii.gz",
i_sub, i_run, c(rep(1:10,n_pics)), c(rep(1:n_pics, each = n_blocks)))))
# output file
out_dir <- file.path(dirs$searchlight, "rel_vols_4D")
if(!dir.exists(out_dir)){dir.create(out_dir, recursive = TRUE)}
fn <- file.path(out_dir, sprintf("%s_run%02d_rel_vols.nii.gz",
i_sub, i_run))
# merge the files
fslmerge(infiles = files, direction = "t", outfile = fn,
retimg = FALSE, verbose = FALSE)
# change datatype to short to save space and memory for the cluster jobs
fslmaths(file = fn, outfile = fn, retimg = FALSE,
opts = "-odt short", opts_after_outfile = TRUE, verbose = FALSE)
}
}7.2 First-level RSA searchlight
We further probed how temporal distances between events shaped representational change using searchlight analyses. Using the procedures described above, we calculated pattern similarity change values for search spheres with a radius of 3 voxels around the center voxel. Search spheres were centered on all brain voxels within our field of view. Within a given search sphere, only gray matter voxels were analyzed. Search spheres not containing more than 25 gray matter voxels were discarded.
Building on the prepared data, we want to run the searchlight analysis. To speed things up, we split the brain into chunks that we analyze in parallel jobs.
Set analysis parameters
First, let’s set some analysis parameters. These include
- the radius of the search spheres. This is the radius around the center voxel, so that the diameter of the spheres is given by 2*radius + 1
- the minimum number of valid voxels that a sphere has to contain to be analyzed
- the number of chunks into which we split the analysis to speed up the computations using the HPC cluster
Further, we define a function that returns the voxels in the sphere of a given radius around. This canonical sphere definition will then be combined with the actual sphere center coordinates. It is based on the way the sphere definition is implemented in the fmri package for R.
# how large should the searchlights be
radius <- 3 # radius around center voxel
# how many voxels (counting the center) must be in a searchlight for it to be run; smaller ones are skipped.
min_voxels <- 25
# how many chunks to split the analysis up into (for speed reasons)
n_chunks <- 40
# function to get searchsphere of given radius (based on: https://rdrr.io/cran/fmri/src/R/searchlight.r)
searchlight <- function(radius){
rad <- as.integer(radius)
nr <- 2*rad+1
indices <- rbind(rep(-rad:rad,nr*nr),
rep(-rad:rad,rep(nr,nr)),
rep(-rad:rad,rep(nr*nr,nr)))
indices[,apply(indices^2,2,sum)<=radius^2]
}Prepare RSA-predictions
To relate pattern similarity change in each search sphere to the temporal structure of the sequences we need the relationships of the events in virtual time. We calculate these here analogously to the main ROI analyses, but save them for each subject separately.
for (i_sub in subjects){
# set up the dataframe to use for RSA
pics <- which(upper.tri(matrix(TRUE, n_pics, n_pics), diag = FALSE), arr.ind=TRUE)
# load the behavioral data
fn <- file.path(dirs$timeline_dat_dir, sprintf("%s_behavior_tbl_timeline.txt", i_sub))
col_types_list <- cols_only(sub_id = col_character(), day = col_factor(),
event = col_integer(), pic = col_integer(),
virtual_time = col_double(), real_time = col_double(),
memory_time = col_double(), memory_order = col_double(),
sorted_day = col_integer())
beh_dat <- read_csv(fn, col_types = col_types_list)
# sort behavioral data according to picture identity
beh_dat_ordered <- beh_dat[order(beh_dat$pic),]
# find the order of comparisons
pics <- which(upper.tri(matrix(nrow = 20, ncol = 20)), arr.ind=TRUE)
# extract the data for the first and second picture in each pair from the behavioral data
pic1_dat <- beh_dat_ordered[pics[,1],]
pic2_dat <- beh_dat_ordered[pics[,2],]
colnames(pic1_dat) <- paste0(colnames(pic1_dat), "1")
colnames(pic2_dat) <- paste0(colnames(pic2_dat), "2")
# combine the two tibbles
pair_dat <- cbind(pic1_dat, pic2_dat)
# reorder the tibble columns to make a bit more sense
pair_dat <- pair_dat %>%
select(sub_id1,
day1, day2,
pic1, pic2,
event1, event2,
virtual_time1, virtual_time2,
real_time1, real_time2,
memory_order1, memory_order2,
memory_time1, memory_time2,
sorted_day1, sorted_day2) %>%
rename(sub_id = sub_id1)
# prepare RSA distance measures
pair_dat <- pair_dat %>%
mutate(
# pair of events from say or different day (dv = deviation code)
same_day = day1 == day2,
same_day_dv = plyr::mapvalues(same_day, from = c(FALSE, TRUE), to = c(-1, 1)),
# absolute difference in time metrics
vir_time_diff = abs(virtual_time1 - virtual_time2),
order_diff = abs(event1 - event2),
real_time_diff = abs(real_time1 - real_time2)) %>%
# z-score the time metric predictors
mutate_at(
c("vir_time_diff", "order_diff", "real_time_diff"), scale)
# write to file
out_dir <- file.path(dirs$searchlight, "rsa_predictions")
if(!dir.exists(out_dir)){dir.create(out_dir, recursive = TRUE)}
fn <- file.path(out_dir, sprintf("%s_rsa_info.txt", i_sub))
write.table(pair_dat, file = fn)
}Create lookup table for searchlight
As a first step, we prepare a table that - for each individual sphere - holds the coordinates of the voxels making up the search spheres. To find these voxels we rely on both a brain mask and a gray matter mask obtained during preprocessing. The brain mask is combined across the 10 preprocessing blocks and all voxels within this brain mask can be sphere centers. We get the coordinates of voxels around these sphere centers given the desired radius. Then, we exclude the non-gray matter voxels from the spheres. We check which spheres still contain the minimum number of voxels (do this here to speed up the subsequent computations) and discard the spheres that don’t hold enough voxels, e.g. because they are in white matter or on the edge of the brain/field of view.
We split this table into chunks to speed up the calculations and save the chunked lookup-tables for later use. To be able to relate the rows of the table (each row is a sphere) back to the brain, we create a nifti image that numbers all voxels which are sphere centers. These numbers correspond to the row names of the lookup table and will later be used to map the RSA results of each sphere back to the brain.
# Set up dataframe
in_df <- tibble(subject = subjects)
# how many chunks to run the analysis in for each
in_df$n_chunks <- n_chunks
# input files with relevant functional volumes
in_df$rel_vol_pre_fn <- file.path(dirs$searchlight, "rel_vols_4D",
sprintf("%s_run%02d_rel_vols.nii.gz", subjects, 1))
in_df$rel_vol_post_fn <- file.path(dirs$searchlight, "rel_vols_4D",
sprintf("%s_run%02d_rel_vols.nii.gz", subjects, 2))
# which voxels should be center voxels of spheres --> typically brain mask
in_df$ctr_mask_fn <- file.path(dirs$samespace_dir, paste0("VIRTEM_P", subjects),
sprintf("VIRTEM_P%s_common_mask_perBlock.nii.gz",
subjects))
# which voxels should be used for correlations --> typically graymatter mask
in_df$feat_mask_fn <- file.path(dirs$samespace_dir, paste0("VIRTEM_P", subjects),
sprintf("VIRTEM_P%s_common_graymatter_mask_brainmasked.nii.gz",
subjects))
# radius of spheres
in_df$radius <- radius
# where to store outputs
in_df$out_path <- file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius), subjects)
in_df$chunk_path <- file.path(in_df$out_path, "chunks")
prepare_searchlight_lut <- function(in_df = in_df[1,]){
if(!dir.exists(in_df$out_path)){dir.create(in_df$out_path, recursive = TRUE)}
if(!dir.exists(in_df$chunk_path)){dir.create(in_df$chunk_path, recursive = TRUE)}
# load the mask with voxels that will be sphere centers (typically brain mask)
ctr_mask_nii <- readNIfTI(in_df$ctr_mask_fn, reorient=FALSE)
if (!all.equal(unique(c(img_data(ctr_mask_nii))), c(0,1))){
stop("center mask not binary!")}
# load the mask with voxels that will be the features used for correlations
feat_mask_nii <- readNIfTI(in_df$feat_mask_fn, reorient=FALSE)
feat_mask_array <- array(0, dim = dim(feat_mask_nii))
feat_mask_array[feat_mask_nii > 0.7 ] <- 1
if (!all.equal(unique(c(feat_mask_array)), c(0,1))){
stop("make sure feature mask is binary!")}
# make sure mask dimensions match functional images and each other
if (!all.equal(dim(ctr_mask_nii), dim(feat_mask_array))){
stop("mask dimensions don't match each other!")}
# find all the voxel coordinates that will be sphere centers
ctr_inds <- which(ctr_mask_nii != 0, arr.ind=TRUE) # inds has subscripts for 3D coords of center voxels
if (ncol(ctr_inds) != 3) { stop('wrong dims on inds')}
# 3D array with voxel numbers to save for future reference
vox_num_array <- array(0, dim(ctr_mask_nii)) # make a blank 3d matrix 'brain', all zeros
vox_num_array[ctr_inds] <- 1:nrow(ctr_inds) # and put integers into the voxel places
# create nifti with same header info as input data
vox_num_nii <- copyNIfTIHeader(img = ctr_mask_nii, arr = vox_num_array)
write_nifti(vox_num_nii, file.path(in_df$out_path, "vox_num")) # write out as a NIfTI image
# get the search sphere 3D subscripts
sphere_coords <- searchlight(radius)
# actually fill up the lookup table. This goes through every voxel, so can take some time.
# the table will hold the sphere indices for each center voxel
lookup_tbl <- array(NA, c(nrow(ctr_inds), ncol(sphere_coords)))
rownames(lookup_tbl) <- vox_num_nii[ctr_inds]
for (i in 1:nrow(ctr_inds)) { # i <- 1
# add the sphere coordinates to the current voxel coordinates
curr_sphere <- t(ctr_inds[i,] + sphere_coords)
# remove voxels that are out of bounds (i.e. out of image dimensions)
dims <- dim(vox_num_array)
curr_sphere[which(curr_sphere[,1]>dims[1]),] <- NA
curr_sphere[which(curr_sphere[,2]>dims[2]),] <- NA
curr_sphere[which(curr_sphere[,3]>dims[3]),] <- NA
curr_sphere[which(curr_sphere[,1]<1),] <- NA
curr_sphere[which(curr_sphere[,2]<1),] <- NA
curr_sphere[which(curr_sphere[,3]<1),] <- NA
# remove voxels that are not in the feature mask
curr_sphere[!feat_mask_array[curr_sphere],] <- NA
# store the voxel numbers for this sphere
vox.id <- vox_num_array[ctr_inds[i,1], ctr_inds[i,2], ctr_inds[i,3]]
lookup_tbl[vox.id,] <- vox_num_array[curr_sphere]
}
# replace the zeroes with NA so can sort better at runtime
to_remove <- which(lookup_tbl == 0, arr.ind=TRUE)
lookup_tbl[to_remove] <- NA
# remove rows that don't have the minimum amount of voxels (do this here because this speeds up all later computations)
lookup_tbl <- lookup_tbl[rowSums(!is.na(lookup_tbl)) > min_voxels,]
# write the table, gzipped to save file size
write.table(lookup_tbl, gzfile(file.path(in_df$out_path, sprintf("lookup_radius%d", in_df$radius, ".txt.gz"))), row.names = TRUE)
# vector to split lookupt table into chunks (last chunk will be slightly smaller)
r <- rep(1:in_df$n_chunks,each=ceiling(nrow(lookup_tbl)/in_df$n_chunks))[1: nrow(lookup_tbl)]
for (i_chunk in 1:in_df$n_chunks){
# extract current chunk and write to file
curr_chunk <- lookup_tbl[r == i_chunk,]
write.table(curr_chunk,
gzfile(file.path(in_df$chunk_path,
sprintf("lookup_radius%d_chunk_%02d.txt.gz", in_df$radius, i_chunk))),
row.names = TRUE)
}
}
# run for all datasets
for (i in 1:nrow(in_df)){
prepare_searchlight_lut(in_df = in_df[i,])
}We want to create some diagnostic images of the spheres we created to make sure the resulting search spheres look as expected.
For this a random subject is picked and a number of random spheres is plotted on the gray matter mask.
# pick random subject
rand_sub <- sample(1:n_subs)[1]
# read their lookup-table
lookup_tbl<-read.table(gzfile(
file.path(in_df$out_path[rand_sub],
sprintf("lookup_radius%d", in_df$radius[rand_sub], ".txt.gz"))))
# load their graymatter mask and the nifti linking indices
gm_mask_nii <- in_df[rand_sub,]$feat_mask_fn %>% readNIfTI2()
vox_num_nii <- file.path(in_df$out_path[rand_sub], "vox_num.nii.gz") %>% readNIfTI2()
# pick 10 searchlights to plot at random
do.searchlights <- sample(1:nrow(lookup_tbl))[1:10]
# make a blank 3d matrix 'brain' to put the searchlights into
sphere_test_array <- array(0, dim(gm_mask_nii))
for (i in 1:length(do.searchlights)) { # i <- 1
# voxels of this example sphere
voxs <- unlist(lookup_tbl[do.searchlights[i],], use.names=FALSE)
voxs <- voxs[which(!is.na(voxs))] # get rid of NAs
print(sprintf('sphere %02d: %d voxels', i, length(voxs)))
# put integers into the sphere test array for each searchlight
for (j in 1:length(voxs)) { # j <- 1
coords <- which(vox_num_nii == voxs[j], arr.ind=TRUE) # location of this voxel
if (ncol(coords) != 3 | nrow(coords) != 1) { stop("wrong sized coords")}
sphere_test_array[coords[1], coords[2], coords[3]] <- i # assign this voxel the searchlight number
}
# visualize the sphere at the most frequent coordinate
ortho2(gm_mask_nii, y=sphere_test_array==i, col.y = "red", col.crosshairs = "lightgrey",
xyz = c(statip::mfv1(which(sphere_test_array==i, arr.ind = TRUE)[,1]),
statip::mfv1(which(sphere_test_array==i, arr.ind = TRUE)[,2]),
statip::mfv1(which(sphere_test_array==i, arr.ind = TRUE)[,3])))
}
# create nifti with same header info as input data and save
sphere_test_nii = copyNIfTIHeader(img = gm_mask_nii, arr = sphere_test_array)
write_nifti(sphere_test_nii,
file.path(dirs$searchlight,
sprintf("first_level_%dvox_radius", radius),
sprintf("sphere_test_sub%s", rand_sub))) Dataframe with input for analysis
We start by preparing a data frame with one line for each chunk of each subject. The entries define the files to be used for this chunk and some basic analysis parameters. The function to be defined subsequently takes one row of this data frame as input.
in_df <- tibble(subject = rep(subjects, each=n_chunks))
in_df$chunk <- rep(1:n_chunks, length(subjects))
in_df$lut_file <- file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius), rep(subjects, each=n_chunks),
"chunks", sprintf("lookup_radius%d_chunk_%02d.txt.gz",
radius, rep(1:n_chunks, length(subjects))))
in_df$pred_file <- file.path(dirs$searchlight, "rsa_predictions",
sprintf("%s_rsa_info.txt", rep(subjects, each=n_chunks)))
# input files with relevant functional volumes
in_df$rel_vol_pre_fn <- rep(file.path(dirs$searchlight, "rel_vols_4D",
sprintf("%s_run%02d_rel_vols.nii.gz",
subjects, 1)), each = n_chunks)
in_df$rel_vol_post_fn <- rep(file.path(dirs$searchlight, "rel_vols_4D",
sprintf("%s_run%02d_rel_vols.nii.gz",
subjects, 2)), each = n_chunks)
# file to use to bring data back into brain shape
in_df$vox_num_fn <- rep(file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius), subjects, "vox_num.nii.gz"), each = n_chunks)
# output directory
in_df$out_path <- rep(file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius), subjects, "chunks"), each = n_chunks)
# minimum number of features
in_df$min_voxels <- min_voxels
head(in_df)Searchlight implementation
Here, we define the function that implements RSA within each searchlight. The logic of running the searchlight analysis via the look-up table is inspired by and partly based on this blogpost by Joset A. Etzel and the code accompanying it.
For each search sphere, we implemented linear models to quantify the relationship between representational change and the learned temporal structure. Specifically, we assessed the relationship of pattern similarity change and absolute virtual temporal distances, separately for event pairs from the same sequences and from pairs from different sequences. In a third model, we included all event pairs and tested for an interaction effect of a sequence (same or different) predictor and virtual temporal distances. The t-values of the respective regressors of interest were stored at the center voxel of a given search sphere.
The function goes through the same steps as the ROI-based analysis. First, the data of the voxels in each sphere is extracted. The resulting multi-voxel patterns are correlated between picture presentations to yield a trial-by-trial (200-by-200) correlation matrix. For each pair of images, the trial-wise correlations (10-by-10) are averaged such that comparisons of trials from the same block are excluded (i.e. excluding the diagonal of the 10-by-10 squares). This results in the condition-by-condition similarity matrix (20-by-20), which is then subjected to further analysis based on the temporal structure of events. Specifically, linear models are implemented for the different analyses. The resulting t-values are stored at the location of the center voxels of the respective sphere. Thus, for each model we test, we end up with a nifti image that has t-values at the voxel locations of the current chunk.
# FUNCTION DEFINITION
run_searchlight <- function(in_df){ #in_df <- in_df[1,]
########## SET UP A FEW PARAMS
sprintf("working on: %s", in_df$lut_file)
n_pics <- 20
n_blocks <- 10
# for all 10x10 comparisons we will be averaging all comparisons apart from the diagonal
# to exclude same_block comparisons
no_diag <- matrix(data = TRUE, nrow=n_blocks, ncol=n_blocks)
diag(no_diag)<- FALSE
# indices to extract upper triangle excluding diagonal from pic-by-pic matrices
triu_idx <- which(upper.tri(matrix(TRUE, n_pics, n_pics), diag = FALSE), arr.ind = TRUE)
########## LOAD THE DATA TO BE USED
# load the data about the learned temporal relationships
rsa_info <- read.table(in_df$pred_file)
# initialize dataframe to fill (later ps-change data will be added)
pics <- which(upper.tri(matrix(TRUE, n_pics, n_pics), diag = FALSE), arr.ind=TRUE)
rsa_df <- tibble(pic1 = pics[,1], pic2 = pics[,2])
# make sure the data frames have the same order
if(!(all.equal(rsa_df$pic1, rsa_info$pic1) & all.equal(rsa_df$pic2, rsa_info$pic2))){
stop('order of data frames does not match [RSA preparations]')}
# store the columns we need for RSA in the data frame
rsa_df <- cbind(rsa_df, rsa_info[names(rsa_info) %in% c("same_day", "same_day_dv",
"vir_time_diff", "order_diff",
"real_time_diff")])
# read the look up table
lookup_tbl <- read.table(gzfile(in_df$lut_file), row.names = 1)
# figure out which voxels to run in this chunk
n_voxels <- nrow(lookup_tbl)
# read the functional images (takes ~ 2 mins)
print("starting to load functional images")
func_nii_pre <- readNIfTI(in_df$rel_vol_pre_fn, reorient=FALSE)
func_nii_post <- readNIfTI(in_df$rel_vol_post_fn, reorient=FALSE)
if(!all.equal(dim(func_nii_pre),dim(func_nii_post))){
stop("functional images of different size!")}# make sure inputs have same dimensions
# load the image with voxel numbers (linear subscripts) that will be used to sort back into brain shape
vox_num_nii <- readNIfTI(in_df$vox_num_fn, reorient=FALSE)
if (!all.equal(dim(func_nii_pre)[1:3], dim(vox_num_nii))){
stop("voxel number image dimensions don't match functional data!")}
print("finished loading images")
########## BEGIN THE ACTUAL ANALYSIS
# initialize the output as all 0 (because FSL does not like NAs)
rsa_out_array_same_day_vir_time <- array(0, dim(vox_num_nii))
rsa_out_array_diff_day_vir_time <- array(0, dim(vox_num_nii))
rsa_out_array_all_pairs_vir_time <- array(0, dim(vox_num_nii))
rsa_out_array_interaction <- array(0, dim(vox_num_nii))
#rsa_out_array_n_in_sphere <- array(0, dim(vox_num_nii))
# loop over all voxels (in this chunk)
for (v in 1:n_voxels) { # v <- do.centers[500]
# print a message to show progress
if (v%%100 == 0) { print(paste("at", v, "of", n_voxels))}
# find which voxels belong in this center voxel's searchlight
voxs <- unlist(lookup_tbl[v,], use.names=FALSE)
voxs <- voxs[which(!is.na(voxs))]; # get rid of NAs. There will be NA entries if some surrounding voxels not in gray matter or brain
# how many surrounding voxels must be in the searchlight? Smaller ones (edge of brain) will be skipped.
if (length(voxs) > in_df$min_voxels) {
# initialize arrays to store data of current sphere for the pre and post runs
# images in the rows (voxels in this searchlight only), voxels in the columns
curr_dat_pre <- array(NA, c(n_pics*n_blocks, length(voxs)))
curr_dat_post <- array(NA, c(n_pics*n_blocks, length(voxs)))
# put the data into a matrix to run analysis
for (i in 1:length(voxs)) { # i <- 1
# for the current voxel, take the data from all conditions and store it (separately for pre and post)
coords <- which(vox_num_nii == voxs[i], arr.ind=TRUE)
if (ncol(coords) != 3 | nrow(coords) != 1) { stop("wrong sized coords")}
curr_vox <- func_nii_pre[coords[1], coords[2], coords[3],] # pre
if (sd(curr_vox) > 0) {curr_dat_pre[,i] <- curr_vox} else { stop("zero variance voxel")}
curr_vox <- func_nii_post[coords[1], coords[2], coords[3],] # post
if (sd(curr_vox) > 0) {curr_dat_post[,i] <- curr_vox} else { stop("zero variance voxel")}
}
# data is in repetition (row) by voxel (col) format, so we transpose
# to get a voxel x repetition format
curr_dat_pre <- t(curr_dat_pre)
curr_dat_post <- t(curr_dat_post)
# calculate correlation matrix (trial by trial) for pre and post run
cor_mat_pre_trial <- cor(curr_dat_pre, curr_dat_pre)
cor_mat_post_trial <- cor(curr_dat_post, curr_dat_post)
# initialize condition by condition correlation matrix for pre and post run
corr_mat_pre <- matrix(nrow = 20, ncol = 20)
corr_mat_post <- matrix(nrow = 20, ncol = 20)
# loop over all picture comparisons
for(i_pic1 in 1:20){
for(i_pic2 in 1:20){
# extract the current 10x10 correlation matrix
i1 <- (1+(i_pic1-1)*10):(i_pic1*10)
i2 <- (1+(i_pic2-1)*10):(i_pic2*10)
curr_mat_pre <- cor_mat_pre_trial[i1, i2]
curr_mat_post <- cor_mat_post_trial[i1, i2]
# average the correlations while excluding diagonal (same block comparisons)
corr_mat_pre[i_pic1, i_pic2] <- mean(curr_mat_pre[no_diag])
corr_mat_post[i_pic1, i_pic2] <- mean(curr_mat_post[no_diag])
}
}
# calculate pattern similarity change based on FisherZ-transformed upper
# triangles of the correlation matrices
rsa_df$ps_change <- FisherZ(corr_mat_post[triu_idx]) - FisherZ(corr_mat_pre[triu_idx])
# find 3D-coordinates of this searchlight center to save output
coords <- which(vox_num_nii == as.numeric(rownames(lookup_tbl[v,])), arr.ind=TRUE)
if (ncol(coords) != 3 | nrow(coords) != 1) { stop("wrong sized coords"); }
########## RUN RSA FOR THIS SPHERE
# same day virtual time
fit <- lm(ps_change ~ vir_time_diff, rsa_df[rsa_df$same_day,])
rsa_out_array_same_day_vir_time[coords[1], coords[2], coords[3]] <- coef(summary(fit))[, "t value"][2]
# different day virtual time
fit <- lm(ps_change ~ vir_time_diff, rsa_df[!rsa_df$same_day,])
rsa_out_array_diff_day_vir_time[coords[1], coords[2], coords[3]] <- coef(summary(fit))[, "t value"][2]
# all pairs virtual time
fit <- lm(ps_change ~ vir_time_diff, rsa_df)
rsa_out_array_all_pairs_vir_time[coords[1], coords[2], coords[3]] <- coef(summary(fit))[, "t value"][2]
# day*time interaction
fit <- lm(ps_change ~ vir_time_diff*same_day_dv, rsa_df)
rsa_out_array_interaction[coords[1], coords[2], coords[3]] <- coef(summary(fit))[, "t value"][4]
# number of features
#rsa_out_array_n_in_sphere[coords[1], coords[2], coords[3]] <- length(voxs)
} else{
stop("number of features too small!")
}
}
########## SAVE RESULTS
# create nifti with same header info as input data and write to file
# same day virtual time
rsa_out_nii = copyNIfTIHeader(img = vox_num_nii, arr = rsa_out_array_same_day_vir_time)
write_nifti(rsa_out_nii, file.path(in_df$out_path,
sprintf("%s_searchlight_same_day_vir_time_chunk%02d",
in_df$subject, in_df$chunk)))
# different day virtual time
rsa_out_nii = copyNIfTIHeader(img = vox_num_nii, arr = rsa_out_array_diff_day_vir_time)
write_nifti(rsa_out_nii, file.path(in_df$out_path,
sprintf("%s_searchlight_diff_day_vir_time_chunk%02d",
in_df$subject, in_df$chunk)))
# all pairs virtual time
rsa_out_nii = copyNIfTIHeader(img = vox_num_nii, arr = rsa_out_array_all_pairs_vir_time)
write_nifti(rsa_out_nii, file.path(in_df$out_path,
sprintf("%s_searchlight_all_pairs_vir_time_chunk%02d",
in_df$subject, in_df$chunk)))
# day*time interaction
rsa_out_nii = copyNIfTIHeader(img = vox_num_nii, arr = rsa_out_array_interaction)
write_nifti(rsa_out_nii, file.path(in_df$out_path,
sprintf("%s_searchlight_day_time_interaction_chunk%02d",
in_df$subject, in_df$chunk)))
# no. voxels in sphere
#rsa_out_nii = copyNIfTIHeader(img = vox_num_nii, arr = rsa_out_array_n_in_sphere)
#write_nifti(rsa_out_nii, file.path(in_df$out_path, sprintf("%s_searchlight_n_vox_sphere_chunk%02d", in_df$subject, in_df$chunk)))
}Now actually run the function for each chunk; either in parallel or serially. Because the analysis is slow (several hours per chunk when using 50 chunks), it is not feasible to run it serially.
# next step depends on whether we are in parallel or serial mode
if (!run_parallel){ # run serially
print("Running the searchlight analysis serially. Do this only for testing purposes because the code will run forever")
# run the function for each row of the data frame,
# i.e. for each block in each run for each subject
for(i in 1:nrow(in_df)){
tic(sprintf("Subject %s, chunk %d",
in_df$subject[i], in_df$chunk[i]))
run_searchlight(in_df = in_df[i,])
toc()
}
} else if (run_parallel){ # run in parallel, assumes CBS HTCondor is available
# write the data frame to file
fn <- file.path(here("data","mri", "rsa", "searchlight", sprintf("first_level_%dvox_radius", radius)),
"htc_config_run_searchlight.txt")
fn_def <- cat(sprintf('"fn <- "%s"',fn))
write.table(in_df, fn,
append = FALSE, sep = ",", dec = ".", row.names = FALSE, col.names = TRUE)
# store the function definition as text
func_def <- capture.output(run_searchlight)
func_def[1] <- paste0("run_searchlight <- ",func_def[1])
#func_def <- func_def[-length(func_def)]
#write the Rscript that we want to run
rscript_fn <- here("data","mri", "rsa", "searchlight", sprintf("first_level_%dvox_radius", radius), "run_searchlight.R")
con <- file(rscript_fn)
open(con, "w")
writeLines(c(
"\n# handle input",
"args = commandArgs()",
"i <- as.numeric(args[length(args)])",
"\n#load required packages",
noquote(sprintf('lib_dir <- "%s"',"/data/pt_02261/virtem/virtem_code/R3.6.1/library/Linux")),
'.libPaths(c(lib_dir,.libPaths()))',
'lapply(c("oro.nifti", "assertthat", "dplyr", "neurobase", "DescTools"), library, character.only = TRUE)',
"\n# read the data and transform ROI column back to list",
noquote(sprintf('fn <- "%s"',fn)),
'func_dat_df <- read.table(fn, sep = ",", header = TRUE, stringsAsFactors = FALSE)',
"\n#define the function to run motion GLM",
func_def,
"\n# run the function on one line of the data frame",
"run_searchlight(in_df = func_dat_df[i,])"),con)
close(con)
# folder for condor output
htc_dir <- here("htc_logs", "first_level_searchlight")
if(!dir.exists(htc_dir)){dir.create(htc_dir, recursive = TRUE)}
# write the submit script
fn <- here("data","mri", "rsa", "searchlight", sprintf("first_level_%dvox_radius", radius), "run_searchlight.submit")
con <- file(fn)
open(con, "w")
writeLines(c(
"universe = vanilla",
"executable = /afs/cbs.mpg.de/software/scripts/envwrap",
"request_memory = 13500", # "request_memory = 13500", # works for 3 voxel radius & 40 chunks
"notification = error"
),con)
c <- 0
for (i in 1:nrow(in_df)){
if (!file.exists(file.path(in_df[i,]$out_path,
sprintf("%d_searchlight_all_pairs_vir_time_chunk%02d.nii.gz",
as.numeric(in_df[i,]$subject), in_df[i,]$chunk)))){
c <- c + 1
writeLines(c(
sprintf("\narguments = R+ --version 3.6.1 Rscript %s %d", rscript_fn, i),
sprintf("log = %s/%d.log", htc_dir, i),
sprintf("output = %s/%d.out", htc_dir, i),
sprintf("error = %s/%d.err", htc_dir, i),
sprintf("Queue\n")),con)
}
}
close(con)
# submit to condor and play the waiting game
batch_id <- system(paste("condor_submit", fn), intern = TRUE)
batch_id <- regmatches(batch_id[2], gregexpr("[[:digit:]]+", batch_id[2]))[[1]][2]
cat(sprintf("submitted jobs (ID = %s) for 1st-level searchlights. Time to wait...", batch_id))
pause_until_batch_done(batch_id = batch_id, wait_interval = 1800) # check every half hour if done
#system("condor_q")
}Combine the chunks
Lastly, we need to combine the results across chunks to get the whole picture.
# analyses that we ran the searchlights for
searchlights <- c("same_day_vir_time",
"diff_day_vir_time",
"day_time_interaction",
"all_pairs_vir_time"#,
#"n_vox_sphere"
)
for (i_sub in subjects){
# where are the chinks
in_dir <- file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius), i_sub, "chunks")
for (i_srchlght in searchlights){
# file names of the individual chunks
in_files <- file.path(in_dir, sprintf("%d_searchlight_%s_chunk%02d.nii.gz",
as.numeric(i_sub), i_srchlght, 1:n_chunks))
# load the first chunk image, this will serve as the basis for combining all chunks
comb_nii <- readNIfTI(in_files[1], reorient = FALSE)
for (i_chunk in 2:n_chunks){
# load this chunk's nifti
new_nii <- readNIfTI(in_files[i_chunk], reorient = FALSE)
# find where the data is (i.e. where the image is non-zero)
#new_nii[is.na(new_nii)] <- 0
coords <- which(new_nii !=0, arr.ind = TRUE)
# the combined image should be 0 at all these coordinates
if(!(sum(comb_nii[coords]==0) == nrow(coords))){
stop("chunks overlap!")}
# add this chunk's data to the combined nifti
comb_nii[coords] <- new_nii[coords]
}
# make a simple plot to check that the output is shaped like a brain
fn <- file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius),
i_sub, sprintf("%s_%s", i_sub, i_srchlght))
jpeg(file = paste0(fn,".jpg"), width = 1000, height = 1000, units = "px")
ortho2(comb_nii, NA.x=TRUE)
dev.off()
# write nifti to file
write_nifti(comb_nii, fn)
}
}7.3 Group-level searchlight analysis
After running the first level analysis, we want to run group-level statistics. Group-level stats will be run using permutation-based one sample t-tests as implemented by the sign-flipping procedure in FSL Randomise. This test will be run for all voxels in our field of view as well as in a mask comprising the aHPC and alEC to correct for multiple comparisons within our a priori regions of interest.
First, let’s define some parameters and folders for these analyses.
# searchlight radius (used in file names)
radius <- 3
# analyses that we ran the searchlights for
searchlights <- c("same_day_vir_time",
"diff_day_vir_time",
"day_time_interaction",
"all_pairs_vir_time"
)
# what smoothing to apply
FWHM = 3 # in mm
sigma = FWHM/2.3548 # fslmaths needs sigma not FWHM of kernel
# Output folder
invisible(lapply(file.path(dirs$data4analysis, "searchlight_results",
searchlights),
function(x) if(!dir.exists(x)) dir.create(x, recursive=TRUE)))
# folder for condor output
htc_dir <- here("htc_logs", "srchlght")
if(!exists(htc_dir)){dir.create(htc_dir, recursive = TRUE)}Move first-level images to MNI space
The resulting t-maps were registered to MNI space for group level statistics and spatially smoothed (FWHM 3mm).
The first level searchlight output is in each participant’s common functional space. For group-level stats we need to register the searchlight outputs to MNI (1mm) space.
# name of transformation matrix file to move from the shared functional space to MNI 1mm
func2standard <- here("data", "mri", "processed", "wholebrain",
paste0("VIRTEM_P", rep(subjects, each = length(searchlights)), ".feat"),
"reg", "example_func2standard.mat")
# searchlight result in whole-brain functional space
in_nii_fn <- file.path(dirs$searchlight, sprintf("first_level_%dvox_radius", radius),
rep(subjects, each = length(searchlights)),
sprintf("%s_%s.nii.gz", rep(subjects, each = length(searchlights)), searchlights))
# output folder & file name in MNI space
invisible(lapply(file.path(dirs$searchlight, sprintf("mni_%dvox_radius", radius),
searchlights, "3d"),
function(x) if(!dir.exists(x)) dir.create(x, recursive=TRUE)))
out_nii_fn <- file.path(dirs$searchlight, sprintf("mni_%dvox_radius", radius),
rep(searchlights,n_subs), "3d",
sprintf("%s_%s.nii.gz",
rep(subjects, each = length(searchlights)), searchlights))
# apply FSL flirt to move from whole-brain functional space to 1mm MNI space
invisible(mapply(flirt_apply,
infile = in_nii_fn,
reffile = mni_fname("1"),
initmat = func2standard,
outfile = out_nii_fn,
verbose = FALSE, retimg = FALSE))Merge and smooth first-level images
Now that all first-level images are in MNI space, we are ready to concatenate the images across participants. Then, we subject them to Gaussian smoothing.
# field of view mask
fov_mask <- file.path(dirs$mask_dir, "fov", "fov_mask_mni.nii.gz")
# open the task list
fn <- file.path(htc_dir, "smooth_searchlights_tasklist.txt")
con <- file(fn)
open(con, "w")
for (i_srchlght in 1:length(searchlights)){
# file names of files to merge
in_files <- file.path(dirs$searchlight, sprintf("mni_%dvox_radius", radius),
searchlights[i_srchlght], "3d",
sprintf("%s_%s.nii.gz", subjects,searchlights[i_srchlght]))
# file name of 4D output file
fn_4d <- file.path(dirs$searchlight, sprintf("mni_%dvox_radius", radius),
searchlights[i_srchlght], sprintf("%s_4d.nii.gz",
searchlights[i_srchlght]))
# concatenate the images
fslmerge(infiles = in_files, direction = "t", outfile = fn_4d, retimg = FALSE, verbose = FALSE)
# name of smoothed output file
fn_4d_smooth <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s_4d_smooth_fwhm%d.nii.gz",
searchlights[i_srchlght], FWHM))
# write smoothing command to file
writeLines(sprintf("fslmaths %s -s %f -mas %s %s",
fn_4d, sigma, fov_mask, fn_4d_smooth), con)
}
close(con)
# submit to cluster
cmd <- sprintf("fsl_sub -T 30 -t %s -l %s -M bellmund@cbs.mpg.de -N smooth_srchlghts", fn, htc_dir)
batch_id <- system(cmd, intern = TRUE)
pause_until_batch_done(batch_id = batch_id, wait_interval = 300)Run FSL Randomise
Now we are ready to run FSL Randomise.
Group level statistics were carried out using random sign flipping implemented with FSL Randomise and threshold-free cluster enhancement. We corrected for multiple comparisons using a small volume correction mask including our a priori regions of interest, the anterior hippocampus and the anterior-lateral entorhinal cortex.
# masks to use
gm_mask_fn <- file.path(dirs$mask_dir, "gray_matter", "gray_matter_mask.nii.gz")
svc_mask_fn <- file.path(dirs$mask_dir, "svc", "svc_mask.nii.gz")
# open the tasklist
fn <- file.path(htc_dir, "randomise_searchlights_tasklist.txt")
con <- file(fn)
open(con, "w")
for (i_srchlght in 1:length(searchlights)){
# do we want to look at the positive or negative contrast?
if (any(searchlights[i_srchlght] == c("same_day_vir_time", "day_time_interaction"))){
test_side = ""
} else {
# we want to test for a negative effect for these searchlights
test_side = "_neg"
# multiply by -1 to get the negative contrast
orig_file <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s_4d_smooth_fwhm%d.nii.gz",
searchlights[i_srchlght], FWHM))
mul_neg1_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_4d_smooth_fwhm%d.nii.gz",
searchlights[i_srchlght], test_side, FWHM))
fslmaths(file = orig_file, outfile = mul_neg1_fn, opts = "-mul -1")
}
# 4D input image to run randomise on
in_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_4d_smooth_fwhm%d.nii.gz",
searchlights[i_srchlght], test_side, FWHM))
# output file name for FOV
out_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_randomise_fov_fwhm%d",
searchlights[i_srchlght], test_side, FWHM))
# define randomise command for full FOV and write to file
writeLines(sprintf("randomise -i %s -o %s -1 -T --uncorrp -m %s -n 5000",
in_fn, out_fn, gm_mask_fn),con)
# output file name for SVC
out_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_randomise_svc_fwhm%d",
searchlights[i_srchlght], test_side, FWHM))
# define randomise command for small-volume correction and and write to file
writeLines(sprintf("randomise -i %s -o %s -1 -T --uncorrp -m %s -n 5000",
in_fn, out_fn, svc_mask_fn),con)
}
close(con)
# submit to cluster
cmd <- sprintf("fsl_sub -T 300 -t %s -l %s -M bellmund@cbs.mpg.de -N randomise_srchlghts", fn, htc_dir)
batch_id <- system(cmd, intern = TRUE)
pause_until_batch_done(batch_id = batch_id, wait_interval = 600)Find searchlight clusters and atlas labels
Next, we do an automated search in the Randomise results. This is done via FSL cluster and the result is stored in a text file. Then, we add atlas labels based on the Harvard-Oxford Cortical/Subcortial Structural Atlas using FSL atlasquery.
Here is a function that runs the atlasquery command and removes some excess characters around the return string.
find_FSL_atlas_label <- function(atlas = "Harvard-Oxford Cortical Structural Atlas", x=-20, y=-3, z=-26){
# build FSL atlasquery command and send to command line
cmd<- sprintf('atlasquery -a "%s" -c %s', atlas,sprintf("%f,%f,%f",x,y,z))
cmd<-sprintf("%s | sed 's/<b>//g' | sed 's/<\\/b><br>/\\,/g'", cmd) # based on FSL autoaq (l.106)
string <- system(cmd,intern=TRUE)
# remove atlas name
label <- stringr::str_remove(string, sprintf("%s,", atlas))
return(label)
}We extract clusters that are significant at p<0.05 after correcting for multiple comparisons using small volume correction. To explore the searchlight results beyond the hippocampal-entorhinal region, we look for clusters at a threshold of p<0.001 uncorrected with a minimum extent of 30 voxels.
We corrected for multiple comparisons using a small volume correction mask including our a priori regions of interest, the anterior hippocampus and the anterior-lateral entorhinal cortex. Further, we used a liberal threshold of puncorrected<0.001 to explore the data for additional effects within our field of view. Exploratory searchlight results are shown in Supplemental Figure 9 and clusters with a minimum extent of 30 voxels are listed in Supplemental Tables 9, 11 and 12.
for (i_srchlght in 1:length(searchlights)){
# searchlight with positive or negative contrast?
if (any(searchlights[i_srchlght] == c("same_day_vir_time", "day_time_interaction"))){
test_side <- ""
} else {
test_side <- "_neg"
}
# file with t-values
t_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_randomise_fov_fwhm%d_tstat1.nii.gz",
searchlights[i_srchlght], test_side, FWHM))
# SVC-CORRECTED CLUSTERS
# file with corrected p-values based on SVC
corrp_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_randomise_svc_fwhm%d_tfce_corrp_tstat1.nii.gz",
searchlights[i_srchlght], test_side, FWHM))
# output text file for clusters significant after small volume correction
cluster_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("cluster_%s%s_svc_corrp.txt",
searchlights[i_srchlght], test_side, FWHM))
cmd <- sprintf('cluster -i %s -t 0.95 -c %s --mm > %s',
corrp_fn, t_fn, cluster_fn)
system(cmd)
# read the cluster file
cluster_df <- readr::read_tsv(cluster_fn,
col_types = c("dddddddddd____"))
colnames(cluster_df) <- c("cluster_index", "n_voxel", "max_1minusp",
"x", "y", "z", "cog_x", "cog_y", "cog_z", "t_extr")
# add columns for atlas label info
cluster_df <- cluster_df %>%
tibble::add_column(Harvard_Oxford_Cortical = NA, .before = "cluster_index") %>%
tibble::add_column(Harvard_Oxford_Subcortical = NA, .before = "cluster_index")
# get atlas label info
if (nrow(cluster_df)>0){
for (i in 1:nrow(cluster_df)){
cluster_df$Harvard_Oxford_Cortical[i] <- find_FSL_atlas_label(
atlas = "Harvard-Oxford Cortical Structural Atlas",
x=cluster_df$x[i], y=cluster_df$y[i], z=cluster_df$z[i])
cluster_df$Harvard_Oxford_Subcortical[i] <- find_FSL_atlas_label(
atlas = "Harvard-Oxford Subcortical Structural Atlas",
x=cluster_df$x[i], y=cluster_df$y[i], z=cluster_df$z[i])
}
}
# write to new file
fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("cluster_%s%s_svc_corrp_atlas.txt",
searchlights[i_srchlght], test_side, FWHM))
write_tsv(cluster_df, path=fn)
# FOV CLUSTERS AT p<0.001 UNCORRECTED
# file with uncorrected p-values in entire FOV
uncorrp_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("%s%s_randomise_fov_fwhm%d_tfce_p_tstat1.nii.gz",
searchlights[i_srchlght], test_side, FWHM))
# output text file for clusters significant after small volume correction
cluster_fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("cluster_%s%s_fov_uncorrp.txt",
searchlights[i_srchlght], test_side, FWHM))
cmd <- sprintf('cluster -i %s -t 0.999 -c %s --mm --minextent=30 > %s',
uncorrp_fn, t_fn, cluster_fn)
system(cmd)
# read the cluster file
cluster_df <- readr::read_tsv(cluster_fn,
col_types = c("dddddddddd____"))
colnames(cluster_df) <- c("cluster_index", "n_voxel", "max_1minusp",
"x", "y", "z", "cog_x", "cog_y", "cog_z", "t_extr")
# add columns for atlas label info
cluster_df <- cluster_df %>%
tibble::add_column(Harvard_Oxford_Cortical = NA, .before = "cluster_index") %>%
tibble::add_column(Harvard_Oxford_Subcortical = NA, .before = "cluster_index")
# get atlas label info
for (i in 1:nrow(cluster_df)){
cluster_df$Harvard_Oxford_Cortical[i] <- find_FSL_atlas_label(
atlas = "Harvard-Oxford Cortical Structural Atlas",
x=cluster_df$cog_x[i], y=cluster_df$cog_y[i], z=cluster_df$cog_z[i])
cluster_df$Harvard_Oxford_Subcortical[i] <- find_FSL_atlas_label(
atlas = "Harvard-Oxford Subcortical Structural Atlas",
x=cluster_df$cog_x[i], y=cluster_df$cog_y[i], z=cluster_df$cog_z[i])
}
# write to new file
fn <- file.path(dirs$data4analysis, "searchlight_results",
searchlights[i_srchlght],
sprintf("cluster_%s%s_fov_uncorrp_atlas.txt",
searchlights[i_srchlght], test_side, FWHM))
write_tsv(cluster_df, path=fn)
}Create outlines of significant clusters
To later show which voxels survive corrections for multiple comparisons based on our small-volume correction, we create an outline of the clusters surviving at corrected p<0.05. We do this by dilating a binarized mask of this effect with a spherical kernel with a radius of 2mm.
voxels within black outline are significant after correction for multiple comparisons using small volume correction
Same Day Event Pairs
# to create an outline of the significant cluster, we threshold the small-volume corrected p-image
# at 0.95 (i.e. 1-0.05) and binarize it. Then we dilate it using a spherical kernel.
corrpsvc_fn <- file.path(dirs$data4analysis, "searchlight_results", "same_day_vir_time",
"same_day_vir_time_randomise_svc_fwhm3_tfce_corrp_tstat1.nii.gz")
outl_fn <- file.path(dirs$data4analysis, "searchlight_results", "same_day_vir_time",
"same_day_vir_time_randomise_svc_fwhm3_tfce_corrp_outline.nii.gz")
fslmaths(file=corrpsvc_fn,outfile = outl_fn, opts = "-thr 0.95 -bin",
retimg = FALSE, verbose = FALSE)
outl_nii <- fslmaths(file=outl_fn, outfile = outl_fn,
opts = sprintf("-kernel sphere 2 -dilM -sub %s", outl_fn),
verbose = FALSE)Sequence-Time Interaction
# to create an outline of the significant cluster, we threshold the small-volume corrected p-image
# at 0.95 (i.e. 1-0.05) and binarize it. Then we dilate it using a spherical kernel.
corrpsvc_fn <- file.path(dirs$data4analysis, "searchlight_results","day_time_interaction",
"day_time_interaction_randomise_svc_fwhm3_tfce_corrp_tstat1.nii.gz")
outl_fn <- file.path(dirs$data4analysis, "searchlight_results","day_time_interaction",
"day_time_interaction_randomise_svc_fwhm3_tfce_corrp_outline.nii.gz")
fslmaths(file=corrpsvc_fn,outfile = outl_fn, opts = "-thr 0.95 -bin",
retimg = FALSE, verbose = FALSE)
outl_nii <- fslmaths(file=outl_fn, outfile = outl_fn,
opts = sprintf("-kernel sphere 2 -dilM -sub %s", outl_fn),
verbose = FALSE)7.4 Prepare RSA in searchlight peak
To show that overlapping clusters of voxels drive both the within- and across sequence effects, we will run the across-sequence analysis in the cluster of voxels that show the within-day effect. Because these are independent comparisons, we can use the within-day searchlight to define a region of interest for the across-day analysis.
To test whether within- and across-sequence representations overlap, we defined an ROI based on the within-sequence searchlight analysis. Specifically, voxels belonging to the cluster around the peak voxel, thresholded at p<0.01 uncorrected within our small volume correction mask, were included. The analysis of representational change was then carried out as described for the other ROIs above.
Create the ROI mask
The first steps are to prepare the ROI mask we want to use. Thus, we need to threshold the ROI from the main analysis in MNI space and move the resulting mask to each participant’s functional space for the analysis.
# threshold the searchlight results
in_fn <- file.path(dirs$data4analysis, "searchlight_results", "same_day_vir_time",
"same_day_vir_time_randomise_svc_fwhm3_tfce_p_tstat1.nii.gz")
svc_fn <- file.path(dirs$mask_dir, "svc", "svc_mask.nii.gz")
bin_fn <- here("data", "mri", "rois", "mni_masks", "searchlight_same-day_svc.nii.gz")
fsl_thresh(file = in_fn, outfile = bin_fn, thresh = 0.99, opts = sprintf("-bin -mas %s", svc_fn),
verbose = FALSE, retimg = FALSE)
# peak cluster is in left hemisphere, so don't include any voxels in right hemisphere
lh_fn <- here("data", "mri", "rois", "mni_masks", "left_hemi.nii.gz")
fsl_maths(file=bin_fn, outfile = lh_fn, opts = "-mul 0 -add 1 -roi 91 182 0 218 0 182 0 1",
retimg = FALSE, verbose = FALSE)
fsl_mul(file=bin_fn, file2=lh_fn, outfile=bin_fn,
retimg = FALSE, verbose = FALSE)
# let's have a look at the mask we created
mni_nii <- readNIfTI2(mni_fname("1"))
roi_nii <- readNIfTI2(bin_fn)
coords <- c(statip::mfv1(which(roi_nii==1, arr.ind = TRUE)[,1]),
statip::mfv1(which(roi_nii==1, arr.ind = TRUE)[,2]),
statip::mfv1(which(roi_nii==1, arr.ind = TRUE)[,3]))
ortho2(mni_nii, y = roi_nii, xyz = coords, add.orient = TRUE)
# check the number of voxels in this ROI
sum(c(roi_nii))The resulting ROI mask is now coregistered from MNI 1mm space to the analysis space of the wholebrain functional sequence. Finally, it is thresholded at a probability of 0.5.
samespace_dir <- here("data", "mri", "rois", "same-day_searchlight", "samespace")
if (!dir.exists(samespace_dir)){dir.create(samespace_dir, recursive = TRUE)}
# name of transformation matrix file to move from highres to functional space
standard2func <- here("data", "mri", "processed", "wholebrain",
paste0("VIRTEM_P", subjects, ".feat"),
"reg", "standard2example_func.mat")
# use the mean EPI of wholebrain image as a reference
mean_epi <- here("data", "mri", "processed", "samespace", paste0("VIRTEM_P", subjects),
paste0("VIRTEM_P", subjects, "_wholebrain.nii.gz"))
# define output files in samespace (ss) = analysis space based on the wholebrain EPI
roi_ss <- file.path(samespace_dir, sprintf("P%s_%s_ss.nii.gz",
subjects, "same-day_searchlight"))
# apply FSL flirt to move ROI from standard to wholebrain functional space
invisible(mapply(flirt_apply, infile = bin_fn, reffile = mean_epi,
initmat = standard2func, outfile = roi_ss,
verbose = FALSE, retimg = FALSE))
# use fslmaths to binarize the masked ROIs using a threshold of 0.5
out <- mapply(fsl_thresh, file = roi_ss, outfile = roi_ss, thresh = 0.5, opts = "-bin",
verbose = FALSE, retimg = FALSE)Calculating the correlation matrices
The following steps are analogous to the preparation of the functional data in the ROI and searchlight analyses. For the main searchlight analyses we already cleaned the voxel-wise time series and extracted the volumes relevant to RSA. Thus, to calculate the correlation matrices and to calculate pattern similarity changes, we fall back onto the scripts from the main ROI analyses and the searchlight.
# for all 10x10 comparisons we will be averaging all comparisons apart from the diagonal
# to exclude same_block comparisons
no_diag <- matrix(data = TRUE, nrow=10, ncol = 10)
diag(no_diag)<- FALSE
out_dir <- here("data", "mri","rsa","correlation_matrices", "same-day_searchlight")
if (!dir.exists(out_dir)){dir.create(out_dir, recursive = TRUE)}
for (i_sub in subjects){
# load the ROI based on the searchlight
roi_fn <-file.path(samespace_dir, sprintf("P%s_%s_ss.nii.gz",
i_sub, "same-day_searchlight"))
roi_nii <- readNIfTI(roi_fn, reorient=FALSE)
# array indices of ROI voxels
roi_vox <- which(roi_nii == 1, arr.ind=TRUE)
sprintf("%s: %d voxels\n", i_sub, nrow(roi_vox))
for (i_run in 1:n_runs){
# load the relevant functional volumes
rel_vol_fn <- file.path(dirs$searchlight, "rel_vols_4D",
sprintf("%s_run%02d_rel_vols.nii.gz", i_sub, i_run))
func_nii <- readNIfTI(rel_vol_fn, reorient=FALSE)
# get the ROI voxels
rel_dat <- array(NA, c(n_pics*n_blocks, nrow(roi_vox))); # images in rows (ROI voxels), voxels in columns
for (i in 1:nrow(roi_vox)) { # i <- 1
curr_vox <- func_nii[roi_vox[i,1], roi_vox[i,2], roi_vox[i,3],]
if (sd(curr_vox) > 0) {rel_dat[,i] <- curr_vox} else { stop("zero variance voxel")}
}
# data is in repetition (row) by voxel (col) format, so we transpose
# to get a voxel x repetition format
rel_dat <- t(rel_dat)
# calculate correlation matrix (trial by trial) for pre and post run
cor_mat_trial <- cor(rel_dat, rel_dat)
# initialize condition by condition correlation matrix for pre and post run
corr_mat <- matrix(nrow = 20, ncol = 20)
# loop over all picture comparisons
for(i_pic1 in 1:20){
for(i_pic2 in 1:20){
# extract the current 10x10 correlation matrix
i1 <- (1+(i_pic1-1)*10):(i_pic1*10)
i2 <- (1+(i_pic2-1)*10):(i_pic2*10)
curr_mat <- cor_mat_trial[i1, i2]
# average the correlations while excluding diagonal (same block comparisons)
corr_mat[i_pic1, i_pic2] <- mean(curr_mat[no_diag])
}
}
# save the correlation matrix
fn <- file.path(out_dir, sprintf("%s_%s_%s_corr_mat.txt",
i_sub, "same-day_searchlight", runs[i_run]))
write.table(corr_mat, fn, append = FALSE, sep = ",",
dec = ".", row.names = FALSE, col.names = FALSE)
}
}Calculate Pattern Similarity Change
In the next step, we calculate how the correlation between patterns change from the first to the second picture viewing task, i.e. through the day learning task. To do this, we load both correlation matrices and reduce them to the upper triangle, excluding the diagonal. Then, the correlations are Fisher Z-transformed and the similarity values from the pre-learning picture viewing task are subtracted from those of the post-learning picture viewing task to isolate pattern similarity change. The correlations and their changes are then saved together with information about the pictures that were compared.
This part is based on the corresponding script from the main ROI analyses.
in_dir <- here("data", "mri","rsa","correlation_matrices", "same-day_searchlight")
out_dir <- here("data", "mri","rsa","pattern_similarity_change", "same-day_searchlight")
if (!dir.exists(out_dir)){dir.create(out_dir, recursive = TRUE)}
for(i_sub in subjects){
# load pre correlation matrix
fn <- file.path(in_dir, sprintf("%s_%s_pre_corr_mat.txt",
i_sub, "same-day_searchlight"))
corr_mat_pre <- read.table(fn, sep = ",", dec = ".")
# load post correlation matrix
fn <- file.path(in_dir, sprintf("%s_%s_post_corr_mat.txt",
i_sub, "same-day_searchlight"))
corr_mat_post <- read.table(fn, sep = ",", dec = ".")
# reduce to upper triangle of correlations (without diagonal)
pre_corrs <- corr_mat_pre[upper.tri(corr_mat_pre, diag = FALSE)]
post_corrs <- corr_mat_post[upper.tri(corr_mat_post, diag = FALSE)]
# create a tibble with the correlations
pics <- which(upper.tri(corr_mat_post), arr.ind=TRUE)
corr_df <- tibble(pic1 = pics[,1], pic2 = pics[,2], pre_corrs, post_corrs)
# Fisher Z transform correlations and calculate pattern similarity change
# by subtracting the pre-correlations from the post-correlations
corr_df$ps_change <- FisherZ(corr_df$post_corrs) - FisherZ(corr_df$pre_corrs)
# write to file
fn <- file.path(out_dir, sprintf("%s_%s_pattern_similarity_change.txt",
i_sub, "same-day_searchlight"))
write.table(corr_df, fn, append = FALSE, sep = ",",
dec = ".", row.names = FALSE, col.names = TRUE)
}Combine behavioral and fMRI data for RSA
To create input for the RSA the next step is to combine the similarity change data with the behavioral data, so that we can do meaningful analyses. First the behavioral data is loaded and brought into the same pair-wise format as the similarity change data. Both datasets are combined for each subject and then written to disk.
This part is based on the corresponding script from the main ROI analyses.
in_dir <- here("data", "mri","rsa","pattern_similarity_change", "same-day_searchlight")
out_dir <- here("data", "mri","rsa","data_for_rsa")
if (!dir.exists(out_dir)){dir.create(out_dir, recursive = TRUE)}
for (i_sub in subjects){
# load the behavioral data
fn <- file.path(dirs$timeline_dat_dir, sprintf("%s_behavior_tbl_timeline.txt", i_sub))
col_types_list <- cols_only(sub_id = col_character(), day = col_factor(),
event = col_integer(), pic = col_integer(),
virtual_time = col_double(), real_time = col_double(),
memory_time = col_double(), memory_order = col_double(),
sorted_day = col_integer())
beh_dat <- read_csv(fn, col_types = col_types_list)
# sort behavioral data according to picture identity
beh_dat_ordered <- beh_dat[order(beh_dat$pic),]
# find the order of comparisons
pairs <- which(upper.tri(matrix(nrow = 20, ncol = 20)), arr.ind=TRUE)
# extract the data for the first and second picture in each pair from the behavioral data
pic1_dat <- beh_dat_ordered[pairs[,1],]
pic2_dat <- beh_dat_ordered[pairs[,2],]
colnames(pic1_dat) <- paste0(colnames(pic1_dat), "1")
colnames(pic2_dat) <- paste0(colnames(pic2_dat), "2")
# combine the two tibbles
pair_dat <- cbind(pic1_dat, pic2_dat)
# reorder the tibble columns to make a bit more sense
pair_dat <- pair_dat %>%
select(sub_id1,
day1, day2,
pic1, pic2,
event1, event2,
virtual_time1, virtual_time2,
real_time1, real_time2,
memory_order1, memory_order2,
memory_time1, memory_time2,
sorted_day1, sorted_day2) %>%
rename(sub_id = sub_id1)
rsa_dat <- tibble()
# load the pattern similarity change data for this ROI
fn <- file.path(in_dir, sprintf("%s_%s_pattern_similarity_change.txt",
i_sub, "same-day_searchlight"))
ps_change_dat <- read.csv(fn)
# make sure files have the same order
assertthat::are_equal(c(pair_dat$pic1, pair_dat$pic2), c(ps_change_dat$pic1, ps_change_dat$pic2))
# add column with ROI name
ps_change_dat <- add_column(ps_change_dat, roi = "same-day_searchlight")
# collect the data from this ROI and merge into long data frame
roi_dat <- cbind(pair_dat, ps_change_dat[,3:6])
rsa_dat <- rbind(rsa_dat, roi_dat)
# write to file
fn <- file.path(dirs$rsa_dat_dir, sprintf("%s_data_for_rsa_same-day_searchlight.txt",i_sub))
write.table(rsa_dat, fn, append = FALSE, sep = ",",
dec = ".", row.names = FALSE, col.names = TRUE)
}Finally, the datasets are combined across subjects.
# set up a dataframe to collect the data
rsa_dat = tibble()
for (i_sub in subjects){
# load data from CSV
fn <- file.path(out_dir, sprintf("%s_data_for_rsa_same-day_searchlight.txt",i_sub))
col_types_list <- cols_only(
sub_id = col_character(),
day1 = col_integer(), day2 = col_integer(),
event1 = col_integer(), event2 = col_integer(),
pic1 = col_integer(), pic2 = col_integer(),
virtual_time1 = col_double(), virtual_time2 = col_double(),
real_time1 = col_double(), real_time2 = col_double(),
memory_time1 = col_double(), memory_time2 = col_double(),
memory_order1 = col_double(), memory_order2 = col_double(),
sorted_day1 = col_integer(), sorted_day2 = col_integer(),
pre_corrs = col_double(), post_corrs = col_double(),
ps_change = col_double(), roi = col_factor()
)
sub_dat <- as_tibble(read_csv(fn, col_types = col_types_list))
# append to table with data from all subjects
rsa_dat <- bind_rows(sub_dat, rsa_dat)
}
# sort the data
rsa_dat <- rsa_dat[with(rsa_dat, order(sub_id, day1, day2, event1, event2)),]
# write to file
fn <- file.path(dirs$data4analysis, "rsa_data_in_same-seq_searchlight_peak.txt")
write.table(rsa_dat, fn, append = FALSE, sep = ",",
dec = ".", row.names = FALSE, col.names = TRUE)