Simulating complex design scRNA-seq data with muscat
October 24, 2023
Abstract
muscat: multi-sample multi-group scRNA-seq analysis tools (???) provides a straightforward but effective simulation framework that is anchored to a labeled multi-sample multi-subpopulation scRNA-seq reference dataset, uses (non-zero-inflated) negative binomial (NB) as the canonical distribution for droplet scRNA-seq datasets, and exposes various parameters to modulate: the number of subpopulations and samples simulated, the number of cells per subpopulation (and sample), and the type and magnitude of a wide range of patterns of differential expression.
This vignette serves to provide the underlying theoretical background, to thoroughly describe the various input arguments, and to demonstrate the simulation framework’s current capabilities using some illustrative (not necessarily realistic) examples.
Package
muscat 1.16.0
Contents
For details on the concepts presented here, consider having a look at our publication:
Crowell HL, Soneson C*, Germain P-L*, Calini D,
Collin L, Raposo C, Malhotra D, and Robinson MD:
muscat detects subpopulation-specific state transitions from
multi-sample multi-condition single-cell transcriptomics data.
Nature Communications 11, 6077 (2020).
DOI: 10.1038/s41467-020-19894-4
Load packages
library(cowplot)
library(dplyr)
library(reshape2)
library(muscat)
library(purrr)
library(scater)
library(SingleCellExperiment)Data description
To demonstrate muscat’s simulation framework, we will use a SingleCellExperiment (SCE) containing 10x droplet-based scRNA-seq PBCM data from 8 Lupus patients obtained befor and after 6h-treatment with IFN-\(\beta\) (Kang et al. 2018). The complete raw data, as well as gene and cell metadata is available through the NCBI GEO, accession number GSE96583.
1 Simulation framework
muscat’s simulation framework comprises: i) estimation of negative binomial (NB) parameters from a reference multi-subpopulation, multi-sample dataset; ii) sampling of gene and cell parameters to use for simulation; and, iii) simulation of gene expression data as NB distributions of mixtures thereof. See Fig. 1.
Let \(Y = (y_{gc})\in\mathbb{N}_0^{G\times C}\) denote the count matrix of a multi-sample multi-subpopulation reference dataset with genes \(\mathcal{G} = \{ g_1, \ldots, g_G \}\) and sets of cells \(\mathcal{C}_{sk} = \{ c^{sk}_1, ..., c^{sk}_{C_{sk}} \}\) for each sample \(s\) and subpopulation \(k\) (\(C_{sk}\) is the number of cells for sample \(s\), subpopulation \(k\)). For each gene \(g\), we fit a model to estimate sample-specific means \(\beta_g^s\), for each sample \(s\), and dispersion parameters \(\phi_g\) using ’s function with default parameters. Thus, we model the reference count data as NB distributed:
\[Y_{gc} \sim NB(\mu_{gc}, \phi_g)\]
for gene \(g\) and cell \(c\), where the mean \(\mu_{gc} = \exp(\beta_{g}^{s(c)}) \cdot \lambda_c\). Here, \(\beta_{g}^{s(c)}\) is the relative abundance of gene \(g\) in sample \(s(c)\), \(\lambda_c\) is the library size (total number of counts), and \(\phi_g\) is the dispersion.
Figure 1: Schematic overview of muscat’s simulation framework
Given a count matrix of features by cells and, for each cell, pre-determined subpopulation identifiers as well as sample labels (0), dispersion and sample-wise means are estimated from a negative binomial distribution for each gene (for each subpopulation) (1.1); and library sizes are recorded (1.2). From this set of parameters (dispersions, means, library sizes), gene expression is sampled from a negative binomial distribution. Here, genes are selected to be “type” (subpopulation-specifically expressed; e.g., via marker genes), “state” (change in expression in a condition-specific manner) or equally expressed (relatively) across all samples (2). The result is a matrix of synthetic gene expression data (3).
For each subpopulation, we randomly assign each gene to a given differential distribution (DD) category (Korthauer et al. 2016) according to a probability vector p_dd \(=(p_{EE},p_{EP},p_{DE},p_{DP},p_{DM},p_{DB})\). For each gene and subpopulation, we draw a vector of fold changes (FCs) from a Gamma distribution with shape parameter \(\alpha=4\) and rate \(\beta=4/\mu_\text{logFC}\), where \(\mu_\text{logFC}\) is the desired average logFC across all genes and subpopulations specified via argument lfc. The direction of differential expression is randomized for each gene, with equal probability of up- and down-regulation.
Next, we split the cells in a given subpopulations into two sets (representing treatment groups), \(\mathcal{T}_A\) and \(\mathcal{T}_B\), which are in turn split again into two sets each (representing subpopulations within the given treatment group.), \(\mathcal{T}_{A_1}/\mathcal{T}_{A_2}\) and \(\mathcal{T}_{B_1}/\mathcal{T}_{B_2}\).
For EE genes, counts for \(\mathcal{T}_A\) and \(\mathcal{T}_B\) are drawn using identical means.For EP genes, we multiply the effective means for identical fractions of cells per group by the sample FCs, i.e., cells are split such that \(\dim\mathcal{T}_{A_1} = \dim\mathcal{T}_{B_1}\) and \(\dim\mathcal{T}_{A_2} = \dim\mathcal{T}_{B_2}\). For DE genes, the means of one group, \(A\) or \(B\), are multiplied with the samples FCs. DP genes are simulated analogously to EP genes with \(\dim\mathcal{T}_{A_1} = a\cdot\dim\mathcal{T}_A\) and \(\dim\mathcal{T}_{B_1} = b\cdot\dim\mathcal{T}_B\), where \(a+b=1\) and \(a\neq b\). For DM genes, 50% of cells from one group are simulated at \(\mu\cdot\text{logFC}\). For DB genes, all cells from one group are simulated at \(\mu\cdot\text{logFC}/2\), and the second group is split into equal proportions of cells simulated at \(\mu\) and \(\mu\cdot\text{logFC}\), respectively. See Fig. 2.
Figure 2: Schematic of the various types of differential distributions supported by muscat’s simulation framework
Differential distributions are simulated from a NB distribution or mixtures thereof, according to the definitions of random variables X, Y and Z.
1.1 prepSim: Preparing data for simulation
To prepare a reference SingleCellExperiment (SCE) for simulation of multi-sample multi-group scRNA-seq data, prepSim will
- perform basic filtering of genes and cells
- (optionally) filter for subpopulation-sample instances with a threshold number of cells to assure accurate parameter estimation
- estimate cell (library sizes) and gene parameters (dispersions and sample-specific means)
Importantly, we want to introduce known changes in states across conditions; thus, only replicates from a single condition should go into the simulation. The group to be kept for simulation may be specified via group_keep, in which case samples from all other groups (sce$group_id != group_keep) will be dropped. By default (group_keep = NULL), prepSim will select the first group available as reference.
Arguments min_count, min_cells, min_genes and min_size are used to tune the filtering of genes, cells and subpopulation-instances as follows:
- only genes with a count
> min_countin>= min_cellswill be retained - only cells with a count
> 0in>= min_geneswill be retained - only subpopulation-sample instances with
>= min_sizecells will be retained;min_size = NULLwill skip this step
# estimate simulation parameters
data(example_sce)
ref <- prepSim(example_sce, verbose = FALSE)
# only samples from `ctrl` group are retained
table(ref$sample_id)##
## ctrl101 ctrl107
## 200 200# cell parameters: library sizes
sub <- assay(example_sce[rownames(ref), colnames(ref)])
all.equal(exp(ref$offset), as.numeric(colSums(sub)))## [1] "names for target but not for current"
## [2] "Mean relative difference: 0.4099568"# gene parameters: dispersions & sample-specific means
head(rowData(ref))## DataFrame with 6 rows and 4 columns
## ENSEMBL SYMBOL beta disp
## <character> <character> <DataFrame> <numeric>
## ISG15 ENSG00000187608 ISG15 -7.84574:-0.24711310:-1.039480 4.6360532
## AURKAIP1 ENSG00000175756 AURKAIP1 -7.84859: 0.00768812:-1.171896 0.1784345
## MRPL20 ENSG00000242485 MRPL20 -8.31434:-0.58684477:-0.304827 0.6435324
## SSU72 ENSG00000160075 SSU72 -8.05160:-0.17119703:-0.793222 0.0363892
## RER1 ENSG00000157916 RER1 -7.75327: 0.10731331:-1.261821 0.5046166
## RPL22 ENSG00000116251 RPL22 -8.03553:-0.03357193: 0.143506 0.20236321.2 simData: Simulating complex designs
Provided with a reference SCE as returned by prepSim, a variery of simulation scenarios can be generated using the simData function, which will again return an SCE containg the following elements:
assaycountscontaining the simulated count datacolDatacolumnscluster/sample/group_idcontaining each cells cluster, sample, and group ID (A or B).metadata$gene_infocontaining adata.framelisting, for each gene and cluster- the simulationed DD
category - the sampled
logFC; note that this will only approximatelog2(sim_mean.B/sim_mean.A)for genes of thedecategory as other types of state changes use mixtures for NBs, and will consequently not exhibit a shift in means of the same magnitude aslogFC - the reference
sim_genefrom which dispersionsim_dispand sample-specific meansbeta.<sample_id>were used - the simulated expression means
sim_mean.A/Bfor each group
- the simulationed DD
In the code chunk that follows, we run a simple simulation with
p_dd = c(1,0,...0), i.e., 10% of EE genesnk = 3subpopulations andns = 3replicates for each of 2 groupsng = 1000genes andnc = 2000cells, resulting in2000/2/ns/nk\(\approx111\) cells for 2 groups with 3 samples each and 3 subpopulations
# simulated 10% EE genes
sim <- simData(ref, p_dd = diag(6)[1, ],
nk = 3, ns = 3, nc = 2e3,
ng = 1e3, force = TRUE)
# number of cells per sample and subpopulation
table(sim$sample_id, sim$cluster_id)##
## cluster1 cluster2 cluster3
## sample1.A 102 119 109
## sample2.A 131 85 109
## sample3.A 115 94 126
## sample1.B 117 91 112
## sample2.B 98 121 120
## sample3.B 110 126 115By default, we have drawn a random reference sample from levels(ref$sample_id) for every simulated sample in each group, resulting in an unpaired design:
metadata(sim)$ref_sids## A B
## sample1 "ctrl107" "ctrl101"
## sample2 "ctrl101" "ctrl101"
## sample3 "ctrl101" "ctrl107"Alternatively, we can re-run the above simulation with paired = TRUE such that both groups will use the same set of reference samples, resulting in a paired design:
# simulated paired design
sim <- simData(ref, paired = TRUE,
nk = 3, ns = 3, nc = 2e3,
ng = 1e3, force = TRUE)
# same set of reference samples for both groups
ref_sids <- metadata(sim)$ref_sids
all(ref_sids[, 1] == ref_sids[, 2])## [1] TRUE1.2.1 p_dd: Simulating differential distributions
Argument p_dd specifies the fraction of cells to simulate for each DD category. Its values should thus lie in \([0,1]\) and sum to 1. Expression densities for an exemplary set of genes simulated from the code below is shown in Fig. 3.
# simulate genes from all DD categories
sim <- simData(ref, p_dd = c(0.5, rep(0.1, 5)),
nc = 2e3, ng = 1e3, force = TRUE)We can retrieve the category assigned to each gene in each cluster from the gene_info table stored in the output SCE’s metadata:
gi <- metadata(sim)$gene_info
table(gi$category)##
## ee ep de dp dm db
## 994 206 210 212 202 176
Figure 3: Expression densities for an exemplary set of 3 genes per differential distribution category
Each density corresponds to one sample, lines are colored by group ID, and panels are split by gene and subpopulation.
1.2.2 rel_lfc: Simulating cluster-specific state changes
By default, for each gene and subpopulation, we draw a vector of fold changes (FCs) from a Gamma distribution with rate parameter \(\beta\propto\mu_\text{logFC}\), where \(\mu_\text{logFC}\) is the desired average logFC across all genes and subpopulations specified via argument lfc. This results in state changes that are of same magnitute for each subpopulation.
Now, suppose we wanted to have a subpopulation that does not exhibit any state changes across conditions, or vary the magnitute of changes across subpopulations. To this end, argument rel_lfc supplies a subpopulation-specific factor applied to the FCs sampled for subpopulation. Fig. 4 demonstrates how this manifests in in two-dimensional embeddings of the cells: Here, we generate a set of 3 simulations with
- equal magnitute of change for all subpopulations:
rel_lfc=c(1,1,1) - stronger change for one cluster:
rel_lfc=c(1,1,3) - cluster-specific FC factors with no change for one cluster:
rel_lfc=c(0,1,2)
Figure 4: t-SNEs of exemplary simulations demonstrating rel_lfc’s effect to induce cluster-specific state changes
Cells are colored by cluster ID (top-row) and group ID (bottom-row), respectively. From left to right: No cluster-specific changes, stronger change for cluster3, different logFC factors for all clusters with no change for cluster1.
1.2.3 p_type: Simulating type features
The idea underlying differential state (DS) analysis to test for subpopulation-specific changes in expression across experimental conditions is based on the idea that we i) use stable moleculare signatures (i.e., type features) to group cells into meaningful subpopulations; and, ii) perform statistical tests on state features that are more transiently expression and may be subject to changes in expression upon, for example, treatment or during disease.
The fraction of type features introduced into each subpopulation is specified via argument p_type. Note that, without introducing any differential states, a non-zero fraction of type genes will result in separation of cells into clusters. Fig. 5 demonstrates how increasing values for p_type lead to more and more separation of the cells when coloring by cluster ID, but that the lack of state changes leads to homogenous mixing of cells when coloring by group ID.
Figure 5: t-SNEs of exemplary simulations demonstrating p_type’s effect to introduce type features
Cells are colored by cluster ID (top-row) and group ID (bottom-row), respectively. The percentage of type features increases from left to right (1, 5, 10%). Simulations are pure EE, i.e., all genes are non-differential across groups.
1.3 Simulation a hierarchical cluster structure
simData contains three parameters that control how subpopulations relate to and differ from one another:
p_typedetermines the percentage of type genes exclusive to each clusterphylo_treerepresents a phylogenetic tree specifying of clusters relate to one anotherphylo_parscontrols how branch distances are to be interpreted
Note that, when supplied with a cluster phylogeny, argument nk is ignored and simData extracts the number of clusters to be simulated from phylo_tree.
1.3.1 p_type: Introducing type features
To exemplify the effect of the parameter p_type, we simulate a dataset with \(\approx5\%\) of type genes per cluster, and one group only via probs = list(..., c(1, 0) (i.e., \(\text{Prob}(\textit{cell is in group 2}) = 0\)):
# simulate 5% of type genes; one group only
sim <- simData(ref, p_type = 0.1,
nc = 2e3, ng = 1e3, force = TRUE,
probs = list(NULL, NULL, c(1, 0)))
# do log-library size normalization
sim <- logNormCounts(sim)For visualizing the above simulation, we select for genes that are of class type (rowData()$class == "type") and have a decent simulated expression mean. Furthermore, we sample a subset of cells for each cluster. The resulting heatmap (Fig. 6) shows that the 3 clusters separate well from one another, but that type genes aren’t necessarily expressed higher in a single cluster. This is the case because a gene selected as reference for a type gene in a given cluster may indeed have a lower expression than the gene used for the remainder of clusters.
# extract gene metadata & number of clusters
rd <- rowData(sim)
nk <- nlevels(sim$cluster_id)
# filter for type genes with high expression mean
is_type <- rd$class == "type"
is_high <- rowMeans(assay(sim, "logcounts")) > 1
gs <- rownames(sim)[is_type & is_high]
# sample 100 cells per cluster for plotting
cs <- lapply(split(seq_len(ncol(sim)), sim$cluster_id), sample, 100)
plotHeatmap(sim[, unlist(cs)], features = gs, center = TRUE,
colour_columns_by = "cluster_id", cutree_cols = nk)
Exemplary heatmap demonstrating the effect of p_type to introduce cluster-specific type genes. Included are type genes (= rows) with a simulated expression mean > 1, and a random subset of 100 cells (= columns) per cluster; column annotations represent cluster IDs. Bins are colored by expression scaled in row direction, and both genes and cells are hierarchically clustered.
1.3.2 phylo_tree: Introducing a cluster phylogeny
The scenario illustrated above is arguably not very realistic. Instead, in a biology setting, subpopulations don’t differ from one another by a specific subset of genes, but may share some of the genes decisive for their biological role. I.e., the set type features is not exclusive for every given subpopulation, and some subpopulations are more similar to one another than others.
To introduce a more realistic subpopulation structure, simData can be supplied with a phylogenetic tree, phylo_tree, that specifies the relationship and distances between clusters. The tree should be written in Newick format as in the following example:
# specify cluster phylogeny
tree <- "(('cluster1':0.4,'cluster2':0.4):0.4,('cluster3':
0.5,('cluster4':0.2,'cluster5':0.2,'cluster6':0.2):0.4):0.4);"
# visualize cluster tree
library(phylogram)
plot(read.dendrogram(text = tree))
Figure 7: Exemplary phylogeny
The phylogenetic tree specified via phylo relates 3 clusters such that there are 2 main branches, and clusters 1 and 2 should be more similar to one another than cluster 3.
# simulate 5% of type genes; one group only
sim <- simData(ref,
phylo_tree = tree, phylo_pars = c(0.1, 1),
nc = 800, ng = 1e3, dd = FALSE, force = TRUE)
# do log-library size normalization
sim <- logNormCounts(sim)# extract gene metadata & number of clusters
rd <- rowData(sim)
nk <- nlevels(sim$cluster_id)
# filter for type & shared genes with high expression mean
is_type <- rd$class != "state"
is_high <- rowMeans(assay(sim, "logcounts")) > 1
gs <- rownames(sim)[is_type & is_high]
# sample 100 cells per cluster for plotting
cs <- lapply(split(seq_len(ncol(sim)), sim$cluster_id), sample, 50)
plotHeatmap(sim[, unlist(cs)], features = gs,
center = TRUE, show_rownames = FALSE,
colour_columns_by = "cluster_id")
Exemplary heatmap demonstrating the effect of phylo_tree to introduce a hierarchical cluster structure. Included are 100 randomly sampled non-state, i.e. type or shared, genes (= rows) with a simulated expression mean > 1, and a random subset of 100 cells (= columns) per cluster; column annotations represent cluster IDs. Bins are colored by expression scaled in row direction, and both genes and cells are hierarchically clustered.
1.4 Simulating batch effects
under development.
2 Quality control
As is the case with any simulation, it is crutial to verify the qualitation of the simulated data; i.e., how well key characteristics of the reference data are captured in the simulation. While we have demonstrated that muscats simulation framework is capable of reproducing key features of scRNA-seq dataset at both the single-cell and pseudobulk level (???), simulation quality will vary depending on the reference dataset and could suffer from too extreme simulation parameters. Therefore, we advise anyone interested in using the framework presented herein for any type of method evaluation or comparison to generate countsimQC report (Soneson and Robinson 2018) as it is extremly simple to make and very comprehensive.
The code chunk below (not evaluated here) illustrates how to generate a report comparing an exemplary simData simulation with the reference data provided in ref. Runtimes are mainly determined by argument maxNForCorr and maxNForDisp, and computing a full-blown report can be very time intensive. We thus advice using a sufficient but low number of cells/genes for these steps.
# load required packages
library(countsimQC)
library(DESeq2)
# simulate data
sim <- simData(ref,
ng = nrow(ref),
nc = ncol(ref),
dd = FALSE)
# construct 'DESeqDataSet's for both,
# simulated and reference dataset
dds_sim <- DESeqDataSetFromMatrix(
countData = counts(sim),
colData = colData(sim),
design = ~ sample_id)
dds_ref <- DESeqDataSetFromMatrix(
countData = counts(ref),
colData = colData(ref),
design = ~ sample_id)
dds_list <- list(sim = dds_sim, data = dds_ref)
# generate 'countsimQC' report
countsimQCReport(
ddsList = dds_list,
outputFile = "<file_name>.html",
outputDir = "<output_path>",
outputFormat = "html_document",
maxNForCorr = 200,
maxNForDisp = 500)3 Method benchmarking
A variety of functions for calculation and visualizing performance metrics for evaluation of ranking and binary classification (assignment) methods is provided in the iCOBRA package (Soneson and Robinson 2016).
We firstly define a wrapper that takes as input a method passed pbDS and reformats the results as a data.frame in tidy format, which is in turn right_joined with simulation gene metadata. As each methods may return results for different subsets of gene-subpopulation instances, the latter steps assures that the dimensions of all method results will match.
# 'm' is a character string specifying a valid `pbDS` method
.run_method <- function(m) {
res <- pbDS(pb, method = m, verbose = FALSE)
tbl <- resDS(sim, res)
left_join(gi, tbl, by = c("gene", "cluster_id"))
}Having computed result data.frames for a set of methods, we next define a wrapper that prepares the data for evaluation with iCOBRA using the COBRAData constructor, and calculates any performance measures of interest (specified via aspects) with calculate_performance:
# 'x' is a list of result 'data.frame's
.calc_perf <- function(x, facet = NULL) {
cd <- COBRAData(truth = gi,
pval = data.frame(bind_cols(map(x, "p_val"))),
padj = data.frame(bind_cols(map(x, "p_adj.loc"))))
perf <- calculate_performance(cd,
binary_truth = "is_de", maxsplit = 1e6,
splv = ifelse(is.null(facet), "none", facet),
aspects = c("fdrtpr", "fdrtprcurve", "curve"))
}Putting it all together, we can finally simulate some data, run a set of DS analysis methods, calculate their performance, and plot a variety of performance metrics depending on the aspects calculated by .calc_perf:
# simulation with all DD types
sim <- simData(ref,
p_dd = c(rep(0.3, 2), rep(0.1, 4)),
ng = 1e3, nc = 2e3, ns = 3, force = TRUE)
# aggregate to pseudobulks
pb <- aggregateData(sim)
# extract gene metadata
gi <- metadata(sim)$gene_info
# add truth column (must be numeric!)
gi$is_de <- !gi$category %in% c("ee", "ep")
gi$is_de <- as.numeric(gi$is_de)
# specify methods for comparison & run them
# (must set names for methods to show in visualizations!)
names(mids) <- mids <- c("edgeR", "DESeq2", "limma-trend", "limma-voom")
res <- lapply(mids, .run_method)
# calculate performance measures
# and prep. for plotting with 'iCOBRA'
library(iCOBRA)
perf <- .calc_perf(res, "cluster_id")
pd <- prepare_data_for_plot(perf)
# plot FDR-TPR curves by cluster
plot_fdrtprcurve(pd) +
theme(aspect.ratio = 1) +
scale_x_continuous(trans = "sqrt") +
facet_wrap(~ splitval, nrow = 1)
Session info
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] iCOBRA_1.30.0 DESeq2_1.42.0
## [3] phylogram_2.1.0 reshape2_1.4.4
## [5] UpSetR_1.4.0 scater_1.30.0
## [7] scuttle_1.12.0 Matrix_1.6-1.1
## [9] muscData_1.15.0 SingleCellExperiment_1.24.0
## [11] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [13] GenomicRanges_1.54.0 GenomeInfoDb_1.38.0
## [15] IRanges_2.36.0 S4Vectors_0.40.0
## [17] MatrixGenerics_1.14.0 matrixStats_1.0.0
## [19] ExperimentHub_2.10.0 AnnotationHub_3.10.0
## [21] BiocFileCache_2.10.0 dbplyr_2.3.4
## [23] BiocGenerics_0.48.0 purrr_1.0.2
## [25] muscat_1.16.0 limma_3.58.0
## [27] ggplot2_3.4.4 dplyr_1.1.3
## [29] cowplot_1.1.1 BiocStyle_2.30.0
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-7 httr_1.4.7
## [3] RColorBrewer_1.1-3 doParallel_1.0.17
## [5] numDeriv_2016.8-1.1 tools_4.3.1
## [7] sctransform_0.4.1 backports_1.4.1
## [9] DT_0.30 utf8_1.2.4
## [11] R6_2.5.1 uwot_0.1.16
## [13] mgcv_1.9-0 GetoptLong_1.0.5
## [15] withr_2.5.1 prettyunits_1.2.0
## [17] gridExtra_2.3 cli_3.6.1
## [19] Cairo_1.6-1 sandwich_3.0-2
## [21] labeling_0.4.3 sass_0.4.7
## [23] mvtnorm_1.2-3 blme_1.0-5
## [25] parallelly_1.36.0 RSQLite_2.3.1
## [27] generics_0.1.3 shape_1.4.6
## [29] gtools_3.9.4 ggbeeswarm_0.7.2
## [31] fansi_1.0.5 abind_1.4-5
## [33] lifecycle_1.0.3 multcomp_1.4-25
## [35] yaml_2.3.7 edgeR_4.0.0
## [37] gplots_3.1.3 SparseArray_1.2.0
## [39] Rtsne_0.16 grid_4.3.1
## [41] blob_1.2.4 promises_1.2.1
## [43] crayon_1.5.2 shinydashboard_0.7.2
## [45] lattice_0.22-5 beachmat_2.18.0
## [47] KEGGREST_1.42.0 magick_2.8.1
## [49] pillar_1.9.0 knitr_1.44
## [51] ComplexHeatmap_2.18.0 rjson_0.2.21
## [53] boot_1.3-28.1 estimability_1.4.1
## [55] corpcor_1.6.10 future.apply_1.11.0
## [57] codetools_0.2-19 glue_1.6.2
## [59] data.table_1.14.8 vctrs_0.6.4
## [61] png_0.1-8 Rdpack_2.5
## [63] gtable_0.3.4 cachem_1.0.8
## [65] xfun_0.40 rbibutils_2.2.15
## [67] S4Arrays_1.2.0 mime_0.12
## [69] coda_0.19-4 survival_3.5-7
## [71] pheatmap_1.0.12 iterators_1.0.14
## [73] statmod_1.5.0 interactiveDisplayBase_1.40.0
## [75] ellipsis_0.3.2 TH.data_1.1-2
## [77] ROCR_1.0-11 nlme_3.1-163
## [79] pbkrtest_0.5.2 bit64_4.0.5
## [81] progress_1.2.2 EnvStats_2.8.1
## [83] filelock_1.0.2 RcppAnnoy_0.0.21
## [85] bslib_0.5.1 TMB_1.9.6
## [87] irlba_2.3.5.1 vipor_0.4.5
## [89] KernSmooth_2.23-22 colorspace_2.1-0
## [91] DBI_1.1.3 tidyselect_1.2.0
## [93] emmeans_1.8.9 bit_4.0.5
## [95] compiler_4.3.1 curl_5.1.0
## [97] BiocNeighbors_1.20.0 DelayedArray_0.28.0
## [99] bookdown_0.36 scales_1.2.1
## [101] caTools_1.18.2 remaCor_0.0.16
## [103] rappdirs_0.3.3 stringr_1.5.0
## [105] digest_0.6.33 shinyBS_0.61.1
## [107] minqa_1.2.6 variancePartition_1.32.0
## [109] rmarkdown_2.25 aod_1.3.2
## [111] XVector_0.42.0 RhpcBLASctl_0.23-42
## [113] htmltools_0.5.6.1 pkgconfig_2.0.3
## [115] lme4_1.1-34 sparseMatrixStats_1.14.0
## [117] fastmap_1.1.1 htmlwidgets_1.6.2
## [119] rlang_1.1.1 GlobalOptions_0.1.2
## [121] shiny_1.7.5.1 DelayedMatrixStats_1.24.0
## [123] farver_2.1.1 jquerylib_0.1.4
## [125] zoo_1.8-12 jsonlite_1.8.7
## [127] BiocParallel_1.36.0 BiocSingular_1.18.0
## [129] RCurl_1.98-1.12 magrittr_2.0.3
## [131] GenomeInfoDbData_1.2.11 munsell_0.5.0
## [133] Rcpp_1.0.11 ape_5.7-1
## [135] viridis_0.6.4 stringi_1.7.12
## [137] zlibbioc_1.48.0 MASS_7.3-60
## [139] plyr_1.8.9 parallel_4.3.1
## [141] listenv_0.9.0 ggrepel_0.9.4
## [143] Biostrings_2.70.0 splines_4.3.1
## [145] hms_1.1.3 circlize_0.4.15
## [147] locfit_1.5-9.8 ScaledMatrix_1.10.0
## [149] BiocVersion_3.18.0 evaluate_0.22
## [151] BiocManager_1.30.22 nloptr_2.0.3
## [153] foreach_1.5.2 httpuv_1.6.12
## [155] tidyr_1.3.0 future_1.33.0
## [157] clue_0.3-65 rsvd_1.0.5
## [159] broom_1.0.5 xtable_1.8-4
## [161] fANCOVA_0.6-1 later_1.3.1
## [163] viridisLite_0.4.2 tibble_3.2.1
## [165] lmerTest_3.1-3 glmmTMB_1.1.8
## [167] memoise_2.0.1 beeswarm_0.4.0
## [169] AnnotationDbi_1.64.0 cluster_2.1.4
## [171] globals_0.16.2References
Kang, Hyun Min, Meena Subramaniam, Sasha Targ, Michelle Nguyen, Lenka Maliskova, Elizabeth McCarthy, Eunice Wan, et al. 2018. “Multiplexed Droplet Single-Cell RNA-sequencing Using Natural Genetic Variation.” Nature Biotechnology 36 (1): 89–94.
Korthauer, Keegan D, Li-Fang Chu, Michael A Newton, Yuan Li, James Thomson, Ron Stewart, and Christina Kendziorski. 2016. “A Statistical Approach for Identifying Differential Distributions in Single-Cell RNA-seq Experiments.” Genome Biology 17 (1): 222.
Soneson, Charlotte, and Mark D Robinson. 2016. “iCOBRA: Open, Reproducible, Standardized and Live Method Benchmarking.” Nature Methods 13 (4): 283.
———. 2018. “Towards Unified Quality Verification of Synthetic Count Data with countsimQC.” Bioinformatics 34 (4): 691–92.