## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## -----------------------------------------------------------------------------
library(airpart)
suppressPackageStartupMessages(library(SingleCellExperiment))
p.vec <- rep(c(0.2, 0.8, 0.5, 0.5, 0.7, 0.9), each = 2)
set.seed(2021)
sce <- makeSimulatedData(
  mu1 = 2, mu2 = 10, nct = 4, n = 20,
  ngenecl = 25, theta = 20, ncl = 3,
  p.vec = p.vec
)

## -----------------------------------------------------------------------------
unique(rowData(sce)) # the true underlying allelic ratios
table(sce$x) # counts of each cell type
assays(sce)[["a1"]][1:5, 1:5] # allelic counts for the effect allele

## -----------------------------------------------------------------------------
assayNames(sce)
sce$x

## -----------------------------------------------------------------------------
sce <- preprocess(sce)
makeHeatmap(sce)

## -----------------------------------------------------------------------------
cellQCmetrics <- cellQC(sce, mad_detected = 4)
cellQCmetrics

## -----------------------------------------------------------------------------
keep_cell <- (
  cellQCmetrics$filter_sum | # sufficient features (genes)
    cellQCmetrics$filter_detected | # sufficient molecules counted
    # sufficient features expressed compared to spike genes,
    # high quality cells
    cellQCmetrics$filter_spike
)
sce <- sce[, keep_cell]

## ---- eval=FALSE--------------------------------------------------------------
#  featureQCmetric <- featureQC(sce)
#  keep_feature <- (featureQCmetric$filter_celltype &
#    featureQCmetric$filter_sd &
#    featureQCmetric$filter_spike)
#  sce <- sce[keep_feature, ]

## -----------------------------------------------------------------------------
sce <- geneCluster(sce, G = 1:4)
metadata(sce)$geneCluster

## -----------------------------------------------------------------------------
sce.hc <- geneCluster(sce, method = "hierarchical")
metadata(sce.hc)$geneCluster

## -----------------------------------------------------------------------------
summary <- summaryAllelicRatio(sce)
summary

## -----------------------------------------------------------------------------
sapply(1:length(summary), function(i) {
  inst <- summary[[i]]
  inst_order <- inst[order(inst$weighted.mean), ]
  max(diff(inst_order$weighted.mean)) > 0.05
})

## -----------------------------------------------------------------------------
estDisp(sce, genecluster = 1)

## -----------------------------------------------------------------------------
sce_sub <- fusedLasso(sce,
  model = "binomial",
  genecluster = 1, ncores = 1,
  niter = 2
)

## ---- results="asis"----------------------------------------------------------
knitr::kable(metadata(sce_sub)$partition, row.names = FALSE)

## -----------------------------------------------------------------------------
metadata(sce_sub)$lambda

## ---- results='asis', collapse=TRUE-------------------------------------------
sce_sub <- consensusPart(sce_sub)
knitr::kable(metadata(sce_sub)$partition, row.names = FALSE)

## -----------------------------------------------------------------------------
thrs <- 10^seq(from = -2, to = -0.4, by = 0.2)
sce_sub_w <- wilcoxExt(sce, genecluster = 1, threshold = thrs)
knitr::kable(metadata(sce_sub_w)$partition, row.names = FALSE)
metadata(sce_sub_w)$threshold

## ---- warning=FALSE, fig.width=12---------------------------------------------
sce_sub <- allelicRatio(sce_sub, DAItest = TRUE)
makeForest(sce_sub, showtext = TRUE)

## ----results="asis"-----------------------------------------------------------
genepoi <- paste0("gene", seq_len(5))
ar <- extractResult(sce_sub)
knitr::kable(ar[genepoi,])
makeStep(sce_sub[genepoi,])

## -----------------------------------------------------------------------------
s <- extractResult(sce_sub, "svalue")
apply(s[genepoi,],2, function(s){s<0.005})

## -----------------------------------------------------------------------------
adj.p <- mcols(sce_sub)$adj.p.value
adj.p < 0.05

## -----------------------------------------------------------------------------
nct <- 8
p.vec <- (rep(c(
  -3, 0, -3, 3,
  rep(0, nct / 2),
  2, 3, 4, 2
), each = 2) + 5) / 10
sce <- makeSimulatedData(
  mu1 = 2, mu2 = 10, nct = nct, n = 30,
  ngenecl = 50, theta = 20, ncl = 3, p.vec = p.vec
)
sce <- preprocess(sce)

cellQCmetrics <- cellQC(sce, mad_detected = 4)
keep_cell <- (
  cellQCmetrics$filter_sum | # sufficient features (genes)
    cellQCmetrics$filter_detected | # sufficient molecules counted
    # sufficient features expressed compared to spike genes,
    # high quality cells
    cellQCmetrics$filter_spike
)
sce <- sce[, keep_cell]

featureQCmetric <- featureQC(sce)
keep_feature <- (featureQCmetric$filter_celltype &
  featureQCmetric$filter_sd &
  featureQCmetric$filter_spike)
sce <- sce[keep_feature, ]

makeHeatmap(sce)

## -----------------------------------------------------------------------------
sce <- geneCluster(sce, G = 1:4)
table(mcols(sce)$cluster)

## -----------------------------------------------------------------------------
estDisp(sce, genecluster = 1)

## -----------------------------------------------------------------------------
sce_sub <- fusedLasso(sce,
  model = "binomial",
  genecluster = 1, ncores = 1
)

## ---- results="asis"----------------------------------------------------------
knitr::kable(metadata(sce_sub)$partition, row.names = FALSE)

## -----------------------------------------------------------------------------
sce_sub2 <- sce_sub[1:10, ]
sce_sub2 <- allelicRatio(sce_sub2)

## -----------------------------------------------------------------------------
makeForest(sce_sub2)
ar <- extractResult(sce_sub2)
knitr::kable(ar)

## -----------------------------------------------------------------------------
makeViolin(sce_sub2)

## -----------------------------------------------------------------------------
makeHeatmap(sce_sub2)

## -----------------------------------------------------------------------------
makeHeatmap(sce_sub2, order_by_group = FALSE)

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