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