Error in mtm[, "index1"] : incorrect number of dimensions

[LouisLeNezet]

Hi,

I'm trying to run a minimal test dataset for QUILT2 in nf-core.
However I ran into the following error when using a small region.

    QUILT2.R \
        --bamlist=$BAM_LIST \
        --genetic_map_file=map.gz \
         \
         \
        --chr=chr22 \
        --regionStart=16570000 \
        --regionEnd=16610000 \
        --nGen=100 \
        --buffer=10000 \
        --nCores=1 \
        --outputdir="." \
        --reference_vcf_file=1000GP.chr22.vcf.gz \
        --seed=1
    
    
    cat <<-END_VERSIONS > versions.yml
    "QUILT_QUILT2":
        r-base: $(Rscript -e "cat(strsplit(R.version[['version.string']], ' ')[[1]][3])")
        r-quilt2: $(Rscript -e "cat(as.character(utils::packageVersion(\"QUILT2\")))")
    END_VERSIONS
  
  Command exit status:
    1
  
  Command output:
    (empty)
  
  Command error:
    [2025-04-29 14:47:08] Using reference information from:1000GP.chr22.vcf.gz
    [2025-04-29 14:47:08] Using strategy of first imputing using common SNPs and then using all SNPs, with allele frequency threshold:0.001
    [2025-04-29 14:47:08] Begin get sites and haplotypes from reference vcf
    [2025-04-29 14:47:09] End get sites and haplotypes from reference vcf
    [2025-04-29 14:47:09] There were 0 skipped variants when processing the reference VCF (not bi-allelic, not a SNP, not unique position
    [2025-04-29 14:47:09] There are 0 common and 0 rare (0 total) variants in the left buffer region 16560000 <= position < 16570000
    [2025-04-29 14:47:09] There are 419 common and 484 rare (903 total) variants in the central region 16570000 <= position <= 16610000
    [2025-04-29 14:47:09] There are 0 common and 0 rare (0 total) variants in the right buffer region 16610000 < position <= 16620000
    [2025-04-29 14:47:09] Getting and validating genetic map
    [2025-04-29 14:47:14] Match genetic map to desired imputation positions
    [2025-04-29 14:47:14] Done with getting genetic map
    [2025-04-29 14:47:18] There are 0 regions out of 13 below minimum recombination rate, setting them to minimum rate
    [2025-04-29 14:47:18] There are 0 regions out of 13 above maximum recombination rate, setting them to maximum rate
    [2025-04-29 14:47:18] Using nMaxDH = 77
    [2025-04-29 14:47:18] Build mspbwt indices
    [2025-04-29 14:47:18] Done building mspbwt indices
    [2025-04-29 14:47:18] Save converted reference haplotypes
    [2025-04-29 14:47:18] Done converting reference haplotypes
    [2025-04-29 14:47:18] Done processing reference information
    [2025-04-29 14:47:18] Get BAM sample names
    [2025-04-29 14:47:18] Done getting BAM sample names
    [2025-04-29 14:47:18] There are 419 common SNPs and 903 SNPs overall in this region
    [2025-04-29 14:47:19] There are 0 regions out of 28 below minimum recombination rate, setting them to minimum rate
    [2025-04-29 14:47:19] There are 0 regions out of 28 above maximum recombination rate, setting them to maximum rate
    [2025-04-29 14:47:19] Imputing sample: 1
    [2025-04-29 14:47:19] The average depth of this sample is:1.004046141607
    [2025-04-29 14:47:19] There are 137 reads under consideration
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 1 small gibbs
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 1 full
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 2 small gibbs
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 2 full
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 3 small gibbs
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 3 full
    [2025-04-29 14:47:19] i_gibbs=1, i_it = 3 small gibbs (all SNPs)
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 1 small gibbs
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 1 full
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 2 small gibbs
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 2 full
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 3 small gibbs
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 3 full
    [2025-04-29 14:47:19] i_gibbs=2, i_it = 3 small gibbs (all SNPs)
    [2025-04-29 14:47:19] i_gibbs=3, i_it = 1 small gibbs
    [2025-04-29 14:47:19] i_gibbs=3, i_it = 1 full
    [2025-04-29 14:47:19] i_gibbs=3, i_it = 2 small gibbs
    [2025-04-29 14:47:19] i_gibbs=3, i_it = 2 full
    [2025-04-29 14:47:19] i_gibbs=3, i_it = 3 small gibbs
    [2025-04-29 14:47:19] i_gibbs=3, i_it = 3 full
    Error in mtm[, "index1"] : incorrect number of dimensions
    Calls: QUILT ... select_new_haps_mspbwt_v3 -> lapply -> FUN -> lapply -> FUN -> diff
    Execution halted

Here is the link to the PR.
Strangely it does work if I use QUILT.R instead of QUILT2.R so this might be due to a change between the two version.

Do you know what I should do ?

[Zilong-Li]

Hi @LouisLeNezet ,

Thank you for bringing QUILT2 into the nf-core! QUILT2 is tailored for large reference panel, though it should also be able to run for small reference panel. But it seems there are few edge cases it can not handle well, e.g. there are only 903 snps in the genomic chunk from your example.

I'll look into the issue ASAP, and will make a new minor release once we fixed it.

Best,
Zilong.

[Zilong-Li]

Hey @LouisLeNezet ,

This issue has been fixed since version v2.0.2. You can install it via pak::pkg_install("rwdavies/QUILT/QUILT") in R. It will take a while to update the bioconda version.

Let me know if you need more assistance in porting it into Nextflow module.

[LouisLeNezet]

Thanks for fixing this !
I will need for bioconda to be updated and I will update the PR.

[LouisLeNezet]

As the bioconda recipe was failing I've open a new PR to fix it
bioconda/bioconda-recipes#55903

[Zilong-Li]

Thanks! I think we need to fix the bioconda recipe of r-stitch v1.8.1 first, which now depends on the vcfppR package on CRAN that is not in the bioconda yet. I was submitting it to conda forge PR based on the guideline. But now they suggest submitting it to bioconda.

[LouisLeNezet]

Yeah I think it would make sense for stitch to be in bioconda.
Could you ping me, when it will be ready ?

[Zilong-Li]

sure will do. btw, the r-vcfppr is pending review bioconda/bioconda-recipes#55905

[Zilong-Li]

Hi @LouisLeNezet , Sorry for late response. I get r-stitch v1.8.2 updated in bioconda. Will you help fix the recipe of r-quilt v2.0.2 as you previously did? Thanks.

[LouisLeNezet]

Hi @Zilong-Li, I've tried to update it here
bioconda/bioconda-recipes#56700

But the version of r-rcpp is not compatible.
r-mspbwt expect <=1.0.12 and r-quilt >= 1.0.14

Is it possible for r-mspbwt to accept r-rcpp 1.0.14 ?

[Zilong-Li]

uh, annoying dependencies. I just made a new tag release of mspbwt to trigger the autobump in bioconda. It passed the checks bioconda/bioconda-recipes#56708. Will you @LouisLeNezet be able to edit the meta.yaml file there to relax the r-rcpp version before they merge it? Otherwise we need to open another PR.

[LouisLeNezet]

I did and it worked.
I was just wondering if the upper boundery in r-mspwt for r-rcpp was there for a reason.

[Zilong-Li]

I think it is related to a bug in Rcpp 1.0.13 RcppCore/Rcpp#1334. Now quilt2p should work fine with Rcpp 1.0.14

[LouisLeNezet]

Okay, just need now to wait for the recipe to be merged !