## ----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)

## -----------------------------------------------------------------------------
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 = "/")

## -----------------------------------------------------------------------------
unlink(basename(sp_url))
files <- list.files(".")[grepl(counts_file,list.files("."))]
unlink(files)
unlink("spatial", recursive = TRUE)

## -----------------------------------------------------------------------------
download.file(counts_url,basename(counts_url), mode = "wb")
counts_matrix <- load10XExpr(visiumDir = ".", 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,]
files <- list.files(".")[grepl(basename(counts_url),list.files("."))]
unlink(files)

## -----------------------------------------------------------------------------
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, basename(sp_url), mode = "wb")
untar(basename(sp_url))
spCoords <- load10XCoords(visiumDir = ".")
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, ]

## -----------------------------------------------------------------------------
data("optParams")
optParams

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


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

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

## -----------------------------------------------------------------------------
residual_p1_p5<-SpaceMarkers$interacting_genes[[1]]
DE_p1_p5<-SpaceMarkers_DE$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 <- getInteractingGenes(data = counts_matrix,
                                    reconstruction = cogaps_matrix,
                                    optParams = optParams,
                                    spPatterns = spPatterns,
                                    refPattern = "Pattern_1",
                                    mode ="residual",analysis="enrichment")
SpaceMarkers_DE_enrich <- getInteractingGenes(
    data=counts_matrix,reconstruction=NULL,
    optParams = optParams,
    spPatterns = spPatterns,
    refPattern = "Pattern_1",
    mode="DE",analysis="enrichment")
residual_p1_p5_enrichment<-SpaceMarkers_enrich$interacting_genes[[1]]$Gene
DE_p1_p5_enrichment<-SpaceMarkers_DE_enrich$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$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)

## -----------------------------------------------------------------------------
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, ...)
    )
}

## -----------------------------------------------------------------------------
sample_names <- c("BreastCancer")
image_paths <- c("spatial/tissue_lowres_image.png")
scalefactor_paths <- c("spatial/scalefactors_json.json")
tissue_paths <- c("spatial/tissue_positions_list.csv")

images_cl <- list()

for (i in 1:length(sample_names)) {
    images_cl[[i]] <- read.bitmap(image_paths[i])
}

height <- list()

for (i in 1:length(sample_names)) {
    height[[i]] <-  data.frame(height = nrow(images_cl[[i]]))
}

height <- bind_rows(height)

width <- list()

for (i in 1:length(sample_names)) {
    width[[i]] <- data.frame(width = ncol(images_cl[[i]]))
}

width <- bind_rows(width)

grobs <- list()
for (i in 1:length(sample_names)) {
    grobs[[i]]<-rasterGrob(images_cl[[i]],width=unit(1,"npc"),
                            height=unit(1,"npc"))
}

images_tibble <- tibble(sample=factor(sample_names), grob=grobs)
images_tibble$height <- height$height
images_tibble$width <- width$width


scales <- list()

for (i in 1:length(sample_names)) {
    scales[[i]] <- rjson::fromJSON(file = scalefactor_paths[i])
}


## -----------------------------------------------------------------------------
bcs <- list()
for (i in 1:length(sample_names)) {
    bcs[[i]] <- read.csv(tissue_paths[i],col.names=c(
        "barcode","tissue","row","col","imagerow","imagecol"), header = FALSE)
    bcs[[i]]$imagerow <- bcs[[i]]$imagerow * scales[[i]]$tissue_lowres_scalef
    # scale tissue coordinates for lowres image
    bcs[[i]]$imagecol <- bcs[[i]]$imagecol * scales[[i]]$tissue_lowres_scalef
    bcs[[i]]$tissue <- as.factor(bcs[[i]]$tissue)
    bcs[[i]]$height <- height$height[i]
    bcs[[i]]$width <- width$width[i]
}
names(bcs) <- sample_names

## -----------------------------------------------------------------------------
matrix <- list()
for (i in 1:length(sample_names)) {
    matrix[[i]] <- as.data.frame(t(as.matrix(counts_matrix)))
}
umi_sum <- list()
for (i in 1:length(sample_names)) {
    umi_sum[[i]] <- data.frame(barcode =  row.names(matrix[[i]]),
                                sum_umi = Matrix::rowSums(matrix[[i]]))
}
names(umi_sum) <- sample_names

umi_sum <- bind_rows(umi_sum, .id = "sample")
gene_sum <- list()

for (i in 1:length(sample_names)) {
    gene_sum[[i]] <- data.frame(barcode=row.names(
        matrix[[i]]),sum_gene=Matrix::rowSums(matrix[[i]] != 0))

}
names(gene_sum) <- sample_names
gene_sum <- bind_rows(gene_sum, .id = "sample")
bcs_merge <- bind_rows(bcs, .id = "sample")
bcs_merge <- merge(bcs_merge,umi_sum, by = c("barcode", "sample"))
bcs_merge <- merge(bcs_merge,gene_sum, by = c("barcode", "sample"))


## -----------------------------------------------------------------------------
myPalette <- function(numLevels) {
    return(colorRampPalette(c("blue","yellow"))(numLevels))}

