Skip to content

Commit

Permalink
update following BiocCheck warnings
Browse files Browse the repository at this point in the history
  • Loading branch information
josschavezf committed May 14, 2024
1 parent d7365e8 commit 5f386a3
Show file tree
Hide file tree
Showing 40 changed files with 1,307 additions and 770 deletions.
643 changes: 391 additions & 252 deletions R/spatial_enrichment.R

Large diffs are not rendered by default.

770 changes: 498 additions & 272 deletions R/spatial_interaction.R

Large diffs are not rendered by default.

26 changes: 17 additions & 9 deletions R/spdep.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@
#'
#' @param gobject Input a Giotto object.
#' @param method Specify a method name to compute auto correlation.
#' Available methods include \code{"geary.test", "lee.test", "lm.morantest","moran.test"}.
#' Available methods include
#' \code{"geary.test", "lee.test", "lm.morantest","moran.test"}.
#' @param spat_unit spatial unit
#' @param feat_type feature type
#' @param expression_values expression values to use, default = normalized
#' @param spatial_network_to_use spatial network to use, default = spatial_network
#' @param spatial_network_to_use spatial network to use,
#' default = spatial_network
#' @param verbose be verbose
#' @param return_gobject if FALSE, results are returned as data.table.
#' If TRUE, values will be appended to feature metadata
Expand Down Expand Up @@ -36,7 +38,7 @@ spdepAutoCorr <- function(gobject,
feat_type = feat_type
)
} else {
stop("gobject has not been provided\n")
stop("gobject has not been provided")
}

