## -----------------------------------------------------------------------------
#| label: setup_alt_IS_THIS_REALLY_NEEDED
#| echo: false
#| eval: false
# if(! exists("..options_set") || isFALSE(..options_set)){
#   knitr::opts_chunk$set(
#     collapse = TRUE,
#     comment = "#>",
#     dpi = 40,
#     fig_retina = 1,
#     dev = "jpeg"
#     # dev.args = list(quality = 20)
#   )
#   ..options_set <- TRUE
# }

## ----initialize---------------------------------------------------------------
#| label: init
#| echo: false
#| cache: false
knitr::opts_chunk$set(cache = TRUE, autodep = TRUE)
options(width = 100)


## -----------------------------------------------------------------------------
#| label: install1
#| eval: false
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# 
# BiocManager::install("lemur")


## -----------------------------------------------------------------------------
#| label: install2
#| eval: false
# devtools::install_github("const-ae/lemur")


## -----------------------------------------------------------------------------
#| label: preparation
#| echo: false
library("lemur")
set.seed(42)
data("glioblastoma_example_data")
sce <- glioblastoma_example_data[1:50, sample.int(5000, size = 500)]


## -----------------------------------------------------------------------------
#| label: quick_start
#| message: false
#| warning: false
# ... sce is a SingleCellExperiment object with your data 
fit <- lemur(sce, design = ~ patient_id + condition, n_embedding = 15)
fit <- align_harmony(fit)   # This step is optional
fit <- test_de(fit, contrast = cond(condition = "ctrl") - cond(condition = "panobinostat"))
nei <- find_de_neighborhoods(fit, group_by = vars(patient_id, condition))


## -----------------------------------------------------------------------------
#| label: load_packages
#| message: false
#| warning: false
library("tidyverse")
library("SingleCellExperiment")
library("lemur")
set.seed(42)


## -----------------------------------------------------------------------------
#| label: load_data
data("glioblastoma_example_data", package = "lemur")
glioblastoma_example_data


## -----------------------------------------------------------------------------
#| label: fig-raw_umap
#| fig-cap: "UMAP of the example dataset `glioblastoma_example_data`, prior to any between-condition alignment."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
orig_umap <- uwot::umap(as.matrix(t(logcounts(glioblastoma_example_data))))

as_tibble(colData(glioblastoma_example_data)) |>
  mutate(umap = orig_umap) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = patient_id, shape = condition), size = 0.5) +
    labs(title = "UMAP of logcounts") + coord_fixed()


## -----------------------------------------------------------------------------
#| label: fit_lemur
#| message: false
fit <- lemur(glioblastoma_example_data, design = ~ patient_id + condition, 
             n_embedding = 15, test_fraction = 0.5)
fit


## -----------------------------------------------------------------------------
#| label: fitembedding
fit$embedding |> str()


## -----------------------------------------------------------------------------
#| label: align_lemur
fit <- align_harmony(fit)


## -----------------------------------------------------------------------------
#| label: fig-lemur_umap
#| fig-cap: "UMAPs of `fit$embedding`. The points are shown separately for the two conditions, but reside in the same latent space."
#| fig-width: 7
#| fig-height: 3
#| out-width: "100%"
umap <- uwot::umap(t(fit$embedding))

as_tibble(fit$colData) |>
  mutate(umap = umap) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = patient_id), size = 0.5) +
    facet_wrap(vars(condition)) + coord_fixed()


## -----------------------------------------------------------------------------
#| label: lemur_test_de
fit <- test_de(fit, contrast = cond(condition = "panobinostat") - cond(condition = "ctrl"))


## -----------------------------------------------------------------------------
#| label: fig-umap_de
#| fig-subcap: 
#|   - "Superimposed on the same UMAP plot as in @fig-lemur_umap"
#|   - "The histogram shows that the values are negative for most cells."
#| fig-cap: "Differential expression (log fold changes) of GAP43"
#| fig-width: 6
#| fig-height: 4
#| out-width: "100%"
#| layout: "[65, 35]"
df <- tibble(umap = umap) |>
  mutate(de = assay(fit, "DE")["ENSG00000172020", ])
 
ggplot(df, aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = de)) +
    scale_color_gradient2(low = "#FFD800", high= "#0056B9") + coord_fixed()

ggplot(df, aes(x = de)) + geom_histogram(bins = 100)


## -----------------------------------------------------------------------------
#| label: fig-de_neighborhoods
#| message: false 
neighborhoods <- find_de_neighborhoods(fit, group_by = vars(patient_id, condition))

as_tibble(neighborhoods) |>
  left_join(as_tibble(rowData(fit)[,1:2]), by = c("name" = "gene_id")) |>
  relocate(symbol, .before = "name") |>
  arrange(pval) |>
  head(5)


## -----------------------------------------------------------------------------
#| label: fig-umap_de2
#| fig-cap: "Differential expression of CXCL8 superimposed on the same UMAP plot as in @fig-lemur_umap."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
sel_gene <- "ENSG00000169429" # is CXCL8

