BUG: failure: length > 1 in coercion to logical

[HenrikBengtsson]

Running

$ export _R_CHECK_LENGTH_1_LOGIC2_=verbose
$ R CMD check infercnv_1.3.0.tar.gz

reveals a

'length(x) = 2 > 1' in coercion to 'logical(1)'

bug that potentially may cause invalid scientific results. For an explanation on _R_CHECK_LENGTH_1_LOGIC2_, see Section on R 3.6.0 in news() or https://github.com/HenrikBengtsson/Wishlist-for-R/issues/48.

The full output is:

$ R CMD check infercnv_1.3.0.tar.gz 
* using log directory ‘/tmp/infercnv.Rcheck’
* using R version 3.6.1 Patched (2019-10-31 r77348)
* using platform: x86_64-pc-linux-gnu (64-bit)
* using session charset: UTF-8
* checking for file ‘infercnv/DESCRIPTION’ ... OK
* checking extension type ... Package
* this is package ‘infercnv’ version ‘1.3.0’
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for executable files ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking for sufficient/correct file permissions ... OK
* checking whether package ‘infercnv’ can be installed ... OK
* checking installed package size ... NOTE
  installed size is  5.1Mb
  sub-directories of 1Mb or more:
    extdata   3.1Mb
* checking package directory ... OK
* checking ‘build’ directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* checking whether the package can be loaded ... OK
* checking whether the package can be loaded with stated dependencies ... OK
* checking whether the package can be unloaded cleanly ... OK
* checking whether the namespace can be loaded with stated dependencies ... OK
* checking whether the namespace can be unloaded cleanly ... OK
* checking loading without being on the library search path ... OK
* checking dependencies in R code ... NOTE
Unexported object imported by a ':::' call: ‘HiddenMarkov:::makedensity’
  See the note in ?`:::` about the use of this operator.
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... OK
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking contents of ‘data’ directory ... OK
* checking data for non-ASCII characters ... OK
* checking data for ASCII and uncompressed saves ... OK
* checking installed files from ‘inst/doc’ ... OK
* checking files in ‘vignettes’ ... OK
* checking examples ... ERROR
Running examples in ‘infercnv-Ex.R’ failed
The error most likely occurred in:

> ### Name: inferCNVBayesNet
> ### Title: inferCNVBayesNet: Run Bayesian Network Mixture Model To Obtain
> ###   Posterior Probabilities For HMM Predicted States
> ### Aliases: inferCNVBayesNet
> 
> ### ** Examples
> 
> data(data)
> data(annots)
> data(genes)
> data(HMM_states)
> 
> infercnv_obj <- infercnv::CreateInfercnvObject(raw_counts_matrix=data, 
+                                                gene_order_file=genes,
+                                                annotations_file=annots,
+                                                ref_group_names=c("normal"))
INFO [2019-11-07 10:12:22] ::order_reduce:Start.
INFO [2019-11-07 10:12:22] .order_reduce(): expr and order match.
INFO [2019-11-07 10:12:22] ::process_data:order_reduce:Reduction from positional data, new dimensions (r,c) = 8252,20 Total=2817798 Min=0 Max=6927.
INFO [2019-11-07 10:12:22] num genes removed taking into account provided gene ordering list: 2308 = 27.9689772176442% removed.
INFO [2019-11-07 10:12:22] validating infercnv_obj
> 
> out_dir = tempfile()
> infercnv_obj <- infercnv::run(infercnv_obj,
+                               cutoff=1,
+                               out_dir=out_dir, 
+                               cluster_by_groups=TRUE, 
+                               denoise=TRUE,
+                               HMM=TRUE,
+                               num_threads=2,
+                               no_plot=TRUE)
 ----------- FAILURE REPORT -------------- 
 --- failure: length > 1 in coercion to logical ---
 --- srcref --- 
: 
 --- package (from environment) --- 
infercnv
 --- call from context --- 
infercnv::run(infercnv_obj, cutoff = 1, out_dir = out_dir, cluster_by_groups = TRUE, 
    denoise = TRUE, HMM = TRUE, num_threads = 2, no_plot = TRUE)
 --- call from argument --- 
