Merge branch 'glmm' of https://github.com/immunogenomics/presto into glmm (ac17d381) · Commits · github_fork / Presto 20250701

R/glmm.R

+24 −22

Original line number	Diff line number	Diff line
		@@ -303,38 +303,34 @@ presto.presto <- function(


		#' @export
		correct_counts <- function(object, effects_remove, umi_common, verbose=0) {

		correct_counts <- function (object, effects_remove, umi_common, verbose = 0)
		{
		if (missing(umi_common)) {
		umi_common <- mean(log(Matrix::colSums(object$response)))
		# umi_common <- exp(mean(object$meta_data$EXPOSURE))
		}

		## read depth should not be removed but equalized
		effects_remove <- setdiff(effects_remove, 'EXPOSURE')
		effects_remove <- setdiff(effects_remove, "EXPOSURE")
		idx_keep <- object$betanames_df %>%
		tibble::rowid_to_column('idx') %>%
		subset(!grepl(paste(effects_remove, collapse = '\|'), grpvar)) %>%
		tibble::rowid_to_column("idx") %>%
		# subset(!grepl(paste(effects_remove, collapse = "\|"), grpvar)) %>%
		subset(!grpvar %in% effects_remove) %>%
		with(idx)

		if (verbose > 0) {
		message('remove')
		object$betanames_df[-idx_keep, ] %>% with(unique(grpvar)) %>% print()
		message('preserve')
		object$betanames_df[idx_keep, ] %>% with(unique(grpvar)) %>% print()
		message("remove")
		object$betanames_df[-idx_keep, ] %>% with(unique(grpvar)) %>%
		print()
		message("preserve")
		object$betanames_df[idx_keep, ] %>% with(unique(grpvar)) %>%
		print()
		}

		# b1 <- object$beta[which(object$betanames_df$grpvar == 'EXPOSURE'), ]

		design_keep <- object$design[idx_keep, ]
		design_keep['EXPOSURE', ] <- umi_common

		effect_keep <- exp(Matrix::crossprod(design_keep, object$beta[idx_keep, ]))
		design_keep["EXPOSURE", ] <- umi_common
		effect_keep <- exp(Matrix::crossprod(design_keep, object$beta[idx_keep,
		]))
		object$corrected <- Matrix::t(effect_keep + object$epsilon)
		object$corrected <- as(object$corrected, class(object$response)[[1]])
		colnames(object$corrected) <- colnames(object$response)
		row.names(object$corrected) <- row.names(object$response)

		return(object)
		}

		@@ -438,7 +434,8 @@ toptable.presto <- function(
		) %>%
		dplyr::mutate(
		wald = beta / sigma,
		pval = 2 * (1 - pnorm(abs(beta / sigma)))
		pval = -2 * pnorm(abs(beta / sigma), log.p = TRUE)
		# pval = 2 * (1 - pnorm(abs(beta / sigma)))
		) %>%
		dplyr::mutate(fdr = p.adjust(pval, 'BH')) %>%
		subset(pval < max_pval & beta > min_beta & beta < max_beta & fdr < max_fdr)
		@@ -494,7 +491,12 @@ effects.presto <- function(object, effects) {
		suppressMessages({
		object$effect <- dplyr::full_join(beta_tidy, sigma_tidy) %>%
		dplyr::left_join(object$log_mu) %>%
		dplyr::mutate(wald = beta/sigma, pval = 2 * (1 - pnorm(abs(beta/sigma)))) %>%

		dplyr::mutate(
		wald = beta/sigma,
		pval = -2 * pnorm(abs(beta / sigma), log.p = TRUE)
		# pval = 2 * (1 - pnorm(abs(beta/sigma)))
		) %>%
		dplyr::mutate(fdr = p.adjust(pval, "BH"))
		})
		return(object)

Admin message