Question about SNP-based p-values: STAARpipeline VS REGENIE

[Jun-Jun-Xu]

Hi,

I have a question regarding SNP-based association results from STAARpipeline.

I performed GWAS using:

The same genotype data (MAC ≥ 20)

The same phenotype

The same covariates

I ran:

REGENIE (step 1 + step 2)

STAARpipeline_Individual_Analysis.r from STAARpipeline

When applying the conventional genome-wide significance threshold (5 × 10⁻⁸), I observed a large discrepancy:

REGENIE detected 32 significant variants

STAARpipeline detected 25,848 significant variants

Most of the additional variants detected by STAARpipeline are rare variants. Although the significant variants from REGENIE are also mostly rare, STAARpipeline seems to identify many more rare signals.

Summary statistics of p-values from STAARpipeline results (column 6):

min: 3.98e-05
median: 4.93e-05
max: 0.0657

Based on my understanding, STAARpipeline (in the individual SNP analysis) does not apply variant-level weighting in the same way as gene-based STAAR tests. Therefore, I am concerned about whether these additional rare variant signals are reliable.

My questions are:

1. Is it expected that STAARpipeline detects substantially more rare variant signals than REGENIE under the same MAC filter (MAC ≥ 20)?

2. Are these additional rare variant associations likely to be inflated or due to model differences?

3. Would it be more appropriate to restrict interpretation to variants with MAF ≥ 0.01?

Any guidance would be greatly appreciated.

Thank you!

[yuxinyuanqt]

Hi @Jun-Jun-Xu ,

Both REGENIE and STAARpipeline are well-established methods. In principle, single-variant analyses should not produce such dramatic differences under the same data and model settings. The discrepancy is therefore more likely due to differences in the analysis workflow rather than the methods themselves.

To help identify the issue, we would need more details:

1. Genotype data
- Did REGENIE and STAARpipeline use exactly the same genotype data (e.g., BGEN vs aGDS converted from the same source)?

2. Phenotype and sample size
- What is the sample size in each analysis?
- Is the trait continuous or binary?
- If continuous, was inverse normal transformation applied consistently?
- If binary and imbalanced, was use_SPA = TRUE specified when fitting the STAAR null model?

3. Number of variants tested
- Are the total numbers of analyzed variants comparable between the two methods?

4. Covariates
- Were exactly the same covariates included (e.g., age, age2, sex, PCs)?

5. Genetic relatedness control
- How was relatedness handled?
- Was a GRM specified when fitting the STAAR null model?

6. P-value summary

You reported:

 ```
 min: 3.98e-05  
 median: 4.93e-05  
 max: 0.0657
 ```
 This is somewhat confusing — if the minimum p-value is 3.98e-05, then none of the variants reach genome-wide significance (5×10⁻⁸). Could you clarify whether this summary corresponds to all variants or only a subset?

Overall, the reported pattern suggests that something in the workflow may be inconsistent, and further clarification would help pinpoint the cause.

[Jun-Jun-Xu]

Thank you for your comments and questions. Please find clarifications regarding our analysis workflow below:

1. Genotype data: Both REGENIE and STAARpipeline analyses used the same genotype (ukb24308)

2. The trait analyzed is binary, with 2,424 cases and 242,419 controls.

For STAARpipeline, the null model was fit using the following code:

obj_nullmodel <- fit_nullmodel(
cc ~ age + age2 + sex + PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 + PC10,
data = phenotype,
kins = sgrm,
use_sparse = TRUE,
kins_cutoff = 0.022,
id = "name",
family = binomial(link = "logit"),
verbose = TRUE
)

3. Both methods analyzed the same set of variants derived from the above genotype files.

4. For STAARpipeline, the covariates included: age, age2, sex, PC1–PC10. For REGENIE, the covariates included: age, sex, PC1–PC10 (i.e., age² was not included).

5. Individuals with relatedness were excluded.

In the STAARpipeline null model, we also used a GRM matrix (ukb22418_all_b0_v2_maf001_mac100_mis01_hwe_lo38_sort_pruned_sparseGRM.sGRM.RData) to account for residual relatedness.

6. Apologies for the confusion in the previous message. The values reported earlier (min: 3.98e-05, median: 4.93e-05, max: 0.0657) were not p-values but minor allele frequencies (MAFs).

Please let us know if further clarification or additional details are needed.

[yuxinyuanqt]

Hi @Jun-Jun-Xu ,

Thank you for the clarification.

Since your trait is a highly imbalanced case-control phenotype (2,424 cases vs. 242,419 controls), special care is required. For rare variants, score tests relying on normal approximation can inflate type I error in such settings.

Therefore, SPA (Saddlepoint Approximation) or Firth correction is recommended. In STAARpipeline, you should specify:
use_SPA = TRUE
when fitting the null model. The default is use_SPA = FALSE, and without SPA correction, rare variant p-values may be inflated in imbalanced case-control analyses.

Additionally, in Individual_Analysis, consider increasing the minimum MAC threshold (default = 20). Raising it to MAC ≥ 80 or 100 can improve numerical stability, reduce false positives among rare variants, and improve computational efficiency.

Overall, the excess rare variant signals in STAARpipeline are most likely due to the absence of SPA correction under this highly imbalanced design.