tests/testthat/test_downstream.R

In rliger: Linked Inference of Genomic Experimental Relationships

data("pbmc", package = "rliger")

withNewH5Copy <- function(fun) {
    ctrlpath.orig <- system.file("extdata/ctrl.h5", package = "rliger")
    stimpath.orig <- system.file("extdata/stim.h5", package = "rliger")
    if (!file.exists(ctrlpath.orig))
        stop("Cannot find original h5 file at: ", ctrlpath.orig)
    # if (file.exists("ctrltest.h5")) file.remove("ctrltest.h5")
    # if (file.exists("stimtest.h5")) file.remove("stimtest.h5")
    # pwd <- getwd()
    # # Temp setting for GitHub Actions
    # fsep <- ifelse(Sys.info()["sysname"] == "Windows", "\\", "/")
    # if (Sys.info()["sysname"] == "Windows") {
    #     pwd <- file.path("C:\\Users", Sys.info()["user"], "Documents", fsep = fsep)
    # }

    # ctrlpath <- file.path(pwd, "ctrltest.h5", fsep = fsep)
    # stimpath <- file.path(pwd, "stimtest.h5", fsep = fsep)
    ctrlpath <- tempfile(pattern = "ctrltest_", fileext = ".h5")
    stimpath <- tempfile(pattern = "stimtest_", fileext = ".h5")
    cat("Working ctrl H5 file path: ", ctrlpath, "\n")
    cat("Working stim H5 file path: ", stimpath, "\n")
    file.copy(ctrlpath.orig, ctrlpath, copy.mode = TRUE)
    file.copy(stimpath.orig, stimpath, copy.mode = TRUE)
    if (!file.exists(ctrlpath))
        stop("Cannot find copied h5 file at: ", ctrlpath)
    if (!file.exists(stimpath))
        stop("Cannot find copied h5 file at: ", stimpath)

    fun(list(ctrl = ctrlpath, stim = stimpath))

    if (file.exists(ctrlpath)) unlink(ctrlpath)
    if (file.exists(stimpath)) unlink(stimpath)
}

closeH5Liger <- function(object) {
    for (d in names(object)) {
        if (isH5Liger(object, d)) {
            h5file <- getH5File(object, d)
            h5file$close()
        }
    }
}

process <- function(object, f = TRUE, q = TRUE) {
    object <- normalize(object)
    object <- selectGenes(object)
    object <- scaleNotCenter(object)
    if (f) object <- runOnlineINMF(object, k = 20, minibatchSize = 100)
    if (q) object <- quantileNorm(object)
    return(object)
}

is_online <- function() {
    tryCatch({
        readLines("https://cran.r-project.org/", n = 1)
        TRUE
    },
    warning = function(w) invokeRestart("muffleWarning"),
    error = function(e) FALSE)
}

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Clustering
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

