R/basis_set.R

In neuroconductor-devel/spant: MR Spectroscopy Analysis Tools

#' Generate a basis file using TARQUIN.
#' @param basis_file filename of the basis file to be generated.
#' @param metab_data MRS data object to match the generated basis parameters.
#' @param opts list of options to pass to TARQUIN.
#' @examples
#' \dontrun{
#' write_basis_tqn('test.basis',mrs_data,c("--echo","0.04"))
#' }
#' @export
write_basis_tqn <- function(basis_file, metab_data, opts = NULL) {
  # write the metab data to a temp dpt file
  metab_file <- tempfile()
  write_mrs_dpt_v2(metab_file, metab_data)
  
  # run TARQUIN
  cmd = paste(getOption("spant.tqn_cmd"), "--input",metab_file, "--format dpt",
              "--no_fit true", "--water_width 0", "--auto_ref false",
              "--auto_phase false", "--output_basis_lcm", basis_file,
              paste(opts,collapse = " "))
  
  res = system(cmd, intern = TRUE)
  
  if (!file.exists(basis_file)) {
    print(res)
    print(cmd)
    stop("Error loading data with above TARQUIN command.")
  }
}

#' Simulate a basis file using TARQUIN.
#' @param fs sampling frequency
#' @param ft transmitter frequency
#' @param N number of data points
#' @param ref chemical shift reference
#' @param opts list of options to pass to TARQUIN.
#' @examples
#' \dontrun{
#' write_basis_tqn('test.basis',mrs_data,c("--echo","0.04"))
#' }
#' @export
sim_basis_tqn <- function(fs = def_fs(), ft = def_ft(), N = def_N(),
                          ref = def_ref(), opts = NULL) {
  
  metab_data <- sim_zeros(fs = fs, ft = ft, N = N, ref = ref)
  
  # write the metab data to a temp dpt file
  metab_file <- tempfile()
  write_mrs_dpt_v2(metab_file, metab_data)
  
  # specify a basis file 
  basis_file <- tempfile()
  
  # run TARQUIN
  cmd = paste(getOption("spant.tqn_cmd"), "--input", metab_file, "--format dpt",
              "--no_fit true", "--water_width 0", "--auto_ref false",
              "--auto_phase false", "--output_basis_lcm", basis_file,
              paste(opts,collapse = " "))
  
  res = system(cmd, intern = TRUE)
  
  if (!file.exists(basis_file)) {
    print(res)
    print(cmd)
    stop("Error loading data with above TARQUIN command.")
  }
  
  read_basis(basis_file, ref = ref)
}

#' @export
print.basis_set <- function(x, ...) {
  cat("Basis set parameters\n")
  cat("-------------------------------\n")
  cat(paste(c("Trans. freq (MHz)       : ", x$ft * 1e-6, "\n")), sep = "")
  cat(paste(c("Data points             : ", dim(x$data)[1], "\n")), sep = "")
  cat(paste(c("Sampling frequency (Hz) : ", x$fs, "\n")), sep = "")
  cat(paste(c("Elements                : ", dim(x$data)[2], "\n\n")), sep = "")
  cat(paste(c("Names\n")), sep = "")
  cat("-------------------------------\n")
  cat(x$names, sep = ",", fill = 31)
}

#' @export
plot.basis_set <- function(x, n = 1, ...) {
  mrs_basis <- basis2mrs_data(x)
  graphics::plot(mrs_basis, dyn = n, ...)
  #matplot(Re(x$data),type="l",lty=1,col=1)
}

#' @export
stackplot.basis_set <- function(x, ...) {
  mrs_basis <- basis2mrs_data(x)
  stackplot(mrs_basis, ...)
}

