Unverified Commit 968ecbbc authored by ilyakorsunsky's avatar ilyakorsunsky Committed by GitHub
Browse files

in correct_counts, specify exact levels to remove

parent 8ec89c2f
Loading
Loading
Loading
Loading
+16 −20
Original line number Diff line number Diff line
@@ -217,38 +217,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)
}