context("Clustering")
test_that("clustering", {
    skip_if_not_installed("RcppPlanc")
    pbmc <- process(pbmc, f = FALSE, q = FALSE)
    expect_error(runCluster(pbmc), "No cell factor loading available")

    pbmc <- runOnlineINMF(pbmc, k = 20, minibatchSize = 100)
    expect_message(runCluster(pbmc, nRandomStarts = 1),
                   "leiden clustering on unaligned")

    expect_message(runCluster(pbmc, nRandomStarts = 1, method = "louvain"),
                   "louvain clustering on unaligned")

    pbmc <- quantileNorm(pbmc)
    expect_message(pbmc <- runCluster(pbmc, nRandomStarts = 1, saveSNN = TRUE),
                   "leiden clustering on aligned")
    expect_is(defaultCluster(pbmc, droplevels = TRUE), "factor")
    expect_is(pbmc@uns$snn, "dgCMatrix")
    expect_message(pbmc <- runCluster(pbmc, nRandomStarts = 1, method = "louvain"),
                   "louvain clustering on aligned")
    expect_message(defaultCluster(pbmc, name = "louvain_cluster") <- "louvain_cluster",
                   "Cannot have")
    expect_error(defaultCluster(pbmc) <- "notexist", "Selected variable does not exist")
    defaultCluster(pbmc) <- pbmc$leiden_cluster
    expect_identical(pbmc$leiden_cluster, pbmc$defaultCluster)
    expect_error(defaultCluster(pbmc) <- factor(letters), "Length of")
    defaultCluster(pbmc) <- NULL
    defaultCluster(pbmc, name = "leiden") <- unname(pbmc$leiden_cluster)
    expect_identical(pbmc$leiden, pbmc$leiden_cluster)

    fakevar <- pbmc$leiden_cluster
    names(fakevar)[1:26] <- letters
    expect_error(defaultCluster(pbmc) <- fakevar, "Not all `names")



    expect_equal(calcPurity(pbmc, "leiden_cluster", "leiden_cluster"), 1)
    expect_error(calcPurity(pbmc, letters, "leiden_cluster"),
                 "Longer/shorter `trueCluster` than cells considered requires")
    expect_message(calcPurity(pbmc, unname(pbmc$leiden_cluster), "leiden_cluster"), "Assuming unnamed")

    expect_equal(calcARI(pbmc, "leiden_cluster", "leiden_cluster"), 1)
    expect_error(calcARI(pbmc, letters, "leiden_cluster"),
                 "Longer/shorter `trueCluster` than cells considered requires")
    expect_message(calcARI(pbmc, unname(pbmc$leiden_cluster), "leiden_cluster"), "Assuming unnamed")

    expect_equal(calcNMI(pbmc, "leiden_cluster", "leiden_cluster"), 1)
    expect_error(calcNMI(pbmc, letters, "leiden_cluster"),
                 "Longer/shorter `trueCluster` than cells considered requires")
    expect_message(calcNMI(pbmc, unname(pbmc$leiden_cluster), "leiden_cluster"), "Assuming unnamed")
    # Tests for singleton grouping. Need to find the case where there are singletons
    # expect_message(runCluster(pbmc, nRandomStarts = 1,
    #                          partitionType = "CPMVertexPartition"),
    #                "63 singletons identified. 112 final clusters")
    # pbmc <- runCluster(pbmc, nRandomStarts = 1,
    #                          partitionType = "CPMVertexPartition",
    #                          groupSingletons = FALSE)
    # expect_equal(nlevels(pbmc$leiden_cluster), 113)
    # expect_true("singleton" %in% pbmc$leiden_cluster)

    # Test downsampling with the info already here
    pbmc.small1 <- downsample(pbmc, maxCells = 100)
    expect_true(all.equal(sapply(pbmc.small1@datasets, ncol), c(ctrl = 44, stim = 56)))
    idx <- downsample(pbmc, balance = c("dataset"), maxCells = 100,
                              returnIndex = TRUE)
    expect_is(idx, "integer")
    expect_equal(length(idx), 200)

    pbmc <- mapCellMeta(pbmc, from = "dataset", newTo = "pseudo",
                        `ctrl` = "CTRL")
    expect_identical(pbmc$pseudo,
                     setNames(factor(c(rep("CTRL", 300), rep("stim", 300))),
                              colnames(pbmc)))
})

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Dimensionality reduction
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

