sffdr Package Vignette

Introduction

The sffdr package implements the surrogate functional false discovery rate (sfFDR) procedure. This methodology integrates GWAS summary statistics from related traits (i.e., pleiotropy) to increase statistical power within the functional FDR framework. The inputs into sffdr are a set of p-values from a GWAS of interest and a set of p-values from one or many informative GWAS.

The significance quantities estimated by sffdr can be used for a variety of analyses:

  • Functional p-value (fp): Controls the type I error rate and can be used for standard significance analyses.
  • Functional q-value (fq): A measure of significance in terms of the positive FDR (closely related to FDR).
  • Functional local FDR (flfdr): A posterior error probability that is useful for functional fine-mapping and assessing significance of a SNP.

Citing this package

The methods implemented in this package are described in:

Bass AJ, Wallace C. Exploiting pleiotropy to enhance variant discovery with functional false discovery rates. Nature Computational Science; 2025.

Note that this work is an extension of the functional FDR methodology and the software builds on some of the functions in the fFDR package found at https://github.com/StoreyLab/fFDR.

Getting help

To report any bugs or issues related to usage please report it on GitHub at https://github.com/ajbass/sffdr.

Installation

You can install the development version of sffdr from GitHub:

# install development version of package
install.packages("devtools")
devtools::install_github("ajbass/sffdr")

Quick start guide

To demonstrate the package, we will use a sample dataset containing 10,000 SNPs for body mass index (BMI) as the primary trait of interest, with body fat percentage (BFP), cholesterol, and triglycerides as informative traits.

library(sffdr)
data(bmi)

# Define primary p-values and informative p-values
p <- sumstats$bmi
z <- as.matrix(sumstats[,-1])

head(sumstats)
#>         bmi      bfp cholesterol triglycerides
#> 1 0.4959680 0.313475   0.0394766      0.148379
#> 2 0.6937320 0.565400   0.8986590      0.298349
#> 3 0.0843257 0.447650   0.7301960      0.715846
#> 4 0.2623540 0.911920   0.1753410      0.775170
#> 5 0.2150660 0.668330   0.0526082      0.603160
#> 6 0.9501450 0.607349   0.2759850      0.867493

1. Modeling the functional proportion of null tests \(\pi_{0}(z)\)

The first step is to model the relationship between the functional proportion of null tests and the informative traits. We use pi0_model to create a model and fpi0est to estimate \(\pi_{0}(z)\)​.

# Create model
mpi0 <- pi0_model(z)
#> ==================================================
#> Building Pi0 Model (Rank: Signal -> 0.0)
#>   Variables        : 3
#>   Min Discoveries  : 150 (Raw SNPs)
#> --------------------------------------------------
#>   [Auto-Scale]: N = 10,000. Bandwidth set to 200 independent SNPs per knot.
#>   [bfp]: Found 220 discoveries (FDR < 0.25). Calculating knots...
#>   [bfp]: Final Model -> Spline (2 knots: 0.022, 0.1)
#>   [cholesterol]: Fallback -> 1 Knot at 0.1000
#>   [triglycerides]: Found 161 discoveries (FDR < 0.25). Calculating knots...
#>   [triglycerides]: Final Model -> Spline (2 knots: 0.0161, 0.1)
#> --------------------------------------------------
#> ~splines::ns(bfp, knots = c(0.022, 0.1), Boundary.knots = c(0, 1)) + splines::ns(cholesterol, knots = c(0.1), Boundary.knots = c(0, 1)) + splines::ns(triglycerides, knots = c(0.0161, 0.1), Boundary.knots = c(0, 1))
#> ==================================================

# Estimation of functional pi0
fpi0 <- fpi0est(p, mpi0)
#> ==================================================
#> Estimating Functional Pi0
#> ==================================================
#>   Fitting models across 19 lambda values...
#>   |                                                       |                                               |   0%  |                                                       |==                                             |   5%  |                                                       |=====                                          |  11%  |                                                       |=======                                        |  16%  |                                                       |==========                                     |  21%  |                                                       |============                                   |  26%  |                                                       |===============                                |  32%  |                                                       |=================                              |  37%  |                                                       |====================                           |  42%  |                                                       |======================                         |  47%  |                                                       |=========================                      |  53%  |                                                       |===========================                    |  58%  |                                                       |==============================                 |  63%  |                                                       |================================               |  68%  |                                                       |===================================            |  74%  |                                                       |=====================================          |  79%  |                                                       |========================================       |  84%  |                                                       |==========================================     |  89%  |                                                       |=============================================  |  95%  |                                                       |===============================================| 100%                                                         
#>   Selecting optimal lambda via MISE...
#>   Selected lambda: 0.25
#>   Generating predictions...
#>   Done.
#> ==================================================
#>

