R/filterfMRIforNetworkAnalysis.R

In stnava/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Basic pre-processing for BOLD or ASL-based network analysis.
#'
#' This function works for either raw BOLD time-series data, ASL-based BOLD
#' time-series data or ASL-based CBF time series data.  In all 3 cases, this
#' function performs motion-correction, factoring out motion and compcor
#' nuisance paramters, frequency filtering and masking.  The output contains
#' the filtered time series (matrix form), the mask and a vector of temporal
#' signal variance. Some ASL MR sequences allow network analysis of either BOLD
#' or ASL signal.  See "Implementation of Quantitative Perfusion Imaging
#' Techniques for Functional Brain Mapping using Pulsed Arterial Spin Labeling"
#' by Wong et al, 1997 for an overview.  This function employs "surround"
#' techniques for deriving either CBF or BOLD signal from the input ASL.  This
#' is a WIP.
#'
#'
#' @param aslmat The filename to an antsr image or pointer
#' to an antsr image
#' @param tr The sequence's TR value , typically 3 or 4.
#' @param freqLo The lower frequency limit, e.g. 0.01 in band-pass
#' filter
#' @param freqHi The higher frequency limit, e.g. 0.1 in band-pass
#' filter
#' @param cbfnetwork "ASLCBF" A string dictating whether to do nothing special
#' (standard BOLD) or get CBF (ASLCBF) or BOLD (ASLBOLD) signal from ASL
#' @param mask the mask image
#' @param labels the label image
#' @param graphdensity desired density
#' @param seg a segmentation image
#' @param useglasso use sparse inverse covariance for network estimation
#' @param nuisancein nuisance variable data frame
#' @param usesvd bool, to reduce nuisance variables
#' @param robustcorr bool

