Analysis with tacmagic

Data loading

Time-activity curve (TAC) data is loaded via load_tac(), which is a wrapper for format-specific functions. To specify which file format the TAC data is stored as, use the format parameter. Supported formats can be viewed in help(load_tac).

The minimal amount of information required is the TAC data for one or more ROI, including the start and stop times of each frame, the time units and the activity units. This information may be in 1 or more files depending on the format and software that created it.

For example, PMOD’s .tac files contain all of the information, but the TAC .voistat files do not contain start and stop times, but this information could be specified using a .acqtimes file. Support is also available for DFT format, which contains both TAC and volume data.

We processed the PIB PET and T1 MRI data with the PMOD PNEURO software suite to produce a .tac file with TACs for all ROIs in the Hammer’s atlas. The .tac file can be loaded with load_tac():

# Filename is a character string of the file's path on your computer.
filename <- system.file("extdata", "AD06.tac", package="tacmagic")
# Note: This file can also serve as a template if the TAC data is in some other 
# format that is not yet supported.

AD06_tac <- load_tac(filename, format="PMOD")

A TAC object is a data frame with extra attributes including time and activity units. A summary can be printed with the generic print() function.

summary(AD06_tac) 
#> tac object
#>  Activity unit:           kBq/cc 
#>  Time unit:               seconds 
#>  Number of ROIs:          172 
#>  Number of frames:        34 
#>  Time span:               0 - 5400 seconds 
#>  Unique frame durations:  15 30 60 180 300 600 seconds

AD06_tac[1:5,1:5] # the first 5 frames of the first 3 ROIs
#>   start end FL_mid_fr_G_l FL_mid_fr_G_r FL_precen_G_l
#> 1     0  15    0.05156934     0.0830544   -0.05035268
#> 2    15  30    5.76280770     5.9782509    6.15967912
#> 3    30  45   16.07137663    16.3756559   16.69802808
#> 4    45  60   19.11766922    19.4707593   18.52006408
#> 5    60  90   19.77191392    20.3383362   19.55189423

PMOD’s suite also produces .voistat and .acqtimes formats, than can be used to produce the same data if you do not have .tac files:

filename_acq <- system.file("extdata", "AD06.acqtimes", package="tacmagic")
filename_voistat <- system.file("extdata", "AD06_TAC.voistat", package="tacmagic")

tac2 <- load_tac(filename_voistat, format="voistat", acqtimes=filename_acq)

all.equal(AD06_tac, tac2)
#> [1] TRUE

We also used Turku’s magia pipeline to process the same data. It can be loaded similarly, though with units explicitly entered because the information is not encoded in the .mat file:

f_magia <- system.file("extdata", "AD06_tac_magia.mat", package="tacmagic")

AD06_tac_magia <- load_tac(f_magia, format="magia", 
                           time_unit="seconds", activity_unit="kBq/cc")

AD06_tac_magia[1:5,1:5]
#>   start end        CGA         CGP          LOC
#> 1     0  15  0.2070033  0.04986614 -0.001164075
#> 2    15  30  7.8729520  7.12711797  6.058719929
#> 3    30  45 16.6898445 17.37585231 18.097222061
#> 4    45  60 17.7306580 19.33122997 21.146157237
#> 5    60  90 18.9857118 19.57800502 22.009607171

manual <- data.frame(start=c(0:4), end=c(2:6), ROI1=c(10.1:14.2), ROI2=c(11:15))
manual_tac <- as.tac(manual, time_unit="minutes", activity_unit="kBq/cc")

summary(manual_tac)
#> tac object
#>  Activity unit:           kBq/cc 
#>  Time unit:               minutes 
#>  Number of ROIs:          2 
#>  Number of frames:        5 
#>  Time span:               0 - 6 minutes 
#>  Unique frame durations:  2 minutes

Radioactivity unit conversion

Most modern PET tools use kBq/cc as the standard activity units required by the software (e.g. TPCCLIB, PMOD). Most often data will be in this format and not require conversion. However, conversion of TAC objects to these units, or to other radioactivity units for that matter, is possible with change_units(). This is a generic function that works on tac objects as well as numeric objects. For numeric objects, both “to” and “from” units need to be specified. The function works regardless of whether the units are per volume (i.e. kBq is treated the same was as kBq/cc or kBq/mL).

