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

## ----setup, results='hide', message=FALSE, warning=FALSE----------------------
library(scater)
library(ggplot2)
library(BiocParallel)
library(blase)
library(patchwork)

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

## ----concurrency--------------------------------------------------------------
N_CORES <- 2
bpparam <- MulticoreParam(N_CORES)

## ----load_data----------------------------------------------------------------
data(painter_microarray, package = "blase")
data(MCA_PF_SCE, package = "blase")

## ----sc_umaps-----------------------------------------------------------------
plotUMAP(MCA_PF_SCE, colour_by = "STAGE_HR")
#| fig.alt: >
#|   UMAP of Dogga et al. single-cell RNA-seq reference coloured by pseudotime,
#|   starting with Rings, and ending with Schizonts.
plotUMAP(MCA_PF_SCE, colour_by = "slingPseudotime_1")

## ----create_BLASE_object------------------------------------------------------
blase_data <- as.BlaseData(
    MCA_PF_SCE,
    pseudotime_slot = "slingPseudotime_1",
    n_bins = 48,
    split_by = "cells"
)

# Add these bins to the sc metadata
MCA_PF_SCE <- assign_pseudotime_bins(MCA_PF_SCE,
    pseudotime_slot = "slingPseudotime_1",
    n_bins = 48,
    split_by = "cells"
)

## ----calculate_gene_peakedness------------------------------------------------
gene_peakedness_info <- calculate_gene_peakedness(
    MCA_PF_SCE,
    window_pct = 5,
    knots = 18,
    BPPARAM = bpparam
)

genes_to_use <- gene_peakedness_spread_selection(
    MCA_PF_SCE,
    gene_peakedness_info,
    genes_per_bin = 30,
    n_gene_bins = 30
)

head(gene_peakedness_info)

## ----plot_gene_peakedness-----------------------------------------------------
ggplot(gene_peakedness_info, aes(x = peak_pseudotime, y = ratio)) +
    geom_point() +
    ggtitle("All genes")

gene_peakedness_selected_genes <- gene_peakedness_info[
    gene_peakedness_info$gene %in% genes_to_use,
]

## ----plot_gene_peakedness_selection-------------------------------------------
ggplot(gene_peakedness_selected_genes, aes(x = peak_pseudotime, y = ratio)) +
    geom_point() +
    ggtitle("Selected genes")

## ----add_genes_to_blase_data--------------------------------------------------
genes(blase_data) <- genes_to_use

## ----mapping------------------------------------------------------------------
mapping_results <- map_all_best_bins(
    blase_data,
    painter_microarray,
    BPPARAM = bpparam
)

## ----mappingPlot--------------------------------------------------------------
plot_mapping_result_heatmap(mapping_results)

## ----plot_pseudotime_bins-----------------------------------------------------
plotUMAP(MCA_PF_SCE, colour_by = "pseudotime_bin")

## ----transfer_mappings_to_sc--------------------------------------------------
annotated_sce <- annotate_sce(MCA_PF_SCE, mapping_results, include_stats = TRUE)
annotated_sce$BLASE_Annotation <- gsub(
  " hpi", "", annotated_sce$BLASE_Annotation)
annotated_sce$BLASE_Annotation <- as.numeric(annotated_sce$BLASE_Annotation)

annotation_plot <- plotUMAP(annotated_sce, colour_by = "BLASE_Annotation")
corr_plot <- plotUMAP(
  annotated_sce, colour_by = "BLASE_Annotation_Correlation") +
    labs(color = "Correlation")

## ----transfer_mappings_to_scPlot----------------------------------------------
annotation_plot + corr_plot

## ----transfer_painter_stages_to_sc--------------------------------------------
annotated_sce$BLASE_Annotation_transfer <- ""
for (best_bulk in unique(annotated_sce$BLASE_Annotation)) {
    mask <- colData(annotated_sce)[, "BLASE_Annotation"] == best_bulk

    if (best_bulk %in% 0:15) {
        annotated_sce[, mask]$BLASE_Annotation_transfer <- "Ring"
    } else if (best_bulk %in% 16:21) {
        annotated_sce[, mask]$BLASE_Annotation_transfer <- "Ring/Trophozoite"
    } else if (best_bulk %in% 22:32) {
        annotated_sce[, mask]$BLASE_Annotation_transfer <- "Trophozoite"
    } else if (best_bulk %in% 33:48) {
        annotated_sce[, mask]$BLASE_Annotation_transfer <- "Schizont"
    }
}

## ----transfer_painter_stages_to_scPlot----------------------------------------
plotUMAP(annotated_sce, colour_by = "BLASE_Annotation_transfer") +
    plotUMAP(annotated_sce, colour_by = "STAGE_LR")

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