# Evaluate spatial autocorrelation using Giotto
Expand Down Expand Up @@ -72,7 +74,8 @@ spdepAutoCorr <- function(gobject,

result_list <- list()
progressr::with_progress({
if (step_size > 1) pb <- progressr::progressor(steps = nfeats / step_size)
if (step_size > 1) pb <- progressr::progressor(
steps = nfeats / step_size)
result_list <- lapply_flex(
seq_along(feat),
future.packages = c("data.table", "spdep"),
Expand All @@ -84,7 +87,8 @@ spdepAutoCorr <- function(gobject,
)
# Extract the estimated value from the result
result_value <- callSpdepVar$estimate[1]
temp_dt <- data.table(feat_ID = feat[feat_value], value = result_value)
temp_dt <- data.table(
feat_ID = feat[feat_value], value = result_value)
# increment progress
if (exists("pb")) if (feat_value %% step_size == 0) pb()
return(temp_dt)
Expand Down Expand Up @@ -135,10 +139,12 @@ callSpdep <- function(method, ...) {

# Check if 'method' argument is NULL, if so, stop with an error
if (is.null(method)) {
stop("The 'method' argument has not been provided. Please specify a valid method.")
stop("The 'method' argument has not been provided. Please specify a
valid method.")
}

# Check if 'method' exists in the 'spdep' package, if not, stop with an error
# Check if 'method' exists in the 'spdep' package, if not, stop with an
# error
method <- try(eval(get(method, envir = loadNamespace("spdep"))),
silent = TRUE
)
Expand Down Expand Up @@ -178,14 +184,16 @@ callSpdep <- function(method, ...) {
if (all(!(names(methodparam)) %in% allArgs)) {
stop("Invalid or missing parameters.")
}
# A vector of specified arguments that trigger 'spW <- spweights.constants()'
# A vector of specified arguments that trigger
# 'spW <- spweights.constants()'
requiredArgs <- c("n", "n1", "n2", "n3", "nn", "S0", "S1", "S2")

# Check if any of the specified arguments are required by the method
if (any(requiredArgs %in% allArgs)) {
# Obtain arguments from 'spweights.constants'
spW <- spdep::spweights.constants(listw = methodparam$listw)
# Combine user-provided arguments and 'spW', checking only against 'feats' value
# Combine user-provided arguments and 'spW', checking only against
# 'feats' value
combinedParams <- append(methodparam, spW)
} else {
combinedParams <- methodparam
Expand Down
138 changes: 93 additions & 45 deletions R/variable_genes.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,12 +10,14 @@
steps <- 1 / nr_expression_groups
prob_sequence <- seq(0, 1, steps)
prob_sequence[length(prob_sequence)] <- 1
expr_group_breaks <- stats::quantile(feat_in_cells_detected$mean_expr, probs = prob_sequence)
expr_group_breaks <- stats::quantile(
feat_in_cells_detected$mean_expr, probs = prob_sequence)

## remove zero's from cuts if there are too many and make first group zero
if (any(duplicated(expr_group_breaks))) {
m_expr_vector <- feat_in_cells_detected$mean_expr
expr_group_breaks <- stats::quantile(m_expr_vector[m_expr_vector > 0], probs = prob_sequence)
expr_group_breaks <- stats::quantile(
m_expr_vector[m_expr_vector > 0], probs = prob_sequence)
expr_group_breaks[[1]] <- 0
}

Expand All @@ -27,10 +29,12 @@
)
feat_in_cells_detected[, expr_groups := expr_groups]
feat_in_cells_detected[, cov_group_zscore := scale(cov), by = expr_groups]
feat_in_cells_detected[, selected := ifelse(cov_group_zscore > zscore_threshold, "yes", "no")]
feat_in_cells_detected[, selected := ifelse(
cov_group_zscore > zscore_threshold, "yes", "no")]

if (any(isTRUE(show_plot), isTRUE(return_plot), isTRUE(save_plot))) {
pl <- .create_cov_group_hvf_plot(feat_in_cells_detected, nr_expression_groups)
pl <- .create_cov_group_hvf_plot(
feat_in_cells_detected, nr_expression_groups)

return(list(dt = feat_in_cells_detected, pl = pl))
} else {
Expand All @@ -56,14 +60,19 @@
loess_formula <- paste0("cov~log(mean_expr)")
var_col <- "cov"

loess_model_sample <- stats::loess(loess_formula, data = feat_in_cells_detected)
feat_in_cells_detected$pred_cov_feats <- stats::predict(loess_model_sample, newdata = feat_in_cells_detected)
feat_in_cells_detected[, cov_diff := get(var_col) - pred_cov_feats, by = 1:nrow(feat_in_cells_detected)]
loess_model_sample <- stats::loess(
loess_formula, data = feat_in_cells_detected)
feat_in_cells_detected$pred_cov_feats <- stats::predict(
loess_model_sample, newdata = feat_in_cells_detected)
feat_in_cells_detected[, cov_diff := get(var_col) - pred_cov_feats,
by = 1:nrow(feat_in_cells_detected)]
data.table::setorder(feat_in_cells_detected, -cov_diff)
feat_in_cells_detected[, selected := ifelse(cov_diff > difference_in_cov, "yes", "no")]
feat_in_cells_detected[, selected := ifelse(
cov_diff > difference_in_cov, "yes", "no")]

if (any(isTRUE(show_plot), isTRUE(return_plot), isTRUE(save_plot))) {
pl <- .create_cov_loess_hvf_plot(feat_in_cells_detected, difference_in_cov, var_col)
pl <- .create_cov_loess_hvf_plot(
feat_in_cells_detected, difference_in_cov, var_col)

return(list(dt = feat_in_cells_detected, pl = pl))
} else {
Expand All @@ -87,7 +96,8 @@
test <- apply(X = scaled_matrix, MARGIN = 1, FUN = function(x) var(x))
} else {
test <- future.apply::future_apply(
X = scaled_matrix, MARGIN = 1, FUN = function(x) var(x), future.seed = TRUE
X = scaled_matrix, MARGIN = 1, FUN = function(x) var(x),
future.seed = TRUE
)
}

Expand Down Expand Up @@ -209,37 +219,48 @@
#' @param feat_type feature type
#' @param expression_values expression values to use
#' @param method method to calculate highly variable features
#' @param reverse_log_scale reverse log-scale of expression values (default = FALSE)
#' @param reverse_log_scale reverse log-scale of expression values
#' (default = FALSE)
#' @param logbase if `reverse_log_scale` is TRUE, which log base was used?
#' @param expression_threshold expression threshold to consider a gene detected
#' @param nr_expression_groups (cov_groups) number of expression groups for cov_groups
#' @param nr_expression_groups (cov_groups) number of expression groups for
#' cov_groups
#' @param zscore_threshold (cov_groups) zscore to select hvg for cov_groups
#' @param HVFname name for highly variable features in cell metadata
#' @param difference_in_cov (cov_loess) minimum difference in coefficient of variance required
#' @param var_threshold (var_p_resid) variance threshold for features for var_p_resid method
#' @param var_number (var_p_resid) number of top variance features for var_p_resid method
#' @param random_subset random subset to perform HVF detection on. Passing `NULL`
#' runs HVF on all cells.
#' @param difference_in_cov (cov_loess) minimum difference in coefficient of
#' variance required
#' @param var_threshold (var_p_resid) variance threshold for features for
#' var_p_resid method
#' @param var_number (var_p_resid) number of top variance features for
#' var_p_resid method
#' @param random_subset random subset to perform HVF detection on.
#' Passing `NULL` runs HVF on all cells.
#' @param set_seed logical. whether to set a seed when random_subset is used
#' @param seed_number seed number to use when random_subset is used
#' @param show_plot show plot
#' @param return_plot return ggplot object (overridden by `return_gobject`)
#' @param save_plot logical. directly save the plot
#' @param save_param list of saving parameters from [GiottoVisuals::all_plots_save_function()]
#' @param default_save_name default save name for saving, don't change, change save_name in save_param
#' @param save_param list of saving parameters from
#' [GiottoVisuals::all_plots_save_function()]
#' @param default_save_name default save name for saving, don't change, change
#' save_name in save_param
#' @param return_gobject boolean: return giotto object (default = TRUE)
#' @return giotto object highly variable features appended to feature metadata (`fDataDT()`)
#' @return giotto object highly variable features appended to feature metadata
#' (`fDataDT()`)
#' @details
#' Currently we provide 2 ways to calculate highly variable genes:
#'
#' \strong{1. high coeff of variance (COV) within groups: } \cr
#' First genes are binned (\emph{nr_expression_groups}) into average expression groups and
#' the COV for each feature is converted into a z-score within each bin. Features with a z-score
#' higher than the threshold (\emph{zscore_threshold}) are considered highly variable. \cr
#' First genes are binned (\emph{nr_expression_groups}) into average expression
#' groups and the COV for each feature is converted into a z-score within each
#' bin. Features with a z-score higher than the threshold
#' (\emph{zscore_threshold}) are considered highly variable. \cr
#'
#' \strong{2. high COV based on loess regression prediction: } \cr
#' A predicted COV is calculated for each feature using loess regression (COV~log(mean expression))
#' Features that show a higher than predicted COV (\emph{difference_in_cov}) are considered highly variable. \cr
#' A predicted COV is calculated for each feature using loess regression
#' (COV~log(mean expression))
#' Features that show a higher than predicted COV (\emph{difference_in_cov})
#' are considered highly variable. \cr
#'
#' @md
#' @export
Expand Down Expand Up @@ -291,8 +312,10 @@ calculateHVF <- function(gobject,
)

# expression values to be used
values <- match.arg(expression_values, unique(c("normalized", "scaled", "custom", expression_values)))
expr_values <- get_expression_values(
values <- match.arg(
expression_values,
unique(c("normalized", "scaled", "custom", expression_values)))
expr_values <- getExpression(
gobject = gobject,
spat_unit = spat_unit,
feat_type = feat_type,
Expand All @@ -318,13 +341,20 @@ calculateHVF <- function(gobject,


# print, return and save parameters
show_plot <- ifelse(is.na(show_plot), readGiottoInstructions(gobject, param = "show_plot"), show_plot)
save_plot <- ifelse(is.na(save_plot), readGiottoInstructions(gobject, param = "save_plot"), save_plot)
return_plot <- ifelse(is.na(return_plot), readGiottoInstructions(gobject, param = "return_plot"), return_plot)
show_plot <- ifelse(is.na(show_plot),
readGiottoInstructions(gobject, param = "show_plot"),
show_plot)
save_plot <- ifelse(is.na(save_plot),
readGiottoInstructions(gobject, param = "save_plot"),
save_plot)
return_plot <- ifelse(is.na(return_plot),
readGiottoInstructions(gobject, param = "return_plot"),
return_plot)


# method to use
method <- match.arg(method, choices = c("cov_groups", "cov_loess", "var_p_resid"))
method <- match.arg(
method, choices = c("cov_groups", "cov_loess", "var_p_resid"))
# select function to use based on whether future parallelization is planned
calc_cov_fun <- ifelse(
use_parallel,
Expand Down Expand Up @@ -379,14 +409,18 @@ calculateHVF <- function(gobject,

## save plot
if (isTRUE(save_plot)) {
do.call(GiottoVisuals::all_plots_save_function, c(list(gobject = gobject, plot_object = pl, default_save_name = default_save_name), save_param))
do.call(
GiottoVisuals::all_plots_save_function,
c(list(gobject = gobject, plot_object = pl,
default_save_name = default_save_name), save_param))
}

## return plot
if (isTRUE(return_plot)) {
if (isTRUE(return_gobject)) {
cat("return_plot = TRUE and return_gobject = TRUE \n
plot will not be returned to object, but can still be saved with save_plot = TRUE or manually \n")
message("return_plot = TRUE and return_gobject = TRUE \n
plot will not be returned to object, but can still be
saved with save_plot = TRUE or manually")
} else {
return(pl)
}
Expand All @@ -405,16 +439,16 @@ calculateHVF <- function(gobject,
column_names_feat_metadata <- colnames(feat_metadata[])

if (HVFname %in% column_names_feat_metadata) {
cat("\n ", HVFname, " has already been used, will be overwritten \n")
cat(HVFname, " has already been used, will be overwritten")
feat_metadata[][, eval(HVFname) := NULL]

### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
gobject <- setFeatureMetadata(gobject,
x = feat_metadata,
verbose = FALSE,
initialize = FALSE
)
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###
}

if (method == "var_p_resid") {
Expand Down Expand Up @@ -452,22 +486,26 @@ calculateHVF <- function(gobject,


# plot generation ####
.create_cov_group_hvf_plot <- function(feat_in_cells_detected, nr_expression_groups) {
.create_cov_group_hvf_plot <- function(
feat_in_cells_detected, nr_expression_groups) {
pl <- ggplot2::ggplot()
pl <- pl + ggplot2::theme_classic() +
ggplot2::theme(
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 12)
)
pl <- pl + ggplot2::geom_point(data = feat_in_cells_detected, ggplot2::aes_string(x = "mean_expr", y = "cov", color = "selected"))
pl <- pl + ggplot2::geom_point(
data = feat_in_cells_detected,
ggplot2::aes_string(x = "mean_expr", y = "cov", color = "selected"))
pl <- pl + ggplot2::scale_color_manual(
values = c(no = "lightgrey", yes = "orange"),
guide = ggplot2::guide_legend(
title = "HVF",
override.aes = list(size = 5)
)
)
pl <- pl + ggplot2::facet_wrap(~expr_groups, ncol = nr_expression_groups, scales = "free_x")
pl <- pl + ggplot2::facet_wrap(
~expr_groups, ncol = nr_expression_groups, scales = "free_x")
pl <- pl + ggplot2::theme(
axis.text.x = ggplot2::element_blank(),
strip.text = ggplot2::element_text(size = 4)
Expand All @@ -477,17 +515,26 @@ calculateHVF <- function(gobject,
}


.create_cov_loess_hvf_plot <- function(feat_in_cells_detected, difference_in_cov, var_col) {
.create_cov_loess_hvf_plot <- function(
feat_in_cells_detected, difference_in_cov, var_col) {
pl <- ggplot2::ggplot()
pl <- pl + ggplot2::theme_classic() +
ggplot2::theme(
axis.title = ggplot2::element_text(size = 14),
axis.text = ggplot2::element_text(size = 12)
)
pl <- pl + ggplot2::geom_point(data = feat_in_cells_detected, ggplot2::aes_string(x = "log(mean_expr)", y = var_col, color = "selected"))
pl <- pl + ggplot2::geom_line(data = feat_in_cells_detected, ggplot2::aes_string(x = "log(mean_expr)", y = "pred_cov_feats"), color = "blue")
pl <- pl + ggplot2::geom_point(
data = feat_in_cells_detected,
ggplot2::aes_string(x = "log(mean_expr)", y = var_col,
color = "selected"))
pl <- pl + ggplot2::geom_line(
data = feat_in_cells_detected,
ggplot2::aes_string(x = "log(mean_expr)", y = "pred_cov_feats"),
color = "blue")
hvg_line <- paste0("pred_cov_feats+", difference_in_cov)
pl <- pl + ggplot2::geom_line(data = feat_in_cells_detected, ggplot2::aes_string(x = "log(mean_expr)", y = hvg_line), linetype = 2)
pl <- pl + ggplot2::geom_line(
data = feat_in_cells_detected,
ggplot2::aes_string(x = "log(mean_expr)", y = hvg_line), linetype = 2)
pl <- pl + ggplot2::labs(x = "log(mean expression)", y = var_col)
pl <- pl + ggplot2::scale_color_manual(
values = c(no = "lightgrey", yes = "orange"),
Expand All @@ -502,7 +549,8 @@ calculateHVF <- function(gobject,

.create_calc_var_hvf_plot <- function(dt_res) {
pl <- ggplot2::ggplot()
pl <- pl + ggplot2::geom_point(data = dt_res, aes_string(x = "rank", y = "var", color = "selected"))
pl <- pl + ggplot2::geom_point(
data = dt_res, aes_string(x = "rank", y = "var", color = "selected"))
pl <- pl + ggplot2::scale_x_reverse()
pl <- pl + ggplot2::theme_classic() + ggplot2::theme(
axis.title = ggplot2::element_text(size = 14),
Expand Down
Loading

0 comments on commit 5f386a3

Please sign in to comment.