HMM && HMM_type == "i6"
 --- R stacktrace ---
where 1: infercnv::run(infercnv_obj, cutoff = 1, out_dir = out_dir, cluster_by_groups = TRUE, 
    denoise = TRUE, HMM = TRUE, num_threads = 2, no_plot = TRUE)

 --- value of length: 2 type: logical ---
[1]  TRUE FALSE
 --- function from context --- 
function (infercnv_obj, cutoff = 1, min_cells_per_gene = 3, out_dir = NULL, 
    window_length = 101, smooth_method = c("pyramidinal", "runmeans", 
        "coordinates"), num_ref_groups = NULL, ref_subtract_use_mean_bounds = TRUE, 
    cluster_by_groups = FALSE, cluster_references = TRUE, k_obs_groups = 1, 
    hclust_method = "ward.D2", max_centered_threshold = 3, scale_data = FALSE, 
    HMM = FALSE, HMM_transition_prob = 1e-06, HMM_report_by = c("subcluster", 
        "consensus", "cell"), HMM_type = c("i6", "i3"), HMM_i3_pval = 0.05, 
    HMM_i3_use_KS = TRUE, BayesMaxPNormal = 0.5, sim_method = "meanvar", 
    sim_foreground = FALSE, reassignCNVs = TRUE, analysis_mode = c("samples", 
        "subclusters", "cells"), tumor_subcluster_partition_method = c("random_trees", 
        "qnorm", "pheight", "qgamma", "shc"), tumor_subcluster_pval = 0.1, 
    denoise = FALSE, noise_filter = NA, sd_amplifier = 1.5, noise_logistic = FALSE, 
    outlier_method_bound = "average_bound", outlier_lower_bound = NA, 
    outlier_upper_bound = NA, final_scale_limits = NULL, final_center_val = NULL, 
    debug = FALSE, num_threads = 4, plot_steps = FALSE, resume_mode = TRUE, 
    png_res = 300, plot_probabilities = TRUE, diagnostics = FALSE, 
    remove_genes_at_chr_ends = FALSE, prune_outliers = FALSE, 
    mask_nonDE_genes = FALSE, mask_nonDE_pval = 0.05, test.use = "wilcoxon", 
    require_DE_all_normals = "any", hspike_aggregate_normals = FALSE, 
    no_plot = FALSE, no_prelim_plot = FALSE, output_format = "png", 
    useRaster = TRUE, up_to_step = 100) 
{
    smooth_method = match.arg(smooth_method)
    if (HMM && HMM_type == "i6" && smooth_method == "coordinates") {
        flog.error("Cannot use 'coordinate' smoothing method with i6 HMM model at this time.")
        stop("Incompatible HMM mode and smoothing method.")
    }
    if (smooth_method == "coordinates" && window_length < 10000) {
        window_length = 1e+07
        flog.warn(paste0("smooth_method set to 'coordinates', but window_length ", 
            "is less than 10.000, setting it to 10.000.000. Please ", 
            "set a different value > 10.000 if this default does not seem suitable."))
    }
    HMM_report_by = match.arg(HMM_report_by)
    analysis_mode = match.arg(analysis_mode)
    if (analysis_mode == "samples" && HMM_report_by == "subclusters") {
        HMM_report_by = "samples"
        flog.warn(paste0("analysis_mode is \"samples\" but HMM_report_by is \"subclusters\", "), 
            "changing HMM_report_by to \"samples\".")
    }
    else if (analysis_mode == "subclusters" && HMM_report_by == 
        "samples") {
        HMM_report_by = "subclusters"
        flog.warn(paste0("analysis_mode is \"subclusters\" but HMM_report_by is \"samples\", "), 
            "changing HMM_report_by to \"subclusters\".")
    }
    else if (analysis_mode == "cells") {
        flog.warn(paste0("analysis_mode is \"cells\" but HMM_report_by is not, "), 
            "changing HMM_report_by to \"cells\".")
        HMM_report_by = "cells"
    }
    tumor_subcluster_partition_method = match.arg(tumor_subcluster_partition_method)
    HMM_type = match.arg(HMM_type)
    if (debug) {
        flog.threshold(DEBUG)
    }
    else {
        flog.threshold(INFO)
    }
    flog.info(paste("::process_data:Start", sep = ""))
    infercnv.env$GLOBAL_NUM_THREADS <- num_threads
    if (is.null(out_dir)) {
        flog.error("Error, out_dir is NULL, please provide a path.")
        stop("out_dir is NULL")
    }
    call_match = match.call()
    if (Reduce("|", is(call_match[["out_dir"]]) == "call")) {
        out_dir = eval(call_match[["out_dir"]])
        call_match[["out_dir"]] = out_dir
    }
    non_default_args = call_match[2:length(call_match)]
    if (out_dir != "." && !file.exists(out_dir)) {
        flog.info(paste0("Creating output path ", out_dir))
        dir.create(out_dir)
    }
    infercnv_obj@options = c(infercnv_obj@options, lapply(non_default_args, 
        deparse))
    call_match[[1]] <- as.symbol(".get_relevant_args_list")
    reload_info = eval(call_match)
    reload_info$relevant_args
    reload_info$expected_file_names
    skip_hmm = 0
    skip_mcmc = 0
    skip_past = 0
    if (resume_mode) {
        flog.info("Checking for saved results.")
        for (i in rev(seq_along(reload_info$expected_file_names))) {
            if (file.exists(reload_info$expected_file_names[[i]])) {
                flog.info(paste0("Trying to reload from step ", 
                  i))
                if ((i == 17 || i == 18 || i == 19) && skip_hmm == 
                  0) {
                  hmm.infercnv_obj = readRDS(reload_info$expected_file_names[[i]])
                  if (!.compare_args(infercnv_obj@options, unlist(reload_info$relevant_args[1:i]), 
                    hmm.infercnv_obj@options)) {
                    rm(hmm.infercnv_obj)
                    invisible(gc())
                  }
                  else {
                    hmm.infercnv_obj@options = infercnv_obj@options
                    skip_hmm = i - 16
                    flog.info(paste0("Using backup HMM from step ", 
                      i))
                  }
                }
                else {
                  reloaded_infercnv_obj = readRDS(reload_info$expected_file_names[[i]])
                  if (skip_past > i) {
                    if (20 > i) {
                      break
                    }
                    else {
                      next
                    }
                  }
                  if (.compare_args(infercnv_obj@options, unlist(reload_info$relevant_args[1:i]), 
                    reloaded_infercnv_obj@options)) {
                    options_backup = infercnv_obj@options
                    infercnv_obj = reloaded_infercnv_obj
                    rm(reloaded_infercnv_obj)
                    infercnv_obj@options = options_backup
                    skip_past = i
                    flog.info(paste0("Using backup from step ", 
                      i))
                    if (i < 21) {
                      break
                    }
                  }
                  else {
                    rm(reloaded_infercnv_obj)
                    invisible(gc())
                  }
                }
            }
        }
    }
    step_count = 0
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %d: incoming data\n", step_count))
    if (skip_past < step_count) {
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: Removing lowly expressed genes\n", 
        step_count))
    if (skip_past < step_count) {
        infercnv_obj <- require_above_min_mean_expr_cutoff(infercnv_obj, 
            cutoff)
        infercnv_obj <- require_above_min_cells_ref(infercnv_obj, 
            min_cells_per_gene = min_cells_per_gene)
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: normalization by sequencing depth\n", 
        step_count))
    resume_file_token = ifelse((HMM), paste0("HMM", HMM_type), 
        "")
    if (skip_past < step_count) {
        infercnv_obj <- normalize_counts_by_seq_depth(infercnv_obj)
        if (HMM && HMM_type == "i6") {
            infercnv_obj <- .build_and_add_hspike(infercnv_obj, 
                sim_method = sim_method, aggregate_normals = hspike_aggregate_normals)
            if (sim_foreground) {
                infercnv_obj <- .sim_foreground(infercnv_obj, 
                  sim_method = sim_method)
            }
        }
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: log transformation of data\n", 
        step_count))
    if (skip_past < step_count) {
        infercnv_obj <- log2xplus1(infercnv_obj)
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        invisible(gc())
        if (plot_steps) {
            plot_cnv(infercnv_obj = infercnv_obj, k_obs_groups = k_obs_groups, 
                cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                out_dir = out_dir, title = sprintf("%02d_log_transformed_data", 
                  step_count), output_filename = sprintf("infercnv.%02d_log_transformed", 
                  step_count), output_format = output_format, 
                write_expr_matrix = TRUE, png_res = png_res, 
                useRaster = useRaster)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (scale_data) {
        flog.info(sprintf("\n\n\tSTEP %02d: scaling all expression data\n", 
            step_count))
        if (skip_past < step_count) {
            infercnv_obj <- scale_infercnv_expr(infercnv_obj)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_scaled", 
                    step_count), output_filename = sprintf("infercnv.%02d_scaled", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (!is.null(num_ref_groups)) {
        if (!has_reference_cells(infercnv_obj)) {
            stop("Error, no reference cells defined. Cannot split them into groups as requested")
        }
        flog.info(sprintf("\n\n\tSTEP %02d: splitting reference data into %d clusters\n", 
            step_count, num_ref_groups))
        if (skip_past < step_count) {
            infercnv_obj <- split_references(infercnv_obj, num_groups = num_ref_groups, 
                hclust_method = hclust_method)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (analysis_mode == "subclusters" & tumor_subcluster_partition_method == 
        "random_trees") {
        flog.info(sprintf("\n\n\tSTEP %02d: computing tumor subclusters via %s\n", 
            step_count, tumor_subcluster_partition_method))
        resume_file_token = paste0(resume_file_token, ".rand_trees")
        if (skip_past < step_count) {
            infercnv_obj <- define_signif_tumor_subclusters_via_random_smooothed_trees(infercnv_obj, 
                p_val = tumor_subcluster_pval, hclust_method = hclust_method, 
                cluster_by_groups = cluster_by_groups)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_tumor_subclusters.%s", 
                    step_count, tumor_subcluster_partition_method), 
                  output_filename = sprintf("infercnv.%02d_tumor_subclusters.%s", 
                    step_count, tumor_subcluster_partition_method), 
                  output_format = output_format, write_expr_matrix = TRUE, 
                  png_res = png_res, useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: removing average of reference data (before smoothing)\n", 
        step_count))
    if (skip_past < step_count) {
        infercnv_obj <- subtract_ref_expr_from_obs(infercnv_obj, 
            inv_log = FALSE, use_bounds = ref_subtract_use_mean_bounds)
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        invisible(gc())
        if (plot_steps) {
            plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                out_dir = out_dir, title = sprintf("%02d_remove_average", 
                  step_count), output_filename = sprintf("infercnv.%02d_remove_average", 
                  step_count), output_format = output_format, 
                write_expr_matrix = TRUE, png_res = png_res, 
                useRaster = useRaster)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (!is.na(max_centered_threshold)) {
        flog.info(sprintf("\n\n\tSTEP %02d: apply max centered expression threshold: %s\n", 
            step_count, max_centered_threshold))
        if (skip_past < step_count) {
            threshold = max_centered_threshold
            if (is.character(max_centered_threshold) && max_centered_threshold == 
                "auto") {
                threshold = mean(abs(get_average_bounds(infercnv_obj)))
                flog.info(sprintf("Setting max centered threshoolds via auto to: +- %g", 
                  threshold))
            }
            infercnv_obj <- apply_max_threshold_bounds(infercnv_obj, 
                threshold = threshold)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_apply_max_centered_expr_threshold", 
                    step_count), output_filename = sprintf("infercnv.%02d_apply_max_centred_expr_threshold", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: Smoothing data per cell by chromosome\n", 
        step_count))
    if (skip_past < step_count) {
        if (smooth_method == "runmeans") {
            infercnv_obj <- smooth_by_chromosome_runmeans(infercnv_obj, 
                window_length)
        }
        else if (smooth_method == "pyramidinal") {
            infercnv_obj <- smooth_by_chromosome(infercnv_obj, 
                window_length = window_length, smooth_ends = TRUE)
        }
        else if (smooth_method == "coordinates") {
            infercnv_obj <- smooth_by_chromosome_coordinates(infercnv_obj, 
                window_length = window_length)
        }
        else {
            stop(sprintf("Error, don't recognize smoothing method: %s", 
                smooth_method))
        }
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        invisible(gc())
        if (plot_steps) {
            plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                out_dir = out_dir, title = sprintf("%02d_smoothed_by_chr", 
                  step_count), output_filename = sprintf("infercnv.%02d_smoothed_by_chr", 
                  step_count), output_format = output_format, 
                write_expr_matrix = TRUE, png_res = png_res, 
                useRaster = useRaster)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: re-centering data across chromosome after smoothing\n", 
        step_count))
    if (skip_past < step_count) {
        infercnv_obj <- center_cell_expr_across_chromosome(infercnv_obj, 
            method = "median")
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        invisible(gc())
        if (plot_steps) {
            plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                out_dir = out_dir, title = sprintf("%02d_centering_of_smoothed", 
                  step_count), output_filename = sprintf("infercnv.%02d_centering_of_smoothed", 
                  step_count), output_format = output_format, 
                write_expr_matrix = TRUE, png_res = png_res, 
                useRaster = useRaster)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: removing average of reference data (after smoothing)\n", 
        step_count))
    if (skip_past < step_count) {
        infercnv_obj <- subtract_ref_expr_from_obs(infercnv_obj, 
            inv_log = FALSE, use_bounds = ref_subtract_use_mean_bounds)
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        invisible(gc())
        if (plot_steps) {
            plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                out_dir = out_dir, title = sprintf("%02d_remove_average", 
                  step_count), output_filename = sprintf("infercnv.%02d_remove_average", 
                  step_count), output_format = output_format, 
                write_expr_matrix = TRUE, png_res = png_res, 
                useRaster = useRaster)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (remove_genes_at_chr_ends == TRUE && smooth_method != 
        "coordinates") {
        flog.info(sprintf("\n\n\tSTEP %02d: removing genes at chr ends\n", 
            step_count))
        if (skip_past < step_count) {
            infercnv_obj <- remove_genes_at_ends_of_chromosomes(infercnv_obj, 
                window_length)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_remove_genes_at_chr_ends", 
                    step_count), output_filename = sprintf("infercnv.%02d_remove_genes_at_chr_ends", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    flog.info(sprintf("\n\n\tSTEP %02d: invert log2(FC) to FC\n", 
        step_count))
    if (skip_past < step_count) {
        infercnv_obj <- invert_log2(infercnv_obj)
        saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        invisible(gc())
        if (plot_steps) {
            plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                out_dir = out_dir, title = sprintf("%02d_invert_log_transform log(FC)->FC", 
                  step_count), output_filename = sprintf("infercnv.%02d_invert_log_FC", 
                  step_count), output_format = output_format, 
                write_expr_matrix = TRUE, png_res = png_res, 
                useRaster = useRaster)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (analysis_mode == "subclusters" & tumor_subcluster_partition_method != 
        "random_trees") {
        resume_file_token = paste0(resume_file_token, ".", tumor_subcluster_partition_method)
        flog.info(sprintf("\n\n\tSTEP %02d: computing tumor subclusters via %s\n", 
            step_count, tumor_subcluster_partition_method))
        if (skip_past < step_count) {
            infercnv_obj <- define_signif_tumor_subclusters(infercnv_obj, 
                p_val = tumor_subcluster_pval, hclust_method = hclust_method, 
                cluster_by_groups = cluster_by_groups, partition_method = tumor_subcluster_partition_method)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_tumor_subclusters", 
                    step_count), output_filename = sprintf("infercnv.%02d_tumor_subclusters", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    else if (analysis_mode != "subclusters") {
        flog.info(sprintf("\n\n\tSTEP %02d: Clustering samples (not defining tumor subclusters)\n", 
            step_count))
        if (skip_past < step_count) {
            infercnv_obj <- define_signif_tumor_subclusters(infercnv_obj, 
                p_val = tumor_subcluster_pval, hclust_method = hclust_method, 
                cluster_by_groups = cluster_by_groups, partition_method = "none")
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
        }
    }
    if (skip_past < step_count) {
        infercnv_obj_file = file.path(out_dir, "preliminary.infercnv_obj")
        saveRDS(infercnv_obj, file = infercnv_obj_file)
        invisible(gc())
        if (!(no_prelim_plot | no_plot)) {
            prelim_heatmap_png = "infercnv.preliminary.png"
            if (!file.exists(file.path(out_dir, prelim_heatmap_png))) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = "Preliminary infercnv (pre-noise filtering)", 
                  output_filename = "infercnv.preliminary", output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (prune_outliers) {
        flog.info(sprintf("\n\n\tSTEP %02d: Removing outliers\n", 
            step_count))
        if (skip_past < step_count) {
            infercnv_obj = remove_outliers_norm(infercnv_obj, 
                out_method = outlier_method_bound, lower_bound = outlier_lower_bound, 
                upper_bound = outlier_upper_bound)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_removed_outliers", 
                    step_count), output_filename = sprintf("infercnv.%02d_removed_outliers", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    hmm_resume_file_token = paste0(resume_file_token, ".hmm_mode-", 
        analysis_mode)
    if (skip_hmm < 1) {
        if (HMM) {
            flog.info(sprintf("\n\n\tSTEP %02d: HMM-based CNV prediction\n", 
                step_count))
            if (HMM_type == "i6") {
                hmm_center = 3
                hmm_state_range = c(0, 6)
            }
            else {
                hmm_center = 2
                hmm_state_range = c(1, 3)
            }
            if (analysis_mode == "subclusters") {
                if (HMM_type == "i6") {
                  hmm.infercnv_obj <- predict_CNV_via_HMM_on_tumor_subclusters(infercnv_obj, 
                    t = HMM_transition_prob)
                }
                else if (HMM_type == "i3") {
                  hmm.infercnv_obj <- i3HMM_predict_CNV_via_HMM_on_tumor_subclusters(infercnv_obj, 
                    i3_p_val = HMM_i3_pval, t = HMM_transition_prob, 
                    use_KS = HMM_i3_use_KS)
                }
                else {
                  stop("Error, not recognizing HMM_type")
                }
                if (tumor_subcluster_partition_method == "random_trees") {
                  hmm.infercnv_obj <- .redo_hierarchical_clustering(hmm.infercnv_obj, 
                    hclust_method = hclust_method)
                }
            }
            else if (analysis_mode == "cells") {
                if (HMM_type == "i6") {
                  hmm.infercnv_obj <- predict_CNV_via_HMM_on_indiv_cells(infercnv_obj, 
                    t = HMM_transition_prob)
                }
                else if (HMM_type == "i3") {
                  hmm.infercnv_obj <- i3HMM_predict_CNV_via_HMM_on_indiv_cells(infercnv_obj, 
                    i3_p_val = HMM_i3_pval, t = HMM_transition_prob, 
                    use_KS = HMM_i3_use_KS)
                }
                else {
                  stop("Error, not recognizing HMM_type")
                }
            }
            else {
                if (HMM_type == "i6") {
                  hmm.infercnv_obj <- predict_CNV_via_HMM_on_whole_tumor_samples(infercnv_obj, 
                    t = HMM_transition_prob)
                }
                else if (HMM_type == "i3") {
                  hmm.infercnv_obj <- i3HMM_predict_CNV_via_HMM_on_tumor_subclusters(infercnv_obj, 
                    i3_p_val = HMM_i3_pval, t = HMM_transition_prob, 
                    use_KS = HMM_i3_use_KS)
                }
                else {
                  stop("Error, not recognizing HMM_type")
                }
            }
            saveRDS(hmm.infercnv_obj, reload_info$expected_file_names[[step_count]])
            generate_cnv_region_reports(hmm.infercnv_obj, output_filename_prefix = sprintf("%02d_HMM_pred%s", 
                step_count, hmm_resume_file_token), out_dir = out_dir, 
                ignore_neutral_state = hmm_center, by = HMM_report_by)
            invisible(gc())
            if (!no_plot) {
                plot_cnv(infercnv_obj = hmm.infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_HMM_preds", 
                    step_count), output_filename = sprintf("infercnv.%02d_HMM_pred%s", 
                    step_count, hmm_resume_file_token), output_format = output_format, 
                  write_expr_matrix = TRUE, x.center = hmm_center, 
                  x.range = hmm_state_range, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (skip_hmm < 2) {
        if (HMM == TRUE && BayesMaxPNormal > 0 && length(unique(apply(hmm.infercnv_obj@expr.data, 
            2, unique))) != 1) {
            flog.info(sprintf("\n\n\tSTEP %02d: Run Bayesian Network Model on HMM predicted CNV's\n", 
                step_count))
            mcmc_obj <- infercnv::inferCNVBayesNet(infercnv_obj = infercnv_obj, 
                HMM_states = hmm.infercnv_obj@expr.data, file_dir = out_dir, 
                no_plot = no_plot, postMcmcMethod = "removeCNV", 
                out_dir = file.path(out_dir, sprintf("BayesNetOutput.%s", 
                  hmm_resume_file_token)), resume_file_token = hmm_resume_file_token, 
                quietly = TRUE, CORES = num_threads, plotingProbs = plot_probabilities, 
                diagnostics = diagnostics, HMM_type = HMM_type, 
                k_obs_groups = k_obs_groups, cluster_by_groups = cluster_by_groups, 
                reassignCNVs = reassignCNVs)
            mcmc_obj_file = file.path(out_dir, sprintf("%02d_HMM_pred.Bayes_Net%s.mcmc_obj", 
                step_count, hmm_resume_file_token))
            saveRDS(mcmc_obj, file = mcmc_obj_file)
            mcmc_obj_hmm_states_list <- infercnv::filterHighPNormals(MCMC_inferCNV_obj = mcmc_obj, 
                HMM_states = hmm.infercnv_obj@expr.data, BayesMaxPNormal = BayesMaxPNormal)
            hmm_states_highPnormCNVsRemoved.matrix <- mcmc_obj_hmm_states_list[[2]]
            hmm.infercnv_obj@expr.data <- hmm_states_highPnormCNVsRemoved.matrix
            saveRDS(hmm.infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (!no_plot) {
                plot_cnv(infercnv_obj = hmm.infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_HMM_preds_Bayes_Net", 
                    step_count), output_filename = sprintf("infercnv.%02d_HMM_pred.Bayes_Net.Pnorm_%g", 
                    step_count, BayesMaxPNormal), output_format = output_format, 
                  write_expr_matrix = TRUE, x.center = 3, x.range = c(0, 
                    6), png_res = png_res, useRaster = useRaster)
            }
            adjust_genes_regions_report(mcmc_obj_hmm_states_list[[1]], 
                input_filename_prefix = sprintf("%02d_HMM_pred%s", 
                  (step_count - 1), hmm_resume_file_token), output_filename_prefix = sprintf("HMM_CNV_predictions.%s.Pnorm_%g", 
                  hmm_resume_file_token, BayesMaxPNormal), out_dir = out_dir)
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (skip_hmm < 3) {
        if (HMM) {
            flog.info(sprintf("\n\n\tSTEP %02d: Converting HMM-based CNV states to repr expr vals\n", 
                step_count))
            if (HMM_type == "i6") {
                hmm.infercnv_obj <- assign_HMM_states_to_proxy_expr_vals(hmm.infercnv_obj)
            }
            else if (HMM_type == "i3") {
                hmm.infercnv_obj <- i3HMM_assign_HMM_states_to_proxy_expr_vals(hmm.infercnv_obj)
            }
            saveRDS(hmm.infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (!no_plot) {
                plot_cnv(infercnv_obj = hmm.infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_HMM_preds.repr_intensities", 
                    step_count), output_filename = sprintf("infercnv.%02d_HMM_pred%s.Pnorm_%g.repr_intensities", 
                    step_count, hmm_resume_file_token, BayesMaxPNormal), 
                  output_format = output_format, write_expr_matrix = TRUE, 
                  x.center = 1, x.range = c(-1, 3), png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (mask_nonDE_genes) {
        if (!has_reference_cells(infercnv_obj)) {
            stop("Error, cannot mask non-DE genes when there are no normal references set")
        }
        flog.info(sprintf("\n\n\tSTEP %02d: Identify and mask non-DE genes\n", 
            step_count))
        if (skip_past < step_count) {
            infercnv_obj <- mask_non_DE_genes_basic(infercnv_obj, 
                p_val_thresh = mask_nonDE_pval, test.use = test.use, 
                center_val = mean(infercnv_obj@expr.data), require_DE_all_normals = require_DE_all_normals)
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (plot_steps) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, title = sprintf("%02d_mask_nonDE", 
                    step_count), output_filename = sprintf("infercnv.%02d_mask_nonDE", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    if (up_to_step == step_count) {
        flog.info("Reached up_to_step")
        return(infercnv_obj)
    }
    step_count = step_count + 1
    if (denoise) {
        flog.info(sprintf("\n\n\tSTEP %02d: Denoising\n", step_count))
        if (skip_past < step_count) {
            if (!is.na(noise_filter)) {
                if (noise_filter > 0) {
                  flog.info(paste("::process_data:Remove noise, noise threshold at: ", 
                    noise_filter))
                  infercnv_obj <- clear_noise(infercnv_obj, threshold = noise_filter, 
                    noise_logistic = noise_logistic)
                }
                else {
                }
            }
            else {
                flog.info(paste("::process_data:Remove noise, noise threshold defined via ref mean sd_amplifier: ", 
                  sd_amplifier))
                infercnv_obj <- clear_noise_via_ref_mean_sd(infercnv_obj, 
                  sd_amplifier = sd_amplifier, noise_logistic = noise_logistic)
            }
            saveRDS(infercnv_obj, reload_info$expected_file_names[[step_count]])
            invisible(gc())
            if (!no_plot) {
                plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, 
                  cluster_by_groups = cluster_by_groups, cluster_references = cluster_references, 
                  out_dir = out_dir, color_safe_pal = FALSE, 
                  title = sprintf("%02d_denoised", step_count), 
                  output_filename = sprintf("infercnv.%02d_denoised", 
                    step_count), output_format = output_format, 
                  write_expr_matrix = TRUE, png_res = png_res, 
                  useRaster = useRaster)
            }
        }
    }
    saveRDS(infercnv_obj, file = file.path(out_dir, "run.final.infercnv_obj"))
    if (!no_plot) {
        if (is.null(final_scale_limits)) {
            final_scale_limits = "auto"
        }
        if (is.null(final_center_val)) {
            final_center_val = 1
        }
        flog.info("\n\n## Making the final infercnv heatmap ##")
        invisible(gc())
        plot_cnv(infercnv_obj, k_obs_groups = k_obs_groups, cluster_by_groups = cluster_by_groups, 
            cluster_references = cluster_references, out_dir = out_dir, 
            x.center = final_center_val, x.range = final_scale_limits, 
            title = "inferCNV", output_filename = "infercnv", 
            output_format = output_format, write_expr_matrix = TRUE, 
            png_res = png_res, useRaster = useRaster)
    }
    return(infercnv_obj)
}
<bytecode: 0x55f2a5119b30>
<environment: namespace:infercnv>
 --- function search by body ---
Function run in namespace infercnv has this body.
 ----------- END OF FAILURE REPORT -------------- 
Error in HMM && HMM_type == "i6" : 
  'length(x) = 2 > 1' in coercion to 'logical(1)'
Calls: <Anonymous>
Execution halted
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ...
  Running ‘testthat.R’
 OK
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in ‘inst/doc’ ... OK
* checking running R code from vignettes ...
  ‘inferCNV.Rmd’using ‘UTF-8’... OK
 NONE
* checking re-building of vignette outputs ... OK
* checking PDF version of manual ... OK
* DONE

Status: 1 ERROR, 2 NOTEs
See
  ‘/tmp/infercnv.Rcheck/00check.log’
for details.

[GeorgescuC]

Hi @HenrikBengtsson ,

Thanks for reporting this. This issue is fixed in the latest commit and I will be pushing an update of the package when another update is also ready as it is not in this case a real issue. The error was simply because the HMM_type=c('i6', 'i3') argument was compared in HMM && HMM_type == 'i6' && smooth_method == 'coordinate' before running HMM_type = match.arg(HMM_type), so when the argument was not set, HMM_type was not collapsed yet to the first value "i6" (which is the one being evaluated when the flag is not set).

Regards, Christophe.