Cepo

# install.packages("remotes")
remotes::install_github("Zaoqu-Liu/Cepo")

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Cepo")

library(Cepo)
library(SingleCellExperiment)

# Load example data
data("cellbench", package = "Cepo")

# Run Cepo analysis
result <- Cepo(
    exprsMat = logcounts(cellbench),
    cellTypes = cellbench$celltype
)

# View top identity genes for each cell type
topGenes(result, n = 5)

# Standard analysis with default parameters
result <- Cepo(
    exprsMat = logcounts(sce),
    cellTypes = sce$celltype
)

# Filter lowly expressed genes (expressed in >5% cells)
result <- Cepo(
    exprsMat = logcounts(sce),
    cellTypes = sce$celltype,
    exprsPct = 0.05
)

# Filter by log fold change
result <- Cepo(
    exprsMat = logcounts(sce),
    cellTypes = sce$celltype,
    logfc = 1
)

# Fast p-value computation using normal approximation
result <- Cepo(
    exprsMat = logcounts(sce),
    cellTypes = sce$celltype,
    computePvalue = 200,        # Number of bootstrap iterations
    computeFastPvalue = TRUE    # Use fast approximation
)

# Access p-values
head(result$pvalues)

# Coefficient of Variation (default)
result_cv <- Cepo(exprsMat, cellTypes, variability = "CV")

# Standard Deviation
result_sd <- Cepo(exprsMat, cellTypes, variability = "SD")

# Median Absolute Deviation
result_mad <- Cepo(exprsMat, cellTypes, variability = "MAD")

# Interquartile Range
result_iqr <- Cepo(exprsMat, cellTypes, variability = "IQR")

# Run analysis accounting for batch effects
result <- Cepo(
    exprsMat = logcounts(sce),
    cellTypes = sce$celltype,
    block = sce$batch,
    minCelltype = 2
)

# Access batch-specific and averaged results
names(result)  # Individual batches + "average"

# Use multiple cores for faster computation
result <- Cepo(
    exprsMat = logcounts(sce),
    cellTypes = sce$celltype,
    computePvalue = 200,
    workers = 4  # Number of cores
)

# Plot gene expression densities
plotDensities(
    x = sce,
    cepoOutput = result,
    nGenes = 2,
    assay = "logcounts",
    celltypeColumn = "celltype"
)

# Plot specific genes
plotDensities(
    x = sce,
    cepoOutput = result,
    genes = c("GENE1", "GENE2", "GENE3"),
    assay = "logcounts",
    celltypeColumn = "celltype"
)

library(DelayedArray)
library(HDF5Array)

# Convert to HDF5-backed DelayedArray
da_matrix <- DelayedArray(realize(logcounts(sce), "HDF5Array"))

# Run Cepo (automatically uses R backend for DelayedArray)
result <- Cepo(
    exprsMat = da_matrix,
    cellTypes = sce$celltype
)