#' Apply water reference scaling to a fitting results object to yield metabolite
#' quantities in millimolar (mM) units (mol / kg of tissue water).
#'
#' Details of this method can be found in "Use of tissue water as a
#' concentration reference for proton spectroscopic imaging" by Gasparovic et al
#' MRM 2006 55(6):1219-26. 1.5 Tesla relaxation assumptions are taken from this
#' paper. For 3 Tesla data, relaxation assumptions are taken from "NMR
#' relaxation times in the human brain at 3.0 Tesla" by Wansapura et al J Magn
#' Reson Imaging 1999 9(4):531-8.
#'
#' @param fit_result result object generated from fitting.
#' @param ref_data water reference MRS data object.
#' @param p_vols a numeric vector of partial volumes expressed as percentages.
#' For example, a voxel containing 100% white matter tissue would use :
#' p_vols = c(WM = 100, GM = 0, CSF = 0).
#' @param te the MRS TE in seconds.
#' @param tr the MRS TR in seconds.
#' @param ... additional arguments to get_td_amp function.
#' @return A \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_molal_pvc <- function(fit_result, ref_data, p_vols, te, tr, ...){
# if (!identical(dim(fit_result$data$data)[2:6], dim(ref_data$data)[2:6])) {
# stop("Mismatch between fit result and reference data dimensions.")
# }
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
B0 <- round(fit_result$data$ft / 42.58e6, 1)
corr_factor <- get_corr_factor(te, tr, B0, p_vols[["GM"]], p_vols[["WM"]],
p_vols[["CSF"]])
amp_cols <- fit_result$amp_cols
w_amp <- as.numeric(get_td_amp(ref_data, ...))
# if there is only one water ref amp, extend if needed, eg for fMRS
if ((length(w_amp) == 1) & (nrow(fit_result$res_tab) > 1)) {
w_amp <- rep(w_amp, nrow(fit_result$res_tab))
}
if (length(w_amp) != nrow(fit_result$res_tab)) {
stop("Mismatch between fit result and reference data.")
}
fit_result$res_tab$w_amp <- w_amp
fit_result$res_tab$GM_vol <- p_vols[["GM"]]
fit_result$res_tab$WM_vol <- p_vols[["WM"]]
fit_result$res_tab$CSF_vol <- p_vols[["CSF"]]
if ("Other" %in% names(p_vols)) {
fit_result$res_tab$Other_vol <- p_vols[["Other"]]
}
fit_result$res_tab$GM_frac <- p_vols[["GM"]] /
(p_vols[["GM"]] + p_vols[["WM"]])
fit_result$res_tab_unscaled$GM_vol <- p_vols[["GM"]]
fit_result$res_tab_unscaled$WM_vol <- p_vols[["WM"]]
fit_result$res_tab_unscaled$CSF_vol <- p_vols[["CSF"]]
if ("Other" %in% names(p_vols)) {
fit_result$res_tab_unscaled$Other_vol <- p_vols[["Other"]]
}
fit_result$res_tab_unscaled$GM_frac <- p_vols[["GM"]] /
(p_vols[["GM"]] + p_vols[["WM"]])
# append tables with %GM, %WM, %CSF and %Other
pvc_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, pvc_cols] <- fit_result$res_tab[, pvc_cols] *
corr_factor / w_amp
return(fit_result)
}
#' Apply water reference scaling to a fitting results object to yield metabolite
#' quantities in millimolar (mM) units (mol / kg of tissue water).
#'
#' Note, this function assumes the volume contains a homogeneous voxel, eg pure
#' WM, GM or CSF. Also note that in the case of a homogeneous voxel the
#' relative densities of MR-visible water (eg GM=0.78, WM=0.65, and CSF=0.97)
#' cancel out and don't need to be considered. Use scale_amp_molal_pvc for
#' volumes containing multiple compartments. Details of this method can be
#' found in "Use of tissue water as a concentration reference for proton
#' spectroscopic imaging" by Gasparovic et al MRM 2006 55(6):1219-26.
#'
#' @param fit_result result object generated from fitting.
#' @param ref_data water reference MRS data object.
#' @param te the MRS TE in seconds.
#' @param tr the MRS TR in seconds.
#' @param water_t1 assumed water T1 value.
#' @param water_t2 assumed water T2 value.
#' @param metab_t1 assumed metabolite T1 value.
#' @param metab_t2 assumed metabolite T2 value.
#' @param ... additional arguments to get_td_amp function.
#' @return A \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_molal <- function(fit_result, ref_data, te, tr, water_t1, water_t2,
metab_t1, metab_t2, ...){
# if (!identical(dim(fit_result$data$data)[2:6], dim(ref_data$data)[2:6])) {
# stop("Mismatch between fit result and reference data dimensions.")
# }
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
R_water <- exp(-te / water_t2) * (1.0 - exp(-tr / water_t1))
R_metab <- exp(-te / metab_t2) * (1.0 - exp(-tr / metab_t1))
water_conc <- 55510.0
corr_factor <- R_water / R_metab * water_conc
amp_cols <- fit_result$amp_cols
w_amp <- as.numeric(get_td_amp(ref_data, ...))
# if there is only one water ref amp, extend if needed, eg for fMRS
if ((length(w_amp) == 1) & (nrow(fit_result$res_tab) > 1)) {
w_amp <- rep(w_amp, nrow(fit_result$res_tab))
}
if (length(w_amp) != nrow(fit_result$res_tab)) {
stop("Mismatch between fit result and reference data.")
}
fit_result$res_tab$w_amp <- w_amp
pvc_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, pvc_cols] <- fit_result$res_tab[, pvc_cols] *
corr_factor / w_amp
return(fit_result)
}
#' Apply water reference scaling to a fitting results object to yield metabolite
#' quantities in millimolar (mM) units (mol / Litre of tissue).
#'
#' See the LCModel manual (section 10.2) on water-scaling for details on the
#' assumptions and relevant references. Use this type of concentration scaling
#' to compare fit results with LCModel and TARQUIN defaults. Otherwise
#' scale_amp_molal_pvc is generally the preferred method.
#'
#' @param fit_result a result object generated from fitting.
#' @param ref_data water reference MRS data object.
#' @param w_att water attenuation factor (default = 0.7). Assumes water T2 of
#' 80ms and a TE = 30 ms. exp(-30ms / 80ms) ~ 0.7.
#' @param w_conc assumed water concentration (default = 35880). Default value
#' corresponds to typical white matter. Set to 43300 for gray matter, and 55556
#' for phantom measurements.
#' @param ... additional arguments to get_td_amp function.
#' @return a \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_molar <- function(fit_result, ref_data, w_att = 0.7, w_conc = 35880,
...) {
# if (!identical(dim(fit_result$data$data)[2:6], dim(ref_data$data)[2:6])) {
# stop("Mismatch between fit result and reference data dimensions.")
# }
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
w_amp <- as.numeric(get_td_amp(ref_data, ...))
# if there is only one water ref amp, extend if needed, eg for fMRS
if ((length(w_amp) == 1) & (nrow(fit_result$res_tab) > 1)) {
w_amp <- rep(w_amp, nrow(fit_result$res_tab))
}
if (length(w_amp) != nrow(fit_result$res_tab)) {
stop("Mismatch between fit result and reference data.")
}
fit_result$res_tab$w_amp <- w_amp
amp_cols <- fit_result$amp_cols
ws_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, ws_cols] <- (fit_result$res_tab[, ws_cols] * w_att *
w_conc / w_amp)
fit_result
}
#' Scale metabolite amplitudes as a ratio to the unsuppressed water amplitude.
#' @param fit_result a result object generated from fitting.
#' @param ref_data a water reference MRS data object.
#' @param ... additional arguments to get_td_amp function.
#' @return a \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_water_ratio <- function(fit_result, ref_data, ...) {
# if (!identical(dim(fit_result$data$data)[2:6], dim(ref_data$data)[2:6])) {
# stop("Mismatch between fit result and reference data dimensions.")
# }
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
w_amp <- as.numeric(get_td_amp(ref_data, ...))
if (length(w_amp) != nrow(fit_result$res_tab)) {
stop("Mismatch between fit result and reference data.")
}
# if there is only one water ref amp, extend if needed, eg for fMRS
if ((length(w_amp) == 1) & (nrow(fit_result$res_tab) > 1)) {
w_amp <- rep(w_amp, nrow(fit_result$res_tab))
}
fit_result$res_tab$w_amp <- w_amp
amp_cols <- fit_result$amp_cols
ws_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, ws_cols] <- fit_result$res_tab[, ws_cols] / w_amp
fit_result
}
#' Scale fitted amplitudes to a ratio of signal amplitude.
#' @param fit_result a result object generated from fitting.
#' @param name the signal name to use as a denominator (usually, "tCr" or
#' "tNAA").
#' @param use_mean_value scales the result by the mean of the signal when set to
#' TRUE.
#' @return a \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_ratio <- function(fit_result, name, use_mean_value = FALSE) {
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
ratio_amp <- as.numeric(fit_result$res_tab[[name]])
if (use_mean_value) ratio_amp <- mean(ratio_amp)
amp_cols <- fit_result$amp_cols
ws_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, ws_cols] <- fit_result$res_tab[, ws_cols] / ratio_amp
fit_result
}
#' Scale fitted amplitudes to a ratio of signal amplitude.
#' @param fit_result a result object generated from fitting.
#' @param value the number use as a denominator.
#' @return a \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_ratio_value <- function(fit_result, value) {
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
ratio_amp <- value
amp_cols <- fit_result$amp_cols
ws_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, ws_cols] <- fit_result$res_tab[, ws_cols] / ratio_amp
fit_result
}
get_corr_factor <- function(te, tr, B0, gm_vol, wm_vol, csf_vol) {
# Correction factor calcualted according to the method of Gasparovic et al (MRM 55:1219-1226 2006)
# NOTE - gives concs as Mol/kg of water NOT Mol/liter of tissue like default LCM/TQN analysis.
if ((B0 == 3.0) | (B0 == 2.9)) {
# Wanasapura values given in Harris paper
t1_gm <- 1.331
t2_gm <- 0.110
t1_wm <- 0.832
t2_wm <- 0.0792
t1_csf <- 3.817
t2_csf <- 0.503
t1_metab <- 1.15
t2_metab <- 0.3
} else if (B0 == 1.5) {
# values from Gasparovic 2006 MRM paper
t1_gm <- 1.304
t2_gm <- 0.093
t1_wm <- 0.660
t2_wm <- 0.073
t1_csf <- 2.93
t2_csf <- 0.23
t1_metab <- 1.15
t2_metab <- 0.3
} else {
stop("Error. Relaxation values not available for this field strength.")
}
# MR-visible water densities
gm_vis <- 0.78
wm_vis <- 0.65
csf_vis <- 0.97
# molal concentration (moles/gram) of MR-visible water
water_conc <- 55510.0
# fractions of water attributable to GM, WM and CSF
f_gm <- gm_vol * gm_vis /
(gm_vol * gm_vis + wm_vol * wm_vis + csf_vol * csf_vis)
f_wm <- wm_vol * wm_vis /
(gm_vol * gm_vis + wm_vol * wm_vis + csf_vol * csf_vis)
f_csf <- csf_vol * csf_vis /
(gm_vol * gm_vis + wm_vol * wm_vis + csf_vol * csf_vis)
# This might give the result in Mol/kg?
# f_gm <- gm_vol * gm_vis / ( gm_vol + wm_vol + csf_vol )
# f_wm <- wm_vol * wm_vis / ( gm_vol + wm_vol + csf_vol )
# f_csf <- csf_vol * csf_vis / ( gm_vol + wm_vol + csf_vol )
# Relaxtion attenuation factors
R_h2o_gm <- exp(-te / t2_gm) * (1.0 - exp(-tr / t1_gm))
R_h2o_wm <- exp(-te / t2_wm) * (1.0 - exp(-tr / t1_wm))
R_h2o_csf <- exp(-te / t2_csf) * (1.0 - exp(-tr / t1_csf))
R_metab <- exp(-te / t2_metab) * (1.0 - exp(-tr / t1_metab))
corr_factor <- ((f_gm * R_h2o_gm + f_wm * R_h2o_wm + f_csf * R_h2o_csf) /
((1 - f_csf) * R_metab)) * water_conc
return(corr_factor)
}
#' Convert default LCM/TARQUIN concentration scaling to molal units with partial
#' volume correction.
#' @param fit_result a \code{fit_result} object to apply partial volume
#' correction.
#' @param p_vols a numeric vector of partial volumes expressed as percentages.
#' For example, a voxel containing 100% white matter tissue would use :
#' p_vols = c(WM = 100, GM = 0, CSF = 0).
#' @param te the MRS TE.
#' @param tr the MRS TR.
#' @return a \code{fit_result} object with a rescaled results table.
#' @export
scale_amp_molar2molal_pvc <- function(fit_result, p_vols, te, tr){
# check if res_tab_unscaled exists, and if not create it
if (is.null(fit_result$res_tab_unscaled)) {
fit_result$res_tab_unscaled <- fit_result$res_tab
} else {
fit_result$res_tab <- fit_result$res_tab_unscaled
}
B0 <- round(fit_result$data$ft / 42.58e6,1)
corr_factor <- get_corr_factor(te, tr, B0, p_vols[["GM"]], p_vols[["WM"]],
p_vols[["CSF"]])
amp_cols <- fit_result$amp_cols
default_factor <- 35880 * 0.7
fit_result$res_tab$GM_vol <- p_vols[["GM"]]
fit_result$res_tab$WM_vol <- p_vols[["WM"]]
fit_result$res_tab$CSF_vol <- p_vols[["CSF"]]
if ("Other" %in% names(p_vols)) {
fit_result$res_tab$Other_vol <- p_vols[["Other"]]
}
# append tables with %GM, %WM, %CSF and %Other
pvc_cols <- 6:(5 + amp_cols * 2)
fit_result$res_tab[, pvc_cols] <- fit_result$res_tab[, pvc_cols] /
default_factor * corr_factor
return(fit_result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.