R/aba.R

In hypercompetent/xstitch: Cross-stitch pattern generation.

#' Query the Allen Brain Atlas API for experiment IDs based on gene symbol.
#'
#' @param gene_symbol A character vector containing a gene symbol.
#' @param plane The plane of the ISH experiment to search for. Options are "coronal" and "saggital". Default is "coronal".
#'
#' @return a character vector of valid experiment IDs
#' @export
#'
#' @examples
#' ids <- aba_find_gene_ids("Pvalb", "coronal"
#' )
aba_find_gene_ids <- function(gene_symbol,
                             plane = "coronal") {
  query <- paste0(
    "http://api.brain-map.org/api/v2/data/query.xml?criteria=
     model::SectionDataSet,
     rma::criteria,[failed$eq'false'],
     products[abbreviation$eq'Mouse'],
     plane_of_section[name$eq'", plane, "'],
     genes[acronym$eq'", gene_symbol, "']"
  )
  query <- gsub("[ \n]+", "", query)

  id_url <- xml2::read_xml(query)
  id_nodes <- xml2::xml_find_all(id_url, xpath = "//Response/section-data-sets/section-data-set/id")
  ids <- as.character(xml_contents(id_nodes))

  ids
}

#' Query the Allen Brain Atlas API to retrieve an ISH energy array for an experiment ID
#'
#' @param id ISH experiment id, like that generated by aba_find_gene_ids()
#'
#' @return A 3-dimensional array of energy values with dims = c(67, 41, 58).
#'
#' @export
#'
#' @examples
#' library(dplyr)
#'
#' pvalb_arr <- aba_find_gene_ids("Pvalb") %>%
#'   aba_get_energy_array()
#'
aba_get_energy_array <- function(id) {
  if(length(id) > 1) {
    warning("Length of id > 1. Will use the first id value")
    id <- id[1]
  }

  array_dims <- c(67, 41, 58)

  query <- paste0(
    "http://api.brain-map.org/grid_data/download/",
    id,
    "?include=energy"
  )
  query <- gsub("[ \n]+", "", query)

  temp <- tempfile(fileext = ".zip")

  download.file(query, temp, mode = "wb")

  raw_file <- unz(temp, "energy.raw", "rb")
  vol_raw <- readBin(raw_file, "double", size = 4, n = prod(array_dims))
  close(raw_file)

  file.remove(temp)

  array(vol_raw, dim = array_dims)
}

#' Slice an Allen Brain Atlas ISH Array to return a matrix from a single plane.
#'
#' @param arr The 3-dimensional ISH array to slice
#' @param slice A numeric value for the slice to use.
#' @param plane Which plane to return. Options are "coronal", "saggital", and "horizontal". Default is "coronal"
#'
#' @return A 2-dimensional array of ISH values
#' @export
#'
#' @examples
#' library(dplyr)
#'
#' pvalb_mat <- aba_find_gene_ids("Pvalb") %>%
#'   aba_get_energy_array() %>%
#'   aba_slice_array(42, "coronal")
#'
aba_slice_array <- function(arr,
                            slice,
                            plane = "coronal") {

  if(plane == "coronal") {
    arr[slice, , ]
  } else if(plane == "saggital") {
    arr[, slice, ]
  } else if(plane == "horizontal") {
    arr[, , slice]
  }

}

#' Retrieve the Allen Brain Atlas Annotation ID array for Gridded expression data.
#'
#' @return a 3-dimensional array with dims = c(67, 41, 58)
#' @export
#'
#' @examples
#' anno_arr <- aba_get_annotation_array()
#'
aba_get_annotation_array <- function() {
  array_dims <- c(67, 41, 58)

  temp <- tempfile()

  query <- "http://download.alleninstitute.org/informatics-archive/current-release/mouse_annotation/P56_Mouse_gridAnnotation.zip"
  download.file(query, temp)

  raw_file <- unz(temp, "gridAnnotation.raw", "rb")
  vol_raw <- readBin(raw_file, "integer", size = 4, n = prod(array_dims))
  close(raw_file)

  file.remove(temp)

  array(vol_raw, dim = array_dims)
}


