Merge pull request #17 from ropensci/devel

Update to v0.3.0

Author: eebrown

.Rbuildignore

@@ -5,6 +5,7 @@
 ^\.travis\.yml$
 ^Meta$
 ^doc$
+.buildconfig
 .git
 scratch
 .gitignore

.gitignore

@@ -3,5 +3,6 @@ Meta
 analysis
 scratch
 .Rhistory
+.buildconfig
 ..Rcheck
 man/.Rhistory

CRAN-RELEASE

@@ -1,7 +1,7 @@
 Package: tacmagic
 Type: Package
 Title: Positron Emission Tomography Time-Activity Curve Analysis
-Version: 0.2.1
+Version: 0.3.0
 Authors@R: c(person("Eric", "Brown", 
                     role = c("aut", "cre"),
                     email = "[email protected]",
@@ -44,5 +44,6 @@ Suggests:
     testthat,
     knitr,
     rmarkdown,
-    covr
+    covr,
+    vdiffr
 VignetteBuilder: knitr

DESCRIPTION

@@ -1,6 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(as.data.frame,tac)
+S3method(change_units,numeric)
+S3method(change_units,tac)
 S3method(plot,ref_Logan)
 S3method(plot,tac)
 S3method(summary,tac)
@@ -10,6 +12,7 @@ export(as.tac)
 export(batch_load)
 export(batch_tm)
 export(batch_voistat)
+export(change_units)
 export(cutoff_aiz)
 export(dvr)
 export(load_tac)
@@ -21,9 +24,11 @@ export(roi_ham_pib)
 export(roi_ham_stand)
 export(save_tac)
 export(split_pvc)
+export(suv)
 export(suvr)
 export(suvr_auc)
 export(tac_roi)
+export(tac_suv)
 importFrom(R.matlab,readMat)
 importFrom(grDevices,rainbow)
 importFrom(graphics,abline)

NAMESPACE

@@ -1,11 +1,21 @@
-# tacmagic 0.2.1
+# tacmagic News
 
-## Minor changes
+## tacmagic 0.3.0 (2019-06-06)
+
+### Major changes
+
+- Radioactivity unit conversion can now be done for tac objects and numeric objects
+- Conversion of a tac object to standardized uptake value (SUV) by controlling for weight and injected dose
+- Vignette updated with new features
+
+## tacmagic 0.2.1 (2019-03-06)
+
+### Minor changes
 
 This version has no changes to features but addresses minor, mostly stylistic, issues for the initial CRAN release.
 
-# tacmagic 0.2.0
+## tacmagic 0.2.0 (2019-02-26)
 
-## Major changes
+### Major changes
 
 This is the first release as a complete R package, and the version that has undergone open review with rOpenSci. The main features of tacmagic 0.2.0 are implemented and tested: loading PET time activity curve (tac) data from multiple formats, merging ROIs weighted for volume, calculating binding potential models including SUVR and DVR, basic plotting, and calculation of cut-off values. As tacmagic is still a new package and testing to date is based on few example files, it is essential to verify that all functions are behaving as expected when applied to your own data.

NEWS.md

@@ -0,0 +1,143 @@
+##################################
+## tacmagic - PET Analysis in R ##
+## SUV.R                        ##
+## (C) Eric E. Brown  2018-2019 ##
+## Beta version--check all work ##
+##################################
+
+#' Calculate average SUV over time window, or maximum SUV
+#'
+#' Calculate the standardized uptake value (SUV) from a tac object, the
+#' participant's weight, and the tracer dose. These values may be in the tac
+#' object or manually supplied. The weight must be in kg, and the tracer units
+#' must be specified. The dose is converted to MBq, the tac is converted to
+#' kBq/cc, and the final SUV units are thus in g/cc. Aside from the tac object,
+#' the remaining parameters should be left NULL if the required data is in the
+#' tac object attributes (as can be done with batch_load()).
+#'
+#'@export
+#'@param tac time-activity curve object (decay-corrected)
+#'@param SUV_def vector of start times for window for SUV weighted average, or
+#'       alternatively, "max" for the maximum ROI SUV value
+#'@param dose the injected tracer dose
+#'@param dose_unit unit of tracer dose (e.g. "MBq", "kBq", "mCi"...)
+#'@param weight_kg the participant's weight in kg
+#'@param ... When called from tm_batch, unused parameters may be supplied
+#'@return table of SUV values
+#'@family SUV functions
+#'@examples
+#' f <- system.file("extdata", "AD06.tac", package="tacmagic")
+#' fv <- system.file("extdata", "AD06_TAC.voistat", package="tacmagic")
+#' AD06_tac <- load_tac(f, format="PMOD")
+#' AD06_volume <- load_vol(fv, format="voistat")
+#' AD06 <- tac_roi(tac=AD06_tac, volumes=AD06_volume, ROI_def=roi_ham_pib(),
+#'                 merge=FALSE, PVC=FALSE)
+#' # dose and weight are fabricated for the example
+#' AD06_suvmax <- suv(AD06, "max", dose = 9.0, dose_unit = "mCi",
+#'                      weight_kg = 70)
+#' AD06_suv <- suv(AD06, c(3000, 3300, 3600), dose = 9.0, dose_unit = "mCi",
+#'                   weight_kg = 70)
+suv <- function(tac, SUV_def, dose=NULL, dose_unit=NULL, weight_kg=NULL, ...) {
+
+  if (!(any(is.null(attributes(tac)$tracer_dose),
+           is.null(attributes(tac)$weight_kg),
+           is.null(attributes(tac)$dose_unit)))) {
+
+    dose <- attributes(tac)$tracer_dose
+    weight_kg <- attributes(tac)$weight_kg
+    dose_unit <- attributes(tac)$dose_unit
+  }
+
+  validate_suv_params(tac=tac, dose=dose, dose_unit=dose_unit,
+                      weight_kg=weight_kg)
+
+  if ( (SUV_def[1] != "max") && (!all(SUV_def %in% tac$start)) ) {
+    stop("The SUV definition must be valid tac start times or \"max\"")
+  }
+
+  tac <- tac_suv(tac=tac, dose=dose, dose_unit=dose_unit, weight_kg=weight_kg)
+
+  SUVtable <- new_table(tac, "SUV")
+
+  if (SUV_def[1] == "max") {
+    for (ROI in names(tac)[-(1:2)]) {
+      SUVtable[ROI, "SUV"] <- max(tac[,ROI])
+    }
+  } else {
+
+    frames <- match(SUV_def, tac$start)
+    frame_weights <- tac$end[frames] - tac$start[frames]
+
+    for (ROI in names(tac)[-(1:2)]) {
+      SUVtable[ROI, "SUV"] <- weighted.mean(tac[frames,][,ROI], frame_weights)
+    }
+
+  }
+
+ return(SUVtable)
+}
+
+#' Calculate SUV from TAC
+#'
+#' Calculate the standardized uptake value (SUV) time-activity curve from a tac
+#' object, the participant's weight, and the tracer dose. The weight must be in
+#' kg, and the tracer dose must be specified. The dose is converted to MBq, the
+#' tac is converted to kBq/cc, and the final SUV units are thus in g/cc. Aside
+#' from the tac object, the remaining parameters should be left NULL if the
+#' required data is in the tac object attributes (as can be done with
+#' batch_load().
+#'
+#'@export
+#'@param tac time-activity curve object (decay-corrected)
+#'@param dose the injected tracer dose
+#'@param dose_unit unit of tracer dose (e.g. "MBq", "kBq", "mCi"...)
+#'@param weight_kg the participant's weight in kg
+#'@return tac object with SUV values
+#'@family SUV functions
+#'@examples
+#' f <- system.file("extdata", "AD06.tac", package="tacmagic")
+#' fv <- system.file("extdata", "AD06_TAC.voistat", package="tacmagic")
+#' AD06_tac <- load_tac(f, format="PMOD")
+#' AD06_volume <- load_vol(fv, format="voistat")
+#' AD06 <- tac_roi(tac=AD06_tac, volumes=AD06_volume, ROI_def=roi_ham_pib(),
+#'                 merge=FALSE, PVC=FALSE)
+#' # dose and weight are fabricated for the example
+#' AD06_suv <- tac_suv(AD06, dose = 9.0, dose_unit = "mCi", weight_kg = 70)
+tac_suv <- function(tac, dose=NULL, dose_unit=NULL, weight_kg=NULL) {
+
+  if (!(any(is.null(attributes(tac)$tracer_dose),
+           is.null(attributes(tac)$weight_kg),
+           is.null(attributes(tac)$dose_unit)))) {
+
+    dose <- attributes(tac)$tracer_dose
+    weight_kg <- attributes(tac)$weight_kg
+    dose_unit <- attributes(tac)$dose_unit
+  }
+
+  validate_suv_params(tac, dose, dose_unit, weight_kg)
+
+  if (attributes(tac)$activity_unit != "kBq/cc") {
+    tac <- change_units(tac, to_unit = "kBq/cc")
+  }
+
+  if (dose_unit != "MBq") dose <- change_units(dose, to_unit = "MBq", dose_unit)
+
+  tac[,-(1:2)] <- tac[,-(1:2)] / (dose / weight_kg)
+
+  attributes(tac)$activity_unit <- "g/mL"
+
+  if(!validate_tac(tac)) stop("The resulting tac object did not validate.")
+
+  return(tac)
+}
+
+#'@noRd
+validate_suv_params <- function(tac, dose, dose_unit, weight_kg) {
+
+  if (!validate_tac(tac)) stop("Supplied tac file did not validate.")
+  if (!is.numeric(dose)) stop("Dose must be numeric.")
+  if (!(dose_unit %in% names(get_activity_unit_index()))) stop("Bad dose unit.")
+  if (!is.numeric(weight_kg)) stop("Weight must be numeric.")
+
+  return(TRUE)
+}

R/SUV.R

@@ -13,7 +13,8 @@
 #' @noRd
 model_definitions <- function() {
   return(c(SUVR=suvr, 
-           Logan=DVR_all_ref_Logan
+           Logan=DVR_all_ref_Logan,
+           SUV=suv
          ))
 }
 

R/batch_models.R

@@ -82,6 +82,10 @@ batch_tm <- function(all_tacs, models, custom_model=NULL, ...) {
 #'@param ROI_def Object that defines combined ROIs, see ROI_definitions.R
 #'@param PVC For PVC, true where the data is stored as _C in same tac file
 #'@param merge Passes value to tac_roi(); T to also incl. original atomic ROIs
+#'@param tracer_dose optionally, a vector of tracer doses (in the same order as
+#'       participants), for SUV
+#'@param dose_unit if tracer_dose is specified, note the unit (e.g "MBq")
+#'@param weight_kg optionally, a vector of participant weights in kg, for SUV
 #'@return A list of data.frames, each is a participant's TACs
 #'@family Batch functions
 #'@examples
@@ -94,7 +98,8 @@ batch_tm <- function(all_tacs, models, custom_model=NULL, ...) {
 batch_load <- function(participants, dir="", tac_file_suffix=".tac",
                        tac_format="PMOD", roi_m=FALSE, PVC=NULL, 
                        vol_file_suffix=NULL, vol_format=NULL, 
-                       merge=NULL, ROI_def=NULL) {
+                       merge=NULL, ROI_def=NULL, tracer_dose=NULL,
+                       dose_unit=NULL, weight_kg=NULL) {
 
   if (!roi_m) {
     if (!all(c(is.null(vol_format), is.null(vol_file_suffix), is.null(ROI_def), 
@@ -111,9 +116,34 @@ batch_load <- function(participants, dir="", tac_file_suffix=".tac",
 
   names(r) <- participants
 
+  if (!is.null(c(tracer_dose, dose_unit, weight_kg))) {
+    validate_batch_suv_data(participants, tracer_dose, dose_unit, weight_kg)
+    for (i in seq_along(r)) {
+      attributes(r[[i]])$tracer_dose <- tracer_dose[[i]]
+      attributes(r[[i]])$dose_unit <- dose_unit
+      attributes(r[[i]])$weight_kg <- weight_kg[[i]]
+    }
+  }
+
   return(r)
 }
 
+#'@noRd
+validate_batch_suv_data <- function(participants, tracer_dose, dose_unit,
+                                    weight_kg) {
+  if (length(participants) != length(tracer_dose)) stop("Need tracer dose for
+                                                         each participant")
+  if (!is.numeric(tracer_dose)) stop("Dose must be numeric.")
+
+  if (length(participants) != length(weight_kg)) stop("Need weight for each
+                                                       participant")
+  if (!is.numeric(weight_kg)) stop("Weight must be numeric.")
+
+  if (length(dose_unit) != 1) stop("Only provide 1 dose unit.")
+  if (!(dose_unit %in% names(get_activity_unit_index()))) stop("Bad dose unit.")
+
+}
+
 #' Obtain values from voistat files (using load_voistat() for a batch.
 #'
 #' For a vector of participant IDs and correspondingly named .voistat files,

R/batches.R

@@ -188,3 +188,4 @@ load_voistat <- function(filename, ROI_def=NULL, model="VALUE") {
     names(VALUEtable) <- model
     return(VALUEtable)
 }
+

R/loading.R

@@ -33,7 +33,8 @@ validate_tac <- function(tac) {
     status <- FALSE
   }
 
-  if (!(attributes(tac)$activity_unit %in% c("kBq/cc", "nCi/cc", "Bq/cc"))) {
+  if (!(attributes(tac)$activity_unit %in% c("kBq/cc", "nCi/cc", "Bq/cc", 
+                                             "g/mL"))) {
     message("TAC data missing attribute activity_unit")
     status <- FALSE
   }

R/loading_utilities.R

@@ -0,0 +1,81 @@
+##################################
+## tacmagic - PET Analysis in R ##
+## units.R                      ##
+## (C) Eric E. Brown  2018-2019 ##
+## Beta version--check all work ##
+##################################
+
+
+#' @noRd
+get_activity_unit_index <- function() {
+  # The equivalent to 1 kBq
+  unit_index <- list(`Bq`  = 0.001,
+                     `kBq` = 1,
+                     `MBq` = 1000,
+                     `nCi` = 0.037,
+                     `uCi` = 37,
+                     `mCi` = 37000,
+                     `Ci`  = 37000000)
+  return(unit_index)
+}
+
+#' Convert radioactivity units
+#'
+#' Change the radioactivity units of a tac or numeric object to the specified
+#' desired units (e.g. Bq, kBq, MBq, nCi, uCi, mCi, Ci). For convenience, if the
+#' unit is per volume ("x/cc" or "x/mL"), the "/cc" part is ignored for the 
+#' conversion.
+#'
+#'@export 
+#'@param x time-activity curve or numeric object
+#'@param to_unit the desired unit (e.g. "kBq")
+#'@param from_unit not used for tac object (it is in the tac object), but for 
+#' numeric objects, must be specified (e.g. "nCi")
+#'@return the converted object, same type as x
+#'@family unit functions
+#'@examples
+#' f <- system.file("extdata", "AD06.tac", package="tacmagic")
+#' AD06_tac <- load_tac(f, format="PMOD")
+#' AD06_tac_nCicc <- change_units(AD06_tac, to_unit = "nCi/cc")
+#'
+#' change_units(5, to_unit = "kBq", from_unit = "nCi")
+#' change_units(0.185, to_unit = "nCi", from_unit = "kBq")
+change_units <- function(x, to_unit, from_unit) UseMethod("change_units")
+
+#'@noRd
+#'@export
+change_units.tac <- function(x, to_unit, from_unit = NULL) {
+
+  if(!validate_tac(x)) stop("Supplied tac file did not validate.")
+  if(!is.null(from_unit)) stop("The from_unit is found within the tac object.")
+
+  unit_index <- get_activity_unit_index()
+
+  from <- sub("([A-Za-z]+).*", "\\1", attributes(x)$activity_unit)
+  to <- sub("([A-Za-z]+).*", "\\1", to_unit)
+
+  x[,-(1:2)] <- x[,-(1:2)] * ( unit_index[[from]] / unit_index[[to]] )
+
+  attributes(x)$activity_unit <- to_unit
+
+  if(!validate_tac(x)) stop("Converted tac file did not validate.") 
+
+  return(x)
+  
+}
+
+#'@noRd
+#'@export
+change_units.numeric <- function(x, to_unit, from_unit) {
+
+  unit_index <- get_activity_unit_index()
+
+  from <- sub("([A-Za-z]+).*", "\\1", from_unit)
+  to <- sub("([A-Za-z]+).*", "\\1", to_unit)
+
+  x <- x * ( unit_index[[from]] / unit_index[[to]] )
+
+  return(x)
+
+}
+