## ----setup, include = FALSE---------------------------------------------------
knitr::knit_hooks$set(pngquant = knitr::hook_pngquant)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "ragg_png",
  dpi = 72,
  fig.retina = 2,
  fig.align = "center",
  out.width = "100%",
  pngquant = "--speed=1 --quality=1-5"
)

## ----message=FALSE------------------------------------------------------------
library(scDiagnostics)
library(SingleCellExperiment)
library(SpatialExperiment)

set.seed(100)

## ----message=FALSE------------------------------------------------------------
data("merfish_reference_data")
data("merfish_query_data")

class(merfish_reference_data)
table(merfish_reference_data$cell_type_merged)
table(merfish_query_data$tier2_merged[merfish_query_data$cell_type_merged == "Fibroblast"])

## ----fig.height=6, fig.width=10-----------------------------------------------
plotCellTypePCA(
    query_data = merfish_query_data,
    reference_data = merfish_reference_data,
    cell_types = unique(merfish_reference_data$cell_type_merged),
    query_cell_type_col = "cell_type_merged",
    ref_cell_type_col = "cell_type_merged",
    pc_subset = 1:3)

## ----fig.height=5, fig.width=10-----------------------------------------------
anomaly_output <- detectAnomaly(
    reference_data = merfish_reference_data,
    query_data = merfish_query_data,
    ref_cell_type_col = "cell_type_merged",
    query_cell_type_col = "cell_type_merged",
    cell_types = "Fibroblast",
    pc_subset = 1:5,
    n_tree = 500)

is_anomalous <- anomaly_output[["Fibroblast"]]$query_anomaly
mean(is_anomalous)

## ----fig.height=5, fig.width=10-----------------------------------------------
plot(anomaly_output, cell_type = "Fibroblast", pc_subset = 1:3, data_type = "query")

## ----message=FALSE------------------------------------------------------------
fibro_query <- merfish_query_data[, merfish_query_data$cell_type_merged == "Fibroblast"]
tapply(is_anomalous, fibro_query$tier2_merged, mean)

## ----message=FALSE------------------------------------------------------------
ecm_signature <- c("Col1a2", "Timp2", "Col6a1", "Sparc", "Dpt")

gene_shifts <- calculateGeneShifts(
    query_data = merfish_query_data,
    reference_data = merfish_reference_data,
    query_cell_type_col = "cell_type_merged",
    ref_cell_type_col = "cell_type_merged",
    cell_types = "Fibroblast",
    pc_subset = 1:5,
    genes_to_analyze = ecm_signature,
    detect_anomalies = TRUE,
    anomaly_comparison = TRUE)

gene_shifts$PC1[order(gene_shifts$PC1$p_adjusted), ]

## ----fig.height=6, fig.width=8------------------------------------------------
plot(gene_shifts, cell_type = "Fibroblast", pc_subset = 1:5,
    plot_type = "heatmap", plot_by = "p_adjusted", n_genes = 5,
    show_anomalies = TRUE)

## ----fig.height=5, fig.width=8------------------------------------------------
gene_shifts_focused <- calculateGeneShifts(
    query_data = merfish_query_data,
    reference_data = merfish_reference_data,
    query_cell_type_col = "cell_type_merged",
    ref_cell_type_col = "cell_type_merged",
    cell_types = "Fibroblast",
    pc_subset = 1:5,
    genes_to_analyze = c("Col1a2", "Sparc"),
    detect_anomalies = TRUE,
    anomaly_comparison = TRUE)

plot(gene_shifts_focused, cell_type = "Fibroblast", pc_subset = 1:5,
    plot_type = "barplot", plot_by = "p_adjusted", n_genes = 2,
    show_anomalies = TRUE, pseudo_bulk = TRUE)

## ----SessionInfo, echo=FALSE, message=FALSE, warning=FALSE, comment=NA--------
options(width = 80)
sessionInfo()