aba_filter_anno_df <- function(anno,
                              slice_num,
                              direction = "coronal") {
  library(dplyr)

  if(direction == "coronal") {
    anno <- anno[anno$x == slice_num,]
    anno <- select(anno, -x)
    names(anno)[1:2] <- c("y","x")
  } else if(direction == "horizontal") {
    anno <- anno[anno$y == slice_num,]
    anno <- select(anno, -y)
    names(anno)[1:2] <- c("y","x")
  } else if(direction == "saggital") {
    anno <- anno[anno$z == slice_num,]
    anno <- select(anno, -z)
    names(anno)[1:2] <- c("y","x")
  }
  anno
}

#' Generate array annotation boundary segments for ggplot2
#'
#' @param anno A 3-d annotation array like that generated by aba_get_annotation_array()
#' @param slice_num A number indicating which slice to use.
#' @param plane Which plane to return. Options are "coronal", "saggital", and "horizontal". Default is "coronal"
#'
#' @return a data.frame of segment coordinates with x, xend, y, and yend.
#' @export
#'
#' @examples
#' library(dplyr)
#' library(ggplot2)
#'
#' boundary_df <- aba_get_annotation_array() %>%
#'   aba_symmetrize_matrix("min") %>%
#'   aba_melt_3d("atlas_id") %>%
#'   filter(atlas_id != 0) %>%
#'   aba_slice_anno_boundaries(30, "coronal")
#'
#' ggplot(boundary_df) +
#'   geom_segment(aes(x = x, xend = xend,
#'                    y = y, yend = yend)) +
#'   scale_y_reverse()
#'
aba_slice_anno_boundaries <- function(anno,
                                      slice_num,
                                      plane = "coronal") {

  slice_anno <- aba_filter_anno_df(anno, slice_num, plane) %>%
    rename(atlas_x = x,
           atlas_y = y)

  slice_anno <- slice_anno %>%
    mutate(atlas_id = ifelse(atlas_x < 30, paste0(atlas_id, "_r"), paste0(atlas_id, "_l")))

  min_x_lines <- slice_anno %>%
    arrange(atlas_x, atlas_y) %>%
    mutate(border = lag(atlas_id) != atlas_id | lag(atlas_y, default = -100) != atlas_y - 1) %>%
    filter(border) %>%
    mutate(y = atlas_y - 0.5,
           yend = atlas_y - 0.5) %>%
    mutate(x = atlas_x - 0.5,
           xend = atlas_x + 0.5) %>%
    select(x, xend, y, yend) %>%
    unique()

  max_x_lines <- slice_anno %>%
    arrange(atlas_x, atlas_y) %>%
    mutate(border = lead(atlas_id) != atlas_id | lead(atlas_y, default = -100) != atlas_y + 1) %>%
    filter(border) %>%
    mutate(y = atlas_y + 0.5,
           yend = atlas_y + 0.5) %>%
    mutate(x = atlas_x - 0.5,
           xend = atlas_x + 0.5) %>%
    select(x, xend, y, yend) %>%
    unique()


  min_y_lines <- slice_anno %>%
    arrange(atlas_y, atlas_x) %>%
    mutate(border = lag(atlas_id) != atlas_id | lag(atlas_x, default = -100) != atlas_x - 1) %>%
    filter(border) %>%
    mutate(x = atlas_x - 0.5,
           xend = atlas_x - 0.5) %>%
    mutate(y = atlas_y - 0.5,
           yend = atlas_y + 0.5) %>%
    select(x, xend, y, yend) %>%
    unique()

  max_y_lines <- slice_anno %>%
    arrange(atlas_y, atlas_x) %>%
    mutate(border = lead(atlas_id) != atlas_id | lead(atlas_x, default = -100) != atlas_x + 1) %>%
    filter(border) %>%
    mutate(x = atlas_x + 0.5,
           xend = atlas_x + 0.5) %>%
    mutate(y = atlas_y - 0.5,
           yend = atlas_y + 0.5) %>%
    select(x, xend, y, yend) %>%
    unique()

  all_lines <- rbind(min_x_lines,
                     max_x_lines,
                     min_y_lines,
                     max_y_lines)

}

#' Get the Allen Brain Atlas structure graph ontology from the Allen API
#'
#' @return a nested list containing the allen brain atlas ontology.
#' @export
#'
#' @examples
#' ont <- aba_get_ontology()
#'
aba_get_ontology <- function() {

  # Download the ontology JSON file
  temp <- tempfile()
  download.file("http://api.brain-map.org/api/v2/structure_graph_download/1.json", temp)

  # Read the JSON
  raw_ontology <- jsonlite::fromJSON(temp)[["msg"]]

  return(raw_ontology)
}

