## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, message = FALSE---------------------------------------------------
# MISTy
library(mistyR)
library(future)

# SpatialExperiment
library(SpatialExperiment)
library(SingleCellExperiment)
library(SummarizedExperiment)

# data manipulation
library(Matrix)
library(tibble)
library(dplyr)
library(purrr)

# normalization
library(sctransform)

# resource
library(progeny)

# plotting
library(ggplot2)

# setup parallel execution
plan(multisession)

## -----------------------------------------------------------------------------
run_misty_spe <- function(slide,
                          # SpatialExperiment object with spatial omics data.
                          view.assays,
                          # Named list of assays for each view.
                          view.features = NULL,
                          # Named list of features/markers to use.
                          # Use all by default.
                          view.types,
                          # Named list of the type of view to construct
                          # from the assay.
                          view.params,
                          # Named list with parameters (NULL or value)
                          # for each view.
                          spot.ids = NULL,
                          # spot IDs to use. Use all by default.
                          out.alias = "results"
                          # folder name for output
) {

  # Extracting geometry
  geometry <- as.data.frame(spatialData(slide)) %>%
    select(array_row, array_col)

  # Extracting data
  view.data <- map(view.assays,
    extract_spe_data,
    geometry = geometry,
    slide = slide
  )

  # Constructing and running a workflow
  build_misty_pipeline(
    view.data = view.data,
    view.features = view.features,
    view.types = view.types,
    view.params = view.params,
    geometry = geometry,
    spot.ids = spot.ids,
    out.alias = out.alias
  )
}

## -----------------------------------------------------------------------------
# Extracts data from an specific assay from a SpatialExperiment object
# and aligns the IDs to the geometry
extract_spe_data <- function(slide,
                             assay,
                             geometry) {
  data <- altExp(slide, assay) %>%
    assay() %>%
    t() %>%
    as_tibble(rownames = NA)

  return(data %>% dplyr::slice(match(rownames(.), rownames(geometry))))
}


# Filters data to contain only features of interest
filter_data_features <- function(data,
                                 features) {
  if (is.null(features)) features <- colnames(data)

  return(data %>% rownames_to_column() %>%
    select(rowname, all_of(features)) %>% rename_with(make.names) %>%
    column_to_rownames())
}

## -----------------------------------------------------------------------------
# Builds views depending on the paramaters defined
create_default_views <- function(data,
                                 view.type,
                                 view.param,
                                 view.name,
                                 spot.ids,
                                 geometry) {
  view.data.init <- create_initial_view(data)

  if (!(view.type %in% c("intra", "para", "juxta"))) {
    view.type <- "intra"
  }

  if (view.type == "intra") {
    data.red <- view.data.tmp$data %>%
      rownames_to_column() %>%
      filter(rowname %in% spot.ids) %>%
      select(-rowname)
  } else if (view.type == "para") {
    view.data.tmp <- view.data.init %>%
      add_paraview(geometry, l = view.param)

    data.ix <- paste0("paraview.", view.param)
    data.red <- view.data.tmp[[data.ix]]$data %>%
      mutate(rowname = rownames(data)) %>%
      filter(rowname %in% spot.ids) %>%
      select(-rowname)
  } else if (view.type == "juxta") {
    view.data.tmp <- view.data.init %>%
      add_juxtaview(
        positions = geometry,
        neighbor.thr = view.param
      )

    data.ix <- paste0("juxtaview.", view.param)
    data.red <- view.data.tmp[[data.ix]]$data %>%
      mutate(rowname = rownames(data)) %>%
      filter(rowname %in% spot.ids) %>%
      select(-rowname)
  }

  if (is.null(view.param) == TRUE) {
    misty.view <- create_view(
      paste0(view.name),
      data.red
    )
  } else {
    misty.view <- create_view(
      paste0(view.name, "_", view.param),
      data.red
    )
  }

  return(misty.view)
}