change_units(5, to_unit = "kBq", from_unit = "nCi")
#> [1] 0.185
change_units(0.5, to_unit = "nCi/cc", from_unit = "kBq/cc")
#> [1] 13.51351

With tac objects, as the activity units are stored in the object, they should not be provided to change_units():

AD06_nCi <- change_units(AD06_tac, to_unit = "nCi/cc")
summary(AD06_nCi)
#> tac object
#>  Activity unit:           nCi/cc 
#>  Time unit:               seconds 
#>  Number of ROIs:          172 
#>  Number of frames:        34 
#>  Time span:               0 - 5400 seconds 
#>  Unique frame durations:  15 30 60 180 300 600 seconds

ROI merging

Often it is desirable to combine TAC ROIs into larger ROIs. For example, if the PET analysis pipeline created TACs for each atlas ROI, your analysis may call for merging these atomic ROIs into larger regions, such as merging all of the atlas ROIs that make up the frontal lobe into a single frontal lobe ROI.

If this is done, the means should be weighted for the relative volumes of the atomic ROIs. If volume information is available, tac_roi() provides this functionality.

In PMOD’s software, volume information is available in .voistat files. Units do not matter because it is the relative volume information that is needed.

In addition to TAC and volume information, we must specify which atomic ROIs make up the merged ROI. This is done by providing a named list, where the names are the merged ROIs and the list items are themselves lists of the atomic ROIs that make up each merged ROI. For the Hammer’s atlas, and as an example, typical data is provided in roi_ham_stand(), roi_ham_full(), or roi_ham_pib().

AD06_volume <- load_vol(filename_voistat, format="voistat")

roi_ham_pib()[1:2] # The first 2 definitions of merged ROIs, as an example.
#> $leftfrontal
#>  [1] "FL_mid_fr_G_l"       "FL_precen_G_l"       "FL_strai_G_l"       
#>  [4] "FL_OFC_AOG_l"        "FL_inf_fr_G_l"       "FL_sup_fr_G_l"      
#>  [7] "FL_OFC_MOG_l"        "FL_OFC_LOG_l"        "FL_OFC_POG_l"       
#> [10] "Subgen_antCing_l"    "Subcall_area_l"      "Presubgen_antCing_l"
#> 
#> $rightfrontal
#>  [1] "FL_mid_fr_G_r"       "FL_precen_G_r"       "FL_strai_G_r"       
#>  [4] "FL_OFC_AOG_r"        "FL_inf_fr_G_r"       "FL_sup_fr_G_r"      
#>  [7] "FL_OFC_MOG_r"        "FL_OFC_LOG_r"        "FL_OFC_POG_r"       
#> [10] "Subgen_antCing_r"    "Subcall_area_r"      "Presubgen_antCing_r"

AD06 <- tac_roi(tac=AD06_tac,           # The TAC file we loaded above.
                volumes=AD06_volume,    # Volume information loaded.
                ROI_def=roi_ham_pib(),  # ROI definitions for the Hammers atlas
                merge=F,                # T to also return atomic ROIs
                PVC=F                   # to use _C ROIs (PMOD convention)            
                )

AD06[1:5,1:5]
#>   start end leftfrontal rightfrontal lefttemporal
#> 1     0  15  0.02912641   0.05734614    0.0612965
#> 2    15  30  6.12655470   6.23215512    6.3573101
#> 3    30  45 16.06737144  16.00035525   15.4122525
#> 4    45  60 18.31769407  18.45778499   17.3612789
#> 5    60  90 19.24138550  19.18263937   17.9358282

Plotting

Basic TAC plotting can be done by calling plot, which accepts two TAC objects, e.g. from 2 participants or group means. The ROIs to plot are specified as a vector of ROI names as they appear in the TAC object. As the TAC object contains time unit information, the plot can convert to desired units, which can be specified with the time argument.

plot(AD06,                                                    # TAC data
     ROIs=c("frontal", "temporal", "parietal", "cerebellum"), # ROIs to plot
     time="minutes",                   # Convert x axis from seconds to minutes
     title="PIB time activity curves for AD06"        # A title for the plot
     )

Standardized uptake value (SUV)

As the activity in the TAC is impacted by the dose of the radiotracer administered and the participant’s body weight, a value adjusted for these factors is sometimes used, the SUV:

$$SUV = \frac{Ct}{\frac{Dose}{Weight}}$$

Where activity is measured in kBq/mL (kBq/cc), the dose is in MBq, and the weight is in kg, the radioactivity units cancel and the SUV units are g/mL.