Note on Knots: Ideally, knots should cover the distribution of non-null p-values in your informative studies. Our algorithm uses an FDR threshold to choose the location and number of knots. Depending on the study, you may want to adjust these values (see ?pi0_model). Note that this toy example contains only 10,000 p-values, so the knot locations are not optimal since the function was designed for GWAS-scale data.

2. Estimating the functional significance quantities

With the functional \(\pi_{0}\)​ estimated, we can now calculate the FDR quantities and functional p-values.

# Apply sfFDR
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)   
#> ==================================================
#> Running sfFDR Analysis (v1.1.4)
#> ==================================================
#>   Transforming surrogate variable...
#>   Fitting kernel density...
#>     High-Res Tcub (nn=0.05, h=0.10772)... ✓
#>   Enforcing conditional monotonicity using the 'min' method...
#>     Dividing data into 10 conditional surrogate bins...
#>   Computing functional local FDRs...
#>   Computing functional p-values and q-values...
#>   Done.
#> ==================================================

# Extract results
results <- data.frame(
  p_value = p,
  fp_value = sffdr_out$fpvalues,
  fq_value = sffdr_out$fqvalues,
  flfdr = sffdr_out$flfdr
)

# Plot significance results (focusing on high-significance SNPs)
plot(sffdr_out, rng = c(0, 1e-6))

3. Incorporating linkage disequilibrium (LD)

In standard GWAS, SNPs are often correlated due to LD. To account for this, you should provide a vector indicating which SNPs are LD-independent (e.g., via pruning). sffdr uses these independent SNPs to fit the model, then computes the significance quantities for the entire dataset.

# logical vector: TRUE if SNP is LD-independent
# (In this example, all SNPs are already independent)
indep_snps <- rep(TRUE, length(p))

# Fit model using LD-independent subset
mpi0 <- pi0_model(z = z, indep_snps = indep_snps)
#> ==================================================
#> Building Pi0 Model (Rank: Signal -> 0.0)
#>   Variables        : 3
#>   Min Discoveries  : 150 (Effective Independent)
#> --------------------------------------------------
#>   [Auto-Scale]: N = 10,000. Bandwidth set to 200 independent SNPs per knot.
#>   [bfp]: Found 220 independent discoveries (FDR < 0.25). Calculating knots...
#>   [bfp]: Final Model -> Spline (2 knots: 0.022, 0.1)
#>   [cholesterol]: Fallback -> 1 Knot at 0.1000
#>   [triglycerides]: Found 161 independent discoveries (FDR < 0.25). Calculating knots...
#>   [triglycerides]: Final Model -> Spline (2 knots: 0.0161, 0.1)
#> --------------------------------------------------
#> ~splines::ns(bfp, knots = c(0.022, 0.1), Boundary.knots = c(0, 1)) + splines::ns(cholesterol, knots = c(0.1), Boundary.knots = c(0, 1)) + splines::ns(triglycerides, knots = c(0.0161, 0.1), Boundary.knots = c(0, 1))
#> ==================================================

fpi0 <- fpi0est(p, mpi0, indep_snps = indep_snps)
#> ==================================================
#> Estimating Functional Pi0
#> ==================================================
#>   Fitting models across 19 lambda values...
#>   |                                                       |                                               |   0%  |                                                       |==                                             |   5%  |                                                       |=====                                          |  11%  |                                                       |=======                                        |  16%  |                                                       |==========                                     |  21%  |                                                       |============                                   |  26%  |                                                       |===============                                |  32%  |                                                       |=================                              |  37%  |                                                       |====================                           |  42%  |                                                       |======================                         |  47%  |                                                       |=========================                      |  53%  |                                                       |===========================                    |  58%  |                                                       |==============================                 |  63%  |                                                       |================================               |  68%  |                                                       |===================================            |  74%  |                                                       |=====================================          |  79%  |                                                       |========================================       |  84%  |                                                       |==========================================     |  89%  |                                                       |=============================================  |  95%  |                                                       |===============================================| 100%                                                         
#>   Selecting optimal lambda via MISE...
#>   Selected lambda: 0.25
#>   Generating predictions...
#>   Done.
#> ==================================================
#> 

