Error in mutate(., state = run_joint_hmm(pAD = pAD, DP = DP, p_s = p_s,

[pcream]

Hello, I am having some issues in performing the run_numbat() function once finished the initial SNP calling/phasing. This is the output after running it using the v 1.3.0 Docker file pulled and run using apptainer on a computing cluster:

INFO [2023-04-14 14:31:39] Numbat version: 1.3.0 Running under parameters: t = 1e-05 alpha = 1e-04 gamma = 20 min_cells = 50 init_k = 3 max_cost = 3358.8 n_cut = 0 max_iter = 2 max_nni = 100 min_depth = 0 use_loh = auto segs_loh = None call_clonal_loh = FALSE segs_consensus_fix = None multi_allelic = TRUE min_LLR = 5 min_overlap = 0.45 max_entropy = 0.5 skip_nj = TRUE diploid_chroms = None ncores = 12 ncores_nni = 12 common_diploid = TRUE tau = 0.3 check_convergence = FALSE plot = TRUE genome = hg38 Input metrics: 11196 cells INFO [2023-04-14 14:31:40] Mem used: 4.61Gb INFO [2023-04-14 14:31:51] Approximating initial clusters using smoothed expression .. INFO [2023-04-14 14:31:52] Mem used: 4.61Gb INFO [2023-04-14 14:39:24] running hclust... INFO [2023-04-14 14:50:34] Iteration 1 INFO [2023-04-14 14:50:35] Mem used: 16.1Gb INFO [2023-04-14 14:50:45] Expression noise level (MSE): medium (1). INFO [2023-04-14 14:50:45] Running HMMs on 5 cell groups.. INFO [2023-04-14 14:50:54] quadruploid state enabled INFO [2023-04-14 14:50:54] diploid regions: 1a,5a,5c,9a,11a,12a,13a,15a,15c,16a,17a,19a,20a,21a,22a ERROR [2023-04-14 14:50:55] job 1,2,3,4,5 failed ERROR [2023-04-14 14:50:55] Error in mutate(., state = run_joint_hmm(pAD = pAD, DP = DP, p_s = p_s, : ℹ In argument: state = run_joint_hmm(...). ℹ In group 1: CHROM = 1. Caused by error in h(): ! error in evaluating the argument 'x' in selecting a method for function 'rowSums': all x must be integers

The stdout on the slurm scheduler is different interestingly, but similar error:

Numbat version: 1.3.0 Running under parameters: t = 1e-05 alpha = 1e-04 gamma = 20 min_cells = 50 init_k = 3 max_cost = 3358.8 n_cut = 0 max_iter = 2 max_nni = 100 min_depth = 0 use_loh = auto segs_loh = None call_clonal_loh = FALSE segs_consensus_fix = None multi_allelic = TRUE min_LLR = 5 min_overlap = 0.45 max_entropy = 0.5 skip_nj = TRUE diploid_chroms = None ncores = 12 ncores_nni = 1 common_diploid = TRUE tau = 0.3 check_convergence = FALSE plot = TRUE genome = hg38 Input metrics: 11196 cells Mem used: 4.61Gb Approximating initial clusters using smoothed expression .. Mem used: 4.61Gb number of genes left: 10823 running hclust... Iteration 1 Mem used: 16.1Gb Expression noise level (MSE): medium (1). Running HMMs on 5 cell groups.. Error in mutate(., state = run_joint_hmm(pAD = pAD, DP = DP, p_s = p_s, : ℹ In argument: state = run_joint_hmm(...). ℹ In group 1: CHROM = 1. Caused by error in h(): ! error in evaluating the argument 'x' in selecting a method for function 'rowSums': all x must be integers

Error in group_by(): ! Must group by variables found in .data. Column seg is not found. Column sample is not found. Backtrace: ▆

1. ├─numbat::run_numbat(...)
2. │ └─bulk_subtrees %>% ...
3. ├─numbat:::run_group_hmms(...)
4. │ └─... %>% ungroup()
5. ├─dplyr::ungroup(.)
6. ├─dplyr::mutate(., seg_start_index = min(snp_index), seg_end_index = max(snp_index))
7. ├─dplyr::group_by(., seg, sample)
8. └─dplyr:::group_by.data.frame(., seg, sample)
9. └─dplyr::group_by_prepare(.data, ..., .add = .add, error_call = current_env())
   10. └─rlang::abort(bullets, call = error_call) Warning message: In mclapply(bulks %>% split(.$sample), mc.cores = ncores, function(bulk) { : all scheduled cores encountered errors in user code Execution halted

I assume there is some issue in rowSums attempting to perform on non-integers (the 'CHROM = 1" maybe?), but I can't tell when or where it is occurring besides during the run_joint_hmm() function, but looking through the source I can't find what that function is doing in exactly. I have double checked that the denseMatrix, vcf allele file, and reference are all formatted correctly. Is it also possible this is a multithread issue? I assume the mclapply is failing because of the earlier error, but I am not sure. I am also confused on the quadraploid callout, as that's not expected with this dataset. I have also tried this on my own computer using both Windows and Unix R installations and run into the same problem, so I don't think it is a slurm/cluster related problem. I've attached the files used and the simple R script to run the function to see if this can be replicated. Thank you!

Numbat_error.zip

[teng-gao]

Hi @pcream, could you attach the count_mat as well (I think only the allele dat and expression reference are included). Also, are you able to run the ATC2 example?

[pcream]

Apologies, I uploaded the wrong files! Here is a link to the correct ones: link. The expression matrix is too big to upload here.

AGC_allele is the SNP file, test_assay is the dense expression matrix, and numbat_reference is a reference I complied from a publicly available human hematopoietic single cell dataset. I also tried it with the built in ref_hca, but that also failed. I'll run the ATC2 example provided in the github vignette and get back to you.

Thanks!

[pcream]

Update, the ATC2 example was successfully run, I've attached the CNA panel here. Does this match the normal output of Numbat for the example? It looks quite different from the CopyKat CNA analysis run in the original paper, but maybe the Numbat output is quite different from CopyKat. I can link the whole output folder if that would help with figuring out the error.

[teng-gao]

Hi @pcream ,

The plot looks right. The difference in CNV profile is due to diploid baseline estimation - copyKat does not explicitly try to identify the diploid state and always uses the copy number median as baseline. This can be confusing when the tumor is hyper-diploid or hypo-diploid (for more details please check the paper). I was able to replicate the error with your input data, so the problem is caused by the specific input (haven't gotten a chance to look into why yet).

[pcream]

Hello @teng-gao,

Great, makes sense! Just wanted to make sure that was the expected output based on default settings. I tried another sample that was run and sequenced in parallel, using the same processing pipeline as this sample and I get the same error "mutate(., state = run_joint_hmm" etc. I think I will try downloading the ATC2 raw fastqs and processing everything to see if perhaps there is something wrong with my cellranger/pileup processing.

[teng-gao]

Hi @pcream ,

I checked your input again and the problem was that your count_mat was fractional instead of integer. It should be raw counts not normalized expression values.

[pcream]

Hi @teng-gao,

Many thanks!!! I can't believe I didn't notice that, thank you for your help!