Reading raw data and plotting
Load the spant package:
Get the path to a data file included with spant:
fname <- system.file("extdata", "philips_spar_sdat_WS.SDAT", package = "spant")
Read the file and save to the workspace as mrs_data:
mrs_data <- read_mrs(fname)
Output some basic information about the data:
print(mrs_data)
#> MRS Data Parameters
#> ----------------------------------
#> Trans. freq (MHz) : 127.7861
#> FID data points : 1024
#> X,Y,Z dimensions : 1x1x1
#> Dynamics : 1
#> Coils : 1
#> Voxel resolution (mm) : 20x20x20
#> Sampling frequency (Hz) : 2000
#> Repetition time (s) : 2
#> Reference freq. (ppm) : 4.65
#> Nucleus : 1H
#> Spectral domain : FALSE
#> Number of transients : 128
#> Echo time (s) : 0.03
#> Manufacturer : Philips
Plot the spectral region between 5 and 0.5 ppm:
plot(mrs_data, xlim = c(5, 0.5))
Basic preprocessing
Apply a HSVD filter to the residual water region and align the
spectrum to the tNAA resonance at 2.01 ppm:
mrs_proc <- hsvd_filt(mrs_data)
mrs_proc <- align(mrs_proc, 2.01)
plot(mrs_proc, xlim = c(5, 0.5))
Basis simulation
Simulate a typical basis set for short TE brain analysis, print some
basic information and plot:
basis <- sim_basis_1h_brain_press(mrs_proc)
print(basis)
#> Basis set parameters
#> -------------------------------
#> Trans. freq (MHz) : 127.8
#> Data points : 1024
#> Sampling frequency (Hz) : 2000
#> Elements : 27
#>
#> Names
#> -------------------------------
#> -CrCH2,Ala,Asp,Cr,GABA,Glc,Gln,
#> GSH,Glu,GPC,Ins,Lac,Lip09,
#> Lip13a,Lip13b,Lip20,MM09,MM12,
#> MM14,MM17,MM20,NAA,NAAG,PCh,
#> PCr,sIns,Tau
stackplot(basis, xlim = c(4, 0.5), labels = basis$names, y_offset = 5)
Perform ABfit analysis of the processed data
(mrs_proc):
fit_res <- fit_mrs(mrs_proc, basis)
Plot the fit result:
Unscaled amplitudes, CRLB error estimates and other useful fitting
diagnostics, such as SNR, are given in the fit_res results
table:
fit_res$res_tab
#> X Y Z Dynamic Coil X.CrCH2 Ala Asp Cr
#> 1 1 1 1 1 1 0 8.133833e-06 3.547057e-05 4.026724e-05
#> GABA Glc Gln GSH Glu GPC
#> 1 1.697279e-05 2.446714e-06 3.03616e-06 2.22805e-05 6.503224e-05 1.606707e-05
#> Ins Lac Lip09 Lip13a Lip13b Lip20 MM09
#> 1 5.906169e-05 5.802226e-06 2.387703e-05 2.67061e-06 0 0 9.630787e-06
#> MM12 MM14 MM17 MM20 NAA NAAG
#> 1 6.511375e-06 2.603851e-05 2.238405e-05 9.203747e-05 6.011109e-05 1.536527e-05
#> PCh PCr sIns Tau tNAA tCr tCho
#> 1 0 2.101916e-05 6.504113e-06 0 7.547636e-05 6.128639e-05 1.606707e-05
#> Glx tLM09 tLM13 tLM20 X.CrCH2.sd Ala.sd
#> 1 6.80684e-05 3.350782e-05 3.522049e-05 9.203747e-05 2.386862e-06 4.343936e-06
#> Asp.sd Cr.sd GABA.sd Glc.sd Gln.sd GSH.sd
#> 1 9.235332e-06 3.688997e-06 4.574412e-06 4.421685e-06 5.083004e-06 2.020428e-06
#> Glu.sd GPC.sd Ins.sd Lac.sd Lip09.sd Lip13a.sd
#> 1 5.082979e-06 2.602651e-06 2.093228e-06 5.301707e-06 4.119258e-06 1.328304e-05
#> Lip13b.sd Lip20.sd MM09.sd MM12.sd MM14.sd MM17.sd
#> 1 6.47447e-06 7.510912e-06 3.827555e-06 4.59403e-06 7.223428e-06 3.811303e-06
#> MM20.sd NAA.sd NAAG.sd PCh.sd PCr.sd sIns.sd
#> 1 8.595857e-06 1.017712e-06 1.208983e-06 2.236872e-06 3.084592e-06 7.240073e-07
#> Tau.sd tNAA.sd tCr.sd tCho.sd Glx.sd tLM09.sd
#> 1 3.75996e-06 7.031384e-07 5.890483e-07 2.110663e-07 3.169662e-06 1.006794e-06
#> tLM13.sd tLM20.sd phase lw shift asym
#> 1 1.573147e-06 3.016811e-06 11.10961 5.023681 -0.003765049 0.1771065
#> res.deviance res.niter res.info
#> 1 7.307679e-05 28 2
#> res.message bl_ed_pppm
#> 1 Relative error between `par' and the solution is at most `ptol'. 2.364083
#> max_bl_flex_used full_res spec_resid fit_pts ppm_range SNR
#> 1 FALSE 7.754374e-05 7.307679e-05 497 3.8 62.71686
#> SRR FQN tNAA_lw tCr_lw tCho_lw auto_bl_crit_7
#> 1 51.33276 1.492723 0.04562445 0.05189506 0.05438239 -8.900349
#> auto_bl_crit_5.901 auto_bl_crit_4.942 auto_bl_crit_4.12 auto_bl_crit_3.425
#> 1 -8.944159 -8.977354 -9.000064 -9.013366
#> auto_bl_crit_2.844 auto_bl_crit_2.364 auto_bl_crit_1.969 auto_bl_crit_1.647
#> 1 -9.02055 -9.024177 -9.023461 -9.009854
#> auto_bl_crit_1.384 auto_bl_crit_1.17 auto_bl_crit_0.997 auto_bl_crit_0.856
#> 1 -8.958655 -8.844273 -8.690748 -8.562637
#> auto_bl_crit_0.743 auto_bl_crit_0.654 auto_bl_crit_0.593 auto_bl_crit_0.558
#> 1 -8.485478 -8.447144 -8.430182 -8.423061
#> auto_bl_crit_0.54 auto_bl_crit_0.532 auto_bl_crit_0.529
#> 1 -8.420094 -8.418848 -8.418321
Note that signal names appended with “.sd” are the CRLB estimates for
the uncertainty (standard deviation) in the metabolite quantity
estimate. e.g. to calculate the percentage s.d. for tNAA:
fit_res$res_tab$tNAA.sd / fit_res$res_tab$tNAA * 100
#> [1] 0.9316008
Spectral SNR:
fit_res$res_tab$SNR
#> [1] 62.71686
Linewidth of the tNAA resonance in PPM:
fit_res$res_tab$tNAA_lw
#> [1] 0.04562445
Ratios to total-creatine
Amplitude estimates measured by the fitting method are essentially
arbitrary unless scaled to a known reference signal. The simplest
approach for proton-MRS is to simply divide all metabolite values by
total-creatine:
fit_res_tcr_sc <- scale_amp_ratio(fit_res, "tCr")
amps <- fit_amps(fit_res_tcr_sc)
print(t(amps))
#> [,1]
#> X.CrCH2 0.00000000
#> Ala 0.13271841
#> Asp 0.57876749
#> Cr 0.65703389
#> GABA 0.27694224
#> Glc 0.03992264
#> Gln 0.04954051
#> GSH 0.36354731
#> Glu 1.06112038
#> GPC 0.26216377
#> Ins 0.96369978
#> Lac 0.09467397
#> Lip09 0.38959761
#> Lip13a 0.04357590
#> Lip13b 0.00000000
#> Lip20 0.00000000
#> MM09 0.15714396
#> MM12 0.10624503
#> MM14 0.42486603
#> MM17 0.36523691
#> MM20 1.50176027
#> NAA 0.98082276
#> NAAG 0.25071261
#> PCh 0.00000000
#> PCr 0.34296611
#> sIns 0.10612654
#> Tau 0.00000000
#> tNAA 1.23153536
#> tCr 1.00000000
#> tCho 0.26216377
#> Glx 1.11066089
#> tLM09 0.54674157
#> tLM13 0.57468696
#> tLM20 1.50176027
Water reference scaling, AKA “absolute-quantification”
A more sophisticated approach to scaling metabolite values involves
the use of a separate water-reference acquisition - which can be
imported in the standard way:
fname_wref <- system.file("extdata", "philips_spar_sdat_W.SDAT", package = "spant")
mrs_data_wref <- read_mrs(fname_wref)
The following code assumes the voxel contains 100% white matter
tissue and scales the metabolite values into molal (mM) units (mol / kg
tissue water) based on the method described by Gasparovic et al MRM 2006
55(6):1219-26:
p_vols <- c(WM = 100, GM = 0, CSF = 0)
TE = 0.03
TR = 2
fit_res_molal <- scale_amp_molal_pvc(fit_res, mrs_data_wref, p_vols, TE, TR)
fit_res_molal$res_tab$tNAA
#> [1] 12.59788
An alternative method scales the metabolite values into molar (mM)
units (mol / Litre of tissue) based on assumptions outlined in the
LCModel manual and references therein (section 10.2). This approach may
be preferred when comparing results to those obtained LCModel or
TARQUIN.
fit_res_molar <- scale_amp_molar(fit_res, mrs_data_wref)
fit_res_molar$res_tab$tNAA
#> [1] 6.826332
Note, while “absolute” units are attractive, a large number of
assumptions about metabolite and water relaxation rates are necessary to
arrive at these mM estimates. If you’re not confident at being able to
justify these assumptions, scaling to a metabolite reference (eg tCr as
above) is going to be a better option in most cases. Simple metabolite
referenced ratios also have the benefit of being more reproducible due
to the simplicity of the approach.