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

## ----install bioc, eval=FALSE-------------------------------------------------
# if (!requireNamespace("BiocManager")) {
#     install.packages("BiocManager")
# }
# BiocManager::install("BatchSVG")

## ----install github, eval = FALSE---------------------------------------------
# remotes::install("christinehou11/BatchSVG")

## ----library, message=FALSE---------------------------------------------------
library(BatchSVG)
library(WeberDivechaLCdata)
library(cowplot)
library(ggplot2)

## ----load data, comment=NA, warning=FALSE, message=FALSE----------------------
spe <- WeberDivechaLCdata_Visium()
spe

## ----load nnsvg, comment=NA, warning=FALSE, message=FALSE---------------------
svgs <- read.csv("svgs_LC.csv", check.names = FALSE)

## ----feature select, comment = NA, warning=FALSE------------------------------
list_batch_df <- featureSelect(input = spe, 
    batch_effect = "sample_id", VGs = svgs$gene_id)

## ----eval=FALSE---------------------------------------------------------------
# list_batch_df <- featureSelect(input = spe,
#     batch_effect = "sample_id", VGs = svgs$gene_id, verbose = FALSE)

## ----feature select class, comment = NA, warning=FALSE------------------------
class(list_batch_df)

## ----feature select print, comment = NA, warning=FALSE------------------------
head(list_batch_df$sample_id)

## ----svg, comment=NA, warning=FALSE, message=FALSE----------------------------
plots <- svg_nSD(list_batch_df = list_batch_df, 
            sd_interval_dev = 5, sd_interval_rank = 5)

## ----figure 1, warning=FALSE, message=FALSE, fig.width=10, fig.height=8-------
plots$sample_id

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

vec_dev <- list_batch_df$sample_id$nSD_dev_sample_id[list_batch_df$sample_id$nSD_dev_sample_id > 0]
x_values <- 2:round(max(vec_dev))

counts <- sapply(x_values, function(x) sum(vec_dev > x))

df <- data.frame(
  threshold = x_values,
  count = counts
)

p1 <- ggplot(df, aes(x = threshold, y = count)) +
  geom_point(color = "blue", size = 2) +
  labs(x = "Threshold Value",
       y = "Number of SVGs > Threshold",
       title = "Relative Change in Deviance") +
  theme_minimal()

vec_rank <- list_batch_df$sample_id$nSD_rank_sample_id[list_batch_df$sample_id$nSD_rank_sample_id > 0]
x_values <- 2:round(max(vec_rank))

counts <- sapply(x_values, function(x) sum(vec_rank > x))

df <- data.frame(
  threshold = x_values,
  count = counts
)

p2 <- ggplot(df, aes(x = threshold, y = count)) +
  geom_point(color = "blue", size = 2) +
  labs(x = "Threshold Value",
       y = "Number of SVGs > Threshold",
       title = "Rank Difference") +
  theme_minimal()

plot_grid(p1,p2)

## ----bias detect dev, comment = NA, message=FALSE, warning=FALSE--------------
bias_dev <- biasDetect(list_batch_df = list_batch_df, 
    threshold = "dev", nSD_dev = 10)

## ----table 1, comment = NA, message=FALSE, warning=FALSE----------------------
head(bias_dev$sample_id$Table)

## ----size change, message=FALSE, warning=FALSE--------------------------------
# size default = 3
bias_dev_size <- biasDetect(list_batch_df = list_batch_df, 
    threshold = "dev", nSD_dev = 10, plot_point_size = 4)

## ----figure 2, warning=FALSE, message=FALSE, fig.width= 10, fig.height=4------
plot_grid(bias_dev$sample_id$Plot, bias_dev_size$sample_id$Plot)

## ----bias detect rank, comment = NA, message=FALSE, warning=FALSE-------------
bias_rank <- biasDetect(list_batch_df = list_batch_df, 
    threshold = "rank", nSD_rank = 8)

## ----table 2, comment = NA, message=FALSE, warning=FALSE----------------------
head(bias_rank$sample_id$Table)

## ----figure 3, message=FALSE, warning=FALSE, fig.width= 10, fig.height=4------
# shape default = 16
bias_rank_shape <- biasDetect(list_batch_df = list_batch_df, 
    threshold = "rank", nSD_rank = 8, plot_point_shape = 2)

plot_grid(bias_rank$sample_id$Plot, bias_rank_shape$sample_id$Plot)

## ----both, comment = NA, message=FALSE, warning=FALSE-------------------------
bias_both <- biasDetect(list_batch_df = list_batch_df, threshold = "both",
    nSD_dev = 10, nSD_rank = 8)

## ----table 3, comment = NA, message=FALSE, warning=FALSE----------------------
head(bias_both$sample_id$Table)

## ----figure 4, message=FALSE, warning=FALSE, fig.width= 10, fig.height=8------
# color default = "YlOrRd"
bias_both_color <- biasDetect(list_batch_df = list_batch_df, 
    threshold = "both", nSD_dev = 10, nSD_rank = 8, plot_palette = "Greens")

plot_grid(bias_both$sample_id$Plot, bias_both_color$sample_id$Plot,nrow = 2)

## ----figure 5, message=FALSE, warning=FALSE, fig.width= 10, fig.height=8------
# text size default = 3
bias_both_color_text <- biasDetect(list_batch_df = list_batch_df, 
    threshold = "both", nSD_dev = 10, nSD_rank = 8, 
    plot_palette = c("Blues"), plot_text_size = 4)

plot_grid(bias_both$sample_id$Plot, bias_both_color_text$sample_id$Plot,nrow = 2)

## ----new svgs, comment = NA, message=FALSE, warning=FALSE---------------------
bias_both_df <- bias_both$sample_id$Table
svgs_filt <- setdiff(svgs$gene_id, bias_both_df$gene_id)
svgs_filt_spe <- svgs[svgs$gene_id %in% svgs_filt, ]
nrow(svgs_filt_spe)

## ----session info-------------------------------------------------------------
## Session info
sessionInfo()