## -----------------------------------------------------------------------------
gene_list <- c()
sp_genes <- SpaceMarkers$interacting_genes[[1]]
interactions <- unique(sp_genes$`Pattern_1 x Pattern_5`)
n_genes <- 3
for (g in 1:length(interactions)){
    
    df <- sp_genes %>% dplyr::filter(
        sp_genes$`Pattern_1 x Pattern_5` == interactions[g] & abs(
        sp_genes$Dunn.zP2_P1) > 1 )
    df <- df[!is.na(df$Gene),]
    valid_genes <- min(nrow(df),n_genes)
    print(paste0("Top ",valid_genes," genes for ",interactions[g]))
    print(df$Gene[1:valid_genes])
    gene_list <- c(gene_list,df$Gene[1:valid_genes])
    
}


## ----message = FALSE, warning=FALSE-------------------------------------------

plots <- list()
# default size = 1.75, stroke = 0.5
for (g in gene_list){
    for (i in 1:length(sample_names)) {
        plots[[length(plots)+1]] <- bcs_merge %>%dplyr::filter(
            sample ==sample_names[i]) %>% bind_cols(as.data.table(
            matrix[i])[,g, with=FALSE]) %>% ggplot(aes_string(
                x='imagecol', y='imagerow', fill=g, alpha = g)) +geom_spatial(
                data=images_tibble[i,], aes(grob=grob), x=0.5, y=0.5)+
            geom_point(shape = 21, colour = "black", size = 1.1, stroke = 0.2)+
            coord_cartesian(expand=FALSE)+scale_fill_gradientn(
                colours = myPalette(100))+xlim(0,max(bcs_merge %>%dplyr::filter(
                    sample ==sample_names[i]) %>% select(width)))+ylim(max(
                        bcs_merge %>%dplyr::filter(sample ==sample_names[i])%>%
                            select(height)),0)+xlab("") +ylab("") + ggtitle(
                                sample_names[i])+
            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())
    }
}

## -----------------------------------------------------------------------------
region <- SpaceMarkers$hotspots[,1]
region <- ifelse(!is.na(region)&!is.na(SpaceMarkers$hotspots[,2]),
                    "Interacting",ifelse(!is.na(region),region,
                                            SpaceMarkers$hotspots[,2]))
region <- factor(region, levels = c("Pattern_1","Interacting","Pattern_5"))
plist <- list()
mplot2 <- t(as.matrix(counts_matrix[,!is.na(region)]))
mplot2 <- as.data.frame(as.matrix(mplot2))
mplot2 <- cbind(mplot2,region = region[!is.na(region)])
for (ii in 1:length(gene_list)){
    plist[[ii]]<- mplot2 %>% ggplot( aes_string(x='region',y=gene_list[ii],
                                                fill='region'))+geom_boxplot()+
        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_list[ii]," Expression (Log)")) + xlab("") 
}

## -----------------------------------------------------------------------------
head(residual_p1_p5 %>% dplyr::filter(
    sp_genes$`Pattern_1 x Pattern_5` == "vsBoth"),n_genes)

## ----message=FALSE, warning=FALSE, dpi=36, fig.width=12, fig.height=5---------
plot_grid(plotlist = list(plist[[1]],plots[[1]]))
plot_grid(plotlist = list(plist[[2]],plots[[2]]))

## ----message=FALSE, dpi=36, warning=FALSE, fig.width=12, fig.height=5---------
plot_grid(plotlist = list(plist[[3]],plots[[3]]))

## -----------------------------------------------------------------------------
head(sp_genes %>% dplyr::filter(
    sp_genes$`Pattern_1 x Pattern_5` == "vsPattern_1"),n_genes - 1 )

## ----message=FALSE, dpi=36, warning=FALSE, fig.width=12, fig.height=5---------
plot_grid(plotlist = list(plist[[4]],plots[[4]]))
plot_grid(plotlist = list(plist[[5]],plots[[5]]))

## -----------------------------------------------------------------------------
head(sp_genes %>% dplyr::filter(
    sp_genes$`Pattern_1 x Pattern_5`=="vsPattern_5"),n_genes - 1)


## ----message=FALSE, dpi=36, warning=FALSE, fig.width=12, fig.height=5---------
plot_grid(plotlist = list(plist[[6]],plots[[6]]))
plot_grid(plotlist = list(plist[[7]],plots[[7]]))

## -----------------------------------------------------------------------------
unlink(basename(sp_url))
unlink("spatial", recursive = TRUE)
files <- list.files(".")[grepl(basename(counts_url),list.files("."))]
unlink(files)

## ----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')
paramDescript = c('A data frame of spatial coordinates and patterns.')
paramTable = data.frame(parameters, paramDescript)
knitr::kable(paramTable, col.names = c("Argument","Description"))


## ----echo=FALSE---------------------------------------------------------------
parameters = c(
        'data','reconstruction', 'optParams','spPatterns',
            'refPattern','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 vector of patterns to compare to the \'refPattern\'',
        'A string specifying the type of analysis')
paramTable = data.frame(parameters, paramDescript)
knitr::kable(paramTable, col.names = c("Argument","Description"))


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