inst/doc/scHOT.R

In scHOT: single-cell higher order testing

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  message = FALSE,
  warning = FALSE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
library(SingleCellExperiment)
library(ggplot2)
library(scHOT)
library(scater)
library(matrixStats)

## -----------------------------------------------------------------------------
data(liver)

liver_pseudotime_hep <- liver$liver_pseudotime_hep
liver_branch_hep <- liver$liver_branch_hep
first_branch_cells <- liver$first_branch_cells

## -----------------------------------------------------------------------------
gene_to_test <- as.matrix(c("Birc5", "H2afz", "Tacc3"))

## -----------------------------------------------------------------------------
scHOT_traj <- scHOT_buildFromMatrix(
  mat = liver_branch_hep[,first_branch_cells],
  cellData = list(pseudotime = liver_pseudotime_hep[first_branch_cells]),
  positionType = "trajectory",
  positionColData = "pseudotime")
scHOT_traj

## -----------------------------------------------------------------------------
scater::plotExpression(scHOT_traj, c("Sall4"),
                       exprs_values = "expression", x = "pseudotime")

## -----------------------------------------------------------------------------
scHOT_traj_wrap = scHOT(scHOT_traj,
                        testingScaffold = gene_to_test,
                        higherOrderFunction = matrixStats::weightedVar,
                        higherOrderFunctionType = "weighted",
                        numberPermutations = 50)

## -----------------------------------------------------------------------------
slot(scHOT_traj_wrap, "scHOT_output")
scHOT_traj_wrap@scHOT_output

slot(scHOT_traj_wrap, "scHOT_output")$higherOrderSequence

plotHigherOrderSequence(scHOT_traj_wrap, gene_to_test)

plotOrderedExpression(scHOT_traj_wrap, gene_to_test) + 
  facet_wrap(~gene, scales = "free_y")

## -----------------------------------------------------------------------------
scHOT_traj@testingScaffold
scHOT_traj <- scHOT_addTestingScaffold(scHOT_traj, gene_to_test)
scHOT_traj@testingScaffold

## -----------------------------------------------------------------------------
scHOT_traj@weightMatrix
scHOT_traj <- scHOT_setWeightMatrix(scHOT_traj,
                                    positionType = "trajectory",
                                    positionColData = c("pseudotime"),
                                    nrow.out = NULL,
                                    span = 0.25)
dim(scHOT_traj@weightMatrix)
class(scHOT_traj@weightMatrix)
plot(scHOT_traj@weightMatrix[50,])

## -----------------------------------------------------------------------------
scHOT_traj <- scHOT_setWeightMatrix(scHOT_traj,
                                    positionType = "trajectory",
                                    positionColData = c("pseudotime"),
                                    nrow.out = 50,
                                    span = 0.25)
dim(scHOT_traj@weightMatrix)
class(scHOT_traj@weightMatrix)
plot(scHOT_traj@weightMatrix[50,])

## -----------------------------------------------------------------------------
scHOT_traj <- scHOT_calculateGlobalHigherOrderFunction(
  scHOT_traj, 
  higherOrderFunction = matrixStats::weightedVar,
  higherOrderFunctionType = "weighted")

slot(scHOT_traj, "scHOT_output")
apply(assay(scHOT_traj, "expression")[c(gene_to_test),],1,var)

## -----------------------------------------------------------------------------
scHOT_traj <- scHOT_setPermutationScaffold(scHOT_traj, 
                                           numberPermutations = c(20,50,20),
                                           storePermutations = c(TRUE, FALSE, TRUE))
slot(scHOT_traj, "scHOT_output")

## -----------------------------------------------------------------------------
scHOT_traj <- scHOT_calculateHigherOrderTestStatistics(
  scHOT_traj,
  higherOrderSummaryFunction = sd)

slot(scHOT_traj, "scHOT_output")

## -----------------------------------------------------------------------------
system.time(scHOT_traj <- scHOT_performPermutationTest(
  scHOT_traj, 
  verbose = TRUE,
  parallel = FALSE))

slot(scHOT_traj, "scHOT_output")

## -----------------------------------------------------------------------------
scHOT_traj <- scHOT_estimatePvalues(scHOT_traj,
                                    nperm_estimate = 10000,
                                    maxDist = 5)
slot(scHOT_traj, "scHOT_output")

## -----------------------------------------------------------------------------
plotHigherOrderSequence(scHOT_traj, gene_to_test)

## -----------------------------------------------------------------------------
data(MOB_subset)
sce_MOB_subset <- MOB_subset$sce_MOB_subset

sce_MOB_subset

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT_buildFromSCE(sce_MOB_subset,
                                    assayName = "logcounts",
                                    positionType = "spatial",
                                    positionColData = c("x", "y"))

scHOT_spatial