p <- tibble(umap = umap) |>
  mutate(de = assay(fit, "DE")[sel_gene,]) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = de)) +
    scale_color_gradient2(low = "#FFD800", high= "#0056B9") +
    coord_fixed()
p


## -----------------------------------------------------------------------------
#| label: fig-umap_de3
#| fig-cap: "Same as @fig-umap_de2, with an attempt to draw a neighborhood boundary."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
neighborhood_coordinates <- neighborhoods |>
  dplyr::filter(name == sel_gene) |>
  unnest(c(neighborhood)) |>
  dplyr::rename(cell_id = neighborhood) |>
  left_join(tibble(cell_id = rownames(umap), umap), by = "cell_id") |>
  dplyr::select(name, cell_id, umap)

p + geom_density2d(data = neighborhood_coordinates, breaks = seq(0.2, 0.6, by = 0.1), 
                   contour_var = "ndensity", color = "#808080") 


## -----------------------------------------------------------------------------
#| label: fig-volcano_plot
#| fig-cap: "Volcano plot, each point corresponds to one neighborhood."
#| fig-width: 6
#| fig-height: 4
#| out-width: "50%"
neighborhoods |>
  drop_na() |>
  ggplot(aes(x = lfc, y = -log10(pval))) +
    geom_point(aes(col  = adj_pval < 0.1)) 


## -----------------------------------------------------------------------------
#| label: fig-Neighborhood_size_vs_significance
#| fig-cap: "Neighborhood size vs neighborhood significance."
#| fig-width: 6
#| fig-height: 4
#| out-width: "50%"
neighborhoods |>
  drop_na() |>
  ggplot(aes(x = n_cells, y = -log10(pval))) +
    geom_point(aes(color  = adj_pval < 0.1)) 


## -----------------------------------------------------------------------------
#| label: fig-tumor_cell_annotation1
#| fig-cap: "A simple gating strategy to find tumor cells"
#| fig-width: 6
#| fig-height: 4
#| out-width: "60%"
tumor_label_df <- tibble(cell_id = colnames(fit),
       chr7_total_expr  = colMeans(logcounts(fit)[rowData(fit)$chromosome == "7",]),
       chr10_total_expr = colMeans(logcounts(fit)[rowData(fit)$chromosome == "10",])) |>
  mutate(is_tumor = chr7_total_expr > 0.8 & chr10_total_expr < 2.5)

ggplot(tumor_label_df, aes(x = chr10_total_expr, y = chr7_total_expr)) +
    geom_point(aes(color = is_tumor), size = 0.5) +
    geom_hline(yintercept = 0.8) +
    geom_vline(xintercept = 2.5) 

## -----------------------------------------------------------------------------
#| label: fig-tumor_cell_annotation2
#| fig-cap: "The tumor cells are enriched in parts of the big (left) blob."
#| fig-width: 6
#| fig-height: 4
#| out-width: "60%"
tibble(umap = umap) |>
  mutate(is_tumor = tumor_label_df$is_tumor) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(color = is_tumor), size = 0.5) +
    facet_wrap(vars(is_tumor)) + coord_fixed()


## -----------------------------------------------------------------------------
#| label: tumor_de_neighborhood
#| paged.print: false
tumor_fit <- fit[, tumor_label_df$is_tumor]
tum_nei <- find_de_neighborhoods(tumor_fit, group_by = vars(patient_id, condition), verbose = FALSE)

as_tibble(tum_nei) |>
  left_join(as_tibble(rowData(fit)[,1:2]), by = c("name" = "gene_id")) |>
  dplyr::relocate(symbol, .before = "name") |>
  filter(adj_pval < 0.1) |>
  arrange(did_pval)  |>
  dplyr::select(symbol, name, neighborhood, n_cells, adj_pval, lfc, did_pval, did_lfc) |>
  print(n = 10)


## -----------------------------------------------------------------------------
#| label: fig-tumor_de_neighborhood_plot
#| paged.print: false
#| fig-cap: ""
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
sel_gene <- "ENSG00000142534" # is RPS11

as_tibble(colData(fit)) |>
  mutate(expr = assay(fit, "logcounts")[sel_gene,]) |>
  mutate(is_tumor = tumor_label_df$is_tumor) |>
  mutate(in_neighborhood = id %in% filter(tum_nei, name == sel_gene)$neighborhood[[1]]) |>
  ggplot(aes(x = condition, y = expr)) +
    geom_jitter(size = 0.3, stroke = 0) +
    geom_point(data = . %>% summarize(expr = mean(expr), .by = c(condition, patient_id, is_tumor, in_neighborhood)),
               aes(color = patient_id), size = 2) +
    stat_summary(fun.data = mean_se, geom = "crossbar", color = "red") +
    facet_wrap(vars(is_tumor, in_neighborhood), labeller = label_both) 



## -----------------------------------------------------------------------------
#| label: fix_linear_coef
#| message: false
#| warning: false
fit <- lemur(sce, design = ~ patient_id + condition, n_embedding = 15, linear_coefficient_estimator = "zero")


## -----------------------------------------------------------------------------
#| label: session_info
sessionInfo()

