import SIMS imaging data into MALDIquant::MassSpectrum objects

importSIMS.R

importSIMS <- function(path) {
  stopifnot(file.exists(path) && file.info(path)$isdir)

  ## find all txt files 
  files <- list.files(path, pattern="^.*\\.txt$", full.names=TRUE, recursive=TRUE)
  files <- normalizePath(files)

  ## read metadata (first 8 lines) of the SIMS file
  ll <- lapply(files, function(f)readLines(f, n=8))

  ## the 2nd line of the SIMS file looks like the following:
  # Source Interval: 424.42 u (ID: 17)
  ## get the 2nd line of each file
  l2 <- unlist(lapply(ll, "[", 2))

  ## we fetch the mass using regular expressions
  mass <- as.double(
            gsub("^# Source Interval: ([[:digit:]\\.]+) *u.*$", "\\1", l2))

  ## the 8th line of the SIMS file looks like the following:
  # Image Size: 100 x 100 pixels
  ## get the 8th line of each file
  l8 <- unlist(lapply(ll, "[", 8))
  ## we fetch the pixel size using regular expressions
  imageSize <- gsub("^# Image Size: ([[:digit:] x]+) *pixels$", "\\1", l8)

  ## split the strings at the " x " and turn into a matrix
  sizes <- do.call(cbind, lapply(strsplit(imageSize, " ?x ?"), as.numeric))

  ## find max values for x and y
  mc <- max.col(sizes)
  nx <- sizes[1, mc[1]]
  ny <- sizes[2, mc[2]]

  ## preallocate memory for a 3D array (x, y, mass)
  arr <- array(0L, dim=c(nx, ny, length(mass)))

  ## fill array
  for (i in seq(along=files)) {
    ## read current mass slice/layer (skip the first 9 lines with comments)
    slice <- read.table(files[i], skip=9, stringsAsFactors=FALSE,
                        col.names=c("x", "y", "intensity"), row.names=NULL,
                        colClasses=c("integer", "integer", "double"))
    ## R start with index 1 and the file with index 0
    indices <- cbind(x=slice[, "x"]+1, y=slice[, "y"]+1, mass=i)

    ## copy slice data to array
    arr[indices] <- slice[, "intensity"]
  }
  
  ## sort mass increasing
  o <- order(mass)
  mass <- mass[o]
  ## sort array by mass
  arr <- arr[,,o]


  ## turn 3D array into a list of MALDIquant::MassSpectrum objects
  require("MALDIquant")

  ## preallocate memory 
  l <- vector(mode="list", length=nx*ny)

  ## move through all x/y combinations and create a MALDIquant::MassSpectrum
  ## object for them
  for (i in 1:dim(arr)[1]) {
    for (j in 1:dim(arr)[2]) {
      l[[ ((i-1)*ny)+j ]] <- 
        createMassSpectrum(mass=mass, intensity=arr[i, j, ], 
                           metaData=list( imaging=list(pos=c(x=i-1, y=j-1))))
    }
  }
  return(l)
}