R/compare.R
In Oshlack/splatter: Simple Simulation of Single-cell RNA Sequencing Data
Defines functions
tidyStatSumm
summariseKS
summariseStats
summariseDiff
makeOverallPanel
makeDiffPanel
makeCompPanel
diffSCEs
compareSCEs
Documented in
compareSCEs
diffSCEs
makeCompPanel
makeDiffPanel
makeOverallPanel
summariseDiff
#' Compare SingleCellExperiment objects
#'
#' Combine the data from several SingleCellExperiment objects and produce some
#' basic plots comparing them.
#'
#' @param sces named list of SingleCellExperiment objects to combine and
#' compare.
#' @param point.size size of points in scatter plots.
#' @param point.alpha opacity of points in scatter plots.
#' @param fits whether to include fits in scatter plots.
#' @param colours vector of colours to use for each dataset.
#'
#' @details
#' The returned list has three items:
#'
#' \describe{
#' \item{\code{RowData}}{Combined row data from the provided
#' SingleCellExperiments.}
#' \item{\code{ColData}}{Combined column data from the provided
#' SingleCellExperiments.}
#' \item{\code{Plots}}{Comparison plots
#' \describe{
#' \item{\code{Means}}{Boxplot of mean distribution.}
#' \item{\code{Variances}}{Boxplot of variance distribution.}
#' \item{\code{MeanVar}}{Scatter plot with fitted lines showing the
#' mean-variance relationship.}
#' \item{\code{LibrarySizes}}{Boxplot of the library size
#' distribution.}
#' \item{\code{ZerosGene}}{Boxplot of the percentage of each gene
#' that is zero.}
#' \item{\code{ZerosCell}}{Boxplot of the percentage of each cell
#' that is zero.}
#' \item{\code{MeanZeros}}{Scatter plot with fitted lines showing
#' the mean-zeros relationship.}
#' \item{\code{VarGeneCor}}{Heatmap of correlation of the 100 most
#' variable genes.}
#' }
#' }
#' }
#'
#' The plots returned by this function are created using
#' \code{\link[ggplot2]{ggplot}} and are only a sample of the kind of plots you
#' might like to consider. The data used to create these plots is also returned
#' and should be in the correct format to allow you to create further plots
#' using \code{\link[ggplot2]{ggplot}}.
#'
#' @return List containing the combined datasets and plots.
#' @examples
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))
#' names(comparison)
#' names(comparison$Plots)
#' @importFrom S4Vectors metadata<- metadata
#' @importFrom SingleCellExperiment cpm<- cpm
#' @importFrom rlang .data
#' @importFrom stats cor
#' @export
compareSCEs <- function(sces, point.size = 0.1, point.alpha = 0.1,
fits = TRUE, colours = NULL) {
checkDependencies(deps = "ggplot2")
checkmate::assertList(
sces,
types = "SingleCellExperiment",
any.missing = FALSE,
min.len = 1,
names = "unique"
)
checkmate::assertNumber(point.size, finite = TRUE)
checkmate::assertNumber(point.alpha, lower = 0, upper = 1)
checkmate::assertLogical(fits, any.missing = FALSE, len = 1)
if (!is.null(colours)) {
checkmate::assertCharacter(
colours,
any.missing = FALSE,
len = length(sces)
)
} else {
colours <- scales::hue_pal()(length(sces))
}
for (name in names(sces)) {
sce <- sces[[name]]
rowData(sce)$Dataset <- name
colData(sce)$Dataset <- name
sce <- scuttle::addPerCellQC(sce)
sce <- scuttle::addPerFeatureQC(sce)
cpm(sce) <- as.matrix(scuttle::calculateCPM(sce))
sce <- addFeatureStats(sce, "counts")
sce <- addFeatureStats(sce, "cpm")
sce <- addFeatureStats(sce, "cpm", log = TRUE)
n.features <- colData(sce)$detected
colData(sce)$PctZero <- 100 * (1 - n.features / nrow(sce))
rowData(sce)$PctZero <- 100 - rowData(sce)$detected
var.genes <- rev(order(rowData(sce)$VarLogCPM))[seq_len(100)]
var.cpm <- log2(cpm(sce)[var.genes, ] + 1)
var.cors <- as.data.frame.table(cor(t(var.cpm), method = "spearman"))
colnames(var.cors) <- c("GeneA", "GeneB", "Correlation")
var.cors$VarGeneA <- rep(paste0("VarGene", seq_len(100)), 100)
var.cors$VarGeneB <- rep(paste0("VarGene", seq_len(100)), each = 100)
var.cors$Dataset <- name
var.cors <- var.cors[, c(
"Dataset", "GeneA", "GeneB", "VarGeneA",
"VarGeneB", "Correlation"
)]
metadata(sce)$VarGeneCorrelation <- var.cors
sces[[name]] <- sce
}
features <- rowData(sces[[1]])
cells <- colData(sces[[1]])
var.cors <- metadata(sces[[1]])$VarGeneCorrelation
if (length(sces) > 1) {
for (name in names(sces)[-1]) {
sce <- sces[[name]]
features <- rbindMatched(features, rowData(sce))
cells <- rbindMatched(cells, colData(sce))
var.cors <- rbindMatched(var.cors, metadata(sce)$VarGeneCorrelation)
}
}
features$Dataset <- factor(features$Dataset, levels = names(sces))
cells$Dataset <- factor(cells$Dataset, levels = names(sces))
var.cors$Dataset <- factor(var.cors$Dataset, levels = names(sces))
features <- data.frame(features)
cells <- data.frame(cells)
means <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$MeanLogCPM, colour = .data$Dataset
)
) +
ggplot2::geom_violin(ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab(expression(paste("Mean ", log[2], "(CPM + 1)"))) +
ggplot2::ggtitle("Distribution of mean expression") +
ggplot2::theme_minimal()
vars <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$VarLogCPM, colour = .data$Dataset
)
) +
ggplot2::geom_violin(ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab(expression(paste("Variance ", log[2], "(CPM + 1)"))) +
ggplot2::ggtitle("Distribution of variance") +
ggplot2::theme_minimal()
mean.var <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$MeanLogCPM, y = .data$VarLogCPM,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab(expression(paste("Mean ", log[2], "(CPM + 1)"))) +
ggplot2::ylab(expression(paste("Variance ", log[2], "(CPM + 1)"))) +
ggplot2::ggtitle("Mean-variance relationship") +
ggplot2::theme_minimal()
libs <- ggplot2::ggplot(
cells,
ggplot2::aes(x = .data$Dataset, y = .data$sum, colour = .data$Dataset)
) +
ggplot2::geom_violin(
ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_y_continuous(labels = scales::comma) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab("Total counts per cell") +
ggplot2::ggtitle("Distribution of library sizes") +
ggplot2::theme_minimal()
z.gene <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_violin(
ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_y_continuous(limits = c(0, 100)) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab("Percentage zeros per gene") +
ggplot2::ggtitle("Distribution of zeros per gene") +
ggplot2::theme_minimal()
z.cell <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$Dataset, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_violin(
ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_y_continuous(limits = c(0, 100)) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab("Percentage zeros per cell") +
ggplot2::ggtitle("Distribution of zeros per cell") +
ggplot2::theme_minimal()
mean.zeros <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$mean, y = .data$PctZero,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_x_log10(labels = scales::comma) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab("Mean count") +
ggplot2::ylab("Percentage zeros") +
ggplot2::ggtitle("Mean-zeros relationship") +
ggplot2::theme_minimal()
var.correlation <- ggplot2::ggplot(
var.cors,
ggplot2::aes(
x = .data$VarGeneA, y = .data$VarGeneB,
fill = .data$Correlation
)
) +
ggplot2::geom_tile() +
ggplot2::scale_fill_distiller(palette = "RdBu", limits = c(-1, 1)) +
ggplot2::coord_fixed() +
ggplot2::facet_wrap(~Dataset) +
ggplot2::ggtitle("Correlation - 100 variable genes") +
ggplot2::theme_minimal() +
ggplot2::theme(
axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank()
)
if (fits) {
mean.var <- mean.var +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
mean.zeros <- mean.zeros +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
}
comparison <- list(
RowData = features,
ColData = cells,
Plots = list(
Means = means,
Variances = vars,
MeanVar = mean.var,
LibrarySizes = libs,
ZerosGene = z.gene,
ZerosCell = z.cell,
MeanZeros = mean.zeros,
VarGeneCor = var.correlation
)
)
return(comparison)
}
#' Diff SingleCellExperiment objects
#'
#' Combine the data from several SingleCellExperiment objects and produce some
#' basic plots comparing them to a reference.
#'
#' @param sces named list of SingleCellExperiment objects to combine and
#' compare.
#' @param ref string giving the name of the SingleCellExperiment to use as the
#' reference
#' @param point.size size of points in scatter plots.
#' @param point.alpha opacity of points in scatter plots.
#' @param fits whether to include fits in scatter plots.
#' @param colours vector of colours to use for each dataset.
#'
#' @details
#'
#' This function aims to look at the differences between a reference
#' SingleCellExperiment and one or more others. It requires each
#' SingleCellExperiment to have the same dimensions. Properties are compared by
#' ranks, for example when comparing the means the values are ordered and the
#' differences between the reference and another dataset plotted. A series of
#' Q-Q plots are also returned.
#'
#' The returned list has five items:
#'
#' \describe{
#' \item{\code{Reference}}{The SingleCellExperiment used as the reference.}
#' \item{\code{RowData}}{Combined feature data from the provided
#' SingleCellExperiments.}
#' \item{\code{ColData}}{Combined column data from the provided
#' SingleCellExperiments.}
#' \item{\code{Plots}}{Difference plots
#' \describe{
#' \item{\code{Means}}{Boxplot of mean differences.}
#' \item{\code{Variances}}{Boxplot of variance differences.}
#' \item{\code{MeanVar}}{Scatter plot showing the difference from
#' the reference variance across expression ranks.}
#' \item{\code{LibraeySizes}}{Boxplot of the library size
#' differences.}
#' \item{\code{ZerosGene}}{Boxplot of the differences in the
#' percentage of each gene that is zero.}
#' \item{\code{ZerosCell}}{Boxplot of the differences in the
#' percentage of each cell that is zero.}
#' \item{\code{MeanZeros}}{Scatter plot showing the difference from
#' the reference percentage of zeros across expression ranks.}
#' }
#' }
#' \item{\code{QQPlots}}{Quantile-Quantile plots
#' \describe{
#' \item{\code{Means}}{Q-Q plot of the means.}
#' \item{\code{Variances}}{Q-Q plot of the variances.}
#' \item{\code{LibrarySizes}}{Q-Q plot of the library sizes.}
#' \item{\code{ZerosGene}}{Q-Q plot of the percentage of zeros per
#' gene.}
#' \item{\code{ZerosCell}}{Q-Q plot of the percentage of zeros per
#' cell.}
#' }
#' }
#' }
#'
#' The plots returned by this function are created using
#' \code{\link[ggplot2]{ggplot}} and are only a sample of the kind of plots you
#' might like to consider. The data used to create these plots is also returned
#' and should be in the correct format to allow you to create further plots
#' using \code{\link[ggplot2]{ggplot}}.
#'
#' @return List containing the combined datasets and plots.
#' @examples
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' names(difference)
#' names(difference$Plots)
#' @importFrom SingleCellExperiment cpm<-
#' @importFrom rlang .data
#' @export
diffSCEs <- function(sces, ref, point.size = 0.1, point.alpha = 0.1,
fits = TRUE, colours = NULL) {
checkDependencies(deps = "ggplot2")
checkmate::assertList(
sces,
types = "SingleCellExperiment",
any.missing = FALSE,
min.len = 2,
names = "unique"
)
checkmate::assertString(ref)
checkmate::assertNumber(point.size, finite = TRUE)
checkmate::assertNumber(point.alpha, lower = 0, upper = 1)
checkmate::assertLogical(fits, any.missing = FALSE, len = 1)
if (!(ref %in% names(sces))) {
stop("'ref' must be the name of a SingleCellExperiment in 'sces'")
} else {
ref.idx <- which(names(sces) == ref)
}
if (!is.null(colours)) {
checkmate::assertCharacter(
colours,
any.missing = FALSE,
len = length(sces) - 1
)
} else {
colours <- scales::hue_pal()(length(sces))
colours <- colours[-ref.idx]
}
ref.dim <- dim(sces[[ref]])
for (name in names(sces)) {
sce <- sces[[name]]
if (!identical(dim(sce), ref.dim)) {
stop("all datasets in 'sces' must have the same dimensions")
}
rowData(sce)$Dataset <- name
colData(sce)$Dataset <- name
sce <- scuttle::addPerCellQC(sce)
sce <- scuttle::addPerFeatureQC(sce)
cpm(sce) <- as.matrix(scuttle::calculateCPM(sce))
sce <- addFeatureStats(sce, "counts")
sce <- addFeatureStats(sce, "cpm", log = TRUE)
n.features <- colData(sce)$detected
colData(sce)$PctZero <- 100 * (1 - n.features / nrow(sce))
rowData(sce)$RankCounts <- rank(rowData(sce)$mean)
rowData(sce)$PctZero <- 100 - rowData(sce)$detected
sces[[name]] <- sce
}
ref.sce <- sces[[ref]]
ref.means <- sort(rowData(ref.sce)$MeanLogCPM)
ref.vars <- sort(rowData(ref.sce)$VarLogCPM)
ref.libs <- sort(colData(ref.sce)$sum)
ref.z.gene <- sort(rowData(ref.sce)$PctZero)
ref.z.cell <- sort(colData(ref.sce)$PctZero)
ref.rank.ord <- order(rowData(ref.sce)$RankCounts)
ref.vars.rank <- rowData(ref.sce)$VarLogCPM[ref.rank.ord]
ref.z.gene.rank <- rowData(ref.sce)$PctZero[ref.rank.ord]
for (name in names(sces)) {
sce <- sces[[name]]
rowData(sce)$RefRankMeanLogCPM <- ref.means[
rank(rowData(sce)$MeanLogCPM)
]
rowData(sce)$RankDiffMeanLogCPM <- rowData(sce)$MeanLogCPM -
rowData(sce)$RefRankMeanLogCPM
rowData(sce)$RefRankVarLogCPM <- ref.vars[rank(rowData(sce)$VarLogCPM)]
rowData(sce)$RankDiffVarLogCPM <- rowData(sce)$VarLogCPM -
rowData(sce)$RefRankVarLogCPM
colData(sce)$RefRankLibSize <- ref.libs[rank(colData(sce)$sum)]
colData(sce)$RankDiffLibSize <- colData(sce)$sum -
colData(sce)$RefRankLibSize
rowData(sce)$RefRankZeros <- ref.z.gene[rank(rowData(sce)$PctZero)]
rowData(sce)$RankDiffZeros <- rowData(sce)$PctZero -
rowData(sce)$RefRankZeros
colData(sce)$RefRankZeros <- ref.z.cell[rank(colData(sce)$PctZero)]
colData(sce)$RankDiffZeros <- colData(sce)$PctZero -
colData(sce)$RefRankZeros
rowData(sce)$MeanRankVarDiff <- rowData(sce)$VarLogCPM -
ref.vars.rank[rowData(sce)$RankCounts]
rowData(sce)$MeanRankZerosDiff <- rowData(sce)$PctZero -
ref.z.gene.rank[rowData(sce)$RankCounts]
sces[[name]] <- sce
}
ref.sce <- sces[[ref]]
sces[[ref]] <- NULL
features <- rowData(sces[[1]])
cells <- colData(sces[[1]])
if (length(sces) > 1) {
for (name in names(sces)[-1]) {
sce <- sces[[name]]
features <- rbindMatched(features, rowData(sce))
cells <- rbindMatched(cells, colData(sce))
}
}
features$Dataset <- factor(features$Dataset, levels = names(sces))
cells$Dataset <- factor(cells$Dataset, levels = names(sces))
features <- data.frame(features)
cells <- data.frame(cells)
means <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffMeanLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(expression(paste(
"Rank difference mean ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Difference in mean expression") +
ggplot2::theme_minimal()
vars <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffVarLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(expression(paste(
"Rank difference variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Difference in variance") +
ggplot2::theme_minimal()
mean.var <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RankCounts, y = .data$MeanRankVarDiff,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab("Expression rank") +
ggplot2::ylab(expression(paste(
"Difference in variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Difference in mean-variance relationship") +
ggplot2::theme_minimal()
libs <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffLibSize,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(paste("Rank difference library size")) +
ggplot2::ggtitle("Difference in library sizes") +
ggplot2::theme_minimal()
z.gene <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffZeros,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(paste("Rank difference percentage zeros")) +
ggplot2::ggtitle("Difference in zeros per gene") +
ggplot2::theme_minimal()
z.cell <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffZeros,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(paste("Rank difference percentage zeros")) +
ggplot2::ggtitle("Difference in zeros per cell") +
ggplot2::theme_minimal()
mean.zeros <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RankCounts, y = .data$MeanRankZerosDiff,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab("Expression rank") +
ggplot2::ylab("Difference in percentage zeros per gene") +
ggplot2::ggtitle("Difference in mean-zeros relationship") +
ggplot2::theme_minimal()
means.qq <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RefRankMeanLogCPM, y = .data$MeanLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab(expression(paste(
"Reference mean ", log[2],
"(CPM + 1)"
))) +
ggplot2::ylab(expression(paste(
"Alternative mean ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Ranked means") +
ggplot2::theme_minimal()
vars.qq <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RefRankVarLogCPM, y = .data$VarLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab(expression(paste(
"Reference variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ylab(expression(paste(
"Alternative variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Ranked variances") +
ggplot2::theme_minimal()
libs.qq <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$RefRankLibSize, y = .data$sum,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab("Reference library size") +
ggplot2::ylab("Alternative library size") +
ggplot2::ggtitle("Ranked library sizes") +
ggplot2::theme_minimal()
z.gene.qq <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RefRankZeros, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab("Reference percentage zeros") +
ggplot2::ylab("Alternative percentage zeros") +
ggplot2::ggtitle("Ranked percentage zeros per gene") +
ggplot2::theme_minimal()
z.cell.qq <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$RefRankZeros, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab("Reference percentage zeros") +
ggplot2::ylab("Alternative percentage zeros") +
ggplot2::ggtitle("Ranked percentage zeros per cell") +
ggplot2::theme_minimal()
if (fits) {
mean.var <- mean.var +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
mean.zeros <- mean.zeros +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
}
comparison <- list(
Reference = ref.sce,
RowData = features,
ColData = cells,
Plots = list(
Means = means,
Variances = vars,
MeanVar = mean.var,
LibrarySizes = libs,
ZerosGene = z.gene,
ZerosCell = z.cell,
MeanZeros = mean.zeros
),
QQPlots = list(
Means = means.qq,
Variances = vars.qq,
LibrarySizes = libs.qq,
ZerosGene = z.gene.qq,
ZerosCell = z.cell.qq
)
)
return(comparison)
}
#' Make comparison panel
#'
#' Combine the plots from \code{compareSCEs} into a single panel.
#'
#' @param comp list returned by \code{\link{compareSCEs}}.
#' @param title title for the panel.
#' @param labels vector of labels for each of the seven plots.
#'
#' @return Combined panel plot
#'
#' @examples
#' \donttest{
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))
#' panel <- makeCompPanel(comparison)
#' }
#'
#' @export
makeCompPanel <- function(comp, title = "Comparison",
labels = c(
"Means", "Variance",
"Mean-variance relationship",
"Library size", "Zeros per gene",
"Zeros per cell",
"Mean-zeros relationship"
)) {
checkDependencies(deps = c("ggplot2", "cowplot"))
checkmate::assertList(comp, any.missing = FALSE, len = 3)
checkmate::checkString(title)
checkmate::checkCharacter(labels, len = 7)
plots <- list(
p1 = comp$Plots$Means,
p2 = comp$Plots$Variances,
p3 = comp$Plots$MeanVar,
p4 = comp$Plots$LibrarySizes,
p5 = comp$Plots$ZerosGene,
p6 = comp$Plots$ZerosCell,
p7 = comp$Plots$MeanZeros
)
# Remove titles and legends
for (plot in names(plots)) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(
legend.position = "none",
plot.title = ggplot2::element_blank()
)
}
# Remove x-axis title from some plots
for (plot in paste0("p", c(1, 2, 4, 5, 6))) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(axis.title.x = ggplot2::element_blank())
}
plots$leg <- cowplot::get_legend(
plots$p3 + ggplot2::theme(legend.position = "bottom")
)
panel <- cowplot::ggdraw() +
cowplot::draw_label(title, 0.5, 0.98, fontface = "bold", size = 18) +
cowplot::draw_label(
labels[1], 0.01, 0.95,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p1, 0.0, 0.74, 0.5, 0.20) +
cowplot::draw_label(
labels[2], 0.51, 0.95,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p2, 0.5, 0.74, 0.5, 0.20) +
cowplot::draw_label(
labels[3], 0.01, 0.70,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p3, 0.0, 0.49, 0.5, 0.20) +
cowplot::draw_label(
labels[4], 0.51, 0.70,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p4, 0.5, 0.49, 0.5, 0.20) +
cowplot::draw_label(
labels[5], 0.01, 0.45,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p5, 0.0, 0.24, 0.5, 0.20) +
cowplot::draw_label(
labels[6], 0.51, 0.45,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p6, 0.5, 0.24, 0.5, 0.20) +
cowplot::draw_label(
labels[7], 0.01, 0.21,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p7, 0.0, 0.00, 0.5, 0.20) +
cowplot::draw_plot(plots$leg, 0.5, 0.00, 0.5, 0.20)
return(panel)
}
#' Make difference panel
#'
#' Combine the plots from \code{diffSCEs} into a single panel.
#'
#' @param diff list returned by \code{\link{diffSCEs}}.
#' @param title title for the panel.
#' @param labels vector of labels for each of the seven sections.
#'
#' @return Combined panel plot
#'
#' @examples
#' \donttest{
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' panel <- makeDiffPanel(difference)
#' }
#'
#' @export
makeDiffPanel <- function(diff, title = "Difference comparison",
labels = c(
"Means", "Variance", "Library size",
"Zeros per cell", "Zeros per gene",
"Mean-variance relationship",
"Mean-zeros relationship"
)) {
checkDependencies(deps = c("ggplot2", "cowplot"))
checkmate::assertList(diff, any.missing = FALSE, len = 5)
checkmate::checkString(title)
checkmate::checkCharacter(labels, len = 7)
plots <- list(
p1 = diff$Plots$Means,
p2 = diff$QQPlots$Means,
p3 = diff$Plots$Variances,
p4 = diff$QQPlots$Variances,
p5 = diff$Plots$MeanVar,
p6 = diff$Plots$LibrarySizes,
p7 = diff$QQPlots$LibrarySizes,
p8 = diff$Plots$ZerosCell,
p9 = diff$QQPlots$ZerosCell,
p10 = diff$Plots$ZerosGene,
p11 = diff$QQPlots$ZerosGene,
p12 = diff$Plots$MeanZeros
)
# Remove titles and legends
for (plot in names(plots)) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(
legend.position = "none",
plot.title = ggplot2::element_blank()
)
}
# Remove x-axis title from some plots
for (plot in paste0("p", c(1, 3, 6, 8, 10))) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(axis.title.x = ggplot2::element_blank())
}
plots$leg <- cowplot::get_legend(
plots$p5 + ggplot2::theme(legend.position = "bottom")
)
panel <- cowplot::ggdraw() +
cowplot::draw_label(
title, 0.5, 0.98,
fontface = "bold", size = 18
) +
cowplot::draw_label(
labels[1], 0.0, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p1, 0.0, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p2, 0.0, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[2], 0.21, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p3, 0.21, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p4, 0.21, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[6], 0.0, 0.30,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p5, 0.0, 0.0, 0.38, 0.29) +
cowplot::draw_label(
labels[3], 0.41, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p6, 0.41, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p7, 0.41, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[4], 0.61, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p8, 0.61, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p9, 0.61, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[7], 0.41, 0.30,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p12, 0.41, 0.0, 0.38, 0.29) +
cowplot::draw_label(
labels[5], 0.81, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p10, 0.81, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p11, 0.81, 0.32, 0.18, 0.29) +
cowplot::draw_plot(plots$leg, 0.81, 0.0, 0.2, 0.29)
return(panel)
}
#' Make overall panel
#'
#' Combine the plots from \code{compSCEs} and \code{diffSCEs} into a
#' single panel.
#'
#' @param comp list returned by \code{\link{compareSCEs}}.
#' @param diff list returned by \code{\link{diffSCEs}}.
#' @param title title for the panel.
#' @param row.labels vector of labels for each of the seven rows.
#'
#' @return Combined panel plot
#'
#' @examples
#' \donttest{
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' panel <- makeOverallPanel(comparison, difference)
#' }
#'
#' @export
makeOverallPanel <- function(comp, diff, title = "Overall comparison",
row.labels = c(
"Means", "Variance",
"Mean-variance relationship",
"Library size", "Zeros per cell",
"Zeros per gene",
"Mean-zeros relationship"
)) {
checkDependencies(deps = c("ggplot2", "cowplot"))
checkmate::assertList(comp, any.missing = FALSE, len = 3)
checkmate::assertList(diff, any.missing = FALSE, len = 5)
checkmate::checkString(title)
checkmate::checkCharacter(row.labels, len = 7)
plots <- list(
p1 = comp$Plots$Means,
p2 = diff$Plots$Means,
p3 = diff$QQPlots$Means,
p4 = comp$Plots$Variances,
p5 = diff$Plots$Variances,
p6 = diff$QQPlots$Variances,
p7 = comp$Plots$MeanVar,
p8 = diff$Plots$MeanVar,
p9 = comp$Plots$LibrarySizes,
p10 = diff$Plots$LibrarySizes,
p11 = diff$QQPlots$LibrarySizes,
p12 = comp$Plots$ZerosCell,
p13 = diff$Plots$ZerosCell,
p14 = diff$QQPlots$ZerosCell,
p15 = comp$Plots$ZerosGene,
p16 = diff$Plots$ZerosGene,
p17 = diff$QQPlots$ZerosGene,
p18 = comp$Plots$MeanZeros,
p19 = diff$Plots$MeanZeros
)
# Remove titles and legends
for (plot in names(plots)) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(
legend.position = "none",
plot.title = ggplot2::element_blank()
)
}
# Remove x-axis title from some plots
for (plot in paste0("p", c(1, 2, 4, 5, 9, 10, 12, 13, 15, 16))) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(axis.title.x = ggplot2::element_blank())
}
plots$leg <- cowplot::get_legend(
plots$p7 + ggplot2::theme(legend.position = "bottom")
)
panel <- cowplot::ggdraw() +
cowplot::draw_label(
title, 0.5, 0.995,
fontface = "bold", size = 18
) +
cowplot::draw_label(
row.labels[1], 0.01, 0.985,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p1, 0.00, 0.86, 0.32, 0.12) +
cowplot::draw_plot(plots$p2, 0.34, 0.86, 0.32, 0.12) +
cowplot::draw_plot(plots$p3, 0.67, 0.86, 0.32, 0.12) +
cowplot::draw_label(
row.labels[2], 0.01, 0.845,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p4, 0.00, 0.72, 0.32, 0.12) +
cowplot::draw_plot(plots$p5, 0.34, 0.72, 0.32, 0.12) +
cowplot::draw_plot(plots$p6, 0.67, 0.72, 0.32, 0.12) +
cowplot::draw_label(
row.labels[3], 0.01, 0.705,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p7, 0.00, 0.58, 0.49, 0.12) +
cowplot::draw_plot(plots$p8, 0.51, 0.58, 0.49, 0.12) +
cowplot::draw_label(
row.labels[4], 0.01, 0.56,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p9, 0.00, 0.44, 0.32, 0.12) +
cowplot::draw_plot(plots$p10, 0.34, 0.44, 0.32, 0.12) +
cowplot::draw_plot(plots$p11, 0.67, 0.44, 0.32, 0.12) +
cowplot::draw_label(
row.labels[5], 0.01, 0.425,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p12, 0.00, 0.30, 0.32, 0.12) +
cowplot::draw_plot(plots$p13, 0.34, 0.30, 0.32, 0.12) +
cowplot::draw_plot(plots$p14, 0.67, 0.30, 0.32, 0.12) +
cowplot::draw_label(
row.labels[6], 0.01, 0.285,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p15, 0.00, 0.16, 0.32, 0.12) +
cowplot::draw_plot(plots$p16, 0.34, 0.16, 0.32, 0.12) +
cowplot::draw_plot(plots$p17, 0.67, 0.16, 0.32, 0.12) +
cowplot::draw_label(
row.labels[7], 0.01, 0.145,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p18, 0.00, 0.02, 0.49, 0.12) +
cowplot::draw_plot(plots$p19, 0.51, 0.02, 0.49, 0.12) +
cowplot::draw_plot(plots$leg, 0.00, 0.00, 1.00, 0.02)
return(panel)
}
#' Summarise diffSCEs
#'
#' Summarise the results of \code{\link{diffSCEs}}. Calculates the Median
#' Absolute Deviation (MAD), Mean Absolute Error (MAE), Root Mean Squared
#' Error (RMSE) and Kolmogorov-Smirnov (KS) statistics for the various
#' properties and ranks them.
#'
#' @param diff Output from \code{\link{diffSCEs}}
#'
#' @return data.frame with MADs, MAEs, RMSEs, scaled statistics and ranks
#' @examples
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' summary <- summariseDiff(difference)
#' head(summary)
#' @export
#' @importFrom SummarizedExperiment rowData
summariseDiff <- function(diff) {
row.stats <- c(
Mean = "RankDiffMeanLogCPM",
Variance = "RankDiffVarLogCPM",
ZerosGene = "RankDiffZeros",
MeanVar = "MeanRankVarDiff",
MeanZeros = "MeanRankZerosDiff"
)
row.ks.stats <- c(
Mean = "MeanLogCPM",
Variance = "VarLogCPM",
ZerosGene = "PctZero",
MeanVar = NA,
MeanZeros = NA
)
row.mad <- summariseStats(diff$RowData, "Dataset", row.stats, "MAD")
row.mae <- summariseStats(diff$RowData, "Dataset", row.stats, "MAE")
row.rmse <- summariseStats(diff$RowData, "Dataset", row.stats, "RMSE")
row.ks <- summariseKS(
diff$RowData,
SummarizedExperiment::rowData(diff$Reference),
"Dataset",
row.ks.stats
)
row.list <- list(row.mad, row.mae, row.rmse, row.ks)
row.list <- lapply(row.list, function(summ) {
summ[, -c(1, 2)]
})
row.summ <- data.frame(
Dataset = row.mad$Dataset,
Statistic = row.mad$Statistic
)
row.list <- c(row.summ, row.list)
row.summ <- do.call("cbind", row.list)
col.stats <- c(
LibSize = "RankDiffLibSize",
ZerosCell = "RankDiffZeros"
)
col.ks.stats <- c(
LibSize = "sum",
ZerosCell = "PctZero"
)
col.mad <- summariseStats(diff$ColData, "Dataset", col.stats, "MAD")
col.mae <- summariseStats(diff$ColData, "Dataset", col.stats, "MAE")
col.rmse <- summariseStats(diff$ColData, "Dataset", col.stats, "RMSE")
col.ks <- summariseKS(
diff$ColData,
SummarizedExperiment::colData(diff$Reference),
"Dataset",
col.ks.stats
)
col.list <- list(col.mad, col.mae, col.rmse, col.ks)
col.list <- lapply(col.list, function(summ) {
summ[, -c(1, 2)]
})
col.summ <- data.frame(
Dataset = col.mad$Dataset,
Statistic = col.mad$Statistic
)
col.list <- c(col.summ, col.list)
col.summ <- do.call("cbind", col.list)
summary <- rbind(row.summ, col.summ)
return(summary)
}
#' Summarise statistics
#'
#' Summarise columns of a data.frame using a single measure.
#'
#' @param data The data.frame to summarise
#' @param split.col Name of the column used to split the dataset
#' @param stat.cols Names of the columns to summarise. If this vector is named
#' those names will be used in the output.
#' @param measure The measure to use for summarisation.
#'
#' @return data.frame with the summarised measure, scaled and ranked
#'
#' @importFrom stats aggregate
#'
#' @noRd
summariseStats <- function(data, split.col, stat.cols,
measure = c("MAD", "MAE", "RMSE")) {
measure <- match.arg(measure)
if (is.null(names(stat.cols))) {
names(stat.cols) <- stat.cols
}
switch(measure,
"MAD" = {
measure_fun <- function(x) {
median(abs(x))
}
},
"MAE" = {
measure_fun <- function(x) {
mean(abs(x))
}
},
"RMSE" = {
measure_fun <- function(x) {
sqrt(mean(abs(x^2)))
}
}
)
summ <- aggregate(
data[, stat.cols],
list(Dataset = data[[split.col]]),
measure_fun
)
colnames(summ) <- c(split.col, names(stat.cols))
tidy.summ <- tidyStatSumm(summ, measure)
return(tidy.summ)
}
#' Summarise KS
#'
#' Summarise columns of a data.frame compared to a reference using the KS test.
#'
#' @param data The data.frame to summarise
#' @param ref The reference data.frame
#' @param split.col Name of the column used to split the dataset
#' @param stat.cols Names of the columns to summarise. If this vector is named
#' those names will be used in the output.
#'
#' @return data.frame with the summarised measure, scaled and ranked
#' @importFrom stats ks.test
#'
#' @noRd
summariseKS <- function(data, ref, split.col, stat.cols) {
if (is.null(names(stat.cols))) {
names(stat.cols) <- stat.cols
}
splits <- unique(data[[split.col]])
summ <- expand.grid(
Dataset = splits,
Statistic = names(stat.cols),
stringsAsFactors = FALSE
)
ks.res <- mapply(function(split, stat.name) {
stat <- stat.cols[stat.name]
if (!is.na(stat)) {
data.stat <- data[data[[split.col]] == split, stat]
ref.stat <- ref[[stat]]
ks <- suppressWarnings(ks.test(ref.stat, data.stat))
ks.out <- c(KS = unname(ks$statistic), KSPVal = ks$p.value)
} else {
ks.out <- c(KS = NA, KSPVal = NA)
}
return(ks.out)
}, summ$Dataset, summ$Statistic)
summ$KS <- ks.res["KS", ]
summ$KSPVal <- ks.res["KSPVal", ]
ks.ranks <- lapply(split(summ, summ$Statistic), function(x) {
rank(x$KS)
})
ks.ranks <- unlist(ks.ranks)
summ$KSRank <- ks.ranks
summ$KSRank[is.na(summ$KS)] <- NA
return(summ)
}
#' Tidy summarised statistics
#'
#' Convert a statistic summary to tidy format and add ranks and scaled values
#'
#' @param stat.summ The summary to convert
#' @param measure The name of the summarisation measure
#'
#' @return tidy data.frame with the summarised measure, scaled and ranked
#'
#' @noRd
tidyStatSumm <- function(stat.summ, measure = c("MAD", "MAE", "RMSE")) {
measure <- match.arg(measure)
summ.mat <- t(stat.summ[, -1])
colnames(summ.mat) <- stat.summ[, 1]
scale.summ <- apply(summ.mat, 1, scale)
# Check if apply has returned a vector
if (is.vector(scale.summ)) {
scale.summ <- t(as.matrix(scale.summ))
}
rank.summ <- apply(summ.mat, 1, rank)
if (is.vector(rank.summ)) {
rank.summ <- t(as.matrix(rank.summ))
}
tidy.summ <- as.data.frame.table(t(summ.mat))
colnames(tidy.summ) <- c("Dataset", "Statistic", measure)
tidy.scale <- as.data.frame.table(scale.summ)
tidy.rank <- as.data.frame.table(rank.summ)
tidy.summ[[paste0(measure, "Scaled")]] <- tidy.scale[, 3]
tidy.summ[[paste0(measure, "Rank")]] <- tidy.rank[, 3]
return(tidy.summ)
}
Oshlack/splatter documentation built on May 8, 2024, 5:33 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions tidyStatSumm summariseKS summariseStats summariseDiff makeOverallPanel makeDiffPanel makeCompPanel diffSCEs compareSCEs
Documented in compareSCEs diffSCEs makeCompPanel makeDiffPanel makeOverallPanel summariseDiff
#' Compare SingleCellExperiment objects
#'
#' Combine the data from several SingleCellExperiment objects and produce some
#' basic plots comparing them.
#'
#' @param sces named list of SingleCellExperiment objects to combine and
#' compare.
#' @param point.size size of points in scatter plots.
#' @param point.alpha opacity of points in scatter plots.
#' @param fits whether to include fits in scatter plots.
#' @param colours vector of colours to use for each dataset.
#'
#' @details
#' The returned list has three items:
#'
#' \describe{
#' \item{\code{RowData}}{Combined row data from the provided
#' SingleCellExperiments.}
#' \item{\code{ColData}}{Combined column data from the provided
#' SingleCellExperiments.}
#' \item{\code{Plots}}{Comparison plots
#' \describe{
#' \item{\code{Means}}{Boxplot of mean distribution.}
#' \item{\code{Variances}}{Boxplot of variance distribution.}
#' \item{\code{MeanVar}}{Scatter plot with fitted lines showing the
#' mean-variance relationship.}
#' \item{\code{LibrarySizes}}{Boxplot of the library size
#' distribution.}
#' \item{\code{ZerosGene}}{Boxplot of the percentage of each gene
#' that is zero.}
#' \item{\code{ZerosCell}}{Boxplot of the percentage of each cell
#' that is zero.}
#' \item{\code{MeanZeros}}{Scatter plot with fitted lines showing
#' the mean-zeros relationship.}
#' \item{\code{VarGeneCor}}{Heatmap of correlation of the 100 most
#' variable genes.}
#' }
#' }
#' }
#'
#' The plots returned by this function are created using
#' \code{\link[ggplot2]{ggplot}} and are only a sample of the kind of plots you
#' might like to consider. The data used to create these plots is also returned
#' and should be in the correct format to allow you to create further plots
#' using \code{\link[ggplot2]{ggplot}}.
#'
#' @return List containing the combined datasets and plots.
#' @examples
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))
#' names(comparison)
#' names(comparison$Plots)
#' @importFrom S4Vectors metadata<- metadata
#' @importFrom SingleCellExperiment cpm<- cpm
#' @importFrom rlang .data
#' @importFrom stats cor
#' @export
compareSCEs <- function(sces, point.size = 0.1, point.alpha = 0.1,
fits = TRUE, colours = NULL) {
checkDependencies(deps = "ggplot2")
checkmate::assertList(
sces,
types = "SingleCellExperiment",
any.missing = FALSE,
min.len = 1,
names = "unique"
)
checkmate::assertNumber(point.size, finite = TRUE)
checkmate::assertNumber(point.alpha, lower = 0, upper = 1)
checkmate::assertLogical(fits, any.missing = FALSE, len = 1)
if (!is.null(colours)) {
checkmate::assertCharacter(
colours,
any.missing = FALSE,
len = length(sces)
)
} else {
colours <- scales::hue_pal()(length(sces))
}
for (name in names(sces)) {
sce <- sces[[name]]
rowData(sce)$Dataset <- name
colData(sce)$Dataset <- name
sce <- scuttle::addPerCellQC(sce)
sce <- scuttle::addPerFeatureQC(sce)
cpm(sce) <- as.matrix(scuttle::calculateCPM(sce))
sce <- addFeatureStats(sce, "counts")
sce <- addFeatureStats(sce, "cpm")
sce <- addFeatureStats(sce, "cpm", log = TRUE)
n.features <- colData(sce)$detected
colData(sce)$PctZero <- 100 * (1 - n.features / nrow(sce))
rowData(sce)$PctZero <- 100 - rowData(sce)$detected
var.genes <- rev(order(rowData(sce)$VarLogCPM))[seq_len(100)]
var.cpm <- log2(cpm(sce)[var.genes, ] + 1)
var.cors <- as.data.frame.table(cor(t(var.cpm), method = "spearman"))
colnames(var.cors) <- c("GeneA", "GeneB", "Correlation")
var.cors$VarGeneA <- rep(paste0("VarGene", seq_len(100)), 100)
var.cors$VarGeneB <- rep(paste0("VarGene", seq_len(100)), each = 100)
var.cors$Dataset <- name
var.cors <- var.cors[, c(
"Dataset", "GeneA", "GeneB", "VarGeneA",
"VarGeneB", "Correlation"
)]
metadata(sce)$VarGeneCorrelation <- var.cors
sces[[name]] <- sce
}
features <- rowData(sces[[1]])
cells <- colData(sces[[1]])
var.cors <- metadata(sces[[1]])$VarGeneCorrelation
if (length(sces) > 1) {
for (name in names(sces)[-1]) {
sce <- sces[[name]]
features <- rbindMatched(features, rowData(sce))
cells <- rbindMatched(cells, colData(sce))
var.cors <- rbindMatched(var.cors, metadata(sce)$VarGeneCorrelation)
}
}
features$Dataset <- factor(features$Dataset, levels = names(sces))
cells$Dataset <- factor(cells$Dataset, levels = names(sces))
var.cors$Dataset <- factor(var.cors$Dataset, levels = names(sces))
features <- data.frame(features)
cells <- data.frame(cells)
means <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$MeanLogCPM, colour = .data$Dataset
)
) +
ggplot2::geom_violin(ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab(expression(paste("Mean ", log[2], "(CPM + 1)"))) +
ggplot2::ggtitle("Distribution of mean expression") +
ggplot2::theme_minimal()
vars <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$VarLogCPM, colour = .data$Dataset
)
) +
ggplot2::geom_violin(ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab(expression(paste("Variance ", log[2], "(CPM + 1)"))) +
ggplot2::ggtitle("Distribution of variance") +
ggplot2::theme_minimal()
mean.var <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$MeanLogCPM, y = .data$VarLogCPM,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab(expression(paste("Mean ", log[2], "(CPM + 1)"))) +
ggplot2::ylab(expression(paste("Variance ", log[2], "(CPM + 1)"))) +
ggplot2::ggtitle("Mean-variance relationship") +
ggplot2::theme_minimal()
libs <- ggplot2::ggplot(
cells,
ggplot2::aes(x = .data$Dataset, y = .data$sum, colour = .data$Dataset)
) +
ggplot2::geom_violin(
ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_y_continuous(labels = scales::comma) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab("Total counts per cell") +
ggplot2::ggtitle("Distribution of library sizes") +
ggplot2::theme_minimal()
z.gene <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_violin(
ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_y_continuous(limits = c(0, 100)) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab("Percentage zeros per gene") +
ggplot2::ggtitle("Distribution of zeros per gene") +
ggplot2::theme_minimal()
z.cell <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$Dataset, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_violin(
ggplot2::aes(fill = .data$Dataset),
draw_quantiles = c(0.25, 0.5, 0.75),
colour = "white", alpha = 0.3, linewidth = 0.8
) +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_y_continuous(limits = c(0, 100)) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::ylab("Percentage zeros per cell") +
ggplot2::ggtitle("Distribution of zeros per cell") +
ggplot2::theme_minimal()
mean.zeros <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$mean, y = .data$PctZero,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_x_log10(labels = scales::comma) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab("Mean count") +
ggplot2::ylab("Percentage zeros") +
ggplot2::ggtitle("Mean-zeros relationship") +
ggplot2::theme_minimal()
var.correlation <- ggplot2::ggplot(
var.cors,
ggplot2::aes(
x = .data$VarGeneA, y = .data$VarGeneB,
fill = .data$Correlation
)
) +
ggplot2::geom_tile() +
ggplot2::scale_fill_distiller(palette = "RdBu", limits = c(-1, 1)) +
ggplot2::coord_fixed() +
ggplot2::facet_wrap(~Dataset) +
ggplot2::ggtitle("Correlation - 100 variable genes") +
ggplot2::theme_minimal() +
ggplot2::theme(
axis.title = ggplot2::element_blank(),
axis.text = ggplot2::element_blank()
)
if (fits) {
mean.var <- mean.var +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
mean.zeros <- mean.zeros +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
}
comparison <- list(
RowData = features,
ColData = cells,
Plots = list(
Means = means,
Variances = vars,
MeanVar = mean.var,
LibrarySizes = libs,
ZerosGene = z.gene,
ZerosCell = z.cell,
MeanZeros = mean.zeros,
VarGeneCor = var.correlation
)
)
return(comparison)
}
#' Diff SingleCellExperiment objects
#'
#' Combine the data from several SingleCellExperiment objects and produce some
#' basic plots comparing them to a reference.
#'
#' @param sces named list of SingleCellExperiment objects to combine and
#' compare.
#' @param ref string giving the name of the SingleCellExperiment to use as the
#' reference
#' @param point.size size of points in scatter plots.
#' @param point.alpha opacity of points in scatter plots.
#' @param fits whether to include fits in scatter plots.
#' @param colours vector of colours to use for each dataset.
#'
#' @details
#'
#' This function aims to look at the differences between a reference
#' SingleCellExperiment and one or more others. It requires each
#' SingleCellExperiment to have the same dimensions. Properties are compared by
#' ranks, for example when comparing the means the values are ordered and the
#' differences between the reference and another dataset plotted. A series of
#' Q-Q plots are also returned.
#'
#' The returned list has five items:
#'
#' \describe{
#' \item{\code{Reference}}{The SingleCellExperiment used as the reference.}
#' \item{\code{RowData}}{Combined feature data from the provided
#' SingleCellExperiments.}
#' \item{\code{ColData}}{Combined column data from the provided
#' SingleCellExperiments.}
#' \item{\code{Plots}}{Difference plots
#' \describe{
#' \item{\code{Means}}{Boxplot of mean differences.}
#' \item{\code{Variances}}{Boxplot of variance differences.}
#' \item{\code{MeanVar}}{Scatter plot showing the difference from
#' the reference variance across expression ranks.}
#' \item{\code{LibraeySizes}}{Boxplot of the library size
#' differences.}
#' \item{\code{ZerosGene}}{Boxplot of the differences in the
#' percentage of each gene that is zero.}
#' \item{\code{ZerosCell}}{Boxplot of the differences in the
#' percentage of each cell that is zero.}
#' \item{\code{MeanZeros}}{Scatter plot showing the difference from
#' the reference percentage of zeros across expression ranks.}
#' }
#' }
#' \item{\code{QQPlots}}{Quantile-Quantile plots
#' \describe{
#' \item{\code{Means}}{Q-Q plot of the means.}
#' \item{\code{Variances}}{Q-Q plot of the variances.}
#' \item{\code{LibrarySizes}}{Q-Q plot of the library sizes.}
#' \item{\code{ZerosGene}}{Q-Q plot of the percentage of zeros per
#' gene.}
#' \item{\code{ZerosCell}}{Q-Q plot of the percentage of zeros per
#' cell.}
#' }
#' }
#' }
#'
#' The plots returned by this function are created using
#' \code{\link[ggplot2]{ggplot}} and are only a sample of the kind of plots you
#' might like to consider. The data used to create these plots is also returned
#' and should be in the correct format to allow you to create further plots
#' using \code{\link[ggplot2]{ggplot}}.
#'
#' @return List containing the combined datasets and plots.
#' @examples
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' names(difference)
#' names(difference$Plots)
#' @importFrom SingleCellExperiment cpm<-
#' @importFrom rlang .data
#' @export
diffSCEs <- function(sces, ref, point.size = 0.1, point.alpha = 0.1,
fits = TRUE, colours = NULL) {
checkDependencies(deps = "ggplot2")
checkmate::assertList(
sces,
types = "SingleCellExperiment",
any.missing = FALSE,
min.len = 2,
names = "unique"
)
checkmate::assertString(ref)
checkmate::assertNumber(point.size, finite = TRUE)
checkmate::assertNumber(point.alpha, lower = 0, upper = 1)
checkmate::assertLogical(fits, any.missing = FALSE, len = 1)
if (!(ref %in% names(sces))) {
stop("'ref' must be the name of a SingleCellExperiment in 'sces'")
} else {
ref.idx <- which(names(sces) == ref)
}
if (!is.null(colours)) {
checkmate::assertCharacter(
colours,
any.missing = FALSE,
len = length(sces) - 1
)
} else {
colours <- scales::hue_pal()(length(sces))
colours <- colours[-ref.idx]
}
ref.dim <- dim(sces[[ref]])
for (name in names(sces)) {
sce <- sces[[name]]
if (!identical(dim(sce), ref.dim)) {
stop("all datasets in 'sces' must have the same dimensions")
}
rowData(sce)$Dataset <- name
colData(sce)$Dataset <- name
sce <- scuttle::addPerCellQC(sce)
sce <- scuttle::addPerFeatureQC(sce)
cpm(sce) <- as.matrix(scuttle::calculateCPM(sce))
sce <- addFeatureStats(sce, "counts")
sce <- addFeatureStats(sce, "cpm", log = TRUE)
n.features <- colData(sce)$detected
colData(sce)$PctZero <- 100 * (1 - n.features / nrow(sce))
rowData(sce)$RankCounts <- rank(rowData(sce)$mean)
rowData(sce)$PctZero <- 100 - rowData(sce)$detected
sces[[name]] <- sce
}
ref.sce <- sces[[ref]]
ref.means <- sort(rowData(ref.sce)$MeanLogCPM)
ref.vars <- sort(rowData(ref.sce)$VarLogCPM)
ref.libs <- sort(colData(ref.sce)$sum)
ref.z.gene <- sort(rowData(ref.sce)$PctZero)
ref.z.cell <- sort(colData(ref.sce)$PctZero)
ref.rank.ord <- order(rowData(ref.sce)$RankCounts)
ref.vars.rank <- rowData(ref.sce)$VarLogCPM[ref.rank.ord]
ref.z.gene.rank <- rowData(ref.sce)$PctZero[ref.rank.ord]
for (name in names(sces)) {
sce <- sces[[name]]
rowData(sce)$RefRankMeanLogCPM <- ref.means[
rank(rowData(sce)$MeanLogCPM)
]
rowData(sce)$RankDiffMeanLogCPM <- rowData(sce)$MeanLogCPM -
rowData(sce)$RefRankMeanLogCPM
rowData(sce)$RefRankVarLogCPM <- ref.vars[rank(rowData(sce)$VarLogCPM)]
rowData(sce)$RankDiffVarLogCPM <- rowData(sce)$VarLogCPM -
rowData(sce)$RefRankVarLogCPM
colData(sce)$RefRankLibSize <- ref.libs[rank(colData(sce)$sum)]
colData(sce)$RankDiffLibSize <- colData(sce)$sum -
colData(sce)$RefRankLibSize
rowData(sce)$RefRankZeros <- ref.z.gene[rank(rowData(sce)$PctZero)]
rowData(sce)$RankDiffZeros <- rowData(sce)$PctZero -
rowData(sce)$RefRankZeros
colData(sce)$RefRankZeros <- ref.z.cell[rank(colData(sce)$PctZero)]
colData(sce)$RankDiffZeros <- colData(sce)$PctZero -
colData(sce)$RefRankZeros
rowData(sce)$MeanRankVarDiff <- rowData(sce)$VarLogCPM -
ref.vars.rank[rowData(sce)$RankCounts]
rowData(sce)$MeanRankZerosDiff <- rowData(sce)$PctZero -
ref.z.gene.rank[rowData(sce)$RankCounts]
sces[[name]] <- sce
}
ref.sce <- sces[[ref]]
sces[[ref]] <- NULL
features <- rowData(sces[[1]])
cells <- colData(sces[[1]])
if (length(sces) > 1) {
for (name in names(sces)[-1]) {
sce <- sces[[name]]
features <- rbindMatched(features, rowData(sce))
cells <- rbindMatched(cells, colData(sce))
}
}
features$Dataset <- factor(features$Dataset, levels = names(sces))
cells$Dataset <- factor(cells$Dataset, levels = names(sces))
features <- data.frame(features)
cells <- data.frame(cells)
means <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffMeanLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(expression(paste(
"Rank difference mean ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Difference in mean expression") +
ggplot2::theme_minimal()
vars <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffVarLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(expression(paste(
"Rank difference variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Difference in variance") +
ggplot2::theme_minimal()
mean.var <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RankCounts, y = .data$MeanRankVarDiff,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab("Expression rank") +
ggplot2::ylab(expression(paste(
"Difference in variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Difference in mean-variance relationship") +
ggplot2::theme_minimal()
libs <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffLibSize,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(paste("Rank difference library size")) +
ggplot2::ggtitle("Difference in library sizes") +
ggplot2::theme_minimal()
z.gene <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffZeros,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(paste("Rank difference percentage zeros")) +
ggplot2::ggtitle("Difference in zeros per gene") +
ggplot2::theme_minimal()
z.cell <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$Dataset, y = .data$RankDiffZeros,
colour = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_boxplot(notch = TRUE, width = 0.1, size = 0.8) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::ylab(paste("Rank difference percentage zeros")) +
ggplot2::ggtitle("Difference in zeros per cell") +
ggplot2::theme_minimal()
mean.zeros <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RankCounts, y = .data$MeanRankZerosDiff,
colour = .data$Dataset, fill = .data$Dataset
)
) +
ggplot2::geom_hline(yintercept = 0, colour = "red") +
ggplot2::geom_point(size = point.size, alpha = point.alpha) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::scale_fill_manual(values = colours) +
ggplot2::xlab("Expression rank") +
ggplot2::ylab("Difference in percentage zeros per gene") +
ggplot2::ggtitle("Difference in mean-zeros relationship") +
ggplot2::theme_minimal()
means.qq <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RefRankMeanLogCPM, y = .data$MeanLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab(expression(paste(
"Reference mean ", log[2],
"(CPM + 1)"
))) +
ggplot2::ylab(expression(paste(
"Alternative mean ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Ranked means") +
ggplot2::theme_minimal()
vars.qq <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RefRankVarLogCPM, y = .data$VarLogCPM,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab(expression(paste(
"Reference variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ylab(expression(paste(
"Alternative variance ", log[2],
"(CPM + 1)"
))) +
ggplot2::ggtitle("Ranked variances") +
ggplot2::theme_minimal()
libs.qq <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$RefRankLibSize, y = .data$sum,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab("Reference library size") +
ggplot2::ylab("Alternative library size") +
ggplot2::ggtitle("Ranked library sizes") +
ggplot2::theme_minimal()
z.gene.qq <- ggplot2::ggplot(
features,
ggplot2::aes(
x = .data$RefRankZeros, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab("Reference percentage zeros") +
ggplot2::ylab("Alternative percentage zeros") +
ggplot2::ggtitle("Ranked percentage zeros per gene") +
ggplot2::theme_minimal()
z.cell.qq <- ggplot2::ggplot(
cells,
ggplot2::aes(
x = .data$RefRankZeros, y = .data$PctZero,
colour = .data$Dataset
)
) +
ggplot2::geom_abline(intercept = 0, slope = 1, colour = "red") +
ggplot2::geom_point(size = point.size) +
ggplot2::scale_colour_manual(values = colours) +
ggplot2::xlab("Reference percentage zeros") +
ggplot2::ylab("Alternative percentage zeros") +
ggplot2::ggtitle("Ranked percentage zeros per cell") +
ggplot2::theme_minimal()
if (fits) {
mean.var <- mean.var +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
mean.zeros <- mean.zeros +
ggplot2::geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"))
}
comparison <- list(
Reference = ref.sce,
RowData = features,
ColData = cells,
Plots = list(
Means = means,
Variances = vars,
MeanVar = mean.var,
LibrarySizes = libs,
ZerosGene = z.gene,
ZerosCell = z.cell,
MeanZeros = mean.zeros
),
QQPlots = list(
Means = means.qq,
Variances = vars.qq,
LibrarySizes = libs.qq,
ZerosGene = z.gene.qq,
ZerosCell = z.cell.qq
)
)
return(comparison)
}
#' Make comparison panel
#'
#' Combine the plots from \code{compareSCEs} into a single panel.
#'
#' @param comp list returned by \code{\link{compareSCEs}}.
#' @param title title for the panel.
#' @param labels vector of labels for each of the seven plots.
#'
#' @return Combined panel plot
#'
#' @examples
#' \donttest{
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))
#' panel <- makeCompPanel(comparison)
#' }
#'
#' @export
makeCompPanel <- function(comp, title = "Comparison",
labels = c(
"Means", "Variance",
"Mean-variance relationship",
"Library size", "Zeros per gene",
"Zeros per cell",
"Mean-zeros relationship"
)) {
checkDependencies(deps = c("ggplot2", "cowplot"))
checkmate::assertList(comp, any.missing = FALSE, len = 3)
checkmate::checkString(title)
checkmate::checkCharacter(labels, len = 7)
plots <- list(
p1 = comp$Plots$Means,
p2 = comp$Plots$Variances,
p3 = comp$Plots$MeanVar,
p4 = comp$Plots$LibrarySizes,
p5 = comp$Plots$ZerosGene,
p6 = comp$Plots$ZerosCell,
p7 = comp$Plots$MeanZeros
)
# Remove titles and legends
for (plot in names(plots)) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(
legend.position = "none",
plot.title = ggplot2::element_blank()
)
}
# Remove x-axis title from some plots
for (plot in paste0("p", c(1, 2, 4, 5, 6))) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(axis.title.x = ggplot2::element_blank())
}
plots$leg <- cowplot::get_legend(
plots$p3 + ggplot2::theme(legend.position = "bottom")
)
panel <- cowplot::ggdraw() +
cowplot::draw_label(title, 0.5, 0.98, fontface = "bold", size = 18) +
cowplot::draw_label(
labels[1], 0.01, 0.95,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p1, 0.0, 0.74, 0.5, 0.20) +
cowplot::draw_label(
labels[2], 0.51, 0.95,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p2, 0.5, 0.74, 0.5, 0.20) +
cowplot::draw_label(
labels[3], 0.01, 0.70,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p3, 0.0, 0.49, 0.5, 0.20) +
cowplot::draw_label(
labels[4], 0.51, 0.70,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p4, 0.5, 0.49, 0.5, 0.20) +
cowplot::draw_label(
labels[5], 0.01, 0.45,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p5, 0.0, 0.24, 0.5, 0.20) +
cowplot::draw_label(
labels[6], 0.51, 0.45,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p6, 0.5, 0.24, 0.5, 0.20) +
cowplot::draw_label(
labels[7], 0.01, 0.21,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p7, 0.0, 0.00, 0.5, 0.20) +
cowplot::draw_plot(plots$leg, 0.5, 0.00, 0.5, 0.20)
return(panel)
}
#' Make difference panel
#'
#' Combine the plots from \code{diffSCEs} into a single panel.
#'
#' @param diff list returned by \code{\link{diffSCEs}}.
#' @param title title for the panel.
#' @param labels vector of labels for each of the seven sections.
#'
#' @return Combined panel plot
#'
#' @examples
#' \donttest{
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' panel <- makeDiffPanel(difference)
#' }
#'
#' @export
makeDiffPanel <- function(diff, title = "Difference comparison",
labels = c(
"Means", "Variance", "Library size",
"Zeros per cell", "Zeros per gene",
"Mean-variance relationship",
"Mean-zeros relationship"
)) {
checkDependencies(deps = c("ggplot2", "cowplot"))
checkmate::assertList(diff, any.missing = FALSE, len = 5)
checkmate::checkString(title)
checkmate::checkCharacter(labels, len = 7)
plots <- list(
p1 = diff$Plots$Means,
p2 = diff$QQPlots$Means,
p3 = diff$Plots$Variances,
p4 = diff$QQPlots$Variances,
p5 = diff$Plots$MeanVar,
p6 = diff$Plots$LibrarySizes,
p7 = diff$QQPlots$LibrarySizes,
p8 = diff$Plots$ZerosCell,
p9 = diff$QQPlots$ZerosCell,
p10 = diff$Plots$ZerosGene,
p11 = diff$QQPlots$ZerosGene,
p12 = diff$Plots$MeanZeros
)
# Remove titles and legends
for (plot in names(plots)) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(
legend.position = "none",
plot.title = ggplot2::element_blank()
)
}
# Remove x-axis title from some plots
for (plot in paste0("p", c(1, 3, 6, 8, 10))) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(axis.title.x = ggplot2::element_blank())
}
plots$leg <- cowplot::get_legend(
plots$p5 + ggplot2::theme(legend.position = "bottom")
)
panel <- cowplot::ggdraw() +
cowplot::draw_label(
title, 0.5, 0.98,
fontface = "bold", size = 18
) +
cowplot::draw_label(
labels[1], 0.0, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p1, 0.0, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p2, 0.0, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[2], 0.21, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p3, 0.21, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p4, 0.21, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[6], 0.0, 0.30,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p5, 0.0, 0.0, 0.38, 0.29) +
cowplot::draw_label(
labels[3], 0.41, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p6, 0.41, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p7, 0.41, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[4], 0.61, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p8, 0.61, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p9, 0.61, 0.32, 0.18, 0.29) +
cowplot::draw_label(
labels[7], 0.41, 0.30,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p12, 0.41, 0.0, 0.38, 0.29) +
cowplot::draw_label(
labels[5], 0.81, 0.94,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p10, 0.81, 0.64, 0.18, 0.29) +
cowplot::draw_plot(plots$p11, 0.81, 0.32, 0.18, 0.29) +
cowplot::draw_plot(plots$leg, 0.81, 0.0, 0.2, 0.29)
return(panel)
}
#' Make overall panel
#'
#' Combine the plots from \code{compSCEs} and \code{diffSCEs} into a
#' single panel.
#'
#' @param comp list returned by \code{\link{compareSCEs}}.
#' @param diff list returned by \code{\link{diffSCEs}}.
#' @param title title for the panel.
#' @param row.labels vector of labels for each of the seven rows.
#'
#' @return Combined panel plot
#'
#' @examples
#' \donttest{
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' comparison <- compareSCEs(list(Splat = sim1, Simple = sim2))
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' panel <- makeOverallPanel(comparison, difference)
#' }
#'
#' @export
makeOverallPanel <- function(comp, diff, title = "Overall comparison",
row.labels = c(
"Means", "Variance",
"Mean-variance relationship",
"Library size", "Zeros per cell",
"Zeros per gene",
"Mean-zeros relationship"
)) {
checkDependencies(deps = c("ggplot2", "cowplot"))
checkmate::assertList(comp, any.missing = FALSE, len = 3)
checkmate::assertList(diff, any.missing = FALSE, len = 5)
checkmate::checkString(title)
checkmate::checkCharacter(row.labels, len = 7)
plots <- list(
p1 = comp$Plots$Means,
p2 = diff$Plots$Means,
p3 = diff$QQPlots$Means,
p4 = comp$Plots$Variances,
p5 = diff$Plots$Variances,
p6 = diff$QQPlots$Variances,
p7 = comp$Plots$MeanVar,
p8 = diff$Plots$MeanVar,
p9 = comp$Plots$LibrarySizes,
p10 = diff$Plots$LibrarySizes,
p11 = diff$QQPlots$LibrarySizes,
p12 = comp$Plots$ZerosCell,
p13 = diff$Plots$ZerosCell,
p14 = diff$QQPlots$ZerosCell,
p15 = comp$Plots$ZerosGene,
p16 = diff$Plots$ZerosGene,
p17 = diff$QQPlots$ZerosGene,
p18 = comp$Plots$MeanZeros,
p19 = diff$Plots$MeanZeros
)
# Remove titles and legends
for (plot in names(plots)) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(
legend.position = "none",
plot.title = ggplot2::element_blank()
)
}
# Remove x-axis title from some plots
for (plot in paste0("p", c(1, 2, 4, 5, 9, 10, 12, 13, 15, 16))) {
plots[[plot]] <- plots[[plot]] +
ggplot2::theme(axis.title.x = ggplot2::element_blank())
}
plots$leg <- cowplot::get_legend(
plots$p7 + ggplot2::theme(legend.position = "bottom")
)
panel <- cowplot::ggdraw() +
cowplot::draw_label(
title, 0.5, 0.995,
fontface = "bold", size = 18
) +
cowplot::draw_label(
row.labels[1], 0.01, 0.985,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p1, 0.00, 0.86, 0.32, 0.12) +
cowplot::draw_plot(plots$p2, 0.34, 0.86, 0.32, 0.12) +
cowplot::draw_plot(plots$p3, 0.67, 0.86, 0.32, 0.12) +
cowplot::draw_label(
row.labels[2], 0.01, 0.845,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p4, 0.00, 0.72, 0.32, 0.12) +
cowplot::draw_plot(plots$p5, 0.34, 0.72, 0.32, 0.12) +
cowplot::draw_plot(plots$p6, 0.67, 0.72, 0.32, 0.12) +
cowplot::draw_label(
row.labels[3], 0.01, 0.705,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p7, 0.00, 0.58, 0.49, 0.12) +
cowplot::draw_plot(plots$p8, 0.51, 0.58, 0.49, 0.12) +
cowplot::draw_label(
row.labels[4], 0.01, 0.56,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p9, 0.00, 0.44, 0.32, 0.12) +
cowplot::draw_plot(plots$p10, 0.34, 0.44, 0.32, 0.12) +
cowplot::draw_plot(plots$p11, 0.67, 0.44, 0.32, 0.12) +
cowplot::draw_label(
row.labels[5], 0.01, 0.425,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p12, 0.00, 0.30, 0.32, 0.12) +
cowplot::draw_plot(plots$p13, 0.34, 0.30, 0.32, 0.12) +
cowplot::draw_plot(plots$p14, 0.67, 0.30, 0.32, 0.12) +
cowplot::draw_label(
row.labels[6], 0.01, 0.285,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p15, 0.00, 0.16, 0.32, 0.12) +
cowplot::draw_plot(plots$p16, 0.34, 0.16, 0.32, 0.12) +
cowplot::draw_plot(plots$p17, 0.67, 0.16, 0.32, 0.12) +
cowplot::draw_label(
row.labels[7], 0.01, 0.145,
fontface = "bold", hjust = 0, vjust = 0
) +
cowplot::draw_plot(plots$p18, 0.00, 0.02, 0.49, 0.12) +
cowplot::draw_plot(plots$p19, 0.51, 0.02, 0.49, 0.12) +
cowplot::draw_plot(plots$leg, 0.00, 0.00, 1.00, 0.02)
return(panel)
}
#' Summarise diffSCEs
#'
#' Summarise the results of \code{\link{diffSCEs}}. Calculates the Median
#' Absolute Deviation (MAD), Mean Absolute Error (MAE), Root Mean Squared
#' Error (RMSE) and Kolmogorov-Smirnov (KS) statistics for the various
#' properties and ranks them.
#'
#' @param diff Output from \code{\link{diffSCEs}}
#'
#' @return data.frame with MADs, MAEs, RMSEs, scaled statistics and ranks
#' @examples
#' sim1 <- splatSimulate(nGenes = 1000, batchCells = 20)
#' sim2 <- simpleSimulate(nGenes = 1000, nCells = 20)
#' difference <- diffSCEs(list(Splat = sim1, Simple = sim2), ref = "Simple")
#' summary <- summariseDiff(difference)
#' head(summary)
#' @export
#' @importFrom SummarizedExperiment rowData
summariseDiff <- function(diff) {
row.stats <- c(
Mean = "RankDiffMeanLogCPM",
Variance = "RankDiffVarLogCPM",
ZerosGene = "RankDiffZeros",
MeanVar = "MeanRankVarDiff",
MeanZeros = "MeanRankZerosDiff"
)
row.ks.stats <- c(
Mean = "MeanLogCPM",
Variance = "VarLogCPM",
ZerosGene = "PctZero",
MeanVar = NA,
MeanZeros = NA
)
row.mad <- summariseStats(diff$RowData, "Dataset", row.stats, "MAD")
row.mae <- summariseStats(diff$RowData, "Dataset", row.stats, "MAE")
row.rmse <- summariseStats(diff$RowData, "Dataset", row.stats, "RMSE")
row.ks <- summariseKS(
diff$RowData,
SummarizedExperiment::rowData(diff$Reference),
"Dataset",
row.ks.stats
)
row.list <- list(row.mad, row.mae, row.rmse, row.ks)
row.list <- lapply(row.list, function(summ) {
summ[, -c(1, 2)]
})
row.summ <- data.frame(
Dataset = row.mad$Dataset,
Statistic = row.mad$Statistic
)
row.list <- c(row.summ, row.list)
row.summ <- do.call("cbind", row.list)
col.stats <- c(
LibSize = "RankDiffLibSize",
ZerosCell = "RankDiffZeros"
)
col.ks.stats <- c(
LibSize = "sum",
ZerosCell = "PctZero"
)
col.mad <- summariseStats(diff$ColData, "Dataset", col.stats, "MAD")
col.mae <- summariseStats(diff$ColData, "Dataset", col.stats, "MAE")
col.rmse <- summariseStats(diff$ColData, "Dataset", col.stats, "RMSE")
col.ks <- summariseKS(
diff$ColData,
SummarizedExperiment::colData(diff$Reference),
"Dataset",
col.ks.stats
)
col.list <- list(col.mad, col.mae, col.rmse, col.ks)
col.list <- lapply(col.list, function(summ) {
summ[, -c(1, 2)]
})
col.summ <- data.frame(
Dataset = col.mad$Dataset,
Statistic = col.mad$Statistic
)
col.list <- c(col.summ, col.list)
col.summ <- do.call("cbind", col.list)
summary <- rbind(row.summ, col.summ)
return(summary)
}
#' Summarise statistics
#'
#' Summarise columns of a data.frame using a single measure.
#'
#' @param data The data.frame to summarise
#' @param split.col Name of the column used to split the dataset
#' @param stat.cols Names of the columns to summarise. If this vector is named
#' those names will be used in the output.
#' @param measure The measure to use for summarisation.
#'
#' @return data.frame with the summarised measure, scaled and ranked
#'
#' @importFrom stats aggregate
#'
#' @noRd
summariseStats <- function(data, split.col, stat.cols,
measure = c("MAD", "MAE", "RMSE")) {
measure <- match.arg(measure)
if (is.null(names(stat.cols))) {
names(stat.cols) <- stat.cols
}
switch(measure,
"MAD" = {
measure_fun <- function(x) {
median(abs(x))
}
},
"MAE" = {
measure_fun <- function(x) {
mean(abs(x))
}
},
"RMSE" = {
measure_fun <- function(x) {
sqrt(mean(abs(x^2)))
}
}
)
summ <- aggregate(
data[, stat.cols],
list(Dataset = data[[split.col]]),
measure_fun
)
colnames(summ) <- c(split.col, names(stat.cols))
tidy.summ <- tidyStatSumm(summ, measure)
return(tidy.summ)
}
#' Summarise KS
#'
#' Summarise columns of a data.frame compared to a reference using the KS test.
#'
#' @param data The data.frame to summarise
#' @param ref The reference data.frame
#' @param split.col Name of the column used to split the dataset
#' @param stat.cols Names of the columns to summarise. If this vector is named
#' those names will be used in the output.
#'
#' @return data.frame with the summarised measure, scaled and ranked
#' @importFrom stats ks.test
#'
#' @noRd
summariseKS <- function(data, ref, split.col, stat.cols) {
if (is.null(names(stat.cols))) {
names(stat.cols) <- stat.cols
}
splits <- unique(data[[split.col]])
summ <- expand.grid(
Dataset = splits,
Statistic = names(stat.cols),
stringsAsFactors = FALSE
)
ks.res <- mapply(function(split, stat.name) {
stat <- stat.cols[stat.name]
if (!is.na(stat)) {
data.stat <- data[data[[split.col]] == split, stat]
ref.stat <- ref[[stat]]
ks <- suppressWarnings(ks.test(ref.stat, data.stat))
ks.out <- c(KS = unname(ks$statistic), KSPVal = ks$p.value)
} else {
ks.out <- c(KS = NA, KSPVal = NA)
}
return(ks.out)
}, summ$Dataset, summ$Statistic)
summ$KS <- ks.res["KS", ]
summ$KSPVal <- ks.res["KSPVal", ]
ks.ranks <- lapply(split(summ, summ$Statistic), function(x) {
rank(x$KS)
})
ks.ranks <- unlist(ks.ranks)
summ$KSRank <- ks.ranks
summ$KSRank[is.na(summ$KS)] <- NA
return(summ)
}
#' Tidy summarised statistics
#'
#' Convert a statistic summary to tidy format and add ranks and scaled values
#'
#' @param stat.summ The summary to convert
#' @param measure The name of the summarisation measure
#'
#' @return tidy data.frame with the summarised measure, scaled and ranked
#'
#' @noRd
tidyStatSumm <- function(stat.summ, measure = c("MAD", "MAE", "RMSE")) {
measure <- match.arg(measure)
summ.mat <- t(stat.summ[, -1])
colnames(summ.mat) <- stat.summ[, 1]
scale.summ <- apply(summ.mat, 1, scale)
# Check if apply has returned a vector
if (is.vector(scale.summ)) {
scale.summ <- t(as.matrix(scale.summ))
}
rank.summ <- apply(summ.mat, 1, rank)
if (is.vector(rank.summ)) {
rank.summ <- t(as.matrix(rank.summ))
}
tidy.summ <- as.data.frame.table(t(summ.mat))
colnames(tidy.summ) <- c("Dataset", "Statistic", measure)
tidy.scale <- as.data.frame.table(scale.summ)
tidy.rank <- as.data.frame.table(rank.summ)
tidy.summ[[paste0(measure, "Scaled")]] <- tidy.scale[, 3]
tidy.summ[[paste0(measure, "Rank")]] <- tidy.rank[, 3]
return(tidy.summ)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.