## -----------------------------------------------------------------------------
pairs <- t(combn(rownames(sce_MOB_subset),2))
rownames(pairs) <- apply(pairs,1,paste0,collapse = "_")
head(pairs)

set.seed(2020)
pairs <- pairs[sample(nrow(pairs), 20), ]
if (!"Arrb1_Mtor" %in% rownames(pairs)) {
pairs <- rbind(pairs, "Arrb1_Mtor" = c("Arrb1", "Mtor"))
}
if (!"Dnm1l_Fam63b" %in% rownames(pairs)) {
pairs <- rbind(pairs, "Dnm1l_Fam63b" = c("Dnm1l", "Fam63b"))
}

scHOT_spatial <- scHOT_addTestingScaffold(scHOT_spatial, pairs)
scHOT_spatial@testingScaffold

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT_setWeightMatrix(scHOT_spatial,
                                       positionColData = c("x","y"),
                                       positionType = "spatial",
                                       nrow.out = NULL,
                                       span = 0.05)

dim(slot(scHOT_spatial, "weightMatrix"))

## -----------------------------------------------------------------------------
cellID = 75
ggplot(as.data.frame(colData(scHOT_spatial)), aes(x = -x, y = y)) + 
  geom_point(aes(colour = slot(scHOT_spatial, "weightMatrix")[cellID,],
                 size = slot(scHOT_spatial, "weightMatrix")[cellID,])) + 
  scale_colour_gradient(low = "black", high = "purple") + 
  scale_size_continuous(range = c(1,5)) +
  theme_classic() +
  guides(colour = guide_legend(title = "Spatial Weight"),
         size = guide_legend(title = "Spatial Weight")) + 
    ggtitle(paste0("Central cell: ", cellID))

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT_calculateGlobalHigherOrderFunction(
  scHOT_spatial, 
  higherOrderFunction = weightedSpearman,
  higherOrderFunctionType = "weighted")

slot(scHOT_spatial, "scHOT_output")

head(diag(cor(t(assay(scHOT_spatial, "expression")[pairs[,1],]),
         t(assay(scHOT_spatial, "expression")[pairs[,2],]),
         method = "spearman")))

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT_setPermutationScaffold(scHOT_spatial, 
                                              numberPermutations = 50,
                                              numberScaffold = 10)

slot(scHOT_spatial, "scHOT_output")

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT_calculateHigherOrderTestStatistics(scHOT_spatial)

slot(scHOT_spatial, "scHOT_output")

## -----------------------------------------------------------------------------
system.time(scHOT_spatial <- scHOT_performPermutationTest(
  scHOT_spatial, 
  verbose = TRUE,
  parallel = FALSE))

slot(scHOT_spatial, "scHOT_output")

## -----------------------------------------------------------------------------
scHOT_plotPermutationDistributions(scHOT_spatial)

scHOT_spatial <- scHOT_estimatePvalues(scHOT_spatial,
                                       nperm_estimate = 100,
                                       maxDist = 0.1)
slot(scHOT_spatial, "scHOT_output")

ggplot(as.data.frame(slot(scHOT_spatial, "scHOT_output")),
       aes(x = -log10(pvalPermutations), y = -log10(pvalEstimated))) + 
  geom_point() + 
  theme_classic() + 
  geom_abline(slope = 1, intercept = 0) +
  xlab("Permutation -log10(p-value)") +
  ylab("Estimated -log10(p-value)") +
  NULL

## -----------------------------------------------------------------------------
colData(scHOT_spatial)[, "-x"] <- -colData(scHOT_spatial)[, "x"]
plotHigherOrderSequence(scHOT_spatial, c("Dnm1l_Fam63b"),
                        positionColData = c("-x", "y"))

plotOrderedExpression(scHOT_spatial, c("Dnm1l", "Fam63b"),
                      positionColData = c("-x", "y"),
                      assayName = "expression")

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT_stripOutput(scHOT_spatial, force = TRUE)

scHOT_spatial
slot(scHOT_spatial, "scHOT_output")

## -----------------------------------------------------------------------------
scHOT_spatial <- scHOT(scHOT_spatial,
                       testingScaffold = pairs,
                       positionType = "spatial",
                       positionColData = c("x", "y"),
                       nrow.out = NULL,
                       higherOrderFunction = weightedSpearman,
                       higherOrderFunctionType = "weighted",
                       numberPermutations = 50,
                       numberScaffold = 10,
                       higherOrderSummaryFunction = sd,
                       parallel = FALSE,
                       verbose = FALSE,
                       span = 0.05)

slot(scHOT_spatial, "scHOT_output")

## -----------------------------------------------------------------------------
plotOrderedExpression(scHOT_spatial, c("Arrb1", "Mtor"),
                      positionColData = c("-x", "y"),
                      assayName = "expression")

plotHigherOrderSequence(scHOT_spatial, "Arrb1_Mtor",
                        positionColData = c("-x", "y"))

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