## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----setup, results='hide', message=FALSE, warning=FALSE----------------------
library(blase)
library(SingleCellExperiment)
library(tradeSeq)
library(slingshot)
library(scran)
library(scater)
library(BiocParallel)
library(ami)
library(patchwork)
library(reshape2)

## ----randomseed---------------------------------------------------------------
RNGversion("3.5.0")
SEED <- 7
set.seed(SEED)

## ----concurrency--------------------------------------------------------------
N_CORES <- 4
if (ami::using_ci()) {
    N_CORES <- 2
}

## ----expressionOverPseudotime-------------------------------------------------
# Adapted from
# https://codepal.ai/code-generator/query/n4dEA6I9/plot-sine-wave-ggplot2-r
x <- seq(0, 2 * pi, length.out = 500)
gene1_expr <- 10 * (sin(x + 0.1) + 1)
gene2_expr <- 11 * (sin(x - pi) + 1)
gene3_expr <- 7 * (sin(x + pi / 2) + 1)

genes <- data.frame(gene1 = gene1_expr, gene2 = gene2_expr, gene3 = gene3_expr)
genes_melt <- melt(t(genes), varnames = c("gene", "x"))
genes_melt$x <- as.numeric(genes_melt$x)

## ----expressionOverPseudotimePlot---------------------------------------------
ggplot(genes_melt, aes = aes()) +
    geom_line(aes(x = x, y = value, color = gene)) +
    geom_vline(xintercept = 0, linetype = "dashed") +
    geom_vline(xintercept = 24, linetype = "dashed") +
    geom_vline(xintercept = 181, linetype = "dashed") +
    geom_vline(xintercept = 243, linetype = "dashed") +
    geom_vline(xintercept = 315, linetype = "dashed") +
    geom_vline(xintercept = 479, linetype = "dashed") +
    geom_vline(xintercept = 500, linetype = "dashed") +
    labs(
        x = "pseudotime",
        y = "Gene Expression",
        title = "Variable Genes Over Pseudotime"
    )

## ----installBlase, eval=FALSE-------------------------------------------------
# if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager")
# 
# BiocManager::install("blase")

## ----loadSCE------------------------------------------------------------------
data(tradeSeq_BLASE_example_sce, package = "blase")
sce <- tradeSeq_BLASE_example_sce

## ----setupSCEPlot-------------------------------------------------------------
plotUMAP(sce, text_by = "celltype", colour_by = "celltype") +
    plotUMAP(sce, colour_by = "pseudotime") +
    patchwork::plot_layout(ncol = 1, axis_title = "collect")

## ----TradeSeqGeneSelection----------------------------------------------------
# Use consecutive for genes that change over time
assoRes <- associationTest(
    sce,
    lineages = TRUE,
    global = FALSE,
    contrastType = "consecutive"
)
genelist <- blase::get_top_n_genes(assoRes, n_genes = 200, lineage = 1)

## ----findBestParams-----------------------------------------------------------
res <- blase::find_best_params(
    sce,
    genelist,
    split_by = "pseudotime_range",
    pseudotime_slot = "pseudotime",
    bins_count_range = c(5, 10),
    gene_count_range = c(40, 80)
)

## ----findBestParamsPlot-------------------------------------------------------
blase::plot_find_best_params_results(res)

## ----checkBinsChoice----------------------------------------------------------
blase_data <- as.BlaseData(sce, pseudotime_slot = "pseudotime", n_bins = 10)
genes(blase_data) <- genelist[1:80]

## ----checkBinsChoicePlot------------------------------------------------------
evaluate_parameters(blase_data, make_plot = TRUE)

## ----evaluateTopGenes---------------------------------------------------------
evaluate_top_n_genes(blase_data)

## ----doBLASEMapping-----------------------------------------------------------
bulks_df <- DataFrame(row.names = rownames(counts(sce)))
for (type in unique(sce$celltype)) {
    bulks_df <- cbind(
        bulks_df,
        rowSums(normcounts(subset(sce, , celltype == type)))
    )
}

colnames(bulks_df) <- unique(sce$celltype)

multipotent_progenitors_result <- map_best_bin(
    blase_data, "Multipotent progenitors", bulks_df
)
multipotent_progenitors_result
basophils_result <- map_best_bin(blase_data, "Basophils", bulks_df)
basophils_result
neutrophils_result <- map_best_bin(blase_data, "Neutrophils", bulks_df)
neutrophils_result

## ----fastMapping--------------------------------------------------------------
map_all_best_bins(blase_data, bulks_df)

## ----plotMappingResults-------------------------------------------------------
plot_mapping_result_heatmap(list(
    multipotent_progenitors_result,
    basophils_result,
    neutrophils_result
), annotate_correlation = TRUE, text_background = TRUE)

## ----plotMappingResultCorr----------------------------------------------------
plot_mapping_result_corr(basophils_result)

## ----assignPseudotimeBins-----------------------------------------------------
binned_sce <- assign_pseudotime_bins(
    sce,
    split_by = "pseudotime_range",
    n_bins = 10,
    pseudotime_slot = "pseudotime"
)

## ----assignPseudotimeBinsPlot-------------------------------------------------
plot_mapping_result(
    binned_sce,
    multipotent_progenitors_result,
    group_by_slot = "celltype"
)

## ----plotBinPopulation--------------------------------------------------------
plot_bin_population(binned_sce, 1, group_by_slot = "celltype")

## ----sessionInfo--------------------------------------------------------------
sessionInfo()