## -----------------------------------------------------------------------------
# Builds automatic MISTy workflow and runs it
build_misty_pipeline <- function(view.data,
                                 view.features,
                                 view.types,
                                 view.params,
                                 geometry,
                                 spot.ids = NULL,
                                 out.alias = "default") {

  # Adding all spots ids in case they are not defined
  if (is.null(spot.ids)) {
    spot.ids <- rownames(view.data[[1]])
  }

  # First filter the features from the data
  view.data.filt <- map2(view.data, view.features, filter_data_features)

  # Create initial view
  views.main <- create_initial_view(view.data.filt[[1]] %>%
    rownames_to_column() %>%
    filter(rowname %in% spot.ids) %>% select(-rowname))

  # Create other views
  view.names <- names(view.data.filt)

  all.views <- pmap(list(
    view.data.filt[-1],
    view.types[-1],
    view.params[-1],
    view.names[-1]
  ),
  create_default_views,
  spot.ids = spot.ids,
  geometry = geometry
  )

  pline.views <- add_views(
    views.main,
    unlist(all.views, recursive = FALSE)
  )


  # Run MISTy
  run_misty(pline.views, out.alias)
}

## -----------------------------------------------------------------------------
plotMolecules_adapted <- function(spe,
                                  molecule = NULL,
                                  x_coord = "array_col",
                                  y_coord = "array_row",
                                  palette = NULL,
                                  alt_assay = "lognorm") {
  if (is.null(palette)) {
    palette <- "yellow"
  }
  if (!is.null(palette) && length(palette) == 1) {
    palette <- c("black", palette)
  }
  df_plot <- spatialData(spe)[, c(x_coord, y_coord), drop = FALSE]
  mRNA_counts <-
    as.numeric(assay(altExp(spe, alt_assay))[molecule, , drop = FALSE])
  stopifnot(length(mRNA_counts) == nrow(df_plot))

  df_plot <- cbind(df_plot, expression = mRNA_counts)

  df_plot <- as.data.frame(df_plot) %>% mutate(array_row = array_row * -1)

  p <- ggplot(
    df_plot,
    aes_string(x = x_coord, y = y_coord, color = "expression")
  ) +
    geom_point(size = 2.5) +
    scale_color_gradient(
      low = palette[1],
      high = palette[2], trans = "sqrt"
    ) +
    coord_fixed() +
    ggtitle(molecule) +
    theme_void()
  p
}

## ----get_data, include=FALSE--------------------------------------------------
cleanup <- FALSE
if (!dir.exists("breast_A_1")) {
  download.file("https://www.dropbox.com/s/igby4csbt9u4uuf/breast_A_1.tgz?dl=1",
    destfile = "breast_A_1.tgz", mode ="wb", quiet = TRUE
  )
  untar("breast_A_1.tgz", tar = "internal")
  cleanup <- TRUE
}

## ---- warning=FALSE-----------------------------------------------------------
# Load and normalize using SCT
folder <- "breast_A_1"

# rename to create the required file structure
file.rename(
  "breast_A_1/V1_Breast_Cancer_Block_A_Section_1_filtered_feature_bc_matrix.h5",
  "breast_A_1/filtered_feature_bc_matrix.h5"
)

spe <- read10xVisium(folder, type = "HDF5", data = "filtered", 
                     images = "lowres")


# normalize data
# counts are of class DelayedMatrix which is incompatible with vst
sct.data <- vst(as(counts(spe), "dgCMatrix"), verbosity = 0)$y

# Dropping duplicates
spe <- spe[!duplicated(rowData(spe)), ]

# Getting relevant genes
gene.dict <- as_tibble(rowData(spe), rownames = NA) %>%
  rownames_to_column("ENSEMBL") %>%
  filter(ENSEMBL %in% rownames(sct.data))

# Re-naming normalized data with gene symbols
sct.data <- sct.data[gene.dict %>% pull("ENSEMBL"), colnames(spe)]
rownames(sct.data) <- gene.dict %>% pull(symbol)

# undo file rename
file.rename(
  "breast_A_1/filtered_feature_bc_matrix.h5",
  "breast_A_1/V1_Breast_Cancer_Block_A_Section_1_filtered_feature_bc_matrix.h5"
)

## -----------------------------------------------------------------------------
coverage <- rowSums(sct.data > 0) / ncol(sct.data)
slide.markers <- names(coverage[coverage >= 0.05])

## -----------------------------------------------------------------------------
estrogen.footprints <- getModel(top = 15) %>%
  rownames_to_column("gene") %>%
  filter(Estrogen != 0, gene %in% slide.markers) %>%
  pull(gene)

