## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse  = TRUE,
    comment   = "#>",
    fig.width = 7,
    fig.height = 5,
    message   = FALSE,
    warning   = FALSE
)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# BiocManager::install("scFastDE")

## ----quickstart---------------------------------------------------------------
library(scFastDE)
library(SingleCellExperiment)

# Simulate a multi-donor, two-condition dataset (unpaired)
set.seed(2024)
n_genes <- 500
n_cells <- 120

counts <- matrix(rpois(n_genes * n_cells, 8L),
                 nrow = n_genes, ncol = n_cells)
rownames(counts) <- paste0("Gene", seq_len(n_genes))
colnames(counts) <- paste0("Cell", seq_len(n_cells))

# Inject DE signal: first 20 genes upregulated in treatment
counts[1:20, 61:120] <- counts[1:20, 61:120] * 4L

sce <- SingleCellExperiment(assays = list(counts = counts))
sce$donor     <- rep(paste0("Donor", 1:12), each = 10)
sce$cell_type <- "CD4_Tcell"
sce$condition <- rep(c("ctrl", "treat"), each = 60)

sce

## ----filter-------------------------------------------------------------------
sce <- filterSparseDonors(sce,
                           donor     = "donor",
                           cell_type = "cell_type",
                           min_cells = 5)
ncol(sce)

## ----pseudobulk---------------------------------------------------------------
pb <- fastPseudobulk(sce,
                      donor       = "donor",
                      cell_type   = "cell_type",
                      target_type = "CD4_Tcell")

cat("Pseudo-bulk matrix:", nrow(pb$pseudobulk),
    "genes x", ncol(pb$pseudobulk), "donors\n")

cat("\nDonor weights (sqrt of cell count):\n")
print(round(pb$donor_weights, 2))

## ----fastde-------------------------------------------------------------------
result <- fastDE(sce,
                  donor       = "donor",
                  cell_type   = "cell_type",
                  condition   = "condition",
                  target_type = "CD4_Tcell",
                  min_cells   = 5,
                  min_cpm     = 1,
                  min_donors  = 2)

result

## ----de_results---------------------------------------------------------------
# Top significant genes
dt <- as.data.frame(deTable(result))
dt_sorted <- dt[order(dt$adj.P.Val), ]
head(dt_sorted[, c("gene", "logFC", "P.Value", "adj.P.Val")], 15)

## ----de_summary---------------------------------------------------------------
# Summary
sig <- sum(dt$adj.P.Val < 0.05, na.rm = TRUE)
cat(sprintf("\n%d / %d genes significant (FDR < 0.05)\n",
            sig, nrow(dt)))
cat(sprintf("Injected genes recovered: %d / 20\n",
            sum(paste0("Gene", 1:20) %in%
                rownames(dt)[dt$adj.P.Val < 0.05])))

## ----volcano, fig.cap = "Volcano plot of DE results. Red = upregulated, blue = downregulated at FDR < 0.05 and |logFC| > 1."----
plotDEResults(result, fdr_thresh = 0.05, lfc_thresh = 1, top_n = 10)

## ----weighting_demo-----------------------------------------------------------
# Create an unbalanced dataset
set.seed(99)
sce_unbal <- sce
# Simulate donor 1 having only 2 cells by keeping just 2
keep <- c(which(sce$donor != "Donor1"),
          which(sce$donor == "Donor1")[1:2])
sce_unbal <- sce_unbal[, keep]

pb_unbal <- fastPseudobulk(sce_unbal,
                             donor       = "donor",
                             cell_type   = "cell_type",
                             target_type = "CD4_Tcell")

cat("Donor weights (unbalanced):\n")
print(round(pb_unbal$donor_weights, 2))

## ----paired_setup-------------------------------------------------------------
# Simulate paired data: SAME 6 donors in BOTH conditions
set.seed(2025)
n_genes <- 500
n_donors <- 6
cells_per_donor_cond <- 15
n_cells <- n_donors * 2 * cells_per_donor_cond  # 180 cells

counts <- matrix(rpois(n_genes * n_cells, 8L),
                 nrow = n_genes, ncol = n_cells)
rownames(counts) <- paste0("Gene", seq_len(n_genes))
colnames(counts) <- paste0("Cell", seq_len(n_cells))

# Inject DE signal: first 25 genes upregulated in stimulated condition
stim_cols <- seq(n_donors * cells_per_donor_cond + 1, n_cells)
counts[1:25, stim_cols] <- counts[1:25, stim_cols] * 4L

sce_paired <- SingleCellExperiment(assays = list(counts = counts))

# KEY: same donors (D1–D6) appear in BOTH conditions
sce_paired$donor <- rep(
    rep(paste0("D", seq_len(n_donors)), each = cells_per_donor_cond),
    2
)
sce_paired$cell_type <- "Monocyte"
sce_paired$condition <- rep(c("ctrl", "stim"),
                             each = n_donors * cells_per_donor_cond)

# Verify: each donor appears in both conditions
cat("Cells per donor and condition:\n")
print(table(sce_paired$donor, sce_paired$condition))

## ----paired_pseudobulk--------------------------------------------------------
pb_paired <- fastPseudobulk(sce_paired,
                              donor       = "donor",
                              cell_type   = "cell_type",
                              target_type = "Monocyte",
                              condition   = "condition")

cat(sprintf("Pseudo-bulk: %d genes × %d samples\n",
            nrow(pb_paired$pseudobulk),
            ncol(pb_paired$pseudobulk)))
cat("  (6 donors × 2 conditions = 12 samples)\n\n")

cat("Sample info:\n")
print(pb_paired$sample_info)

## ----paired_de----------------------------------------------------------------
result_paired <- fastDE(sce_paired,
                         donor       = "donor",
                         cell_type   = "cell_type",
                         condition   = "condition",
                         target_type = "Monocyte",
                         min_cells   = 5,
                         min_cpm     = 1,
                         min_donors  = 2)

result_paired

## ----paired_results-----------------------------------------------------------
dt_paired <- as.data.frame(deTable(result_paired))
dt_paired_sorted <- dt_paired[order(dt_paired$adj.P.Val), ]

sig_paired <- sum(dt_paired$adj.P.Val < 0.05, na.rm = TRUE)
recovered  <- sum(paste0("Gene", 1:25) %in%
                  rownames(dt_paired)[dt_paired$adj.P.Val < 0.05])

cat(sprintf("Significant genes (FDR < 0.05): %d / %d\n",
            sig_paired, nrow(dt_paired)))
cat(sprintf("Injected DE genes recovered:    %d / 25\n", recovered))

## ----paired_volcano, fig.cap = "Volcano plot for paired design. The paired model correctly blocks on donor, recovering the injected DE signal."----
plotDEResults(result_paired, fdr_thresh = 0.05, lfc_thresh = 1,
             top_n = 10)

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