airpart demo

airpart.R

load_all("~/bioc/airpart/airpart")
suppressPackageStartupMessages(library(SingleCellExperiment))
suppressPackageStartupMessages(library(SummarizedExperiment))
p.vec <- c(rep(c(0.5, 0.8), each = 2), 0.65)
set.seed(1)
sce <- makeSimulatedData(
  mu1 = 2, mu2 = 10, nct = 5, n = 20,
  ngenecl = 20, theta = 20, ncl = 1,
  p.vec = p.vec
)
sce <- sce[,c(1:10,31:40,41:50,71:80,81:90)]
sce <- preprocess(sce)
makeRatioHeatmap(sce)

cols <- rep(RColorBrewer::brewer.pal(4, "Paired"), each=10)
boxplot(counts(sce), col=cols, xaxt="n", ylab="total")
boxplot(assay(sce, "ratio"), col=cols, xaxt="n", ylab="ratio")

sce <- geneCluster(sce, G=1:4)
smm <- summaryAllelicRatio(sce)
metadata(smm)$summary
sce <- fusedLasso(sce, model = "binomial", genecluster = 1, ncores = 1, se.rule.mult = 1)
metadata(sce)$partition
sce <- allelicRatio(sce)
makeViolin(sce)
makeRatioHeatmap(sce, order_by_group = FALSE)

### now with gene and donor effects

load_all("~/bioc/airpart/airpart")
suppressPackageStartupMessages(library(SingleCellExperiment))
suppressPackageStartupMessages(library(SummarizedExperiment))

true <- c(1,1,2,2,2,3,3,3,4,4)
log.odds <- c(-.2,-.2,-.1,-.1,-.1,.1,.1,.1,.2,.2) * 5
odds <- exp(c(log.odds, log.odds + 2))
p.vec <- odds/(1 + odds)
plot(p.vec, ylim=c(0,1))

library(pbapply)
res <- pbreplicate(100, {
  sce <- makeSimulatedData(mu1 = 20, mu2 = 20, nct = 10, n = 10, ngenecl = 2, theta = 30, ncl = 2, p.vec = p.vec)
  sce <- sce[1:3,]
  mcols(sce)$cluster <- 1
  sce <- preprocess(sce)
  sce <- fusedLasso(sce, model = "binomial", genecluster = 1, ncores = 1, se.rule.mult = 1)
  part1 <- as.numeric(metadata(sce)$partition[,"part"])
  f <- ratio ~ p(x, pen = "gflasso") + gene
  sce <- fusedLasso(sce, formula = f, model = "binomial", genecluster = 1, ncores = 1, se.rule.mult = 1)
  part2 <- as.numeric(metadata(sce)$partition[,"part"])
  c(mclust::adjustedRandIndex(part1, true),
    mclust::adjustedRandIndex(part2, true))
})
df <- data.frame(adjustedRand=c(res[1,],res[2,]),
                 method=rep(c("no-gene-term","with-gene-term"),
                            each=ncol(res)))
library(ggplot2)
ggplot(df, aes(method, adjustedRand)) +
  geom_violin() + ggforce::geom_sina(alpha=.5)
### now spatial pattern

library(airpart)
suppressPackageStartupMessages(library(SummarizedExperiment))
library(spatialDmelxsim)
se <- spatialDmelxsim()
assays(se) <- assays(se)[1:2]
se <- se[mcols(se)$svASE,]
plot(se$normSlice)
se$x <- factor(floor(se$normSlice/3))
se <- se[,order(se$x)]
rownames(se) <- mcols(se)$paper_symbol
table(se$x)

assay(se, "a1")[is.na(assay(se, "a1"))] <- 0
assay(se, "a2")[is.na(assay(se, "a2"))] <- 0

se <- preprocess(se)
makeHeatmap(se)
se <- geneCluster(se, G=12, method="GMM")
makeHeatmap(se, genecluster=11, show_row_names=TRUE)

nct <- 10
adj_m <- makeOffByOneAdjMat(nct)
rownames(adj_m) <- colnames(adj_m) <- 0:9

library(smurf)
f <- ratio ~ p(x, pen = "ggflasso")
se_sub <- fusedLasso(se,
  formula = f, adj.matrix = list(x = adj_m),
  model = "binomial", lambda.max=100,
  genecluster = 11, ncores = 1, se.rule.mult = 1
)

makeHeatmap(se_sub, order_by_group=FALSE, show_row_names=TRUE)
se_sub <- allelicRatio(se_sub)

library(ggplot2)
makeViolin(se_sub) + ylim(0,1)
makeStep(se_sub, xlab = "grouped slices")