## ----global.options, include = FALSE------------------------------------------
knitr::opts_knit$set(
    collapse = TRUE,
    comment = "#>",
    fig.align   = 'center'
)

knitr::opts_chunk$set(out.extra = 'style="display:block; margin:auto;"')


## ----eval = FALSE-------------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("SpaceMarkers")

## ----message = FALSE, warning = FALSE-----------------------------------------
library(SpaceMarkers)

## -----------------------------------------------------------------------------
data_dir <- "visiumBrCA"
unlink(file.path(data_dir), recursive = TRUE)
dir.create(data_dir,showWarnings = FALSE)
main_10xlink <- "https://cf.10xgenomics.com/samples/spatial-exp/1.3.0"
counts_folder <- "Visium_Human_Breast_Cancer"
counts_file <- "Visium_Human_Breast_Cancer_filtered_feature_bc_matrix.h5"
counts_url<-paste(c(main_10xlink,counts_folder,counts_file), collapse = "/")
sp_folder <- "Visium_Human_Breast_Cancer"
sp_file <- "Visium_Human_Breast_Cancer_spatial.tar.gz"
sp_url<-paste(c(main_10xlink,sp_folder,sp_file),collapse = "/")

## -----------------------------------------------------------------------------
download.file(counts_url,file.path(data_dir,basename(counts_url)), mode = "wb")
counts_matrix <- load10XExpr(visiumDir = data_dir, h5filename = counts_file)
good_gene_threshold <- 3
goodGenes <- rownames(counts_matrix)[
    apply(counts_matrix,1,function(x) sum(x>0)>=good_gene_threshold)]
counts_matrix <- counts_matrix[goodGenes,]

## -----------------------------------------------------------------------------
cogaps_result <- readRDS(system.file("extdata","CoGAPS_result.rds",
    package="SpaceMarkers",mustWork = TRUE))
features <- intersect(rownames(counts_matrix),rownames(
    slot(cogaps_result,"featureLoadings")))
barcodes <- intersect(colnames(counts_matrix),rownames(
    slot(cogaps_result,"sampleFactors")))
counts_matrix <- counts_matrix[features,barcodes]
cogaps_matrix<-slot(cogaps_result,"featureLoadings")[features,]%*%
    t(slot(cogaps_result,"sampleFactors")[barcodes,])

## -----------------------------------------------------------------------------
download.file(sp_url, file.path(data_dir,basename(sp_url)), mode = "wb")
untar(file.path(data_dir,basename(sp_url)), exdir = file.path(data_dir))
spCoords <- load10XCoords(visiumDir = data_dir, version = "1.0")
rownames(spCoords) <- spCoords$barcode
spCoords <- spCoords[barcodes,]
spPatterns <- cbind(spCoords,slot(cogaps_result,"sampleFactors")[barcodes,])
head(spPatterns)

## -----------------------------------------------------------------------------
data("curated_genes")
spPatterns <- spPatterns[c("barcode","y","x","Pattern_1","Pattern_5")]
counts_matrix <- counts_matrix[curated_genes,]
cogaps_matrix <- cogaps_matrix[curated_genes, ]

## -----------------------------------------------------------------------------
optParams <- getSpatialParameters(spPatterns,visiumDir = data_dir,
                                          resolution = "lowres")

## -----------------------------------------------------------------------------
SpaceMarkers <- getPairwiseInteractingGenes(data = counts_matrix,
                                    reconstruction = cogaps_matrix,
                                    optParams = optParams,
                                    spPatterns = spPatterns,
                                    mode ="residual",analysis="overlap")


## -----------------------------------------------------------------------------
print(head(SpaceMarkers[[1]]$interacting_genes[[1]]))
print(head(SpaceMarkers[[1]]$hotspots))

## -----------------------------------------------------------------------------
SpaceMarkers_DE <- getPairwiseInteractingGenes(
    data=counts_matrix,reconstruction=NULL,
    optParams = optParams,
    spPatterns = spPatterns,
    mode="DE",analysis="overlap")

## -----------------------------------------------------------------------------
residual_p1_p5<-SpaceMarkers[[1]]$interacting_genes[[1]]
DE_p1_p5<-SpaceMarkers_DE[[1]]$interacting_genes[[1]]