#' @return output is a list containing "filteredTimeSeries" "mask"
#' "temporalvar"
#'
#' or
#'
#' 1 -- Failure
#' @author Avants BB
#' @examples
#' \dontrun{
#' if (!exists("fn")) fn <- getANTsRData("pcasl")
#' img <- antsImageRead(fn)
#' mask <- getMask(getAverageOfTimeSeries(img))
#' fmat <- timeseries2matrix(img, mask)
#' myres <- filterfMRIforNetworkAnalysis(fmat, tr = 4, 0.01, 0.1, cbfnetwork = "BOLD", mask = mask)
#' }
#'
#' @export filterfMRIforNetworkAnalysis
filterfMRIforNetworkAnalysis <- function(
    aslmat,
    tr, freqLo = 0.01, freqHi = 0.1,
    cbfnetwork = "ASLCBF",
    mask = NULL,
    labels = NULL,
    graphdensity = 0.5,
    seg = NULL,
    useglasso = NA,
    nuisancein = NA,
    usesvd = FALSE,
    robustcorr = FALSE) {
  pixtype <- "float"
  myusage <- "usage: filterfMRIforNetworkAnalysis( timeSeriesMatrix, tr, freqLo=0.01, freqHi = 0.1, cbfnetwork=c(\"BOLD,ASLCBF,ASLBOLD\") , mask = NULL,  graphdensity = 0.5 )"
  if (nargs() == 0) {
    print(myusage)
    return(NULL)
  }
  if (!is.null(mask)) {
    mask <- check_ants(mask)
  }
  if (!is.numeric(tr) | missing(tr)) {
    print("TR parameter is missing or not numeric type - is typically between 2 and 4 , depending on your fMRI acquisition")
    print(myusage)
    return(NULL)
  }
  if (!is.numeric(freqLo) | !is.numeric(freqHi)) {
    print("freqLo/Hi is not numeric type")
    print(myusage)
    return(NULL)
  }
  if (missing(aslmat)) {
    print("Missing first (image) parameter")
    print(myusage)
    return(NULL)
  }
  freqLo <- freqLo * tr
  freqHi <- freqHi * tr
  # network analysis
  # wb <- (mask > 0)  # whole brain
  if (!usePkg("magic")) {
    print("Need magic package")
    return(NULL)
  }
  leftinds <- magic::shift(c(1:nrow(aslmat)), 1)
  rightinds <- magic::shift(c(1:nrow(aslmat)), -1)
  oaslmat <- aslmat
  if (cbfnetwork == "ASLCBF") {
    # surround subtraction for cbf networks
    aslmat <- oaslmat - 0.5 * (oaslmat[leftinds, ] + oaslmat[rightinds, ])
    taginds <- c(1:(nrow(aslmat) / 2)) * 2
    controlinds <- taginds - 1
    aslmat[controlinds, ] <- aslmat[controlinds, ] * (-1) # ok! done w/asl specific stuff
    # plot( apply( aslmat[controlinds,],1, mean) , type='l') # should be (+) everywhere
  }
  if (cbfnetwork == "ASLBOLD") {
    # surround addition for bold networks
    aslmat <- oaslmat + 0.5 * (oaslmat[leftinds, ] + oaslmat[rightinds, ])
    # plot( apply( aslmat[controlinds,],1, mean) , type='l')
  }
  voxLo <- round((1 / freqLo)) # remove anything below this (high-pass)
  voxHi <- round((1 / freqHi)) # keep anything above this
  myTimeSeries <- ts(aslmat, frequency = 1 / tr)
  filteredTimeSeries <- frequencyFilterfMRI(aslmat, tr = tr, freqLo = freqLo, freqHi = freqHi, opt = "trig")
  vox <- round(ncol(filteredTimeSeries) * 0.5) #  a test voxel
  spec.pgram(filteredTimeSeries[, vox], taper = 0, fast = FALSE, detrend = F, demean = F, log = "n")
  temporalvar <- apply(filteredTimeSeries, 2, var)
  wh <- which(temporalvar == 0)
  for (x in wh) {
    filteredTimeSeries[, x] <- sample(filteredTimeSeries, nrow(filteredTimeSeries))
  }
  if (!is.null(labels)) {
    stopifnot(!is.null(mask))
    # do some network thing here
    oulabels <- sort(unique(labels[labels > 0]))
    whvec <- (mask == 1)
    ulabels <- sort(unique(labels[whvec]))
    if (ulabels[1] == 0) {
      ulabels <- ulabels[2:length(ulabels)]
    }
    labelvec <- labels[whvec]
    if (!is.null(seg)) {
      segvec <- seg[whvec]
    } else {
      segvec <- NA
    }
    labmat <- matrix(data = rep(NA, length(oulabels) * nrow(filteredTimeSeries)), nrow = length(oulabels))
    nrowts <- nrow(filteredTimeSeries)
    for (mylab in ulabels) {
      if (!is.null(seg)) {
        dd <- (labelvec == mylab & segvec == 2)
      } else {
        dd <- labelvec == mylab
      }

      submat <- filteredTimeSeries[, dd]
      # if ( length( c( submat ) ) > nrowts ) myavg<-svd( submat )$u[,1] else myavg<-submat

      if (length(c(submat)) > nrowts) {
        myavg <- apply(submat, MARGIN = 1, FUN = mean)
        if (usesvd) {
          eanat <- svd(submat, nu = 1, nv = 0)
          myavg <- eanat$u[, 1]
        }
      } else {
        myavg <- submat
      }

      if (length(myavg) > 0) {
        labmat[which(ulabels == mylab), ] <- myavg
      } else {
        labmat[which(ulabels == mylab), ] <- NA
      }
      # if ( length(myavg) > 0 ) labmat[ mylab, ]<-myavg else labmat[ mylab, ]<-NA
    }
    nbrainregions <- length(oulabels)
    tlabmat <- t(labmat)
    ocormat <- cor(tlabmat, tlabmat)
    rcormat <- ocormat
    if (robustcorr) {
      for (i in 1:nrow(rcormat)) {
        for (j in i:nrow(rcormat)) {
          if (i != j) {
            a <- scale(tlabmat[, i])
            b <- scale(tlabmat[, j])
            rcormat[i, j] <- rcormat[j, i] <- as.numeric(
              coefficients(MASS::rlm(a ~ b))
            )[2]
          }
        }
      }
    }
    if (!is.na(nuisancein)) {
      tlabmat <- cbind(tlabmat, nuisancein)
      ocormat <- cor(tlabmat, tlabmat)
    }
    cormat <- ocormat
    pcormat <- solve(cormat + (diag(ncol(cormat)) + 0.01))
    diagmag <- sqrt(diag(pcormat) %o% diag(pcormat)) * (-1)
    pcormat <- pcormat / diagmag
    diag(pcormat) <- 1
    gcormat <- pcormat
    if (!is.na(useglasso)) {
      if (useglasso > 0) # treat useglasso as rho
        {
          gcormat <- glasso::glasso(cormat, useglasso)$wi
          diagmag <- sqrt(diag(gcormat) %o% diag(gcormat)) * (-1)
          gcormat <- gcormat / diagmag
          myinds <- (abs(gcormat) < 1.e-4)
          gcormat[myinds] <- 0
          diag(gcormat) <- 1
          rgdens <- (0.5 * sum(gcormat > 0 & gcormat < 1) / (0.5 * ncol(gcormat) * (ncol(gcormat) - 1)))
        }
    }
    if (!is.na(nuisancein)) pcormat <- pcormat[1:nbrainregions, 1:nbrainregions]
    if (!is.na(nuisancein)) gcormat <- pcormat[1:nbrainregions, 1:nbrainregions]
    if (!is.na(nuisancein)) cormat <- cormat[1:nbrainregions, 1:nbrainregions]
    if (!is.na(nuisancein)) ocormat <- ocormat[1:nbrainregions, 1:nbrainregions]
    gmet <- makeGraph(cormat, graphdensity = graphdensity)
    L <- list(
      filteredTimeSeries = filteredTimeSeries,
      temporalvar = temporalvar, network = labmat,
      graph = gmet, corrmat = ocormat, partialcorrmat = pcormat, glassocormat = gcormat, rcormat = rcormat
    )
    L$mask <- mask
    return(L)
  } else {
    L <- list(filteredTimeSeries = filteredTimeSeries, mask = mask, temporalvar = temporalvar)
    L$mask <- mask
    return(L)
  }
}