R/fsl_generate_fsf_lvl1_ev_syntax.R

In UNCDEPENdLab/fmri.pipeline: Analysis of fMRI Data on HPC clusters

#' This function generates syntax for FSL Feat .fsf files for the EVs tab of a first-level fMRI analysis.
#' It accepts a numeric design matrix whose colum names correspond to individual EVs in the model
#'
#' This allows you to generate a design matrix dynamically based on a numeric matrix that has been
#' setup in R using a function such as lm() to handle more complex designs appropriately. The syntax
#' generated by this function can be combined with a .fsf template file to implement the entire LVL2 or LVL3
#' analysis setup.
#'
#' @param regressors A list of regressors containing the specification for each one as an EV.
#'
#' @details
#'
#' Each element of the \code{regressors} list is itself a list specifying how to handle a given regressor.
#' Required elements of the list are:
#'
#' \preformatted{
#'   $waveform: The type of input for this regressor. Options are:
#'      "custom_1"    : 1-column text file containing regressor values (can be convolved or unconvolved)
#'      "custom_3"    : 3-column FSL-format text file containing onset time, duration, and regressor value
#'      "empty"       : an all-zero regressor
#'      "square"      : for block design (not currently supported)
#'      "sinusoid"    : for periodic regressors (not currently supported)
#'      "interaction" : unclear what this is, but exists in .fsf files (not currently supported)
#'
#'   $convolution: The type of convolution function to be applied to the $waveform. Options are:
#'      "none"         : no convolution is applied (esp. useful when using custom_1 files that are already convolved)
#'      "double_gamma" : convolve with FSL-variant of double gamma function
#'      "gaussian"     : convolve with Gaussian function (not tested)
#'      "gamma_basis"  : Gamma basis functions (not tested)
#'      "sine_basis"   : sine basis functions (not tested)
#'      "fir_basis"    : finite impulse response (FIR) basis functions (not tested)
#'      "alternate_double_gamma" : what it sounds like (not tested)
#'
#'   $name: The name of the EV (i.e., regressor) in the model specification and output.
#' }
#'
#' Optional elements of the list are:
#'
#' \preformatted{
#'   $timing_file: the name of the text file containing regressor values. Used by "custom_1" and "custom_3".
#'   $tempfilt: whether to apply temporal filtering to the regressor, as with the fMRI data. Default: 1 (TRUE)
#'   $deriv: whether to add a temporal derivative for this regressor to the model. Default 0 (FALSE)
#'   $convolve_phase: not entirely sure what it does but is 0/1. Default: 0 (FALSE)
#' }
#'
#' @return A character vector containing .fsf syntax for the EV portion of a first-level Feat analysis
#' based on \code{regressors}.
#'
#' @examples
#'
#'   regressors <- list(
#'     list(name="cue_onset", waveform="custom_3", timing_file="r1_cue_onset_3col.txt", convolution="double_gamma"),
#'     list(name="feedback",  waveform="custom_3", timing_file="r1_feedback_3col.txt", convolution="double_gamma"),
#'     list(name="RPE", waveform="custom_1", timing_file="r1_RPEs.1D", convolution="none")
#'   )
#'
#'   ev_syntax <- fsl_generate_fsf_lvl1_ev_syntax(regressors)
#'
#' @author Michael Hallquist
#' @export
fsl_generate_fsf_lvl1_ev_syntax <- function(regressors) {

  stopifnot(is.list(regressors))
  #stopifnot(is.matrix(cmat))

  #helper function to validate defaults for each regressor
  validate_regressors <- function(regressors) {
    lapply(regressors, function(r) {
      if (is.null(r$waveform)) { r$waveform <- "custom_1" } #assume 1-column format
      if (is.null(r$convolution)) { r$convolution <- "none" } #assume no convolution applied
      if (is.null(r$convolve_phase)) { r$convolve_phase <- 0 } #not too sure what this does!
      if (is.null(r$tempfilt)) { r$tempfilt <- 1 } #apply temporal filtering by default
      if (is.null(r$deriv)) { r$deriv <- 0 } #do not add temporal derivative by default
      return(r)
    })
  }

  regressors <- validate_regressors(regressors)

  n_evs <- length(regressors)
  n_tempderiv <- sum(sapply(regressors, function(x) { isTRUE(x$deriv) })) #number of temporal derivatives in the model

  fsf_syntax <- c(
    "# Number of EVs",
    paste0("set fmri(evs_orig) ", n_evs),
    paste0("set fmri(evs_real) ", n_evs + n_tempderiv),
    "set fmri(evs_vox) 0",
    ""
  )

  ## "# Number of contrasts",
  ##   paste0("set fmri(ncon_orig) ", nrow(cmat)),
  ##   paste0("set fmri(ncon_real) ", nrow(cmat)),
  ##   "",
  ##   "# Number of F-tests",
  ##   paste0("set fmri(nftests_orig) ", "0"), #nrow(cmat)), ##NOT SUPPORTED YET
  ##   paste0("set fmri(nftests_real) ", "0") #nrow(cmat))
  ## )


  map_waveform_code <- function(waveform) {
    switch(tolower(waveform),
      square = 0,
      sinusoid = 1,
      custom_1 = 2,
      custom_3 = 3,
      interaction = 4,
      empty = 10)
  }

  map_convolution_code <- function(convolution) {
    switch(tolower(convolution),
      none = 0,
      gaussian = 1,
      gamma = 2,
      double_gamma = 3,
      gamma_basis = 4,
      sine_basis = 5,
      fir_basis = 6,
      alternate_double_gamma = 8)
  }

  #loop over EV columns in the design
  for (j in seq_along(regressors)) {
    fsf_syntax <- c(fsf_syntax,
      paste0("# EV ", j, " title"),
      paste0("set fmri(evtitle", j, ") \"", regressors[[j]]$name, "\""),
      "",
      paste0("# Basic waveform shape (EV ", j, ")"),
      "# 0 : Square",
      "# 1 : Sinusoid",
      "# 2 : Custom (1 entry per volume)",
      "# 3 : Custom (3 column format)",
      "# 4 : Interaction",
      "# 10 : Empty (all zeros)",
      paste0("set fmri(shape", j, ") ", map_waveform_code(regressors[[j]]$waveform)),
      "",
      paste0("# Convolution (EV ", j, ")"),
      "# 0 : None",
      "# 1 : Gaussian",
      "# 2 : Gamma",
      "# 3 : Double-Gamma HRF",
      "# 4 : Gamma basis functions",
      "# 5 : Sine basis functions",
      "# 6 : FIR basis functions",
      "# 8 : Alternate Double-Gamma",
      paste0("set fmri(convolve", j, ") ", map_convolution_code(regressors[[j]]$convolution)),
      "",
      paste0("# Convolve phase (EV ", j, ")"),
      paste0("set fmri(convolve_phase", j, ") ", regressors[[j]]$convolve_phase),
      "",
      paste0("# Apply temporal filtering (EV ", j, ")"),
      paste0("set fmri(tempfilt_yn", j, ") ", regressors[[j]]$tempfilt),
      "",
      paste0("# Add temporal derivative (EV ", j, ")"),
      paste0("set fmri(deriv_yn", j, ") ", regressors[[j]]$deriv)
    )

    #for custom waveforms, we need to specify the text file containing the values
    if (regressors[[j]]$waveform %in% c("custom_1", "custom_3")) {
      fsf_syntax <- c(fsf_syntax,
        "",
        paste0("# Custom EV file (EV ", j, ")"),
        paste0("set fmri(custom", j, ") \"", regressors[[j]]$timing_file, "\""),
        ""
      )
    }

    #setup EV orthogonalizations (none supported at present)
    #also, not sure why you need a section for orthogonalization wrt EV 0, but keeping it for consistency with Feat
    for (w in 0:length(regressors)) {
      fsf_syntax <- c(fsf_syntax,
        paste0("# Orthogonalise EV ", j, " wrt EV ", w),
        paste0("set fmri(ortho", j, ".", w, ") 0"),
        "")
    }

  }

  return(fsf_syntax)

}