Nothing
# TODO add option to explicitly specify the type of metabolite scaling to perform
# include "auto" version which guesses based on the input
#' Standard SVS 1H brain analysis pipeline.
#' @param metab filepath or mrs_data object containing MRS metabolite data.
#' @param basis basis set object to use for analysis.
#' @param w_ref filepath or mrs_data object containing MRS water reference data.
#' @param mri_seg filepath or nifti object containing segmented MRI data.
#' @param mri filepath or nifti object containing anatomical MRI data.
#' @param output_dir directory path to output fitting results.
#' @param extra data.frame with one row containing additional information to be
#' attached to the fit results table.
#' @param decimate option to decimate the input data by a factor of two. The
#' default value of NULL does not perform decimation unless the spectral width
#' is greater than 20 PPM.
#' @param rats_corr option to perform rats correction, defaults to TRUE.
#' @param ecc option to perform water reference based eddy current correction,
#' defaults to FALSE.
#' @param comb_dyns option to combine dynamic scans, defaults to TRUE.
#' @param hsvd_filt option to apply hsvd water removal, defaults to FALSE.
#' @param scale_amps option to scale metabolite amplitude estimates, defaults to
#' TRUE.
#' @param te metabolite mrs data echo time in seconds.
#' @param tr metabolite mrs data repetition time in seconds.
#' @param preproc_only only perform the preprocessing steps and omit fitting.
#' The preprocessed metabolite data will be returned in this case.
#' @return a fit_result or mrs_data object depending on the preproc_only option.
#' @param method analysis method to use, see fit_mrs help.
#' @param opts options to pass to the analysis method.
#' @export
svs_1h_brain_analysis <- function(metab, basis = NULL, w_ref = NULL,
mri_seg = NULL, mri = NULL, output_dir = NULL,
extra = NULL, decimate = NULL,
rats_corr = TRUE, ecc = FALSE,
comb_dyns = TRUE, hsvd_filt = FALSE,
scale_amps = TRUE, te = NULL, tr = NULL,
preproc_only = FALSE, method = "ABFIT",
opts = NULL) {
if (!preproc_only & is.null(basis)) stop("basis argument not set")
# read the data file if not already an mrs_data object
if (class(metab)[[1]] != "mrs_data") metab <- read_mrs(metab)
# TODO check for combined metab + water ref. data, eg GE
# remove fs_path types as they cause problems with comb_fit_list_fit_tables
extra[] = lapply(extra, as.character)
# create the output dir if it doesn't exist
if (is.def(output_dir)) {
if(!dir.exists(output_dir)) dir.create(output_dir)
}
# read the ref data file if not already an mrs_data object
if (is.def(w_ref) & (class(w_ref)[[1]] != "mrs_data")) {
w_ref <- read_mrs(w_ref)
}
if (is.def(mri) & (!("niftiImage" %in% class(mri)))) {
mri <- readNifti(mri)
}
if (is.def(mri_seg) & (!("niftiImage" %in% class(mri_seg)))) {
mri_seg <- readNifti(mri_seg)
}
if (is.def(mri)) RNifti::orientation(mri) <- "RAS"
if (is.def(mri_seg)) RNifti::orientation(mri_seg) <- "RAS"
# combine coils if needed
if (Ncoils(metab) > 1) {
coil_comb_res <- comb_coils(metab, w_ref)
if (is.null(w_ref)) {
metab <- coil_comb_res
} else {
metab <- coil_comb_res$metab
w_ref <- coil_comb_res$ref
}
}
# if the spectral width exceeds 20 PPM, then it should be ok
# to decimate the signal to improve analysis speed
if (is.null(decimate) & ((fs(metab) / metab$ft * 1e6) > 20)) {
decimate <- TRUE
} else {
decimate <- FALSE
}
if (decimate) {
metab <- decimate_mrs_fd(metab)
if (!is.null(w_ref)) w_ref <- decimate_mrs_fd(w_ref)
}
# rats
if (rats_corr & (Ndyns(metab)> 1)) metab <- rats(metab)$corrected
# TODO plot of shifts?
# eddy current correction
if (ecc & (!is.null(w_ref))) metab <- ecc(metab, w_ref)
# combine dynamic scans
if (comb_dyns) metab <- mean_dyns(metab)
# HSVD residual water removal
if (hsvd_filt) metab <- hsvd_filt(metab)
# end here if we're only interested in the preprocessed data
if (preproc_only) return(metab)
fit_res <- fit_mrs(metab, basis = basis, extra = extra, method = method,
opts = opts)
# plot the fit result and output csv
if (is.def(output_dir)) {
grDevices::pdf(file.path(output_dir, "fit_plot.pdf"))
plot(fit_res)
grDevices::dev.off()
utils::write.csv(fit_res$res_tab, file.path(output_dir, "fit_res.csv"))
}
# plot the voxel location on the mri
if (is.def(mri) & is.def(output_dir)) {
# generate the svs voi in the segmented image space
voi <- get_svs_voi(metab, mri)
plot_voi_overlay(mri, voi, file.path(output_dir, "vox_plot.png"))
}
# the calculate water reference frequency, linewidth and water suppression
# efficiency
if (is.def(w_ref)) {
# phase correct based on the phase of the first data point
w_ref <- fp_phase_correct(w_ref)
w_info_re <- peak_info(w_ref, xlim = c(6, 4), mode = "real")
w_info_mod <- peak_info(w_ref, xlim = c(6, 4), mode = "mod")
fit_res$res_tab$water_freq_mod_ppm <- w_info_mod$freq_ppm
fit_res$res_tab$water_freq_real_ppm <- w_info_re$freq_ppm
fit_res$res_tab$water_fwhm_mod_ppm <- w_info_mod$fwhm_ppm
fit_res$res_tab$water_fwhm_real_ppm <- w_info_re$fwhm_ppm
# search for highest peak in the water suppressed data 0.1 ppm either side
# of the water frequency
water_freq <- w_info_mod$freq_ppm
resid_water_peak_mod <- peak_info(metab, xlim = c(water_freq + 0.1,
water_freq - 0.1), mode = "mod")
fit_res$res_tab$water_sup_efficiency_percent <- resid_water_peak_mod$height /
w_info_mod$height * 100
}
if (scale_amps) {
if (is.def(w_ref) & is.def(mri_seg)) {
if (is.null(te)) stop("te not given, amplitude scaling failed")
if (is.null(tr)) stop("tr not given, amplitude scaling failed")
# generate the svs voi in the segmented image space
voi_seg <- get_svs_voi(metab, mri_seg)
# calculate partial volumes
seg <- get_voi_seg(voi_seg, mri_seg)
# do pvc
fit_res <- scale_amp_molal_pvc(fit_res, w_ref, seg, te, tr)
if (is.def(output_dir)) {
plot_voi_overlay_seg(mri_seg, voi_seg, file.path(output_dir,
"vox_seg_plot.png"))
}
} else if (is.def(w_ref) & !is.def(mri_seg)) {
# do straight w scaling default LCM style
fit_res <- scale_amp_molar(fit_res, w_ref)
} else {
# scale to tCr
fit_res <- scale_amp_ratio(fit_res, "tCr")
}
}
return(fit_res)
}
# TODO
# what if metab fname results in metab + ref file, eg GE data?
# options to add, fit_opts, fit_method, format (GE, Siemens etc)
#' Batch interface to the standard SVS 1H brain analysis pipeline.
#' @param metab_list list of file paths or mrs_data objects containing MRS
#' metabolite data.
#' @param w_ref_list list of file paths or mrs_data objects containing MRS
#' water reference data.
#' @param mri_seg_list list of file paths or nifti objects containing segmented
#' MRI data.
#' @param mri_list list of file paths or nifti objects containing anatomical
#' MRI data.
#' @param output_dir_list list of directory paths to output fitting results.
#' @param extra a data frame with the same number of rows as metab_list,
#' containing additional information to be attached to the fit results table.
#' @param ... additional options to be passed to the svs_1h_brain_analysis
#' function.
#' @return a list of fit_result objects.
#' @export
svs_1h_brain_batch_analysis <- function(metab_list, w_ref_list = NULL,
mri_seg_list = NULL, mri_list = NULL,
output_dir_list = NULL, extra = NULL,
...) {
# check input is sensible
metab_n <- length(metab_list)
if (is.def(w_ref_list) & length(w_ref_list) != metab_n) {
stop("Incorrect number of w_ref_list items.")
}
if (is.def(mri_list) & length(mri_list) != metab_n) {
stop("Incorrect number of mri_list items.")
}
if (is.def(mri_seg_list) & length(mri_seg_list) != metab_n) {
stop("Incorrect number of mri_seg_list items.")
}
if (is.def(output_dir_list) & length(output_dir_list) != metab_n) {
stop("Incorrect number of output_dir_list items.")
}
if (is.def(extra) & nrow(extra) != metab_n) {
stop("Incorrect number of rows in extra.")
} else {
# split into a list suitable for mapply
extra <- split(extra, seq(nrow(extra)))
}
# check file paths exist TODO write a function to return bad paths
if (is.null(w_ref_list)) w_ref_list <- vector("list", metab_n)
if (is.null(mri_seg_list)) mri_seg_list <- vector("list", metab_n)
if (is.null(mri_list)) mri_list <- vector("list", metab_n)
if (is.null(output_dir_list)) output_dir_list <- vector("list", metab_n)
if (is.null(extra)) extra <- vector("list", metab_n)
fit_list <- mapply(svs_1h_brain_analysis, metab = metab_list,
w_ref = w_ref_list, mri_seg = mri_seg_list, mri = mri_list,
output_dir = output_dir_list, extra = extra,
MoreArgs = list(...), SIMPLIFY = FALSE, USE.NAMES = FALSE)
return(fit_list)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.