## ----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(limma)
library(ggVennDiagram)

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

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

## ----load_data----------------------------------------------------------------
data(zhang_2021_heat_shock_bulk, 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 = 8,
    split_by = "pseudotime_range"
)

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

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

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

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

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

## ----setUpDE------------------------------------------------------------------
metadata <- data.frame(row.names = seq_len(12))
metadata$strain <- c(
    rep("NF54", 4),
    rep("PB4", 4),
    rep("PB31", 4)
)
metadata$growth_conditions <- rep(c("Normal", "Normal", "HS", "HS"), 3)
metadata$group <- paste0(metadata$strain, "_", metadata$growth_conditions)
rownames(metadata) <- c(
    "NF54_37_1",
    "NF54_37_2",
    "NF54_41_1",
    "NF54_41_2",
    "PB4_37_1",
    "PB4_37_2",
    "PB4_41_1",
    "PB4_41_2",
    "PB31_37_1",
    "PB31_37_2",
    "PB31_41_1",
    "PB31_41_2"
)

design <- model.matrix(~ 0 + group, metadata)
colnames(design) <- gsub("group", "", colnames(design))

contr.matrix <- makeContrasts(
    WT_normal_v_HS = NF54_Normal - NF54_HS,
    PB31_normal_v_HS = PB31_Normal - PB31_HS,
    normal_NF54_v_PB31 = NF54_Normal - PB31_Normal,
    levels = colnames(design)
)

## ----calculateBulkDE----------------------------------------------------------
vfit <- lmFit(zhang_2021_heat_shock_bulk, design)
vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
efit <- eBayes(vfit)

summary(decideTests(efit))

PB31_normal_v_HS_BULK_DE <- topTable(efit, n = Inf, coef = 2)
PB31_normal_v_HS_BULK_DE <- PB31_normal_v_HS_BULK_DE[PB31_normal_v_HS_BULK_DE$adj.P.Val < 0.05 & PB31_normal_v_HS_BULK_DE$logFC > 0, ]
PB31_normal_v_HS_BULK_DE <- PB31_normal_v_HS_BULK_DE[order(-PB31_normal_v_HS_BULK_DE$logFC), ]

## ----pseudobulk---------------------------------------------------------------
pseudobulked_MCA_PF_SCE <- data.frame(row.names = rownames(MCA_PF_SCE))
for (r in mapping_results) {
    pseudobulked_MCA_PF_SCE[bulk_name(r)] <- rowSums(counts(MCA_PF_SCE[, MCA_PF_SCE$pseudotime_bin == best_bin(r)]))
}
pseudobulked_MCA_PF_SCE <- as.matrix(pseudobulked_MCA_PF_SCE)

## -----------------------------------------------------------------------------
## Normalise, not needed for true bulks
par(mfrow = c(1, 2))
v <- voom(pseudobulked_MCA_PF_SCE, design, plot = FALSE)

vfit_sc <- lmFit(v, design)
vfit_sc <- contrasts.fit(vfit_sc, contrasts = contr.matrix)
efit_sc <- eBayes(vfit_sc)

summary(decideTests(efit_sc))

## -----------------------------------------------------------------------------
PB31_normal_v_HS_SC_DE <- topTable(efit_sc, n = Inf, coef = 2)
PB31_normal_v_HS_SC_DE <- PB31_normal_v_HS_SC_DE[PB31_normal_v_HS_SC_DE$adj.P.Val < 0.05 & PB31_normal_v_HS_SC_DE$logFC > 0, ]
PB31_normal_v_HS_SC_DE <- PB31_normal_v_HS_SC_DE[order(-PB31_normal_v_HS_SC_DE$logFC), ]

## -----------------------------------------------------------------------------
ggVennDiagram(list(
    Phenotype = rownames(PB31_normal_v_HS_BULK_DE),
    Development = rownames(PB31_normal_v_HS_SC_DE)
)) + ggtitle("PB31 Normal vs Heat Shock")

## -----------------------------------------------------------------------------
PB31_normal_v_HS_corrected_DE <- rownames(PB31_normal_v_HS_BULK_DE[!(rownames(PB31_normal_v_HS_BULK_DE) %in% rownames(PB31_normal_v_HS_SC_DE)), ])
head(PB31_normal_v_HS_corrected_DE)

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