context("dimensionality reduction")
test_that("dimensionality reduction", {
    skip_if_not_installed("RcppPlanc")
    pbmc <- process(pbmc)
    expect_message(runUMAP(pbmc, useRaw = TRUE),
                   "Generating UMAP on unaligned")
    expect_error(dimRed(pbmc), "available in this")
    expect_message(pbmc <- runUMAP(pbmc, useRaw = FALSE),
                   "Generating UMAP on aligned")
    pbmc@uns$defaultDimRed <- NULL
    expect_message(dimRed(pbmc), "No default")
    defaultDimRed(pbmc) <- "UMAP"
    expect_error(defaultDimRed(pbmc) <- letters, "Can only set one")
    expect_identical(defaultDimRed(pbmc), dimRed(pbmc, "UMAP"))
    expect_equal(dim(dimRed(pbmc)), c(ncol(pbmc), 2))
    expect_no_error(dimRed(pbmc, "UMAP2") <- dimRed(pbmc, "UMAP"))
    expect_equal(nrow(dimRed(pbmc, name = 1, cellIdx = 1:10)), 10)
    expect_equal(nrow(dimRed(pbmc, name = 1, useDatasets = names(pbmc))), ncol(pbmc))
    expect_equal(nrow(dimRed(pbmc, name = "UMAP", cellIdx = 1:10)), 10)
    expect_equal(nrow(dimRed(pbmc, name = "UMAP", useDatasets = names(pbmc))), ncol(pbmc))
    expect_no_error(dimRed(pbmc, 2) <- NULL)

    expect_message(runTSNE(pbmc, useRaw = TRUE),
                   "Generating TSNE \\(Rtsne\\) on unaligned")
    expect_message(pbmc <- runTSNE(pbmc, useRaw = FALSE),
                   "Generating TSNE \\(Rtsne\\) on aligned")
    expect_equal(dim(dimRed(pbmc, "TSNE")), c(ncol(pbmc), 2))

    expect_error(runTSNE(pbmc, method = "fft"),
                 "Please pass in path to FIt-SNE directory as fitsne.path.")
})

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Differential Expression
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

context("Differential Expression")
test_that("wilcoxon", {
    skip_if_not_installed("RcppPlanc")
    expect_error(runMarkerDEG(pbmc),
                 "No `conditionBy` given or default cluster not set")
    pbmc <- process(pbmc)
    pbmc <- runCluster(pbmc, nRandomStarts = 1)
    res0 <- runMarkerDEG(pbmc, conditionBy = "dataset", useDatasets = 1, method = "wilcox")
    expect_true(all(is.nan(res0$pval)))

    res1 <- runMarkerDEG(pbmc, method = "wilcox")
    expect_equal(dim(res1), c(249 * nlevels(pbmc$leiden_cluster), 10))
    res2 <- runMarkerDEG(pbmc, conditionBy = "dataset", splitBy = "leiden_cluster", method = "wilcox")
    expect_is(res2, "data.frame")
    hm1 <- plotMarkerHeatmap(pbmc, res1, dedupBy = "l")
    hm2 <- plotMarkerHeatmap(pbmc, res1, dedupBy = "p")
    expect_is(hm1, "HeatmapList")
    expect_is(hm2, "HeatmapList")
    expect_is(plotVolcano(res1, 0), "ggplot")
    if (requireNamespace("EnhancedVolcano", quietly = TRUE)) {
        expect_is(plotEnhancedVolcano(res1, 0), "ggplot")
    }

    expect_error(getFactorMarkers(pbmc, "ctrl", "stim", factorShareThresh = 0),
                 "No factor passed the dataset specificity threshold")
    expect_message(
        res3 <- getFactorMarkers(pbmc, "ctrl", "stim", printGenes = TRUE)
    )
    expect_is(res3, "list")
    expect_identical(names(res3), c("ctrl", "shared", "stim", "num_factors_V1",
                                   "num_factors_V2"))

    skip_if_not_installed("gprofiler2")
    expect_error(runGOEnrich(res1, group = "a"),
                 "Selected groups not available")
    expect_error(runGOEnrich(res1, group = 0, orderBy = c("logFC", "pval")),
                 "Only one `orderBy`")
    expect_error(runGOEnrich(res1, orderBy = "score"),
                 "`orderBy` should be one of")
    if (is_online()) {
        go1 <- runGOEnrich(res1, group = 0, orderBy = "logFC", significant = FALSE)
        expect_is(go1, "list")
        expect_is(go1$`0`$result, "data.frame")
        go2 <- runGOEnrich(res1, group = 0, orderBy = "pval", significant = FALSE)
        expect_is(go2, "list")
        expect_is(go2$`0`$result, "data.frame")
        go3 <- runGOEnrich(res1, group = c(0, 1), orderBy = "pval", significant = FALSE)

        expect_is(plotGODot(go1, pvalThresh = 1), "ggplot")
        expect_error(plotGODot(go1, group = "ctrl"), "Specified group not available")
        expect_message(plotGODot(go1, group = '0'), "No enough matching")
        expect_is(plotGODot(go3, pvalThresh = 1), "list")
    }
})