## -----------------------------------------------------------------------------
paste(
    "Residual mode identified",dim(residual_p1_p5)[1],
        "interacting genes,while DE mode identified",dim(DE_p1_p5)[1],
        "interacting genes",collapse = NULL)

## -----------------------------------------------------------------------------

compare_genes <- function(ref_list, list2,ref_name = "mode1",
                            list2_name = "mode2", sub_slice = NULL){
    ref_rank <- seq(1,length(ref_list),1)
    list2_ref_rank <- which(list2 %in% ref_list)
    list2_ref_genes <- list2[which(list2 %in% ref_list)]
    ref_genes_only <- ref_list[ !ref_list  %in% list2_ref_genes ]
    mode1 <- data.frame("Gene" = ref_list,"Rank" = ref_rank,"mode"= ref_name)
    mode2 <- data.frame("Gene" = c(list2_ref_genes, ref_genes_only),"Rank" = c(
        list2_ref_rank,rep(NA,length(ref_genes_only))),"mode"= list2_name)
    mode1_mode2 <- merge(mode1, mode2, by = "Gene", all = TRUE) 
    mode1_mode2 <- mode1_mode2[order(mode1_mode2$Rank.x),]
    mode1_mode2 <- subset(mode1_mode2,select = c("Gene","Rank.x","Rank.y"))
    colnames(mode1_mode2) <- c("Gene",paste0(ref_name,"_Rank"),
                                paste0(list2_name,"_Rank"))
    return(mode1_mode2)
}

## -----------------------------------------------------------------------------
res_to_DE <- compare_genes(head(residual_p1_p5$Gene, n = 20),DE_p1_p5$Gene,
                            ref_name="residual",list2_name="DE")
DE_to_res <- compare_genes(head(DE_p1_p5$Gene, n = 20),residual_p1_p5$Gene,
                            ref_name = "DE",list2_name = "residual")

## -----------------------------------------------------------------------------
res_to_DE

## -----------------------------------------------------------------------------
DE_to_res

## -----------------------------------------------------------------------------
SpaceMarkers_enrich <- getPairwiseInteractingGenes(data = counts_matrix,
                                    reconstruction = cogaps_matrix,
                                    optParams = optParams,
                                    spPatterns = spPatterns,
                                    mode ="residual",analysis="enrichment")
SpaceMarkers_DE_enrich <- getPairwiseInteractingGenes(
    data=counts_matrix,reconstruction=NULL,
    optParams = optParams,
    spPatterns = spPatterns,
    mode="DE",analysis="enrichment")
residual_p1_p5_enrichment<-SpaceMarkers_enrich[[1]]$interacting_genes[[1]]$Gene
DE_p1_p5_enrichment<-SpaceMarkers_DE_enrich[[1]]$interacting_genes[[1]]$Gene


## -----------------------------------------------------------------------------
enrich_res_to_de<-compare_genes(
    head(DE_p1_p5_enrichment, 20),
    residual_p1_p5_enrichment,
    ref_name="DE_Enrich",list2_name = "res_Enrich")
enrich_res_to_de

## -----------------------------------------------------------------------------
overlap_enrich_de<-compare_genes(
    head(DE_p1_p5_enrichment,20),
    DE_p1_p5$Gene,
    ref_name="DE_Enrich",
    list2_name="DE_Overlap")
overlap_enrich_de

## -----------------------------------------------------------------------------
tail(SpaceMarkers_DE_enrich[[1]]$interacting_genes[[1]])

## -----------------------------------------------------------------------------
overlap_enrich_res<-compare_genes(
    head(residual_p1_p5$Gene, 20),
    residual_p1_p5_enrichment,
    ref_name ="res_overlap",list2_name="res_enrich")
overlap_enrich_res

## ----message = FALSE, warning=FALSE-------------------------------------------
library(Matrix)
library(rjson)
library(cowplot)
library(RColorBrewer)
library(grid)
library(readbitmap)
library(dplyr)
library(data.table)
library(viridis)
library(hrbrthemes)
library(ggplot2)

## -----------------------------------------------------------------------------
res_enrich <- SpaceMarkers_enrich[[1]]$interacting_genes[[1]]
hotspots <- SpaceMarkers_enrich[[1]]$hotspots
top <- res_enrich %>% arrange(-SpaceMarkersMetric)
print(head(top))

