diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index beb6c1e18a8c82d3eadda642f743e168c46b1dc8..8c03b2c090d45d9bb6449b6c3cff262a03098487 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -6,7 +6,7 @@ stages: before_script: - apt-get update - - apt-get --yes install libglpk-dev + - apt-get --yes install libglpk-dev libicu-dev - mkdir -p rlib - echo 'R_LIBS="rlib"' > .Renviron - echo 'R_LIBS_USER="rlib"' >> .Renviron @@ -14,17 +14,16 @@ before_script: cache: key: r-library - paths: + paths: - rlib prep: stage: .pre script: - R --version - - R -e "install.packages(c('devtools', 'usethis'))" - - R -e "devtools::install_gitlab('choh/littlehelpers')" - - R -e "devtools::install_gitlab('choh/rsfmri')" - - R -e "devtools::install_deps(dependencies = TRUE)" + - R -e "install.packages(c('remotes', 'stringi', 'devtools', 'usethis', 'stringr', 'lintr'), type = 'source', repos = 'https://cran.uni-muenster.de/')" + - R -e "remotes::install_gitlab('choh/littlehelpers')" + - R -e "remotes::install_git('https://git.rwth-aachen.de/christian.hohenfeld/rsfmri.git')" - R -e "devtools::install()" lintr: diff --git a/DESCRIPTION b/DESCRIPTION index f6d0a926487fabeb2e228faad9c04a570162e4a1..ee7976627fa54503f0d517c2203fe62fe2322524 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: rsAnalysis Title: High-Level Tools for rsfMRI analysis -Version: 0.4.5.9000 +Version: 0.4.6 Authors@R: person(given = "Christian", family = "Hohenfeld", @@ -41,6 +41,6 @@ Suggests: withr Remotes: gitlab::choh/littlehelpers, - gitlab::choh/rsfmri + git::https://git.rwth-aachen.de/christian.hohenfeld/rsfmri.git Depends: R (>= 2.10) diff --git a/R/check_motion.R b/R/check_motion.R index 76f1d7329a4ecaec5ad921474ae9e5b5da3f2197..a881816ae6b107b9c0089bfd857b749b07261511 100644 --- a/R/check_motion.R +++ b/R/check_motion.R @@ -58,8 +58,8 @@ check_motion <- function(rawdir = NULL, pattern = NULL, parfile = NULL, spikes <- apply(motion_data, 2, cdiff) - spikes[4:6] <- 100 * pi * spikes[4:6] / 360 - has_spike <- rowSums(abs(spikes)) > 0.5 + spikes[, 1:3] <- 100 * pi * spikes[, 1:3] / 360 + has_spike <- rowSums(abs(spikes)) > spike_cutoff motion_data <- motion_data |> dplyr::mutate( diff --git a/R/compare_measure_group.R b/R/compare_measure_group.R index 6fd54d6083a7adea3c9b27865586b3f5366b405e..2a23b3b4ee0bf86550b8e1f817d67fa5cf7ac6c3 100644 --- a/R/compare_measure_group.R +++ b/R/compare_measure_group.R @@ -5,44 +5,51 @@ #' @param ... Strings giving variables to include in the model formula. #' @param add_name Whether to add interaction term group x name, default FALSE #' @param remove Name of a group to remove +#' @param anova Logical value specifying whether to compute an ANOVA. #' @importFrom rlang .data #' @export compare_measure_group <- function(rawlist, dv, ..., add_name = FALSE, - remove = NULL) { - grplist <- lapply(seq_along(rawlist), - function(x) { - rawlist[[x]]$group <- names(rawlist)[x] - rawlist[[x]] - }) - grplist <- do.call("rbind", grplist) + remove = NULL, anova = TRUE) { + grplist <- lapply(seq_along(rawlist), function(x) { + rawlist[[x]]$group <- names(rawlist)[x] + rawlist[[x]] + }) + grplist <- do.call("rbind", grplist) - if (!is.null(remove)) { - grplist <- grplist %>% dplyr::filter(.data$group != remove) - } + if (!is.null(remove)) { + grplist <- grplist %>% dplyr::filter(.data$group != remove) + } - additional_formula <- paste(..., sep = "+") - additional_formula <- ifelse(length(additional_formula) == 0, "", - paste0("+", additional_formula)) + additional_formula <- paste(..., sep = "+") + additional_formula <- ifelse(length(additional_formula) == 0, "", + paste0("+", additional_formula)) - if (add_name) { - frml <- stats::as.formula(paste(dv, "~ group * name", additional_formula)) - } else { - frml <- stats::as.formula(paste(dv, " ~ group", additional_formula)) - } - mod <- stats::lm(frml, data = grplist) - modtbl <- mod %>% - broom::tidy() %>% - dplyr::mutate(s.value = -log2(.data$p.value)) + if (add_name) { + frml <- stats::as.formula( + paste(dv, "~ group * name", additional_formula)) + } else { + frml <- stats::as.formula(paste(dv, " ~ group", additional_formula)) + } + mod <- stats::lm(frml, data = grplist) + modtbl <- mod %>% + broom::tidy() %>% + dplyr::mutate(s.value = -log2(.data$p.value)) - av <- car::Anova(mod, type = 3) - tbl <- av %>% - broom::tidy() %>% - dplyr::mutate(s.value = -log2(.data$p.value)) + if (anova) { + av <- car::Anova(mod, type = 3) + tbl <- av %>% + broom::tidy() %>% + dplyr::mutate(s.value = -log2(.data$p.value)) + } else { + av <- NA + tbl <- NA + } - bayes <- rstanarm::stan_glm(frml, data = grplist, iter = 2000, refresh = 0) - bayes_tbl <- parameters::parameters(bayes) + bayes <- rstanarm::stan_glm(frml, data = grplist, iter = 2000, refresh = 0) + bayes_tbl <- parameters::model_parameters( + bayes, test = c("pd", "rope", "bf")) - list(lm = mod, lm_tbl = modtbl, - anova = av, anova_tbl = tbl, - bayesian = bayes, bayesian_tbl = bayes_tbl) + list(lm = mod, lm_tbl = modtbl, + anova = av, anova_tbl = tbl, + bayesian = bayes, bayesian_tbl = bayes_tbl) } diff --git a/R/rs_pipeline.R b/R/rs_pipeline.R index a5ec4616f62d9bf2de2e5eb719d6c215ec910cfd..1eacfc204cc4c2c882b62326a42f12468d02193c 100644 --- a/R/rs_pipeline.R +++ b/R/rs_pipeline.R @@ -1,129 +1,235 @@ #' Run a complete resting-state fMRI pipeline. #' #' @param functional Path to a functional dataset. -#' @param anatomy Path to a corresponing anatomical dataset. -#' @param physfile Path to corresponding pyhsiological data. +#' @param anatomy Path to a corresponding anatomical dataset. #' @param std_anat Path to standard anatomical file. -#' @param std_white_matter Path to standard white matter mask. -#' @param std_csf Path to standard CSF mask. +#' @param std_brain Path to skull stripped standard file. #' @param labels_img Path to atlas image. #' @param labels_text Path to atlas label info file. +#' @param magnitude_map Path to magnitude B0 map. +#' @param phase_map Path to phase B0 map. +#' @param dwell_time Dwell time in s. +#' @param ap_map Phase map in AP direction. +#' @param pa_map Phase map in PA direction. +#' @param readout_time Readout time for AP/PA fieldmaps. +#' @param echo_spacing Echo spacing time in ms. +#' @param fnirt_conf Name of the FSL fnirt config file to use. #' @param gsr Boolean whether to use global signal regression. #' @param tr The TR of the data im ms. #' @param odd Bool indicating whether slices were acquired interleaved. +#' @param scale_anat_to Optional value to scale the anatomical image resolution +#' to. +#' @param motion_cutoff Motion correction cutoff value in deg/mm. +#' @param keep_intermediate Should itermediate files be kept on disk? This is +#' useful for debugging. #' #' @return A list containing time courses and a pairwise correlation matrix. #' @export +# nolint start: cyclocomp_linter rs_pipeline <- function( functional, anatomy, - physfile = NULL, std_anat, - std_white_matter, - std_csf, + std_brain, labels_img, labels_text, + magnitude_map = NULL, + phase_map = NULL, + dwell_time = NULL, + ap_map = NULL, + pa_map = NULL, + readout_time = NULL, + echo_spacing = NULL, + fnirt_conf, gsr = TRUE, tr = NULL, - odd = NULL) { + odd = NULL, + scale_anat_to = NULL, + motion_cutoff = 1, + keep_intermediate = FALSE) { + old_scipen <- options(scipen = 200) + on.exit(options(old_scipen)) + + functional <- normalizePath(functional) + input_dir <- dirname(functional) + temp_name <- paste0("proc", as.numeric(Sys.time())) + temp_dir <- file.path(input_dir, temp_name) + dir.create(temp_dir) + + copy_to_temp <- function(in_file, temp_dir) { + old_path <- normalizePath(in_file) + base <- basename(in_file) + new <- file.path(temp_dir, base) + copy_success <- file.copy(old_path, new) + if (!copy_success) { + stop("could not copy") + } + new + } + + functional <- copy_to_temp(functional, temp_dir) + anatomy <- copy_to_temp(anatomy, temp_dir) + + mag_path <- NULL + pha_path <- NULL + if (!is.null(magnitude_map) || !is.null(phase_map)) { + mag_path <- copy_to_temp(magnitude_map, temp_dir) + pha_path <- copy_to_temp(phase_map, temp_dir) + } + + ap_path <- NULL + pa_path <- NULL + if (!is.null(ap_map) || !is.null(pa_map)) { + ap_path <- copy_to_temp(ap_map, temp_dir) + pa_path <- copy_to_temp(pa_map, temp_dir) + } + std_anat <- normalizePath(std_anat) - std_white_matter <- normalizePath(std_white_matter) - std_csf <- normalizePath(std_csf) + std_brain <- normalizePath(std_brain) labels_img <- normalizePath(labels_img) labels_text <- normalizePath(labels_text) - if (!is.null(physfile)) { - input_files <- normalizePath(c(functional, anatomy, physfile)) - } else { - input_files <- normalizePath(c(functional, anatomy)) + print_msg <- function(txt) { + print(paste(date(), txt)) } - input_dir <- dirname(functional) - print("Standardising orientation...") - ori_func <- rsfmri::fsl_reorient2std(input_files[1]) - ori_anat <- rsfmri::fsl_reorient2std(input_files[2]) + print_msg("Standardising orientation...") + ori_func <- rsfmri::fsl_reorient2std(functional) + ori_anat <- rsfmri::fsl_reorient2std(anatomy) - print("Skull stripping...") + if (!is.null(scale_anat_to)) { + print_msg("Rescaling anatomical data") + ori_anat <- rsfmri::rescale_image(ori_anat, res = 2) + } + + print_msg("Skull stripping...") skull_stripped <- rsfmri::skull_stripping(ori_func, anatomy = ori_anat) - print("Co-registering functional to standard...") - normalise <- rsfmri::register_functonal_to_standard( - skull_stripped[1], anatfile = skull_stripped[2], - standard = std_anat) - - print("Extracting WM and CSF time course...") - wm_timecourse <- rsfmri::fsl_meants(normalise[1], - outfile = file.path(input_dir, "wmtc.txt"), - maskfile = std_white_matter) - csf_timecourse <- rsfmri::fsl_meants(normalise[1], - outfile = file.path(input_dir, "csftc.txt"), - maskfile = std_csf) - if (gsr) { - gs_timecourse <- rsfmri::fsl_meants(normalise[1], - outfile = file.path(input_dir, "gs.txt"), - maskfile = file.path(input_dir, "mean_mask.nii.gz")) + print_msg("Correcting motion...") + motion_cor <- rsfmri::fsl_mcflirt(ori_func, meanvol = TRUE) + mean_motion_reg <- gsub(".nii.gz", "_mean_reg.nii.gz", motion_cor) + motion_file <- gsub(".nii.gz", ".par", motion_cor) + + motion_check <- check_motion(parfile = motion_file, cutoff = motion_cutoff, + spike_cutoff = 1.5, plot = FALSE) + + if (motion_check$exceeds_motion > 3 || + motion_check$exceeds_spike > min( + c(20, nrow(motion_check$data) * 0.05))) { + stop("Excessive motion.") } - if (!is.null(physfile) && !is.null(tr)) { - print("Generating physiological confounds...") - tr_s <- tr / 1000 - # need to sort out the other options to be more generalisable - physio_confounds <- rsfmri::fsl_pnm( - input_files[3], normalise[1], tr = tr_s) - physio_confounds <- file.path(input_dir, physio_confounds) + print_msg("Segmentation...") + pve_wm <- rsfmri::fsl_fast(skull_stripped[2], smooth_extent = 10) + binary_wm <- rsfmri::fsl_maths_binary_segmentation(pve_wm) + + map_mag_pha <- all(c(!is.null(magnitude_map), !is.null(phase_map))) + map_ap_pa <- all(c(!is.null(ap_path), !is.null(pa_path))) + + if (map_mag_pha && map_ap_pa) { + stop("Either use Mag/Pha map or AP/PA, not both.") } - print("Correcting motion...") - motion_cor <- rsfmri::fsl_mcflirt(normalise[1]) + if (map_mag_pha || map_ap_pa) { + print_msg("Generating Field Map...") + if (map_mag_pha) { + fm <- rsfmri::fsl_make_fieldmap(mag_path, pha_path, echo_spacing) + } else if (map_ap_pa) { + fm <- rsfmri::fsl_make_fieldmap_pe(ap_path, pa_path, readout_time) + } + + print_msg("EPI B0 Field Map Correction...") + epi_reg <- rsfmri::fsl_epi_reg( + mean_motion_reg, + anatomical = ori_anat, + anatomical_brain = skull_stripped[2], + anatomical_wm_seg = binary_wm, + fmap = fm$fmap, + magnitude_map = fm$mag_norm, + magnitude_brain = fm$mag, + pedir = "-y", + dwell_time = dwell_time + ) + + print_msg("Co-registering functional to standard...") + reg <- rsfmri::register_functonal_to_standard( + functional = motion_cor, bbr_warp = epi_reg, anatfile = ori_anat, + anat_extracted = skull_stripped[2], standard = std_anat, + standard_extracted = std_brain, + fnirt_config = fnirt_conf) + } else { + print_msg("Co-registering functional to standard (1/2)...") + epi_reg <- rsfmri::fsl_epi_reg( + motion_cor, + anatomical = ori_anat, + anatomical_brain = skull_stripped[2], + anatomical_wm_seg = binary_wm + ) + + print_msg("Co-registering functional to standard (2/2)...") + reg <- rsfmri::register_functonal_to_standard( + functional = epi_reg, anatfile = ori_anat, + anat_extracted = skull_stripped[2], standard = std_anat, + standard_extracted = std_brain, + fnirt_config = fnirt_conf) + } + + print_msg("Extracting WM and CSF time course...") + pve_use <- gsub("_\\d.nii.gz$", "", pve_wm) + timecourses <- rsfmri::extract_csf_and_wm( + reg, pve_basename = pve_use, std_brain = std_anat, + anat_to_mni_warp = file.path(temp_dir, "t1_to_mni_warp.nii.gz") + ) + + if (gsr) { + gs_timecourse <- rsfmri::fsl_meants(motion_cor, + outfile = file.path(temp_dir, "gs.txt"), + maskfile = file.path(temp_dir, "mean_mask.nii.gz")) + } if (!is.null(tr) && !is.null(odd)) { - print("Correcting slice scan time...") + print_msg("Correcting slice scan time...") slice_timed <- rsfmri::fsl_slicetimer( - motion_cor[1], tr = as.character(tr), odd = odd) + reg, tr = as.character(tr), odd = odd) } else { # TODO: Make these a list or so, so we can check the last element. - slice_timed <- motion_cor + slice_timed <- reg } # Build file containing noise to "regress out". - motion_file <- paste0(motion_cor, ".par") if (gsr) { - noise <- noise_file(motion_file, wm_timecourse, csf_timecourse, - gs_timecourse) + noise <- noise_file(motion_file, timecourses[1], timecourses[2], + gs_timecourse, first_deriv = TRUE, + second_deriv = TRUE, square = FALSE) } else { - noise <- noise_file(motion_file, wm_timecourse, csf_timecourse, + noise <- noise_file(motion_file, timecourses[1], timecourses[2], first_deriv = TRUE, second_deriv = TRUE, square = FALSE) } noise <- tibble::as_tibble(noise) - noise_file <- file.path(input_dir, "noise.txt") - utils::write.table(noise, file = file.path(input_dir, "noise.txt"), + noise_file <- file.path(temp_dir, "noise.txt") + utils::write.table(noise, file = file.path(temp_dir, "noise.txt"), row.names = FALSE, col.names = FALSE) # call glm to get rid of all the noise - print("Removing noise...") - if (!is.null(physfile)) { - glm_out <- rsfmri::fsl_glm(slice_timed[1], outtype = "both", - predictors = noise_file, demean = TRUE, - physio_list = physio_confounds) - } else { - glm_out <- rsfmri::fsl_glm(slice_timed[1], outtype = "both", - predictors = noise_file, demean = TRUE) - } - print("Extracting AALv4 time courses") + print_msg("Removing noise...") + glm_out <- rsfmri::fsl_glm(slice_timed, outtype = "both", + predictors = noise_file, demean = TRUE) + + print_msg("Extracting region time courses") labels <- rsfmri::fsl_meants(glm_out[1], label = labels_img) - print(labels[1]) rs_tc <- readr::read_table2(labels[1], col_names = FALSE) rs_tc <- rs_tc[, apply(rs_tc, 2, function(x) !(all(is.na(x))))] - info <- readr::read_table(labels_text, - col_names = c("no", "name"), skip = 1) - colnames(rs_tc) <- info[["name"]] + info <- readr::read_table(labels_text, col_names = c("no", "name")) + colnames(rs_tc)[info$no] <- info$name + null_cols <- apply(rs_tc, 2, function(x) all(x == 0)) + rs_tc[null_cols] <- NULL - print("Filtering...") + print_msg("Filtering...") bandpass <- signal::butter(2, c(0.01, 0.15), "pass") rs_filtered <- scale(apply( @@ -132,16 +238,10 @@ rs_pipeline <- function( rs_filtered <- tibble::as_tibble(rs_filtered) - print("Scrub time series...") - motion_check <- check_motion(parfile = motion_file, cutoff = 1, - plot = FALSE) - to_remove <- c(1:3, motion_check$flag_motion, motion_check$flag_spike) - - if (motion_check$exceeds_motion > 3 || - motion_check$exceeds_spike > min(c(20, nrow(rs_filtered) * 0.05))) { - stop("Excessive motion.") - } - rs_filtered <- rs_filtered[-to_remove, ] + print_msg("Scrub time series...") + rs_filtered[c( + motion_check$flag_motion, motion_check$flag_spike)] <- NA_real_ + rs_filtered <- rs_filtered[-c(1:3), ] rs_list <- list() rs_list$time_course <- rs_filtered @@ -150,13 +250,22 @@ rs_pipeline <- function( x2sym <- rlang::sym("x2") corsym <- rlang::sym("cor") - print("Correlating") - rs_cor <- stats::cor(rs_filtered) + print_msg("Correlating...") + rs_cor <- stats::cor(rs_filtered, use = "complete.obs") rs_cor <- tibble::as_tibble(rs_cor, rownames = "x1") %>% tidyr::gather(!!x2sym, !!corsym, -!!x1sym) rs_list$correlation <- rs_cor - print("Done.") + print_msg("Tidying up...") + file.copy(labels[1], input_dir) + file.copy(glm_out[1], input_dir) + file.copy(motion_file, input_dir) + if (!keep_intermediate) { + unlink(temp_dir, recursive = TRUE) + } + + print_msg("Done.") rs_list } +# nolint end diff --git a/man/compare_measure_group.Rd b/man/compare_measure_group.Rd index 4e2ca80057d58f37874ba98625de3ceaec01ad0b..51aa3146df28ad303df7b7155bfeb46b8060b800 100644 --- a/man/compare_measure_group.Rd +++ b/man/compare_measure_group.Rd @@ -4,7 +4,14 @@ \alias{compare_measure_group} \title{Compare a measure between groups.} \usage{ -compare_measure_group(rawlist, dv, ..., add_name = FALSE, remove = NULL) +compare_measure_group( + rawlist, + dv, + ..., + add_name = FALSE, + remove = NULL, + anova = TRUE +) } \arguments{ \item{rawlist}{A list of graphs split by groups.} @@ -16,6 +23,8 @@ compare_measure_group(rawlist, dv, ..., add_name = FALSE, remove = NULL) \item{add_name}{Whether to add interaction term group x name, default FALSE} \item{remove}{Name of a group to remove} + +\item{anova}{Logical value specifying whether to compute an ANOVA.} } \description{ Compare a measure between groups. diff --git a/man/rs_pipeline.Rd b/man/rs_pipeline.Rd index 96d7ce939a97d14f5ad29f91d74a5ed34ef01a7b..5d7de3184397ec0ba0a56abf7012188e60939a42 100644 --- a/man/rs_pipeline.Rd +++ b/man/rs_pipeline.Rd @@ -7,39 +7,68 @@ rs_pipeline( functional, anatomy, - physfile = NULL, std_anat, - std_white_matter, - std_csf, + std_brain, labels_img, labels_text, + magnitude_map = NULL, + phase_map = NULL, + dwell_time = NULL, + ap_map = NULL, + pa_map = NULL, + readout_time = NULL, + echo_spacing = NULL, + fnirt_conf, gsr = TRUE, tr = NULL, - odd = NULL + odd = NULL, + scale_anat_to = NULL, + motion_cutoff = 1, + keep_intermediate = FALSE ) } \arguments{ \item{functional}{Path to a functional dataset.} -\item{anatomy}{Path to a corresponing anatomical dataset.} - -\item{physfile}{Path to corresponding pyhsiological data.} +\item{anatomy}{Path to a corresponding anatomical dataset.} \item{std_anat}{Path to standard anatomical file.} -\item{std_white_matter}{Path to standard white matter mask.} - -\item{std_csf}{Path to standard CSF mask.} +\item{std_brain}{Path to skull stripped standard file.} \item{labels_img}{Path to atlas image.} \item{labels_text}{Path to atlas label info file.} +\item{magnitude_map}{Path to magnitude B0 map.} + +\item{phase_map}{Path to phase B0 map.} + +\item{dwell_time}{Dwell time in s.} + +\item{ap_map}{Phase map in AP direction.} + +\item{pa_map}{Phase map in PA direction.} + +\item{readout_time}{Readout time for AP/PA fieldmaps.} + +\item{echo_spacing}{Echo spacing time in ms.} + +\item{fnirt_conf}{Name of the FSL fnirt config file to use.} + \item{gsr}{Boolean whether to use global signal regression.} \item{tr}{The TR of the data im ms.} \item{odd}{Bool indicating whether slices were acquired interleaved.} + +\item{scale_anat_to}{Optional value to scale the anatomical image resolution +to.} + +\item{motion_cutoff}{Motion correction cutoff value in deg/mm.} + +\item{keep_intermediate}{Should itermediate files be kept on disk? This is +useful for debugging.} } \value{ A list containing time courses and a pairwise correlation matrix.