test_that("pseudo bulk", {
    skip_if_not_installed("RcppPlanc")
    skip_if_not_installed("DESeq2")
    nIsecGenes <- length(Reduce(intersect, lapply(rawData(pbmc), rownames)))
    pbmc <- process(pbmc)
    pbmc <- runCluster(pbmc, nRandomStarts = 1)
    res1 <- runPairwiseDEG(pbmc, groupTest = pbmc$leiden_cluster == 1,
                           groupCtrl = pbmc$leiden_cluster == 2,
                           method = "pseudo")
    expect_is(res1, "data.frame")
    expect_true(all.equal(dim(res1), c(nIsecGenes, 7)))
    res2 <- runPairwiseDEG(pbmc, groupTest = 1, groupCtrl = 2,
                           variable1 = "leiden_cluster",
                           method = "pseudo", useReplicate = "dataset")
    expect_is(res2, "data.frame")
    expect_true(all.equal(dim(res2), c(nIsecGenes, 7)))
    res3 <- runPairwiseDEG(pbmc, groupTest = 1, groupCtrl = 2,
                           variable1 = "leiden_cluster",
                           method = "pseudo")
    expect_true(all.equal(res1[,-2], res3[,-2])) # Different in "group" column
    pbmc$leiden2 <- pbmc$leiden_cluster
    res4 <- runPairwiseDEG(pbmc, groupTest = 1, groupCtrl = 2,
                           variable1 = "leiden_cluster", variable2 = "leiden2",
                           method = "pseudo", useReplicate = "dataset")
    expect_true(all.equal(res2[,-2], res4[,-2])) # Different in "group" column

    expect_error(runPairwiseDEG(pbmc, variable2 = "yo"),
                 "Please see")

    expect_message(
        runPairwiseDEG(
            pbmc, groupTest = pbmc$dataset == "ctrl" & pbmc$leiden_cluster == 0,
            groupCtrl = pbmc$dataset == "stim" & pbmc$leiden_cluster == 0,
            method = "pseudo", useReplicate = "dataset", nPsdRep = 1
        ),
        "Too few replicates"
    )

    pbmc@datasets$ctrl@rawData <- NULL
    expect_error(
        runPairwiseDEG(
            pbmc, groupTest = 1, groupCtrl = 2,
            variable1 = "leiden_cluster",　method = "pseudo",
            useReplicate = "dataset"
        ),
        "not all available for involved datasets"
    )
})

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# GSEA
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

context("GSEA")
test_that("gsea", {
    skip_if_not(is_online())
    skip_if_not_installed("org.Hs.eg.db")
    skip_if_not_installed("reactome.db")
    skip_if_not_installed("fgsea")
    skip_if_not_installed("AnnotationDbi")
    custom <- list(
        `Immune System` = c("9636", "2633", "6282", "6280", "6279", "2207", "2214",
                            "6402", "91543", "6233", "10578", "3553", "5473",
                            "3627", "51316", "929", "972")
    )
    expect_warning({
        expect_is(runGSEA(pbmcPlot, genesets = "Immune System"), "list")
        expect_is(runGSEA(pbmcPlot, customGenesets = custom), "list")
    })
})


#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# ATAC
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

context("ATAC")
data("bmmc")
test_that("ATAC", {
    skip_if_not_installed("RcppPlanc")
    skip_if_not_installed("GenomicRanges")
    skip_if_not_installed("IRanges")
    skip_if_not_installed("psych")
    bmmc <- normalize(bmmc)
    bmmc <- selectGenes(bmmc)
    bmmc <- scaleNotCenter(bmmc)
    bmmc <- runOnlineINMF(bmmc, minibatchSize = 80)
    bmmc <- quantileNorm(bmmc)
    bmmc <- normalizePeak(bmmc)
    bmmc <- imputeKNN(bmmc, reference = "atac", queries = "rna")
    expect_is(dataset(bmmc, "rna"), "ligerATACDataset")
    corr <- linkGenesAndPeaks(
        bmmc, useDataset = "rna",
        pathToCoords = system.file("extdata/hg19_genes.bed",
                                   package = "rliger")
    )
    expect_is(corr, "dgCMatrix")
})