#' Read a basis file in LCModel .basis format.
#' @param basis_file path to basis file.
#' @param ref assumed ppm reference value.
#' @return basis object.
#' @export
read_basis <- function(basis_file, ref = def_ref()) {
  con  <- file(basis_file, open = "r")
  names <- vector()
  data <- vector()
  
  while (length(line <- readLines(con, n = 1, warn = FALSE)) > 0) {
    if (startsWith(line, " NDATAB = ")) {
      N <- as.integer(strsplit(trimws(line), "\\s+")[[1]][3])
      data_lines <- ceiling(2 * N / 6) # assume 6 cols
    } else if (startsWith(line, " HZPPPM = ")) {
      bas_ft <- strsplit(trimws(line), "\\s+")[[1]][3]
      bas_ft <- as.double(gsub(",", "", bas_ft))*1e6
    } else if (startsWith(line, " BADELT = ")) {
      bas_fs <- strsplit(trimws(line), "\\s+")[[1]][3]
      bas_fs <- 1/as.double(gsub(",", "",bas_fs))
    } else if (endsWith(line, "$BASIS")) {
      while (length(line <- readLines(con, n = 1, warn = FALSE)) > 0) {
        if (startsWith(line, " ID = ")) {
          id <- (strsplit(trimws(line), " ")[[1]][3])
          id <- gsub(",", "",id)
          id <- gsub("'", "",id)
          names <- c(names, id)
        } else if (endsWith(line, "$END")) {
          x <- utils::read.fortran(con, "6F13.0", n = data_lines)
          data_pts <- as.vector(t(as.matrix(x)))
          data_pts <- data_pts[seq(1, 2 * N, 2)] +
                      1i * data_pts[seq(2, 2 * N, 2)]
          data_pts <- pracma::ifftshift(data_pts)
          data <- cbind(data, data_pts)
          break
        }
      }
    }
  }
  close(con)
  
  basis_set <- list(data = data, N = N, fs = bas_fs, ft = bas_ft, 
                    names = names, ref = ref)
  
  class(basis_set) <- "basis_set"
  sort_basis(basis_set)
}

#' Write a basis object to an LCModel .basis formatted file.
#' @param basis basis object to be exported.
#' @param basis_file path to basis file to be generated.
#' @param fwhmba parameter used by LCModel.
#' @export
write_basis <- function(basis, basis_file, fwhmba = 0.1) {
  mrs_data <- basis2mrs_data(basis)
  N <- Npts(mrs_data)
  
  sink(basis_file)
  cat(" $SEQPAR\n", sep = "")
  cat(" FWHMBA =  ", fwhmba, " ,\n", sep = "")
  cat(" HZPPPM =  ",mrs_data$ft * 1e-6,",\n", sep = "")
  cat(" ECHOT =  ",mrs_data$meta$EchoTime * 1000,",\n", sep = "")
  cat(" SEQ = 'PRESS' $END\n", sep = "")
  cat(" $BASIS1\n", sep = "")
  cat(" IDBASI = 'SPANT',\n", sep = "")
  cat(" FMTBAS = '(6E13.5)',\n", sep = "")
  cat(" BADELT =  ", 1 / fs(mrs_data), ",\n", sep = "")
  cat(" NDATAB = ", Npts(mrs_data), " $END\n",sep = "")
  
  for (n in 1:Ndyns(mrs_data)) {
    name = basis$names[n]
    cat(" $NMUSED\n")
    cat(" FILERAW = '", name,"',\n", sep = "")
    cat(" METABO_CONTAM = ' ',\n")
    cat(" METABO_SINGLET = ' ',\n")
    cat(" AUTOPH = F,\n")
    cat(" AUTOSC = F,\n")
    cat(" DO_CONTAM = F,\n")
    cat(" FLATEN = F,\n")
    cat(" CONCSC = -1.,\n")
    cat(" DEGZER =  0.,\n")
    cat(" DEGPAP =  0.,\n")
    cat(" DEGPPM =  0.,\n")
    cat(" HWDPHA =  0.100000001,\n")
    cat(" HWDSCA =  0.100000001,\n")
    cat(" PPMAPP =  0.5 -0.5,\n")
    cat(" PPMAPP_CONTAM =  0.  0.,\n")
    cat(" PPMBAS =  0.00999999978,\n")
    cat(" PPMPK =  0.00999999978,\n")
    cat(" PPMPK_CONTAM =  0.,\n")
    cat(" PPMPHA =  0.,\n")
    cat(" PPMSCA =  8.43999958,\n")
    cat(" PPM_SPLIT = -999.,\n")
    cat(" RINTEG =  0.,\n")
    cat(" SDPNTS =  1.,\n")
    cat(" XTRASH =  0. $END\n")
    cat(" $BASIS\n")
    cat(" ID = '", name, "',\n", sep = "")
    cat(" METABO = '", name, "',\n", sep = "")
    cat(" CONC = 1.0,\n")
    cat(" TRAMP = 1.0,\n")
    cat(" VOLUME = 1.0,\n")
    cat(" ISHIFT = 0 $END\n")
    sig <- pracma::fftshift(mrs_data$data[1, 1, 1, 1, n, 1,])
    for (n in seq(1, N, 3)) {
      cat(" ")
      cat(noquote(formatC(c(Re(sig[n]), Im(sig[n])), width = 12,format = "E",
                            digits = 5)))
      
      if ((n + 1) < N) {
        cat(" ")
        cat(noquote(formatC(c(Re(sig[n + 1]), Im(sig[n + 1])), width = 12, 
                              format = "E", digits = 5)))
      }
      if ((n + 2) < N) {
        cat(" ")
        cat(noquote(formatC(c(Re(sig[n + 2]), Im(sig[n + 2])), width = 12, 
                              format = "E", digits = 5)))
      }
      cat("\n")
    }
  }
  sink()
}
  
