Crash in step 2 with: mean(): object has no elements

[pettyalex]

Error in mainMarkerInCPP(genoType, traitType, genoIndex_prev, genoIndex,  :
  mean(): object has no elements
Calls: SPAGMMATtest -> SAIGE.Marker -> mainMarkerInCPP
Execution halted

There's only a handful of lines this error could come from, as it's coming from the Armadillo library. We're investigating inputs, but this seems related to ultra-rare variants around line 1369:

    if(indexForNonZero.n_elem > 0){
    //URindVec.push_back(jm+1);
      double altFreq = arma::mean(genoURVec)/2;
      double altCounts = arma::accu(genoURVec);
      double missingRate = 0;
      double imputeInfo = 1;

We recently updated to the latest version, we had been running v1.0.6 before this but we needed to be able to run all of chr1 and chr2 without truncating, so we updated.

[pettyalex]

Oh, this turned out to be coming from the efficient resampling, for a couple of our variants. SKATExactBin_Work was getting called from the efficient resampling path such that p1 was empty: https://github.com/saigegit/SAIGE/blob/main/src/ER_binary_func.cpp#L210

For now we are planning to work around this by setting --max_MAC_for_ER=0 to fully disable the efficient resampling.

[pettyalex]

My logic was backwards when I posted that, to disable ER we actually used --max_MAC_for_ER=10000000.

This is a real problem with efficient resampling where it will fail for certain datasets. What should SKATExactBin_Work return when there's nothing in SAIGEClass.m_mu for an index? I've seen in other places you've hard coded the p-value to return zero. Would that work? I'll open a PR to hopefully find a way to resolve this with you.

[pettyalex]

Huh. Setting that high max_MAC_for_ER seems to make it fail to converge. I wonder what else we could do.

[weizhou0]

Hi @pettyalex,

Thanks for reporting this issue! ER is expected to be used for genetic variants that are very rare with MAC < 4, for which SPA does not work. For single-variant association tests, we usually set --minMAC=10 or 20 to not test those very rare variants. For group test, since SAIGE-GENE+ collapses ultra-rare variants with MAC < 10 together first, there won't be very rare variants that need ER. But I think it would be still nice to figure out why ER failed for some cases. I will investigate and report back.

Thanks, Wei

[pettyalex]

We are running this with --minMAC=5. I haven't been able to create a minimal reproduction case, unfortunately, and have been reproducing it with real data. Are you able to run SAIGE under a debugger? I haven't ever configured Rcpp code to run under a debugger, but that would really help here.

The scientists are asking me if I can produce a version based on 1.0.6, which they had been using until we noticed that it was truncating positions past 200 million, and just update it with this fix: https://github.com/saigegit/SAIGE/commit/94851282a47c91ab707a167094c101d42afdc246

Here's full arguments that we see this failure in, unfortunately the VCFs are not shareable data. We're seeing the failure at different points in each chromosome.

        --vcfFile=/home/alex/src/SAIGE/test_data/chr21_merged.dose.vcf.bgz \
        --sampleFile=/home/alex/src/SAIGE/test_data/DR_sample_file_EUR \
        --chrom=chr21 \
        --minMAC=5 \
        --GMMATmodelFile=/home/alex/src/SAIGE/test_data/DR_EUR_SAIGE.rda \
        --varianceRatioFile=/home/alex/src/SAIGE/test_data/DR_EUR_SAIGE.varianceRatio.txt \
        --SAIGEOutputFile=/home/alex/src/SAIGE/test_data/DR_EUR_SAIGE.chr21.output.txt \
        --is_output_moreDetails=TRUE \
        --LOCO=FALSE

[pettyalex]

When going into the failing path, StdStat is NaN around here: https://github.com/saigegit/SAIGE/blob/main/src/SAIGE_test.cpp#L426

[pettyalex]

It looks like when SPA fails to converge, we end up with:

        t_pval = t_pval_noSPA;
    pval = pval_noadj;

So maybe that's what we should do in the case where ER fails? I see some but not much use of exceptions in the codebase, is this a situation for catching an exception?

[pettyalex]

I'm investigating now, but pulling the latest 1.1.8 version that you updated last night may have fixed this for us? Does that seem possible?

[weizhou0]

Thanks! 1.1.8 actually is not fixing the ER problem and it improves the SPA convergence, while outputting a warning message when the max MAC for ER is set to be greater than 4.

On Mon, Apr 24, 2023 at 1:24 PM Alex Petty @.***> wrote:

I'm investigating now, but pulling the latest 1.1.8 version that you updated last night may have fixed this for us? Does that seem possible?

— Reply to this email directly, view it on GitHub https://github.com/saigegit/SAIGE/issues/88#issuecomment-1520558150, or unsubscribe https://github.com/notifications/unsubscribe-auth/APJSQMRPDBQWRJAF7X6VYJTXC2ZM7ANCNFSM6AAAAAAWUPDXSA . You are receiving this because you commented.Message ID: @.***>

[pettyalex]

I just realized that I had built my own branch that hacks around it by returning a p-value of 0: https://github.com/saigegit/SAIGE/pull/89

[pettyalex]

@weizhou0 It looks like https://github.com/saigegit/SAIGE/commit/97f87f42ece8ea8f594e82b9acff271d7aa661be fixed this? We will test.

[weizhou0]

Hi Alex,

Yes the version fixed this issue. Thanks!

Wei

On Wed, May 17, 2023 at 4:29 PM Alex Petty @.***> wrote:

@weizhou0 https://github.com/weizhou0 It looks like 97f87f4 https://github.com/saigegit/SAIGE/commit/97f87f42ece8ea8f594e82b9acff271d7aa661be fixed this? We will test.

— Reply to this email directly, view it on GitHub https://github.com/saigegit/SAIGE/issues/88#issuecomment-1552032959, or unsubscribe https://github.com/notifications/unsubscribe-auth/APJSQMUP5GBOGOB2XI3GL2DXGUYJVANCNFSM6AAAAAAWUPDXSA . You are receiving this because you were mentioned.Message ID: @.***>

[pettyalex]

Awesome, we have validated that too.

Thanks so much for SAIGE!