Merge pull request #24 from slamballais/dev

0.8.1

Author: slamballais

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: QDECR
 Type: Package
 Title: Vertex-wise statistical analysis
-Version: 0.8.0
+Version: 0.8.1
 Authors@R: as.person(c(
     "Ryan Muetzel <r.muetzel@erasmusmc.nl> [aut, cre]",
     "Sander Lamballais <s.lamballaistessensohn@erasmusmc.nl> [aut]"

NAMESPACE

@@ -5,10 +5,13 @@ S3method(hist,vw)
 S3method(nobs,vw_fastlm)
 S3method(print,vw)
 S3method(summary,vw_fastlm)
+export(as_mgh)
 export(bsfbm2mgh)
+export(fbm_col_is_0)
 export(fbm_col_mean)
 export(fbm_col_sd)
 export(fbm_col_sum)
+export(fbm_row_is_0)
 export(fbm_row_mean)
 export(fbm_row_sd)
 export(fbm_row_sum)
@@ -21,6 +24,7 @@ export(qdecr_load)
 export(qdecr_save)
 export(qdecr_snap)
 export(reload)
+export(save.mgh)
 export(stacks)
 export(unload)
 importFrom(bigstatsr,cols_along)

NEWS.md

@@ -1,3 +1,16 @@
+# QDECR 0.8.1
+
+## Bug fixes
+* If the estimated smoothness is below 1, we now increase it to 1 to avoid problems down the line.
+* Fixed a bug in qdecr_clusters where it assumes that there is always at least 1 cluster significant.
+* "w-g.pct" files can now be used as a measure by specifying "qdecr_w_g.pct" (underscore instead of hyphen).
+* Added `fwhm` argument to `qdecr_fastlm`, which was missing before.
+
+## New (minor) features
+* Added `cwp_thr` argument to `qdecr_fastlm` and `qdecr` to set the further cluster-wise p-value adjustment (default is 0.025 due to having 2 hemispheres, thus 0.05 / 2).
+* Automatically output two extra files: "significant_clusters.txt" (the output of `summary(vw, annot = TRUE)`) and "stack_names.txt" (the output of `stacks(vw)` and the corresponding stack numbers).
+* Modified `freeview` and `qdecr_snap`. The `mask` argument is now called `sig`. Furthermore, the ranges for the overlay colors are determined dynamically. Finally, users can now set any arguments to Freeview for manipulating surface files (see `freeview --help` on the command line).
+
 # QDECR 0.8.0: Momo
 
 Version 0.8.0 is the first update after public release. It fixes a bunch of mistakes, introduces further modularization, improves the speed and also reduces the RAM load.
@@ -24,4 +37,4 @@ Version 0.7.0 is the first version that is publically released. It is also the v
 * Creation of the vw object and associated functions
 * Handling of imputed datasets (via imp2list)
 * Creation of all summary and plot functions
-* Creation of the `stack` concept
+* Creation of the `stack` concept

R/io_mgh.R

@@ -10,6 +10,8 @@ qdecr_prep_mgh <- function(input_path,
                            files_list = list.files(input_path),
                            mask, hemi, measure, fwhmc, target, 
                            n_cores, dir_tmp, project, backing, verbose) {
+  measure2 <- measure
+  if(measure2 == "w_g.pct") measure2 <- "w-g.pct"
   new_files <- file.path(input_path, files_list, "surf", paste(hemi, measure, fwhmc, target, "mgh", sep = "."))
   n <- length(new_files)
 
@@ -115,7 +117,7 @@ load.annot <- function(input.file) {
        )
 }
 
-#' Save out an MGH file from memory
+#' Save out an MGH file from a big matrix
 #'
 #' @param vol MGH object (as from load.mgh)
 #' @param fname file name to be used to save out the data
@@ -173,6 +175,13 @@ bsfbm2mgh <-function(fbm, fname, filter = NULL) {
   NULL
 }
 
+#' Save out an MGH file from memory
+#'
+#' @param vol MGH object (as from load.mgh)
+#' @param fname file name to be used to save out the data
+#'
+#' @export
+#'
 save.mgh <-function(vol,fname) {
 
   # R translation of save_mgh.m
@@ -251,4 +260,46 @@ save.mgh <-function(vol,fname) {
   nelts <- width * height   # bytes per slice
   #writeBin(vol$x, fid, size = 4, endian = "big")
   writeBin(vol$x, fid, size = 4, endian = "big")
+}
+
+
+#' Create an object that is structured like an mgh object
+#'
+#' @param x the vertex-wise values
+#' @param v (to be added)
+#' @param ndim1 (to be added)
+#' @param ndim2 (to be added)
+#' @param ndim3 (to be added)
+#' @param nframes (to be added)
+#' @param type (to be added)
+#' @param dof (to be added)
+#'
+#' @export
+#'
+
+as_mgh <- function(x, v = NULL, ndim1 = NULL, ndim2 = NULL, ndim3 = NULL, nframes, type, dof) {
+  
+  if (is.vector(x)) {
+    v <- 1L
+    ndim1 <- as.integer(length(x))
+    ndim2 <- 1L
+    ndim3 <- 1L
+    type <- 3L
+    nframes <- 1L
+    dof <- 0L
+  } else {
+    stop("as_mgh only support objects of class `vector` right now.")
+  }
+  
+  out <- list(x = x, 
+              v = v, 
+              ndim1 = ndim1, 
+              ndim2 = ndim2, 
+              ndim3 = ndim3, 
+              nframes = nframes, 
+              type = type, 
+              dof = dof)
+  
+  return(out)
+  
 }

R/plotting.R

@@ -8,29 +8,123 @@
 #'
 #'@param vw The output object of a qdecr call (e.g. qdecr_fastlm)
 #'@param stack Either a numeric or a string for your variable (use `stacks(vw)`)
-#'@param mask Logical; if TRUE, only the significant clusters are shown
+#'@param type The type of map to plot (currently for OLS). Default is `coef`, others are `se`, `t` and `p`
+#'@param sig Logical; if TRUE (default), only the significant clusters are shown
+#'@param ... any arguments that freeview (on command line) normally takes for manipulating surface files
 #'@return NULL
 #'@export
 
-freeview <- function(vw, stack = NULL, mask = TRUE) {
+freeview <- function(vw, stack = NULL, type = c("coef", "se", "t", "p"), sig = TRUE, ...) {
   if (is.null(stack)) stop("`stack` not defined. Please choose: ", paste(stacks(vw), collapse = ", "))
   if (length(stack) > 1) stop("More than 1 stack specified.")
   if(is.character(stack)) stack <- which(stacks(vw) == stack)
   if(is.null(stack) || stack > length(stacks(vw))) stop("specified `stack` is not present in this dataset. Please choose: ", paste(stacks(vw), collapse = ", "))
+  type <- match.arg(type)
+  empty_val <- if(type == "p") 1 else 0
+
+  read_fun <- switch(type, 
+                     coef = qdecr_read_coef, 
+                     se = qdecr_read_se, 
+                     t = qdecr_read_t,
+                     p = qdecr_read_p)
 
-  temp_mgh <- qdecr_read_coef(vw, stack)
-  if (mask) temp_mgh$x <- temp_mgh$x %MASK% qdecr_read_ocn_mask(vw, stack)
-  if (all(temp_mgh$x == 0)) stop("Stack does not contain information (e.g. because of no significant findings), aborting plot.")
+  temp_mgh <- read_fun(vw, stack)
+  if (all(temp_mgh$x == empty_val)) stop("Stack does not contain information (e.g. because of no significant findings), aborting plot.")
+  if (sig) {
+    p_mgh <- qdecr_read_ocn_mask(vw, stack)
+    temp_mgh$x[!p_mgh] <- empty_val
+    if (all(temp_mgh$x == empty_val)) stop("No information in the stack passed the threshold (i.e. `p_thr` was set too strict), aborting plot.")
+  }
+
   temp_mgh_file <- paste0(vw$paths$dir_tmp2, ".stack", stack, ".temp_mgh.mgh")
   save.mgh(temp_mgh, temp_mgh_file)
   on.exit(file.remove(temp_mgh_file), add = TRUE)
 
-  cmdStr <- paste0("freeview --surface ", vw$paths$dir_subj, "/fsaverage/surf/", vw$input$hemi, ".inflated:overlay=", temp_mgh_file)
-  system(cmdStr)
+  surface <- paste0("--surface ", vw$paths$dir_subj, "/fsaverage/surf/", vw$input$hemi, ".inflated")
+
+  freeview_args <- list(...)
+  freeview_args$surface <- surface
+  freeview_args$overlay <- temp_mgh_file
+  
+  if (is.null(freeview_args$overlay_method)) {
+    freeview_args$overlay_method <- "linearopaque"
+    temp_vals <- temp_mgh$x[temp_mgh$x != empty_val]
+    freeview_args$overlay_threshold <- c(min(abs(temp_vals)), max(abs(temp_vals)))
+  }
+  
+  surf_cmd <- do.call(freeview_surf_cmd, freeview_args)
+  cmd_str <- paste("freeview", surf_cmd)
+  cat("Opening Freeview on the command line. Pop-up incoming.\nCall:", cmd_str, "\n")
+  system(cmd_str)
 
   invisible(NULL)
 }
 
+freeview_surf_cmd <- function(surface,
+                              curvature = NULL, 
+                              curvature_method = NULL,
+                              overlay = NULL,
+                              overlay_reg = NULL,
+                              overlay_method = NULL,
+                              overlay_color = NULL,
+                              overlay_threshold = NULL,
+                              overlay_frame = NULL,
+                              correlation = NULL,
+                              color = NULL,
+                              edgecolor = NULL,
+                              edgethickness = NULL,
+                              annot = NULL,
+                              annot_outline = NULL,
+                              name = NULL,
+                              offset = NULL,
+                              visible = NULL,
+                              vector = NULL,
+                              target_surf = NULL,
+                              label = NULL,
+                              label_outline = NULL,
+                              label_color = NULL,
+                              label_centroid = NULL,
+                              label_visible = NULL,
+                              spline = NULL,
+                              vertex = NULL,
+                              vertexcolor = NULL,
+                              goto = NULL,
+                              hide_in_3d = NULL,
+                              all = NULL) {
+  string <- surface
+  if(!is.null(curvature)) string <- paste0(string, ":curvature=", curvature)
+  if(!is.null(curvature_method)) string <- paste0(string, ":curvature_method=", curvature_method)
+  if(!is.null(overlay)) string <- paste0(string, ":overlay=", overlay)
+  if(!is.null(overlay_reg)) string <- paste0(string, ":overlay_reg=", overlay_reg)
+  if(!is.null(overlay_method)) string <- paste0(string, ":overlay_method=", overlay_method)
+  if(!is.null(overlay_color)) string <- paste0(string, ":overlay_color=", collapse(overlay_color, collapse = ","))
+  if(!is.null(overlay_threshold)) string <- paste0(string, ":overlay_threshold=", collapse(overlay_threshold, collapse = ","))
+  if(!is.null(overlay_frame)) string <- paste0(string, ":overlay_frame=", overlay_frame)
+  if(!is.null(correlation)) string <- paste0(string, ":correlation=", correlation)
+  if(!is.null(color)) string <- paste0(string, ":color=", collapse(color, collapse = ","))
+  if(!is.null(edgecolor)) string <- paste0(string, ":edgecolor=", collapse(edgecolor, collapse = ","))
+  if(!is.null(edgethickness)) string <- paste0(string, ":edgethickness=", edgethickness)
+  if(!is.null(annot)) string <- paste0(string, ":annot=", annot)
+  if(!is.null(annot_outline)) string <- paste0(string, ":annot_outline=", annot_outline)
+  if(!is.null(name)) string <- paste0(string, ":name=", name)
+  if(!is.null(offset)) string <- paste0(string, ":offset=", collapse(offset, collapse = ","))
+  if(!is.null(visible)) string <- paste0(string, ":visible=", visible)
+  if(!is.null(vector)) string <- paste0(string, ":vector=", vector)
+  if(!is.null(target_surf)) string <- paste0(string, ":target_surf=", target_surf)
+  if(!is.null(label)) string <- paste0(string, ":label=", label)
+  if(!is.null(label_outline)) string <- paste0(string, ":label_outline=", label_outline)
+  if(!is.null(label_color)) string <- paste0(string, ":label_color=", label_color)
+  if(!is.null(label_centroid)) string <- paste0(string, ":label_centroid=", label_centroid)
+  if(!is.null(label_visible)) string <- paste0(string, ":label_visible=", label_visible)
+  if(!is.null(spline)) string <- paste0(string, ":spline=", spline)
+  if(!is.null(vertex)) string <- paste0(string, ":vertex=", vertex)
+  if(!is.null(vertexcolor)) string <- paste0(string, ":vertexcolor=", collapse(vertexcolor, collapse = ","))
+  if(!is.null(goto)) string <- paste0(string, ":goto=", goto)
+  if(!is.null(hide_in_3d)) string <- paste0(string, ":hide_in_3d=", hide_in_3d)
+  if(!is.null(all)) string <- paste0(string, ":all=", all)
+  string
+}
+
 #'Histograms of the vertex measures
 #'
 #'Plots a histogram of the mean vertex measure, either by vertex or by subject
@@ -64,16 +158,18 @@ hist.vw <- function(vw, qtype = c("vertex", "subject"), xlab = NULL, main = NULL
 #'
 #'@param vw The output object of a qdecr call (e.g. qdecr_fastlm)
 #'@param stack Either a numeric or a string for your variable (use `stacks(vw)`)
+#'@param type The type of map to plot (currently for OLS). Default is `coef`, others are `se`, `t` and `p`
 #'@param ext Extension of the image files that will be stored on disk
 #'@param zoom Float that determines how far the brains are zoomed in
 #'@param compose Logical; if TRUE, a single compiled image will be made (requires Magick++)
 #'@param plot_brain Logical; if TRUE, returns a graphical device with the composed images
 #'@param save_plot Logical; if TRUE, saves the composed image
-#'@param mask Logical; if TRUE, only the significant clusters are shown
+#'@param sig Logical; if TRUE, only the significant clusters are shown
+#'@param ... any arguments that freeview (on command line) normally takes for manipulating surface files
 #'@return NULL
 #'@export
 
-qdecr_snap <- function(vw, stack = NULL, ext = ".tiff", zoom = 1, compose = TRUE, plot_brain = TRUE, save_plot = TRUE, mask = TRUE) {
+qdecr_snap <- function(vw, stack = NULL, type = c("coef", "se", "t", "p"), ext = ".tiff", zoom = 1, compose = TRUE, plot_brain = TRUE, save_plot = TRUE, sig = TRUE, ...) {
   if (is.null(stack)) stop("`stack` not defined. Please choose: ", paste(stacks(vw), collapse = ", "))
   if (length(stack) > 1) stop("More than 1 stack specified.")
   if(is.character(stack)) stack <- which(stacks(vw) == stack)
@@ -98,19 +194,48 @@ qdecr_snap <- function(vw, stack = NULL, ext = ".tiff", zoom = 1, compose = TRUE
                 qsnap_e(180),
                 qsnap(snap_names[4]),
                 "--quit")
-  tfile <-  file.path(vw$paths$dir_tmp, "tmp_snapshot_qdecr.txt")
+  tfile <- file.path(vw$paths$dir_tmp, "tmp_snapshot_qdecr.txt")
   write.table(snap_cmd, tfile, quote = F, row.names = F, col.names = F)
   on.exit(file.remove(tfile))
 
-  temp_mgh <- qdecr_read_coef(vw, stack)
-  if (mask) temp_mgh$x <- temp_mgh$x %MASK% qdecr_read_ocn_mask(vw, stack)
-  if (all(temp_mgh$x == 0)) stop("Stack does not contain information (e.g. because of no significant findings), aborting plot.")
+  type <- match.arg(type)
+  empty_val <- if(type == "p") 1 else 0
+  
+  read_fun <- switch(type, 
+                     coef = qdecr_read_coef, 
+                     se = qdecr_read_se, 
+                     t = qdecr_read_t,
+                     p = qdecr_read_p)
+  
+  temp_mgh <- read_fun(vw, stack)
+  if (all(temp_mgh$x == empty_val)) stop("Stack does not contain information (e.g. because of no significant findings), aborting plot.")
+  if (sig) {
+    p_mgh <- qdecr_read_ocn_mask(vw, stack)
+    temp_mgh$x[!p_mgh] <- empty_val
+    if (all(temp_mgh$x == empty_val)) stop("No information in the stack passed the threshold (i.e. `p_thr` was set too strict), aborting plot.")
+  }
+
   temp_mgh_file <- paste0(vw$paths$dir_tmp2, ".stack", stack, ".temp_mgh.mgh")
   save.mgh(temp_mgh, temp_mgh_file)
   on.exit(file.remove(temp_mgh_file), add = TRUE)
 
-  cmdStr <- paste0("freeview --surface ", vw$paths$dir_subj, "/", vw$input$target, "/surf/", vw$input$hemi, ".inflated:overlay=", temp_mgh_file, " -cmd ", tfile)
-  system(cmdStr)
+  surface <- paste0("--surface ", vw$paths$dir_subj, "/fsaverage/surf/", vw$input$hemi, ".inflated")
+  
+  freeview_args <- list(...)
+  freeview_args$surface <- surface
+  freeview_args$overlay <- temp_mgh_file
+  
+  if (is.null(freeview_args$overlay_method)) {
+    freeview_args$overlay_method <- "linearopaque"
+    temp_vals <- temp_mgh$x[temp_mgh$x != empty_val]
+    freeview_args$overlay_threshold <- c(min(abs(temp_vals)), max(abs(temp_vals)))
+  }
+  
+  surf_cmd <- do.call(freeview_surf_cmd, freeview_args)
+  cmd_str <- paste("freeview", surf_cmd, "-cmd", tfile)
+  cat("Opening Freeview on the command line. Pop-up incoming.\nCall:", cmd_str, "\n")
+  cat(tfile, "contains:", snap_cmd, "\n")
+  system(cmd_str)
 
   if (compose) {
     if (!requireNamespace("magick", quietly = TRUE)) {

R/qdecr.R

@@ -33,9 +33,10 @@ qdecr <- function(id,
                   hemi = c("lh", "rh"),
                   measure = c("thickness", "area", "area.pial", "curv", 
                               "jacobian_white", "pial", "pial_lgi", "sulc", "volume", 
-                              "w-g.pct.mgh", "white.H", "white.K"),
+                              "w_g.pct", "white.H", "white.K"),
                   fwhm = ifelse(measure == "pial_lgi", 5, 10),
                   mcz_thr = 30, 
+                  cwp_thr = 0.025,
                   mgh = NULL,
                   mask = NULL,
                   mask_path = system.file("extdata", paste0(hemi, ".fsaverage.cortex.mask.mgh"), package = "QDECR"),
@@ -176,14 +177,18 @@ qdecr <- function(id,
     message2(paste0("Estimated smoothness is ", vw$post$fwhm_est, ", which is really high. Reduced to 30."), verbose = verbose)
     vw$post$fwhm_est <- 30
   }
+  if(vw$post$fwhm_est < 1) {
+    message2(paste0("Estimated smoothness is ", vw$post$fwhm_est, ", which is really low. Increased to 1."), verbose = verbose)
+    vw$post$fwhm_est <- 1
+  }
 
   catprompt("Clusterwise correction", verbose = verbose)
 
   # Run clusterwise and move all files to the right directory
   for ( i in seq_along(vw$stack$names)){
     message2("\n", verbose = verbose)
     runMriSurfCluster(vw$paths$final_path, vw$paths$dir_fshome, vw$input$hemi, vw$stack$p[[i]], vw$post$fwhm_est, 
-                      mask_path = vw$paths$final_mask_path, verbose = verbose, mczThr = mcz_thr, 
+                      mask_path = vw$paths$final_mask_path, verbose = verbose, mcz_thr = mcz_thr, cwp_thr = cwp_thr,
                       stack = i, stack_name = vw$stack$names[i])
   }
   vw$post$final_mask <- load.mgh(vw$paths$final_mask_path)$x