#' Convert a basis object to an mrs_data object - where basis signals are spread
#' across the dynamic dimension.
#' @param basis basis set object.
#' @param sum_elements return the sum of basis elements (logical)
#' @param amp a vector of scaling factors to apply to each basis element.
#' @param shift a vector of frequency shifts (in PPM) to apply to each basis
#' element.
#' @return an mrs_data object with basis signals spread across the dynamic 
#' dimension or summed.
#' @export
basis2mrs_data <- function(basis, sum_elements = FALSE, amp = NULL,
                           shift = NULL) {
  
  res <- mat2mrs_data(t(basis$data), fs = basis$fs, ft = basis$ft,
                      ref = basis$ref, fd = TRUE)
  
  # scale basis elements
  if (!is.null(amp)) {
    n_sigs <- Ndyns(res)
    if (n_sigs != length(amp)) {
      stop(paste("Error, length of amp does not match the number of basis elements :", dim(basis$data)[2]))
    }
    for (n in 1:n_sigs) {
      res$data[,,,,n ,,] <- res$data[,,,,n ,,] * amp[n]
    }
  }
  
  if (!is.null(shift)) res <- shift(res, shift)
  
  if (sum_elements) res <- sum_dyns(res)
  
  res
}

#' Convert an mrs_data object to basis object - where basis signals are spread
#' across the dynamic dimension in the MRS data.
#' @param mrs_data mrs_data object with basis signals spread across the dynamic dimension.
#' @param names list of names corresponding to basis signals.
#' @return basis set object.
#' @export
mrs_data2basis <- function(mrs_data, names) {
  # transform to FD
  if (!is_fd(mrs_data)) {
      mrs_data <- td2fd(mrs_data)
  }
  
  data <- matrix(mrs_data$data[1, 1, 1, 1,, 1,], 
            nrow = Npts(mrs_data), ncol = Ndyns(mrs_data), byrow = TRUE)
  
  basis_set <- list(data = data, N = Npts(mrs_data), fs = fs(mrs_data), 
                    ft = mrs_data$ft, names = names, ref = mrs_data$ref)
  
  class(basis_set) <- "basis_set"
  basis_set
}

sort_basis <- function(basis) {
  names_sorted <- sort(basis$names, index.return = TRUE)
  basis$names <- names_sorted$x
  basis$data <- basis$data[, names_sorted$ix]
  basis
}

#' Combine a pair of basis set objects.
#' @param basis_a first basis.
#' @param basis_b second basis.
#' @return combined basis set object.
#' @export
append_basis <- function(basis_a, basis_b) {
  basis_a_mrs <- basis2mrs_data(basis_a)
  basis_b_mrs <- basis2mrs_data(basis_b)
  basis_out_mrs <- append_dyns(basis_a_mrs, basis_b_mrs)
  basis_out <- mrs_data2basis(basis_out_mrs, c(basis_a$names, basis_b$names))
  return(basis_out)
}