#' Convert the nested Allen Brain Atlas ontology to a data.frame.
#'
#' @param @ontology The list-based ontology returned by the aba_get_ontology() function.
#' @param ontology_df The current ABA ontology data.frame. Used for recursive unpacking of the list ontology.
#'
#' @return a data.frame containing the ontology, along with a parent_id column for retaining hierarchical information.
#' @export
#'
#' @examples
#' library(dplyr)
#'
#' ont_df <- aba_get_ontology() %>%
#'   aba_flatten_ontology()
#'
aba_flatten_ontology <- function(ontology, ontology_df = NULL) {

  l <- ontology

  if(is.null(ontology_df)) {
    ontology_df <- data.frame(l[names(l) != "children"])[0,]
    ontology_df$n_children <- numeric()
  }

  if("children" %in% names(l)) {

    child_df <- data.frame(l[names(l) != "children"])

    n_children_of_children <- map_dbl(l$children,
                                      function(x) {
                                        if("children" %in% names(x)) {
                                          length(x$children)
                                        } else {
                                          0
                                        }
                                      })

    child_df$n_children <- n_children_of_children

    ontology_df <- rbind(ontology_df, child_df)

    for(i in 1:length(l$children)) {

      child_list <- l$children[[i]]

      ontology_df <- aba_flatten_ontology(child_list, ontology_df)
    }
  }

  return(ontology_df)
}

#' Convert an Allen Brain Atlas 3-d array to a data.frame
#'
#' @param arr a 3-d array to melt
#' @param val A name for the value column. Default is "value".
#'
#' @return a data.frame with columns x, y, z, and val (specified as a parameter).
#' @export
#'
#'
aba_melt_3d <- function(arr,
                        val = "value") {
  df <- reshape2::melt(arr)
  names(df) <- c("x","y","z",val)
  df
}

#' Convert an Allen Brain Atlas 2-d array to a data.frame
#'
#' @param arr a 2-d array to melt, like that generated by aba_slice_array()
#' @param val The name of the value column to use.
#' @param plane The plane of the 2-d array. Options are "coronal", "saggital", and "horizontal". Default is "coronal".
#'
#' @return a data.frame with dimension columns based on the plane, and a value column named for val.
#' @export
#'
#'
aba_melt_2d <- function(arr,
                        val = "value",
                        plane = "coronal") {
  if(plane == "coronal") {
    dim_names <- c("x","y")
  } else if(plane == "saggital") {
    dim_names <- c("y","z")
  } else if(plane == "horizontal") {
    dim_names <- c("x","z")
  }

  df <- reshape2::melt(arr)
  names(df) <- c(dim_names, val)

  df
}

#' Apply hemisphere symmetry to an Allen Brain Atlas array
#'
#' @param arr The ISH array to symmetrize
#' @param fun The function to use for symmetry calculations. Options are "mean", "max", "min", and "mult".
#'
#' @return a 3-d array with the same dimensions as the input
#' @export
#'
aba_symmetrize_array <- function(arr,
                                 fun = "none") {
  if(fun == "none") {
    arr
  } else if(fun == "mean") {
    #(arr1 + arr2)/2
    out_arr <- arr
    out_arr[, ,29:1] <- out_arr[, ,30:58] <- (out_arr[, ,29:1] + out_arr[, ,30:58]) / 2
    out_arr
  } else if(fun == "max") {
    #pmax
    out_arr <- arr
    out_arr[, ,29:1] <- out_arr[, ,30:58] <- pmax(out_arr[, ,29:1], out_arr[, ,30:58])
    out_arr
  } else if(fun == "min") {
    #pmin
    out_arr <- arr
    out_arr[, ,29:1] <- out_arr[, ,30:58] <- pmin(out_arr[, ,29:1], out_arr[, ,30:58])
    out_arr
  } else if(fun == "mult") {
    #arr1 * arr2
    out_arr <- arr
    out_arr[, ,29:1] <- out_arr[, ,30:58] <- out_arr[, ,29:1] * out_arr[, ,30:58]
    out_arr
  }
}