saigegit / SAIGE

Development for SAIGE and SAIGE-GENE(+)
GNU General Public License v3.0
64 stars 28 forks source link

Step 2 - get_newPhi_scaleFactor error in SAIGE-GENE+ depending on annotations #12

Closed tammy-dg closed 1 year ago

tammy-dg commented 2 years ago

Hello all,

I am running into an error on certain annotation sets:

Error in if (abs(q - m1)/sqrt(var1) < Cutoff) { : 
  missing value where TRUE/FALSE needed
Calls: SPAGMMATtest ... SAIGE.Region -> get_newPhi_scaleFactor -> <Anonymous>
Execution halted

Specifically, this error arises at different genes depending on what is passed into --annotation_in_groupTest. For example, when I input --annotation_in_groupTest=PTV_f5,SPDMG_f5,MSDMG_f5 the traceback occurs when testing one particular gene:

[1] "Analyzing Region ACTN4 (12/1451)."
Start analyzing chunk 0.....
In chunks 0-0, 98 markers are ultra-rare and 11 markers are not ultra-rare.
Analyzing chunks (0/1, 0/1)........
Analyzing chunks (1/1, 0/1)........
Analyzing chunks (1/1, 1/1)........
Error in if (abs(q - m1)/sqrt(var1) < Cutoff) { : 
  missing value where TRUE/FALSE needed
Calls: SPAGMMATtest ... SAIGE.Region -> get_newPhi_scaleFactor -> <Anonymous>
Execution halted

However ,when I modify the order of the annotations to test, for example: input --annotation_in_groupTest=MSDMG_f5,SPDMG_f5,PTV_f5, the traceback occurs when testing a different gene:

[1] "Analyzing Region ADM5 (21/1451)."
Error in if (abs(q - m1)/sqrt(var1) < Cutoff) { : 
  missing value where TRUE/FALSE needed
Calls: SPAGMMATtest ... SAIGE.Region -> get_newPhi_scaleFactor -> <Anonymous>
Execution halted

I have attached a subset of the group file that is causing this and log files for the two runs.

chr19.MSDMG_SPDMG_PTV.log chr19.PTV_SPDMG_MSDMG.log PTV_f5-SPDMG_f5-MSDMG_f5.chr19.txt

Thank you for your attention!

saigegit commented 2 years ago

Hi @tammy-dg,

Could you please try v1.0.3 to see if this issue persists. I think it might also be due to the annolist issue.

Thanks, Wei

tammy-dg commented 2 years ago

Hi Wei, thanks for the quick response. Unfortunately the same error is still encountered when I pulled v1.0.3 from docker.

saigegit commented 2 years ago

Hi @tammy-dg,

I investigated the issue and found the reason for the error. The error should be only observed when SKAT-O tests are performed, no weights were applied, and one or more groups do not have any marker. All other groups except for the ones having errors are not affected. The issue should have been fixed in v1.0.4. Can you please give it a try?

Thanks, Wei

tammy-dg commented 2 years ago

Hi Wei,

I've pulled v1.0.4 from docker and confirmed that it works now with these examples.

However, I have a question based on your explanation. In the example above where it was encountering the error with ACTN4, I've attached the results from step 2 for this gene. The output shows that there were samples with markers for each group (as denoted by the non-zero MAC columns), so I don't think that meets the criteria of "one or more groups do not have any marker". I've also confirmed separately that at least one variant in each annotation group had a non-zero MAF in the samples included in the analysis for this gene. So based on this, I was wondering if you had any more insight.

But regardless, thank you greatly for your quick work to patch this!

Tammy

ACTN4.PTV_f5-SPDMG_f5-MSDMG_f5.output.txt

saigegit commented 2 years ago

hi @tammy-dg,

I'm so sorry for taking so long to get back to you. I haven't figured out what can potentially cause the problem. Recently, there was a reported issue about reading markers from VCF. We've found that when the markers are not ordered by chromosome position in the group file, genotype reading can be incomplete. I wonder if this is the possible reason for the problem you saw. We have fixed that issue in version 1.0.6.

Thanks, Wei

tammy-dg commented 1 year ago

Hi Wei, sorry for never following up on this. We ended up using SAIGE-GENE+ with PLINK file inputs instead of VCF input so I cannot comment if this ended up fixing it! But thank you again for your time