## -----------------------------------------------------------------------------
geom_spatial <-  function(mapping = NULL,
                          data = NULL,
                          stat = "identity",
                          position = "identity",
                          na.rm = FALSE,
                          show.legend = NA,
                          inherit.aes = FALSE,...) {
  GeomCustom <- ggproto(
    "GeomCustom",
    Geom,
    setup_data = function(self, data, params) {
      data <- ggproto_parent(Geom, self)$setup_data(data, params)
      data
    },
    draw_group = function(data, panel_scales, coord) {
      vp <- grid::viewport(x=data$x, y=data$y)
      g <- grid::editGrob(data$grob[[1]], vp=vp)
      ggplot2:::ggname("geom_spatial", g)
    },
    required_aes = c("grob","x","y")
  )
  layer(geom = GeomCustom,mapping = mapping,data = data,stat = stat,position=
          position,show.legend = show.legend,inherit.aes =
          inherit.aes,params = list(na.rm = na.rm, ...)
  )
}


## -----------------------------------------------------------------------------
plotSpatialData <- function(spatial_data,positions_path,image_path,jsonPath,
                              feature_col,colors = NULL,
                              title = "Spatial Heatmap",
                              alpha = 0.6,scaled = FALSE,res = "lowres",
                              x_col = "x",y_col = "y",sample = "Sample",
                              plot = TRUE) {
  og_positions <- read.csv(positions_path, header = FALSE)
  og_positions <- og_positions[,c(1,5,6)]
  colnames(og_positions) <- c("barcode",y_col,x_col)
  spatial_data <- dplyr::inner_join(
    og_positions,spatial_data[,c("barcode",feature_col)],by = "barcode")
  images_cl <- readbitmap::read.bitmap(image_path)
  height <-  data.frame(height = nrow(images_cl))
  width <- data.frame(width = ncol(images_cl))
  grobs<-grid::rasterGrob(images_cl,width=unit(1,"npc"),
                          height=unit(1,"npc"))
  images_tibble <- dplyr::tibble(sample=factor(sample), grob=list(grobs))
  images_tibble$height <- height$height
  images_tibble$width <- width$width
  scales <- jsonlite::read_json(jsonPath)
  if (scaled == FALSE){
    res <- names(scales)[grepl(pattern = res,x = names(scales))]
    spatial_data[[x_col]]<-  scales[[res]] * spatial_data[[x_col]]
    spatial_data[[y_col]] <- scales[[res]] * spatial_data[[y_col]]
  }
  spatial_data$height <- height$height
  spatial_data$width <- width$width
  if (is.numeric(spatial_data[[feature_col]])){
    p <- spatial_data %>% ggplot(aes(
      x=.data[[x_col]], y=.data[[y_col]], fill=.data[[feature_col]],
      alpha = .data[[feature_col]])) + geom_spatial(data=images_tibble[1,],
                                                    mapping = aes(grob=grob),
                                                    x=0.5, y=0.5) +
      geom_point(shape = 21, colour = "black", size = 1.1, stroke = 0.2)
  } else {
    p <- spatial_data %>% ggplot(aes(
      x=.data[[x_col]], y=.data[[y_col]], fill=.data[[feature_col]]))+
      geom_spatial(data=images_tibble[1,],mapping = aes(grob=grob),
                   x=0.5, y=0.5) +
      geom_point(shape = 21, colour = "black", size = 1.1, stroke = 0.2,
                 alpha = alpha)
  }
  p <- p + coord_cartesian(expand=FALSE) + 
    xlim(0,max(spatial_data$width)) + ylim(max(spatial_data$height),0) +
    xlab("") + ylab("")+
    ggtitle(title) + theme_set(theme_bw(base_size = 10)) +
    theme(
      panel.grid.major=element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      axis.line = element_line(colour="black"),
      axis.text=element_blank(),axis.ticks = element_blank())
  # Check if the feature_col is continuous or discrete and apply 
  #the appropriate color scale
  if (is.numeric(spatial_data[[feature_col]]) & !is.null(colors)) {
    # Use gradient scale for continuous variables when colors are provided
    p <- p + scale_fill_gradientn(colours = colors)
  } else if (is.numeric(spatial_data[[feature_col]]) & is.null(colors)) {
    # Default to a color palette when no colors are provided for 
    #continuous variables
    p <- p + scale_fill_gradientn(colours = rev(
      RColorBrewer::brewer.pal(name = "RdYlBu", n = 11)))
  } else if (!is.numeric(spatial_data[[feature_col]])) {
    # Use manual scale for discrete variables
    if (!is.null(colors)) {
      p <- p + scale_fill_manual(values = colors)
    } else {
      features <- unique(spatial_data[[feature_col]])
      features <- as.character(features[!is.na(features)])
      if (length(features) > 1){
        message("You can specify your colors. Ensure the length is equal to the
                number of unique features in your feature_col")
      } else {
        p <- p + scale_fill_manual(values = "red")
      }
    }
  }
  if (plot == TRUE) {
    print(p)
  } 
  return(p)
}

## -----------------------------------------------------------------------------
createInteractCol <- function(spHotspots, 
                              interaction_cols = c("T.cells","B-cells")){
  col1 <- spHotspots[,interaction_cols[1]]
  col2 <- spHotspots[,interaction_cols[2]]
  one <- col1
  two <- col2
  one[!is.na(col1)] <- "match"
  two[!is.na(col2)] <- "match"
  both_idx <- which(one == two)
  both <- col1
  both[both_idx] <- "interacting"
  one_only <- setdiff(which(!is.na(col1)),unique(c(which(is.na(col1)),
                                                   both_idx)))
  two_only <- setdiff(which(!is.na(col2)),unique(c(which(is.na(col2)),
                                                   both_idx)))
  both[one_only] <- interaction_cols[1]
  both[two_only] <- interaction_cols[2]
  both <- factor(both,levels = c(interaction_cols[1],"interacting",
                                 interaction_cols[2]))
  return(both)
}

#NB: Since we are likely to plot multipe genes, this function assumes an
#already transposed counts matrix. This saves time and memory in the long run
#for larger counts matrices
plotSpatialExpr <- function(data,gene,hotspots,patterns,
                               remove.na = TRUE,
                               title = "Expression (Log)"){
  counts <- data
  interact <- createInteractCol(spHotspots = hotspots,
                                interaction_cols = patterns)
  df <- cbind(counts,hotspots,data.frame("region" = interact))
  if (remove.na){
    df <- df[!is.na(df$region),]
  }
  p <- df %>% ggplot( aes_string(x='region',y=gene,
                                            fill='region')) + geom_violin() +
    scale_fill_viridis(discrete = TRUE,alpha=0.6) +
    geom_jitter(color="black",size=0.4,alpha=0.9) + theme_ipsum()+
    theme(legend.position="none",plot.title = element_text(size=11),
            axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    ggtitle(paste0(gene,": ",title)) + xlab("")
  return(p)
}




## -----------------------------------------------------------------------------
genes <- top$Gene
counts_df <- as.data.frame(as.matrix(
  t(counts_matrix[rownames(counts_matrix) %in% genes,])))


## -----------------------------------------------------------------------------
image_paths <- file.path(data_dir,"spatial","tissue_lowres_image.png")
scalefactor_paths <- file.path(data_dir,"spatial","scalefactors_json.json")
tissue_paths <- file.path(data_dir,"spatial","tissue_positions_list.csv")

spatialMaps <- list()
exprPlots <- list()

for (g in genes){
  spatialMaps[[length(spatialMaps)+1]] <- plotSpatialData(
    spatial_data = cbind(spPatterns, counts_df), 
                positions_path = tissue_paths,
                image_path = image_paths,
                jsonPath = scalefactor_paths,
                feature_col = g,
                colors = NULL,
                title = g,
                scaled = FALSE, plot = FALSE)
  exprPlots[[length(exprPlots)+1]] <- plotSpatialExpr(
    data = counts_df,gene = g,hotspots = hotspots, 
                   patterns = c("Pattern_1","Pattern_5"))
}

## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------
plotSpatialData(
    spatial_data = cbind(spPatterns, counts_df), 
                positions_path = tissue_paths,
                image_path = image_paths,
                jsonPath = scalefactor_paths,
                feature_col = "Pattern_1",
                colors = NULL,
                title = "Pattern_1",
                scaled = FALSE, plot = FALSE)

## ----message=FALSE, warning=FALSE, dpi=60, fig.width=6------------------------
plotSpatialData(
    spatial_data = cbind(spPatterns, counts_df), 
                positions_path = tissue_paths,
                image_path = image_paths,
                jsonPath = scalefactor_paths,
                feature_col = "Pattern_5",
                colors = NULL,
                title = "Pattern_5",
                scaled = FALSE, plot = FALSE)

## ----message=FALSE, warning=FALSE, dpi=36, fig.width=6, fig.height=2.5--------
plot_grid(plotlist = list(exprPlots[[1]],spatialMaps[[1]]))
plot_grid(plotlist = list(exprPlots[[2]],spatialMaps[[2]]))

## ----message=FALSE, dpi=36, warning=FALSE, fig.width=6, fig.height=2.5--------
plot_grid(plotlist = list(exprPlots[[3]],spatialMaps[[3]]))

## -----------------------------------------------------------------------------
bottom <- res_enrich %>% arrange(SpaceMarkersMetric)
print(head(bottom))

## ----warning=FALSE------------------------------------------------------------
g <- bottom$Gene[1]
p1 <- plotSpatialExpr(
    data = counts_df,gene = g,hotspots = hotspots, 
                   patterns = c("Pattern_1","Pattern_5"))
p2 <- plotSpatialData(
    spatial_data = cbind(spPatterns, counts_df), 
                positions_path = tissue_paths,
                image_path = image_paths,
                jsonPath = scalefactor_paths,
                feature_col = g,
                colors = NULL,
                title = g,
                scaled = FALSE, plot = FALSE)

plot_grid(plotlist = list(p1,p2))


## -----------------------------------------------------------------------------
unlink(file.path(data_dir), recursive = TRUE)

## ----echo=FALSE---------------------------------------------------------------
parameters = c('visiumDir', 'h5filename')
paramDescript = c('A string path to the h5 file with expression information',
                    'A string of the name of the h5 file in the directory')
paramTable = data.frame(parameters, paramDescript)
knitr::kable(paramTable, col.names = c("Argument","Description"))


## ----echo=FALSE---------------------------------------------------------------
parameters = c('visiumDir', 'resolution')
paramDescript = c(
'A path to the location of the the spatial coordinates folder.',
'String values to look for in the .json object;lowres or highres.')
paramTable = data.frame(parameters, paramDescript)
knitr::kable(paramTable, col.names = c("Argument","Description"))


## ----echo=FALSE---------------------------------------------------------------
parameters = c('spPatterns','visiumDir','spatialDir','pattern','sigma',
               'threshold','resolution')
paramDescript = c('A data frame of spatial coordinates and patterns.',
                  'A directory with the spatial and expression data for 
                  the tissue sample',
                  'A directory with spatial data for the tissue sample',
                  'A string of the .json filename with the image parameters',
                  'A numeric value specifying the kernel distribution width',
                  'A numeric value specifying the outlier threshold for the 
                  kernel',
                  'A string specifying the image resolution to scale')
paramTable = data.frame(parameters, paramDescript)
knitr::kable(paramTable, col.names = c("Argument","Description"))


## ----echo=FALSE---------------------------------------------------------------
parameters = c(
        'data','reconstruction', 'optParams','spPatterns',
        'mode', 'minOverlap','hotspotRegions','analysis')
paramDescript = c(
        'An expression matrix of genes and columns being the samples.',
        'Latent feature matrix. NULL if \'DE\' mode is specified',
        'A matrix of sigmaOpts (width) and the thresOpt (outlierthreshold)',
        'A data frame that contains of spatial coordinates and patterns.',
        'A string of the reference pattern for comparison to other patterns',
        'A string specifying either \'residual\' or \'DE\' mode.',
        'A value that specifies the minimum pattern overlap. 50 is the default',
        'A string specifying the type of analysis')
paramTable = data.frame(parameters, paramDescript)
knitr::kable(paramTable, col.names = c("Argument","Description"))


## -----------------------------------------------------------------------------
sessionInfo()