#' Reorder a matrix
#'
#' This function reorders a matrix by rows of columns
#'
#' @param matrix1 a matrix (rows=genes x columns=samples) of gene expression data
#' (e.g., scRNA-seq)
#' @param by.rows By rows (TRUE; default) or by columns
#'
#' @return a reordered matrix
#'
#' @export
reorder_matrix <- function(matrix1,
by.rows = TRUE) {
if (by.rows == TRUE) {
conf.order <- order(apply(matrix1, 1, which.max))
matrix1.reordered <- matrix1[conf.order, ]
} else {
conf.order <- order(apply(matrix1, 2, which.max))
matrix1.reordered <- matrix1[, conf.order]
}
matrix1.reordered
}
#' Compare two cluster sets matched to CCA
#'
#' This function takes cluster calls defined in two different data sets and then
#' determines to what extent these cluster calls match up with cluster calls from CCA.
#'
#' @param cl a matrix (rows=genes x columns=samples) of gene expression data
#' (e.g., scRNA-seq)
#' @param by.rows By rows (TRUE; default) or by columns
#'
#' @return a reordered matrix
#'
#' @export
compareClusterCalls <- function(cl,
ref.cl,
cl.anno,
plot.title = NA, plot.silent = TRUE,
heat.colors = colorRampPalette(c("grey99", "orange", "red"))(100),
row.cl.num = min(length(unique(cl)), length(unique(ref.cl)))) {
library(grid)
library(pheatmap)
conf1 <- table(cl, ref.cl)
conf1 <- sweep(conf1, 1, rowSums(conf1), "/")
conf2 <- reorder_matrix(conf1)
# Cluster co-occurence
cl.prop.cocl <- apply(conf1, 2, function(x) {
grid1 <- expand.grid(x, x)
min.prop <- apply(grid1, 1, min)
})
cl.prop.cocl.total <- apply(cl.prop.cocl, 1, sum)
cl.prop.cocl.m <- matrix(cl.prop.cocl.total, nrow(conf1), nrow(conf1),
dimnames = list(rownames(conf1), rownames(conf1))
)
ph1 <- pheatmap(conf2,
cutree_rows = row.cl.num, clustering_method = "ward.D2",
annotation_row = cl.anno, color = heat.colors, fontsize = 6,
main = plot.title, silent = plot.silent
)
list(conf = conf2, cocl = cl.prop.cocl.m, ph = ph1)
}
#' Get some summary statistics
#'
#' This does the summary. For each group return a vector with N, mean, and sd.
#' This is called by qcPlot.
#'
#' @param data Annotation data frame
#' @param measurevar what variables to measure
#' @param groupvars what to group by
#' @param na.rm how to treat NA
#' @param conf.interval confidence interval (default = .95)
#' @param .drop drop something?
#' @param roundall should everything be rounded
#'
#' @return a data frame of statistics
#'
#' @export
summarySE <- function(data = NULL,
measurevar,
groupvars = NULL,
na.rm = FALSE,
conf.interval = .95,
.drop = TRUE,
roundall = F) {
require(dplyr)
#
names(data)[names(data) == measurevar] <- "measurevar"
datac <- data %>%
select(one_of(groupvars, "measurevar")) %>%
filter(ifelse(na.rm == T, !is.na(measurevar), T)) %>%
mutate(measurevar = as.numeric(measurevar)) %>%
group_by_(c(groupvars)) %>%
summarise(
N = n(),
median = median(measurevar),
mean = mean(measurevar),
max = max(measurevar),
sd = ifelse(N == 1, 0, sd(measurevar)),
q25 = as.numeric(quantile(measurevar, 0.25)),
q75 = as.numeric(quantile(measurevar, 0.75))
) %>%
mutate(se = sd / sqrt(N))
# %>% mutate(ci = se * qt(conf.interval/2 + 0.5, N-1))
if (roundall) {
roundcols <- c("median", "mean", "max", "sd", "q25", "q75", "se", "ci")
datac[roundcols] <- round(datac[roundcols], 3)
}
# datac <- datac %>% mutate(xpos = 1:n())
datac
}
#' Make QC plots
#'
#' Makes QC plots
#'
#' @param anno Annotation data frame for patch-seq
#' @param dendcluster_anno what to cluster by
#' @param groupvars what to group by
#' @param scaleLimits scaleLimits
#' @param scaleBreaks scaleBreaks
#' @param scaleLabels scaleLabels
#' @param ylab ylab
#' @param fileName fileName
#'
#' @return the plot is returned
#'
#' @export
qcPlot <- function(anno,
dendcluster_anno,
name,
scaleLimits = c(-5000, 12000),
scaleBreaks = seq(0, 12000, 2000),
scaleLabels = seq(0, 12, 2),
ylab = "value",
fileName = gsub("\\.", "_", gsub("_label", "", name)),
outputFolder) {
# dendcluster_id is the annotation for cluster ordering based on the current, bootstrapped dendrogram
stats <- summarySE(data = anno, measurevar = name, groupvars = "dendcluster_id")
genes_plot <- ggplot() +
# geom_quasirandom from the ggbeeswarm package
# makes violin-shaped jittered point plots
geom_quasirandom(
data = anno,
aes(
x = dendcluster_id,
y = eval(parse(text = name))
),
color = "skyblue",
# Need to set position_jitter height = 0 to prevent
# jitter on the y-axis, which changes data representation
position = position_jitter(width = .3, height = 0), size = 0.1
) +
# Errorbars built using stats values
geom_errorbar(
data = stats,
aes(x = dendcluster_id, ymin = q25, ymax = q75),
size = 0.2
) +
# Median points from stats
geom_point(
data = stats,
aes(x = dendcluster_id, y = median),
color = "red",
size = 0.5
) +
# Cluster labels as text objects
geom_text(
data = dendcluster_anno,
aes(
x = dendcluster_id, y = 0, label = dendcluster_label,
color = dendcluster_color
),
angle = 90,
hjust = 2,
vjust = 0.3,
size = 2 * 5 / 6
) +
scale_color_identity() +
# Expand the y scale so that the labels are visible
scale_y_continuous(ylab,
limits = scaleLimits,
breaks = scaleBreaks,
labels = scaleLabels
) +
# Remove X-axis title
scale_x_continuous("") +
theme_bw() +
# Theme tuning
theme(
axis.text.x = element_blank(),
axis.ticks = element_blank(),
panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
)
ggsave(paste0(outputFolder, fileName, "_QC.pdf"), genes_plot, width = 8, height = 4, useDingbats = F)
genes_plot
}