# Estimate quantities for all SNPs
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)
#> ==================================================
#> Running sfFDR Analysis (v1.1.4)
#> ==================================================
#>   Transforming surrogate variable...
#>   Fitting kernel density...
#>     High-Res Tcub (nn=0.05, h=0.10772)... ✓
#>   Enforcing conditional monotonicity using the 'min' method...
#>     Dividing data into 10 conditional surrogate bins...
#>   Computing functional local FDRs...
#>   Computing functional p-values and q-values...
#>   Done.
#> ==================================================

Note on perfect-LD copies and scale: SNPs with identical p-values are copies of one signal (perfect LD); sffdr() fits the density on one point per set of copies and evaluates it for every SNP (ties from rounded p-values, which sit almost only in the bulk of the distribution, are left alone). The functional p-values are calibrated in the tail when the fitted tail holds tens of thousands of SNPs that are not copies of each other, which is the case at GWAS scale; at 100,000 SNPs in heavy LD the tail estimate is noisy.

Note on choosing independent SNPs: There are two strategies: (1) prune using plink (e.g., --indep-pairwise 200kb 1 0.3) and use the pruned set of SNPs as the independent set, or (2) use informed sampling using the informative traits (e.g., clumping or selecting a representative SNP in an LD block with the smallest informative p-value). The latter is likely better when the primary study has low power.

Advanced usage

Enforcing monotonicity

The kernel density estimator is adaptive. However, sometimes artifacts—such as LD—can induce spurious “bumps” in the null region which can impact density estimation.

By default, sffdr() enforces a monotone constraint on the density within surrogate bins (bins of ~1,000 SNPs) to stabilize the joint density:

sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0, monotone = TRUE)   # the default; monotone = FALSE turns it off

While this provides robustness in the null region, users should note it may result in a loss of power in specific cases.

Incorporating LD scores

If LD scores (\(l_{j}\)) are available (e.g., from LDSC regression), there are two ways to use them.

  1. As a covariate. SNPs with larger LD scores tag more variation and are more often associated, so the LD score is informative about association and belongs in the functional \(\pi_{0}\) model alongside the informative traits:
# LD score (log scale) added to the pi0 model built by pi0_model()
zt <- cbind(mpi0$zt, ld = log(ld_score))
fmod <- update(mpi0$fmod, ~ . + ns(ld, df = 3))
fpi0 <- fpi0est(p, z = zt, pi0_model = fmod, indep_snps = indep_snps)
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)
  1. As inverse weights in the \(\pi_{0}\) fit (\(w_j = 1/l_{j}\)), in place of LD pruning:
w <- 1 / ld_score
mpi0 <- pi0_model(z, weights = w)
fpi0 <- fpi0est(p, mpi0, weights = w)
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)

sffdr() itself does not use weights: its density is estimated over all tested SNPs, which is the set the functional p-values count errors in, and reweighting that density biases the local FDR (SNPs in LD are each tested, so each counts).

Note on incorporating LD scores: As LD scores are available for common reference panels (e.g., 1000 Genomes), they can be easily downloaded to use in sfFDR. If you have LD scores calculated from your data then that’s even better. However, it is problematic to use LD scores when the reference panel is not matched to your study population. In that case, the pruning (or informative sampling) strategy discussed above is the best option.

Sample overlap with the informative studies

If the primary and informative studies share participants (e.g., traits measured in the same biobank), pass the signed informative z-scores through decorrelate_informative() and use its p-values ($p2) as the informative p-values in pi0_model(); the primary p-values are used unchanged (see ?decorrelate_informative).

Manual Knot Selection and using other informative variables

The selection of knots in the B-spline basis is an important tuning parameter in sffdr. If your surrogate GWAS has extremely high power (very small p-values), we recommend placing more knots at the lower end of the p-value distribution to capture the local density of signals.

# Create model (can include other variables (e.g., MAF) or specify more complicated models)
fmod <- "~ns(bfp, knots = c(0.01, 0.025, 0.05, 0.1))"
fpi0_mod <- fpi0est(p = p,
                    z = mpi0$zt,
                    pi0_model = fmod)

Note on incorporating additional variables: Other informative covariates can be included in the model (e.g., MAF, LD score, etc.). While the LD score annotations can be used (e.g., the baselineLD model), we have found that using the expected \(\chi^{2}\) given the annotations performs better in practice and simplifies the model fitting.