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

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


## -----------------------------------------------------------------------------
#| 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 <- compute_contrasts(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 `glioblastoma_example_data` prior to 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(colour = 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
#| message: false
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. In contrast to @fig-raw_umap, the Harmony alignment procedure has been applied with the aim of making corresponding cells from the different patient samples and conditions lie next to each other."
#| 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(colour = patient_id), size = 0.5) +
    facet_wrap(vars(condition)) + coord_fixed()


## -----------------------------------------------------------------------------
#| label: lemur_compute_contrasts
fit <- compute_contrasts(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(colour = de)) +
    scale_colour_gradient2(low = "#FFD800", high= "#0056B9") + coord_fixed()

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


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


## -----------------------------------------------------------------------------
#| label: topgenes
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 _SERPINE1_ superimposed on the same UMAP plot as in @fig-lemur_umap."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
sel_gene <- "ENSG00000106366" # SERPINE1

dplyr::filter(neighborhoods, name == sel_gene) |>
  dplyr::select(name, adj_pval, lfc, did_adj_pval, did_lfc)

p <- tibble(umap = umap) |>
  mutate(de = assay(fit, "DE")[sel_gene,]) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(colour = de)) +
    scale_colour_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 neighbourhood 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", colour = "#808080") 


## -----------------------------------------------------------------------------
#| label: fig-volcano_plot
#| fig-cap: "Volcano plot, each point corresponds to one of the genes in our dataset, and the neighbourhood that was found for it."
#| fig-width: 4
#| fig-height: 2.5
#| out-width: "50%"
neighborhoods |>
  drop_na() |>
  ggplot(aes(x = lfc, y = -log10(pval))) + geom_point(aes(col = adj_pval < 0.1)) + 
    scale_colour_manual(values= c(`TRUE` = "orangered", `FALSE` = "thistle3"))


## -----------------------------------------------------------------------------
#| label: fig-neighbourhood_size_vs_significance
#| fig-cap: "neighbourhood size vs neighbourhood significance."
#| fig-width: 6
#| fig-height: 4
#| out-width: "50%"
#| eval: false
#| echo: false
# neighborhoods |>
#   drop_na() |>
#   ggplot(aes(x = n_cells, y = -log10(pval))) +
#     geom_point(aes(colour = adj_pval < 0.1))


## -----------------------------------------------------------------------------
#| label: fig-tumour_cell_annotation1
#| fig-cap: "A simple gating strategy to putatively annotate the tumour cells in our data"
#| fig-width: 6
#| fig-height: 4
#| out-width: "60%"
thresh = c(chr7 = 0.8, chr10 = 2.5)
tumour_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_tumour = chr7_total_expr > thresh["chr7"] & chr10_total_expr < thresh["chr10"])

ggplot(tumour_label_df, aes(x = chr10_total_expr, y = chr7_total_expr, colour = is_tumour)) +
  geom_point(size = 0.5) +
  geom_hline(yintercept = thresh["chr7"]) +
  geom_vline(xintercept = thresh["chr10"]) +
  scale_colour_manual(values= c(`TRUE` = "orangered", `FALSE` = "#404040"))

## -----------------------------------------------------------------------------
#| label: fig-tumour_cell_annotation2
#| fig-cap: "The tumour cells appear enriched in the larger of the blobs."
#| fig-width: 6
#| fig-height: 4
#| out-width: "60%"
tibble(umap = umap) |>
  mutate(is_tumour = tumour_label_df$is_tumour) |>
  ggplot(aes(x = umap[,1], y = umap[,2], colour = is_tumour)) +
    geom_point(size = 0.5) +
    facet_wrap(vars(is_tumour)) + coord_fixed() + 
    scale_colour_manual(values= c(`TRUE` = "orangered", `FALSE` = "#404040"))


## -----------------------------------------------------------------------------
#| label: tumour_de_neighbourhood
#| paged.print: false
tumour_fit <- fit[, tumour_label_df$is_tumour]
tum_nei <- find_de_neighborhoods(tumour_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-tumour_de_neighbourhood_plot
#| paged.print: false
#| fig-cap: "CXCL8 expression in different cell subsets"
#| fig-width: 7
#| fig-height: 4
#| out-width: "80%"
sel_gene <- "ENSG00000169429" # CXCL8

as_tibble(colData(fit)) |>
  mutate(expr = assay(fit, "logcounts")[sel_gene, ]) |>
  mutate(is_tumour = tumour_label_df$is_tumour) |>
  mutate(in_neighborhood = id %in% filter(tum_nei, name == sel_gene)$neighborhood[[1]]) |>
  ggplot(aes(x = condition, y = ifelse(expr < 10, expr, 10))) +
    geom_jitter(size = pi/4, stroke = 0) +
    geom_point(data = . %>% summarize(expr = mean(expr), .by = c(condition, patient_id, is_tumour, in_neighborhood)),
               aes(colour = patient_id), size = 2) +
    stat_summary(fun.data = mean_se, geom = "crossbar", colour = "red") +
    facet_wrap(vars(is_tumour, 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()