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

set.seed(100)

## ----message=FALSE------------------------------------------------------------
data("covid_reference_data")
data("covid_query_data")

table(covid_reference_data$author_cell_type_merged)
table(covid_query_data$azimuth_celltype_l1_merged)

## ----fig.height=6, fig.width=10-----------------------------------------------
shared_cell_types <- c("CD14 mono", "CD4 T", "CD8 T", "B cell", "NK_16hi")

plotCellTypePCA(
    query_data = covid_query_data,
    reference_data = covid_reference_data,
    cell_types = shared_cell_types,
    query_cell_type_col = "azimuth_celltype_l1_merged",
    ref_cell_type_col = "author_cell_type_merged",
    pc_subset = 1:3)

## ----fig.height=5, fig.width=10-----------------------------------------------
anomaly_output <- detectAnomaly(
    reference_data = covid_reference_data,
    query_data = covid_query_data,
    ref_cell_type_col = "author_cell_type_merged",
    query_cell_type_col = "azimuth_celltype_l1_merged",
    cell_types = "CD14 mono",
    pc_subset = 1:5,
    n_tree = 500)

mean(anomaly_output[["CD14 mono"]]$query_anomaly)

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

## ----message=FALSE------------------------------------------------------------
yoshida_ifn_signature <- c(
    "BST2", "CMPK2", "EIF2AK2", "EPSTI1", "HERC5", "IFI35", "IFI44L",
    "IFI6", "IFIT3", "ISG15", "LY6E", "MX1", "MX2", "OAS1", "OAS2",
    "PARP9", "PLSCR1", "SAMD9", "SAMD9L", "SP110", "STAT1", "TRIM22",
    "UBE2L6", "XAF1", "IRF7")

gene_shifts <- calculateGeneShifts(
    query_data = covid_query_data[yoshida_ifn_signature, ],
    reference_data = covid_reference_data[yoshida_ifn_signature, ],
    query_cell_type_col = "azimuth_celltype_l1_merged",
    ref_cell_type_col = "author_cell_type_merged",
    cell_types = "CD14 mono",
    pc_subset = 1:5,
    n_top_loadings = 25,
    detect_anomalies = TRUE,
    anomaly_comparison = TRUE)

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

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

## ----fig.height=6, fig.width=8------------------------------------------------
plot(gene_shifts, cell_type = "CD14 mono", pc_subset = 1:5,
    plot_type = "barplot", plot_by = "p_adjusted", n_genes = 15,
    show_anomalies = TRUE, pseudo_bulk = TRUE)

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