hypoxia.footprints <- getModel(top = 15) %>%
  rownames_to_column("gene") %>%
  filter(Hypoxia != 0, gene %in% slide.markers) %>%
  pull(gene)

## -----------------------------------------------------------------------------
sct.exp <- SummarizedExperiment(sct.data[slide.markers, ])
assayNames(sct.exp) <- "SCT"
altExp(spe, "SCT") <- sct.exp

## -----------------------------------------------------------------------------
# Define assay for each view
view.assays <- list(
  "main" = "SCT",
  "para.hypoxia" = "SCT",
  "para.estrogen" = "SCT"
)

# Define features for each view
view.features <- list(
  "main" = hypoxia.footprints,
  "para.hypoxia" = hypoxia.footprints,
  "para.estrogen" = estrogen.footprints
)

# Define spatial context for each view
view.types <- list(
  "main" = "intra",
  "para.hypoxia" = "para",
  "para.estrogen" = "para"
)

# Define additional parameters (l in the case of paraview)
view.params <- list(
  "main" = NULL,
  "para.hypoxia" = 10,
  "para.estrogen" = 10
)

misty.out <- "vignette_model_spe"

## ---- warning=FALSE-----------------------------------------------------------
misty.results <- run_misty_spe(
  slide = spe,
  view.assays = view.assays,
  view.features = view.features,
  view.types = view.types,
  view.params = view.params,
  spot.ids = NULL, # Using the whole slide
  out.alias = misty.out
) %>%
  collect_results()

## ---- warning=FALSE-----------------------------------------------------------
misty.results %>%
  plot_improvement_stats("gain.R2") %>%
  plot_improvement_stats("gain.RMSE")

## -----------------------------------------------------------------------------
misty.results$improvements %>%
  filter(measure == "p.R2") %>%
  arrange(value)

## -----------------------------------------------------------------------------
misty.results %>% plot_view_contributions()

misty.results$contributions.stats %>% filter(target == "PGK1")

## -----------------------------------------------------------------------------
misty.results %>% plot_interaction_heatmap(view = "intra")

## -----------------------------------------------------------------------------
misty.results$importances.aggregated %>% 
  filter(view == "intra", Target == "PGK1") %>%
  arrange(-Importance)


## ---- warning=FALSE, dev='jpeg'-----------------------------------------------
plotMolecules_adapted(spe,
  molecule = "PGK1",
  x_coord = "array_col",
  y_coord = "array_row",
  alt_assay = "SCT"
)

## ---- warning=FALSE, dev='jpeg'-----------------------------------------------
plotMolecules_adapted(spe,
  molecule = "NDRG1",
  x_coord = "array_col",
  y_coord = "array_row",
  alt_assay = "SCT"
)

## -----------------------------------------------------------------------------
misty.results %>% plot_interaction_heatmap(view = "para.hypoxia_10")

## ---- warning=FALSE, dev='jpeg'-----------------------------------------------
plotMolecules_adapted(spe,
  molecule = "PGK1",
  x_coord = "array_col",
  y_coord = "array_row",
  alt_assay = "SCT"
)

## ---- warning=FALSE, dev='jpeg'-----------------------------------------------
plotMolecules_adapted(spe,
  molecule = "PFKFB4",
  x_coord = "array_col",
  y_coord = "array_row",
  alt_assay = "SCT"
)

## -----------------------------------------------------------------------------
misty.results %>% plot_interaction_heatmap(view = "para.estrogen_10")

## ---- warning=FALSE, dev='jpeg'-----------------------------------------------
plotMolecules_adapted(spe,
  molecule = "PGK1",
  x_coord = "array_col",
  y_coord = "array_row",
  alt_assay = "SCT"
)

## ---- warning=FALSE, dev='jpeg'-----------------------------------------------
plotMolecules_adapted(spe,
  molecule = "TPD52L1",
  x_coord = "array_col",
  y_coord = "array_row",
  alt_assay = "SCT"
)

## ----info, echo=FALSE---------------------------------------------------------
sessionInfo()

## ----cleanup, include=FALSE---------------------------------------------------
if (cleanup) {
  unlink(c("breast_A_1.tgz", "breast_A_1/"),
    recursive = TRUE
  )
}
unlink("vignette_model_spe", recursive = TRUE)