With the tac_suv() function, a tac object can be converted to SUV values with units g/mL, as demonstrated below (note the weight and dose are fabricated here for the demonstration).

AD06_suv_tac <- tac_suv(AD06, dose = 8.5, dose_unit = "mCi", weight_kg = 70)

More often, an everage value over a certain time period, or a maximum value may be desired. This can be calculated with the suv() function.

AD06_suv_calc <- suv(AD06, SUV_def = c(3000, 3300, 3600), dose = 8.5, dose_unit = "mCi", weight_kg = 70)
AD06_suv_calc["frontal",]
#> [1] 1.713055
AD06_suv_max <- suv(AD06, SUV_def = "max", dose = 8.5, dose_unit = "mCi", weight_kg = 70)
AD06_suv_max["frontal",]
#> [1] 4.338809

SUV ratio (SUVR)

The standardized uptake value ratio (SUVR) is a simple quantification of PET activity that is commonly used from many tracers including PIB. It is the ratio of the tracer activity over a specified time period (Ct) in a target ROI to a reference region. Using a ratio allows factors that are normally required to calculate an SUV to cancel out, namely tracer dose and patient body weight, and therefore SUVR can be calculated from TAC data alone, i.e. without the need to specify tracer dose or body weight:

$$SUVR = \frac{SUV_{TARGET}}{SUV_{REF}} = \frac{Ct_{TARGET}}{Ct_{REF}}$$

In the literature, SUVR is variably described and calculated using the mean of activity for the frames of the specified time period, or the area under the curve. For PIB, the mean/summed activity has been used, and the time windows have varied from starting at 40-50 minutes and ending at 60-90 minutes.³

The suvr() function calculates SUVR for all regions in a TAC file based on the provided time information (as a vector of frame start times) and the specified reference region (a string). If the frames used are of different durations, the weighted mean is used.

AD06_SUVR <- suvr(AD06,                       # TAC data
                  SUVR_def=c(3000,3300,3600), # = 50-70 minute window
                  ref="cerebellum"            # reference region in TAC data
                  )

AD06_SUVR
#>                    SUVR
#> leftfrontal    2.070909
#> rightfrontal   2.200107
#> lefttemporal   2.043683
#> righttemporal  2.176933
#> leftparietal   2.054232
#> rightparietal  2.169382
#> leftoccipital  1.659831
#> rightoccipital 1.827765
#> leftcingulate  2.348674
#> rightcingulate 2.396375
#> frontal        2.136009
#> temporal       2.110140
#> parietal       2.112264
#> occipital      1.744258
#> cingulate      2.371747
#> cerebellum     1.000000
#> totalcortical  2.073093
#> leftdeep       1.788912
#> rightdeep      1.755389
#> deep           1.772272
#> ventricles     0.902195
#> whitematter    2.035768
#> amyloidcomp    2.231989

An alternative method, using the area under the curve with the mid-frame times as the x-axis is available with suvr_auc() and should provide very similar results.

AD06_altSUVR <- suvr_auc(AD06, SUVR_def=c(3000,3300,3600), ref="cerebellum")

all.equal(AD06_SUVR, AD06_altSUVR) # Should be similar but not exact
#> [1] "Component \"SUVR\": Mean relative difference: 0.007765957"

DVR

The Distribution Volume Ratio (DVR) is a method of quantifying tracer uptake that is used as an alternative to the SUVR in PIB studies, for example. Like SUVR, it can be calculated from TAC data without the need for arterial blood sampling, by making use of a reference region. In this case, it is called the non-invasive Logan plot method. It is calculated with a graphical analysis technique described by Logan et al.⁴

In addition to the TAC data, depending on the tracer, a value for k2’ may need to be specified. For PIB, this has limited effect on the result, but can be specified, and a value of 0.2 has been recommended.⁵

The non-invasive Logan plot works by finding the slope of the line of the following equation after time t* where linearity has been reached:

$$\frac{\int_0^{T}C_{roi}(t)dt}{C_{roi}(t)} = DVR[\frac{\int_0^{T}C_{cer}(t)dt + C_{cer}(t) / k2`}{C_{roi}(T)}] + int $$

AD06_DVR_fr <- DVR_ref_Logan(AD06, 
                             target="frontal", # target ROI
                             ref="cerebellum", # reference region
                             k2prime=0.2,      # suitable k2' for tracer
                             t_star=0,        # 0 to find, or can specify frame
                             )
#> t* is 23

AD06_DVR_fr$DVR
#> [1] 1.780057

plot(AD06_DVR_fr)

AD06_DVR <- DVR_all_ref_Logan(AD06, ref="cerebellum", k2prime=0.2, t_star=23)
#> Trying leftfrontal
#> Trying rightfrontal
#> Trying lefttemporal
#> Trying righttemporal
#> Trying leftparietal
#> Trying rightparietal
#> Trying leftoccipital
#> Trying rightoccipital
#> Trying leftcingulate
#> Trying rightcingulate
#> Trying frontal
#> Trying temporal
#> Trying parietal
#> Trying occipital
#> Trying cingulate
#> Trying cerebellum
#> Trying totalcortical
#> Trying leftdeep
#> Trying rightdeep
#> Trying deep
#> Trying ventricles
#> Trying whitematter
#> Trying amyloidcomp

AD06_DVR
#>                      DVR
#> leftfrontal    1.7351136
#> rightfrontal   1.8241000
#> lefttemporal   1.7023599
#> righttemporal  1.7978415
#> leftparietal   1.7148840
#> rightparietal  1.8014491
#> leftoccipital  1.4303466
#> rightoccipital 1.5521520
#> leftcingulate  1.9240049
#> rightcingulate 1.9542491
#> frontal        1.7800568
#> temporal       1.7500941
#> parietal       1.7585448
#> occipital      1.4915376
#> cingulate      1.9389554
#> cerebellum     1.0000000
#> totalcortical  1.7294061
#> leftdeep       1.5304215
#> rightdeep      1.5307336
#> deep           1.5305573
#> ventricles     0.7435322
#> whitematter    1.6615718
#> amyloidcomp    1.8487819

ADO6_frontal_DVR <- dvr(AD06, target="frontal", ref="cerebellum", k2prime=0.2, 
                        t_star=23)

Batch loading

Data loading is done by batch_load(). See help(batch_load) for the required arguments.

The first argument is a vector of participant IDs that corresponds to file names, e.g.:

participants <- c("participant01", "participant02") if the files are located e.g. /mypath/participant01.tac and /mypath/participant01_TAC.voistat. In this case, the function call might look like:

my_data <- batch_load(participants, dir="/mypath/", tac_format="PMOD", roi_m=T, vol_file_suffix="_TAC.voistat", vol_format="voistat", ROI_def=roi_ham_stand(), merge=F)

The above would load the appropriate TAC and voistat files, perform the ROI merging specified by ROI_def, because roi_m = TRUE, and would return a list where each element represents a participants, e.g. the first participant would be my_data$participant1.

To calculate SUV in batch, the participants dose and weight must be specified when loading with batch_load(), as it is then added to the respective tac objects.

Batch analysis

Once the TAC data is loaded, all analyses can be run using batch_tm(). The output from batch_load() is the first argument for batch_tm(). The models implemented in tacmagic can be specified using the models argument, e.g. models = c("SUVR", "Logan") to calculate both SUVR and Logan DVR. The relevant model parameters will also need to be specified, so see help(batch_tm) for all possible arguments.

Batch example

For the purpose of the vignette, the list of participants will be a list of the full TAC filenames (hence tac_file_suffix=““). In real-world data, the participants parameter can be a list of participant IDs that correspond to the actual filenames, i.e. the filename is made up of dir + participant + tac_file_suffix.

We will also choose not to use the roi_m option in batch_load(), which could be used to combine ROIs as outlined above.

participants <- c(system.file("extdata", "AD06.tac", package="tacmagic"),
                   system.file("extdata", "AD07.tac", package="tacmagic"),
                   system.file("extdata", "AD08.tac", package="tacmagic"))

tacs <- batch_load(participants, dir="", tac_file_suffix="")

# Since the PMOD TAC files used here have 2 copies of ROIs, with and without 
# PVC, we can use split_pvc to keep the PVC-corrected verions. If we had used 
# roi_m here to combine ROIs, we could have specified to use the PVC versions 
# in batch_load() with PVC = TRUE.
tacs <- lapply(tacs, split_pvc, PVC=TRUE)
 
batch <- batch_tm(tacs, models=c("SUVR", "Logan"), ref="Cerebellum_r_C",
                  SUVR_def=c(3000,3300,3600), k2prime=0.2, t_star=23)