## ----eval=FALSE---------------------------------------------------------------
# install.packages("BiocManager")
# BiocManager::install("spoon")
## -----------------------------------------------------------------------------
library(nnSVG)
library(STexampleData)
library(SpatialExperiment)
library(BRISC)
library(BiocParallel)
library(scuttle)
library(Matrix)
library(spoon)
spe <- Visium_mouseCoronal()
## -----------------------------------------------------------------------------
# keep spots over tissue
spe <- spe[, colData(spe)$in_tissue == 1]
# filter out low quality genes
spe <- filter_genes(spe)
# calculate logcounts (log-transformed normalized counts) using scran package
spe <- computeLibraryFactors(spe)
spe <- logNormCounts(spe)
# choose a small number of genes for this example to run quickly
set.seed(3)
ix_random <- sample(seq_len(nrow(spe)), 10)
spe <- spe[ix_random, ]
# remove spots with zero counts
spe <- spe[, colSums(logcounts(spe)) > 0]
## -----------------------------------------------------------------------------
weights <- generate_weights(input = spe,
stabilize = TRUE,
BPPARAM = MulticoreParam(workers = 1,
RNGseed = 4))
## -----------------------------------------------------------------------------
spe <- weighted_nnSVG(input = spe,
w = weights,
BPPARAM = MulticoreParam(workers = 1, RNGseed = 5))
## -----------------------------------------------------------------------------
# display results
rowData(spe)
## ----eval=FALSE---------------------------------------------------------------
# assay_name <- "logcounts"
# weighted_logcounts <- t(weights)*assays(spe)[[assay_name]]
# assay(spe, "weighted_logcounts") <- weighted_logcounts
## ----eval=FALSE---------------------------------------------------------------
# library(SpatialExperiment)
# library(STexampleData)
# library(MASS)
# library(scuttle)
#
# set.seed(1)
#
# #4992 spots and 300 genes
#
# n_genes <- 300
# fraction <- 0.5
# max_sigma.sq <- 1
# low_range_beta <- c(0.5,1)
#
# #check if integer
# stopifnot(n_genes*fraction*0.5 == round(n_genes*fraction*0.5))
#
# #all genes have some nonzero sigma.sq
# sigma.sq <- runif(n_genes, 0.2, max_sigma.sq)
# ground_truth_rank <- rank(-sigma.sq)
#
# #all genes will have nonzero beta values
# beta <- runif(n_genes, log(low_range_beta[1]), log(low_range_beta[2]))
#
# #choose fixed length scale parameter (~medium from nnSVG paper)
#
# scale_length <- 200
#
# params <- data.frame(sigma.sq, beta)
#
# plot(beta, sigma.sq)
#
# #sampling from a poisson distribution - mean controls variance, so we don't specify tau.sq:
# #step 1: use ST example distance matrix instead of creating a new one (Euclidean distance)
#
# spe_demo <- Visium_humanDLPFC()
# points_coord <- spatialCoords(spe_demo)
# n_points <- nrow(points_coord)
#
# pair.points <- cbind(
# matrix( rep(points_coord, each = n_points), ncol = 2, byrow = FALSE),
# rep(1, times = n_points) %x% points_coord # Creating the combinations using kronecker product.
# ) |> data.frame()
# colnames(pair.points) <- c("si.x", "si.y", "sj.x", "sj.y")
#
# #step 2: calculate gaussian process/kernel
#
# kernel.fun <- function(si.x, si.y, sj.x, sj.y, l = 0.2){
# exp(-1*sqrt(((si.x-sj.x)^2+(si.y-sj.y)^2))/l)
# }
#
# C_theta <- with(pair.points, kernel.fun(si.x, si.y, sj.x, sj.y, l = scale_length)) |>
# matrix(nrow = n_points, ncol = n_points)
#
# counts <- matrix(NA, nrow = n_genes, ncol = n_points)
# eta_list <- list()
#
# for (i in c(1:n_genes)) {
#
# sigma.sq_i <- sigma.sq[i]
# beta_i <- beta[i]
#
# #step 3: simulate gaussian process per gene
#
# gp_dat <- mvrnorm(n = 1, rep(0,n_points), sigma.sq_i* C_theta)
#
# #step 4: calculate lambda = exp(beta + gaussian process) per gene
#
# lambda_i <- exp(gp_dat + beta_i)
#
# #step 5: use rpois() to simulate 4992 values per gene
#
# counts_i <- rpois(n = n_points, lambda_i)
#
# #put all counts in matrix
# #orientation: genes x spots
#
# counts[i,] <- counts_i
# }
#
# #create spe using counts and distance matrix
#
# spe <- SpatialExperiment(
# assays = list(counts = counts),
# spatialCoords = points_coord)
#
# rowData(spe)$ground_truth <- ground_truth
# rowData(spe)$ground_truth_rank <- ground_truth_rank
# rowData(spe)$ground_truth_sigma.sq <- sigma.sq
# rowData(spe)$ground_truth_beta <- beta
#
## ----eval=FALSE---------------------------------------------------------------
# library(SpatialExperiment)
# library(nnSVG)
# library(BRISC)
# library(BiocParallel)
# library(scuttle)
#
# spe <- spe[, colSums(counts(spe)) > 0]
# spe <- logNormCounts(spe)
# weights <- generate_weights(input = spe, stabilize = TRUE)
# spe_unweighted <- nnSVG(spe, assay_name = "logcounts")
# spe_weighted <- weighted_nnSVG(input = spe, w = weights)
## ----eval=FALSE---------------------------------------------------------------
# library(ggplot2)
# library(SpatialExperiment)
# library(patchwork)
# library(GGally)
# library(dplyr)
# library(ggridges)
#
# #overlay unweighted and weighted ridge plots
# df_unw <- data.frame(
# rank = rowData(spe_unweighted)$rank,
# mean = rowData(spe_unweighted)$mean,
# method = rep("unw", 300)
# ) %>% mutate(quantile = findInterval(mean,
# quantile(mean, probs=0:9/10))) %>%
# tibble::rownames_to_column()
#
# df_w <- data.frame(
# rank = rowData(spe_weighted)$weighted_rank,
# mean = rowData(spe_weighted)$weighted_mean,
# method = rep("w", 300)
# ) %>% mutate(quantile = findInterval(mean,
# quantile(mean, probs=0:9/10))) %>%
# tibble::rownames_to_column()
#
# df <- rbind(df_unw, df_w) %>%
# mutate(quantile = as.factor(quantile))
#
# ridge_overlay <- ggplot(df, aes(x = rank, y = quantile)) +
# geom_density_ridges2(aes(fill = method), rel_min_height = 0.02, alpha = 0.3, scale = 1.3) +
# theme_ridges(grid = TRUE) +
# labs(
# y = "decile - unw & w mean of logcounts",
# x = "rank",
# title = "Ridge plots: effect of weighting on rank"
# ) +
# scale_fill_manual(labels = c("weighted", "unweighted"), values = c("blue", "red")) +
# coord_cartesian(xlim = c(1, 300)) +
# theme_bw()
#
# #ridge plots separated by noise and signal for unweighted and weighted
# frac <- round(dim(spe_unweighted)[1]*0.1)*0.1
#
# df_unw_signal <- df_unw %>%
# mutate(quantile = as.factor(quantile)) %>%
# group_by(quantile) %>%
# slice_min(order_by = rank, n = frac) %>%
# mutate(grp = "signal")
#
# indices <- as.integer(df_unw_signal$rowname)
#
# df_unw_background <- df_unw[-indices,] %>%
# mutate(quantile = as.factor(quantile)) %>%
# mutate(grp = "background")
#
# df <- rbind(df_unw_signal, df_unw_background)
#
# rank_separated_unw <- ggplot(df, aes(x = rank, y = quantile)) +
# geom_density_ridges2(aes(fill = grp), rel_min_height = 0.02, alpha = 0.3) +
# theme_ridges(grid = TRUE) +
# labs(
# y = "decile - unw mean of logcounts",
# x = "rank",
# title = "Signal unweighted"
# ) +
# guides(fill=guide_legend(title="group")) +
# coord_cartesian(xlim = c(1, 300)) +
# theme_bw()
#
# df_w_signal <- df_w %>%
# mutate(quantile = as.factor(quantile)) %>%
# group_by(quantile) %>%
# slice_min(order_by = rank, n = frac) %>%
# mutate(grp = "signal")
#
# indices <- as.integer(df_w_signal$rowname)
#
# df_w_background <- df_w[-indices,] %>%
# mutate(quantile = as.factor(quantile)) %>%
# mutate(grp = "background")
#
# df <- rbind(df_w_signal, df_w_background)
#
# rank_separated_w <- ggplot(df, aes(x = rank, y = quantile)) +
# geom_density_ridges2(aes(fill = grp), rel_min_height = 0.02, alpha = 0.3) +
# theme_ridges(grid = TRUE) +
# labs(
# y = "decile - w mean of logcounts",
# x = "rank",
# title = "Signal weighted"
# ) +
# guides(fill=guide_legend(title="group")) +
# coord_cartesian(xlim = c(1, 300)) +
# theme_bw()
#
# ridge_signal <- wrap_plots(rank_separated_unw, rank_separated_w, nrow=1, guides = "collect")
## -----------------------------------------------------------------------------
sessionInfo()