| Title: | Surrogate Functional False Discovery Rates for Genome-Wide Association Studies |
|---|---|
| Description: | Pleiotropy-informed significance analysis of genome-wide association studies with surrogate functional false discovery rates (sfFDR). The sfFDR framework adapts the fFDR to leverage informative data from multiple sets of GWAS summary statistics to increase power in study while accommodating for linkage disequilibrium. sfFDR provides estimates of key FDR quantities in a significance analysis such as the functional local FDR and $q$-value, and uses these estimates to derive a functional $p$-value for type I error rate control and a functional local Bayes' factor for post-GWAS analyses (e.g., fine mapping and colocalization). |
| Authors: | Andrew Bass [aut, cre], Chris Wallace [aut] |
| Maintainer: | Andrew Bass <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 1.1.4 |
| Built: | 2026-10-20 19:39:31 UTC |
| Source: | https://github.com/ajbass/sffdr |
A dataset containing a subset of p-values from the UK Biobank.
A data frame with 10,000 rows and 4 columns:
Body mass index
Body fat percentage
Cholesterol
Triglycerides
# Import data data(bmi) # Separate main p-values and informative p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1])# Import data data(bmi) # Separate main p-values and informative p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1])
When the primary study shares samples with the informative studies, their
null z-scores are correlated. Used as covariates in pi0_model,
the informative p-values then predict the primary's null tests, which
inflates false discoveries. decorrelate_informative() removes the
primary's noise from each informative study,
where is the null correlation of study with the primary.
Under the primary's null, is independent of whatever the
informative studies' effects, so the primary's p-values are analysed
unchanged with the decorrelated p-values as covariates.
decorrelate_informative( z_x, z_y, p_x = NULL, p_y = NULL, rho_o = NULL, lfdr_x = NULL, indep_snps = NULL )decorrelate_informative( z_x, z_y, p_x = NULL, p_y = NULL, rho_o = NULL, lfdr_x = NULL, indep_snps = NULL )
z_x |
Numeric vector; the primary study's z-scores. |
z_y |
Numeric vector or matrix; signed z-scores of the informative
studies (one column per study), in the same order as |
p_x, p_y
|
Ignored; kept for compatibility with sffdr 1.1.2. |
rho_o |
Optional overlap correlations of the primary with each
informative study (one number per study). If |
lfdr_x |
Ignored. sffdr 1.1.2 scaled the correction by the primary's local FDR, which left part of the overlap in the covariates. |
indep_snps |
Optional logical or integer vector selecting approximately
independent SNPs, used to estimate |
The correction needs only the overlap correlation of each informative study
with the primary: cross-trait LD score regression intercepts are the
recommended source for GWAS; otherwise they are estimated by
estimate_overlap. Errors in them enter at second order: in
simulations, a correlation of 0.5 given as 0.35 or 0.65 left the FDR within
0.005 of its value at the true correlation.
A list with
decorrelated z-scores of the informative studies (vector or matrix, as z_y)
their p-values, for use as covariates
the overlap correlations used
set.seed(1) n <- 20000 e1 <- rnorm(n) e2 <- 0.5 * e1 + sqrt(1 - 0.5^2) * rnorm(n) z_x <- e1 + ifelse(runif(n) < 0.05, 3, 0) # primary z_y <- e2 + ifelse(runif(n) < 0.1, 3, 0) # informative, same samples dec <- decorrelate_informative(z_x, z_y, rho_o = 0.5) p_x <- 2 * pnorm(-abs(z_x)) mpi0 <- pi0_model(as.matrix(dec$p2)) fpi0 <- fpi0est(p_x, pi0_model_obj = mpi0) res <- sffdr(p_x, fpi0 = fpi0$fpi0)set.seed(1) n <- 20000 e1 <- rnorm(n) e2 <- 0.5 * e1 + sqrt(1 - 0.5^2) * rnorm(n) z_x <- e1 + ifelse(runif(n) < 0.05, 3, 0) # primary z_y <- e2 + ifelse(runif(n) < 0.1, 3, 0) # informative, same samples dec <- decorrelate_informative(z_x, z_y, rho_o = 0.5) p_x <- 2 * pnorm(-abs(z_x)) mpi0 <- pi0_model(as.matrix(dec$p2)) fpi0 <- fpi0est(p_x, pi0_model_obj = mpi0) res <- sffdr(p_x, fpi0 = fpi0$fpi0)
Estimates the null correlation of z-scores between studies that share
samples, for every pair of cbind(z, z_informative). For each pair it
uses the mean product of z-scores inside the box ,
which is unaffected by signal present in only one of the two studies, and
inverts it through the exact truncated bivariate normal moment. Cross-trait
LD score regression intercepts are preferable for GWAS when available.
estimate_overlap(z, z_informative, indep_snps = NULL, threshold = 1.5)estimate_overlap(z, z_informative, indep_snps = NULL, threshold = 1.5)
z |
Numeric vector; the primary study's z-scores. |
z_informative |
Numeric vector or matrix of the informative studies' z-scores. |
indep_snps |
Optional logical or integer vector of approximately independent SNPs to use. |
threshold |
Numeric; the box half-width |
The correlation matrix, primary first.
Estimates the functional proportion of null tests (pi0) using a GLM approach with a constrained binomial family. This function fits models across multiple lambda thresholds and selects the optimal estimate via MISE minimization.
fpi0est( p, pi0_model_obj = NULL, z = NULL, pi0_model = NULL, indep_snps = NULL, weights = NULL, lambda = seq(0.05, 0.95, 0.05), constrained.p = TRUE, tol = 1e-09, maxit = 200, ncores = 1L, verbose = TRUE, ... )fpi0est( p, pi0_model_obj = NULL, z = NULL, pi0_model = NULL, indep_snps = NULL, weights = NULL, lambda = seq(0.05, 0.95, 0.05), constrained.p = TRUE, tol = 1e-09, maxit = 200, ncores = 1L, verbose = TRUE, ... )
p |
Numeric vector of p-values. |
pi0_model_obj |
Optional list object returned by |
z |
Optional data frame or matrix of rank-transformed covariates.
Required if |
pi0_model |
Optional formula (as character string or formula object).
Required if |
indep_snps |
Optional independent SNPs for model fitting: a logical
vector with one entry per p-value, row indices, or SNP ids matched against
|
weights |
Optional numeric vector of weights in (0, 1], one per p-value, for the model fit. For GWAS data, this should be inverse LD scores (values above 1 are clamped to 1 with a warning). If these are available then you don't have to use the indep_snps argument, as the weights will effectively prioritize independent SNPs in the fitting process. |
lambda |
Numeric vector of lambda thresholds in (0, 1). Default is |
constrained.p |
Logical; use constrained binomial family. Default is TRUE. |
tol |
Numeric; convergence tolerance. Default is 1e-9. |
maxit |
Integer; maximum iterations. Default is 200. |
ncores |
Integer; number of cores for parallel lambda fitting via
|
verbose |
Logical; print progress messages. Default is TRUE. |
... |
Additional arguments passed to |
Algorithm:
For each lambda threshold, fit a binomial GLM:
Use constrained binomial family to ensure predictions in (0, 1)
Select optimal lambda via MISE minimization
Usage Patterns:
Pattern 1 (Recommended): Use output from pi0_model
mpi0 <- pi0_model(z) # Clean syntax: fpi0_out <- fpi0est(p, mpi0)
Pattern 2: Manually specify formula and covariates
fpi0_out <- fpi0est(p, z = z_matrix, pi0_model = formula_obj)
An object of class fpi0 (a list) containing:
Numeric vector of functional pi0 estimates for each test, floored at 0.01. Tests with a missing p-value or covariate get 1.
A data frame with one row per lambda value: the MISE, the mean fitted pi0 and whether it was chosen.
The Mean Integrated Squared Error (MISE) at each lambda.
The selected optimal lambda value.
The fitted model at the selected lambda (class
fastglm2; predict() gives ).
# Import data data(bmi) # Separate main p-values and conditioning p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1]) # Apply pi0_model to create model (uses adaptive knot selection) fmod <- pi0_model(z) # Estimate functional pi0 fpi0_out <- fpi0est(p, fmod) fpi0 <- fpi0_out$fpi0 # Apply sffdr sffdr_out <- sffdr(p, fpi0)# Import data data(bmi) # Separate main p-values and conditioning p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1]) # Apply pi0_model to create model (uses adaptive knot selection) fmod <- pi0_model(z) # Estimate functional pi0 fpi0_out <- fpi0est(p, fmod) fpi0 <- fpi0_out$fpi0 # Apply sffdr sffdr_out <- sffdr(p, fpi0)
Calculate functional p-values from functional local FDRs.
fpvalues(lfdr, p = NULL)fpvalues(lfdr, p = NULL)
lfdr |
A vector of functional local FDRs in [0, 1] (values above 1 are clamped to 1). |
p |
A vector of p-values for ranking purposes. Default is NULL. |
Tests are ranked by lfdr, with ties broken by p. The
functional q-value of a test is the mean lfdr of the tests ranked at
or above it, and its functional p-value is the lfdr mass ranked at
or above it as a fraction of the total mass. Tests with identical
lfdr and p are one set: each receives the set's values, so
the result does not depend on the order of the input. When p is
NULL, ties are broken at random.
A list containing:
fp |
Functional p-values. |
fq |
Functional q-values. |
Performs kernel density estimation on p-values (univariate) or joint p-value/covariate pairs (bivariate) using local regression on the probit scale. The estimator is optimized for GWAS data with linkage disequilibrium (LD) and uses adaptive downsampling to prioritize signal-rich regions while maintaining computational efficiency.
kernelEstimator( x, eval.points = x, epsilon = .Machine$double.xmin, epsilon.max = 1 - 1e-04, maxk = 5000, maxit = 200, target_null = 1e+05, trim = 0, nn = NULL, tail_threshold = -2, weights = NULL, verbose = TRUE, ... )kernelEstimator( x, eval.points = x, epsilon = .Machine$double.xmin, epsilon.max = 1 - 1e-04, maxk = 5000, maxit = 200, target_null = 1e+05, trim = 0, nn = NULL, tail_threshold = -2, weights = NULL, verbose = TRUE, ... )
x |
Numeric vector of p-values (for 1D density) or a 2-column matrix where the first column contains an informative covariate/surrogate and the second column contains the p-values. All p-values must be in (0, 1) and the covariate/surrogate must be rank-transformed to be (0,1]. |
eval.points |
Points at which to evaluate the density estimate. Defaults to
|
epsilon |
Lower bound for p-values to prevent numerical issues. P-values
below this are clamped to |
epsilon.max |
Upper bound for p-values. P-values above this are clamped to
|
maxk |
Space allocated for |
maxit |
Maximum number of iterations for the bivariate local likelihood
fit. Default: |
target_null |
Maximum number of null SNPs to include in the weighted fit
(bivariate case only). SNPs in the signal-enriched tail (defined by
|
trim |
Numeric in [0, 1); if > 0, flattens the density for p-values in (1 - trim, 1) to reduce boundary artifacts from p-values artificially clumped near 1. Only applies to the p-value dimension. Default is 0 (no trimming). |
nn |
Nearest-neighbor bandwidth parameter for |
tail_threshold |
Threshold on the probit-transformed p-values
(bivariate case only). SNPs with |
weights |
Optional numeric vector of prior weights in (0, 1] for the fit (values above 1 are clamped to 1). |
verbose |
Logical; whether to print progress messages during fitting. Default is TRUE. |
... |
Additional arguments passed to |
The function implements a multi-stage density estimation procedure:
Probit transformation: P-values are transformed to the normal
scale via qnorm to stabilize variance and handle extreme values.
Adaptive downsampling (bivariate only): To handle large GWAS
datasets efficiently, the null region (p-values above the tail threshold)
is downsampled to target_null SNPs, with inverse-probability weighting
to preserve the density. Signal-enriched SNPs (tail) are always retained.
Cascade fitting: Multiple fitting strategies are attempted in sequence, with decreasing resolution and increasing robustness, until a valid fit is obtained.
Jacobian correction: Density estimates are transformed back to the original p-value scale using the Jacobian of the probit transformation.
The nearest-neighbor bandwidth nn controls smoothing and LD robustness.
By default the window spans about 5000 of the rows the density is fitted on
(the tail plus up to target_null null SNPs), so the estimator averages
over multiple LD blocks (~30-50 blocks in European ancestry populations),
reducing spurious local structure while preserving true signal-covariate
relationships. Versions before 1.1.4 took the fraction relative to all SNPs,
which left about 600 neighbours at one million SNPs.
In the univariate case, tied values are collapsed into one point with the corresponding multiplicity before fitting.
A data frame with columns:
Evaluation points (original scale).
Estimated density at evaluation points (original scale).
Evaluation points on probit scale (qnorm(x)).
Estimated density on probit scale.
The returned object has an attribute "lfit" containing the fitted
locfit object for diagnostics.
## Not run: # 1D density estimation p <- runif(10000, 0, 1) dens <- kernelEstimator(p) plot(dens$x, dens$fx, type = "l") # 2D density: the covariate (rank-transformed) first, the p-values second p <- runif(10000, 0, 1) z <- runif(10000) x_mat <- cbind(z, p) dens <- kernelEstimator(x_mat) ## End(Not run)## Not run: # 1D density estimation p <- runif(10000, 0, 1) dens <- kernelEstimator(p) plot(dens$x, dens$fx, type = "l") # 2D density: the covariate (rank-transformed) first, the p-values second p <- runif(10000, 0, 1) z <- runif(10000) x_mat <- cbind(z, p) dens <- kernelEstimator(x_mat) ## End(Not run)
Enforces strictly non-increasing density conditional on surrogate bins. Note: Data must be sorted by (group, pvalue) before passing.
monoSmooth_conditional(pvalue, density, group)monoSmooth_conditional(pvalue, density, group)
pvalue |
Numeric vector of p-values. |
density |
Numeric vector of estimated densities. |
group |
Integer vector of surrogate bins. |
A numeric vector of smoothed densities.
Fast Conditional PAVA (Decreasing)
monoSmooth_pava(pvalue, density, group)monoSmooth_pava(pvalue, density, group)
pvalue |
A numeric vector of p-values. Must be sorted ascendingly within each group. |
density |
A numeric vector of initial density estimates (e.g., from locfit). |
group |
An integer vector denoting the surrogate bin/stratum for each observation. |
A numeric vector of smoothed, monotonically non-increasing densities.
Generates a natural spline model for the functional proportion of null tests (fpi0) by adaptively selecting knots based on an FDR threshold.
pi0_model( z, indep_snps = NULL, weights = NULL, fdr_threshold = 0.25, min_discoveries = 150, min_snps_per_knot = 100, n_knots = 5, verbose = TRUE, seed = 2026 )pi0_model( z, indep_snps = NULL, weights = NULL, fdr_threshold = 0.25, min_discoveries = 150, min_snps_per_knot = 100, n_knots = 5, verbose = TRUE, seed = 2026 )
z |
Numeric matrix/data.frame of p-values (rows=tests, cols=traits). |
indep_snps |
Independent SNPs (the training subset): a logical vector
with one entry per row of |
weights |
Optional numeric vector of weights in (0, 1] (e.g., inverse LD scores) for weighted spline fitting. Default: NULL. |
fdr_threshold |
FDR threshold for signal definition. Default: 0.25. |
min_discoveries |
Min (LD-independent) significant hits required to place a trait's knots from its discoveries; with fewer, the trait keeps a single fallback knot. Default: 150. |
min_snps_per_knot |
Base minimum significant SNPs per knot interval. Default: 100. Dynamically scales up for larger datasets to prevent overfitting. |
n_knots |
Target knot count. Default: 5. Automatically reduced if insufficient
discoveries (via |
verbose |
Print selection details. Default: TRUE. |
seed |
Optional seed for reproducibility. Default: 2026. |
Independent SNPs determine knot placement; all SNPs train the model to capture signal shape.
For smaller datasets (<100K), consider reducing min_snps_per_knot and min_discoveries to improve knot placement.
Every column must hold p-values: the knots of a column are placed at the
quantiles of its discoveries (qvalue at
fdr_threshold), and a column with values outside [0, 1] gets the
fallback knot with a warning. A null anchor knot at rank 0.1 is always
added.
List containing:
fmod |
Model formula (using |
zt |
Data frame of globally rank-transformed p-values |
Graphical display of the sffdr object
## S3 method for class 'sffdr' plot(x, rng = c(0, 5e-08), ...)## S3 method for class 'sffdr' plot(x, rng = c(0, 5e-08), ...)
x |
A sffdr object. |
rng |
Significance region to show (functional p-value range). Optional. If no functional p-value falls in it, the range is widened to the smallest one percent. |
... |
Additional arguments. Currently unused. |
Nothing of interest.
Andrew J. Bass
## Not run: # import data data(bmi) # Separate main p-values and conditioning p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1]) # Apply pi0_model to create model fmod <- pi0_model(z) # Estimate functional pi0 fpi0_out <- fpi0est(p, z = fmod$zt, pi0_model = fmod$fmod) fpi0 <- fpi0_out$fpi0 # Apply sffdr sffdr_out <- sffdr(p, fpi0) # Plot significance results plot(sffdr_out, rng = c(0, 5e-4)) ## End(Not run)## Not run: # import data data(bmi) # Separate main p-values and conditioning p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1]) # Apply pi0_model to create model fmod <- pi0_model(z) # Estimate functional pi0 fpi0_out <- fpi0est(p, z = fmod$zt, pi0_model = fmod$fmod) fpi0 <- fpi0_out$fpi0 # Apply sffdr sffdr_out <- sffdr(p, fpi0) # Plot significance results plot(sffdr_out, rng = c(0, 5e-4)) ## End(Not run)
Estimate functional p-values, q-values, and local false discovery rates (lfdr) for GWAS data leveraging summary statistics from related traits. Functional p-values map from the functional q-value (FDR-based measure) to a p-value for type I error rate control, accounting for pleiotropy that impacts the prior probability of association.
sffdr( p.value, fpi0, surrogate = NULL, weights = NULL, epsilon = .Machine$double.xmin, nn = NULL, monotone = TRUE, monotone_method = c("min", "pava"), fp_ties = TRUE, seed = 2026, verbose = TRUE, ... )sffdr( p.value, fpi0, surrogate = NULL, weights = NULL, epsilon = .Machine$double.xmin, nn = NULL, monotone = TRUE, monotone_method = c("min", "pava"), fp_ties = TRUE, seed = 2026, verbose = TRUE, ... )
p.value |
Numeric vector of p-values to analyze, in [0, 1]. Exact zeros are clamped to the smallest positive double with a warning; negative values are an error. |
fpi0 |
Numeric vector of functional pi0 estimates, obtained
from |
surrogate |
Optional numeric vector (same length as |
weights |
Ignored, with a warning (kept for compatibility with sffdr 1.1.2); see Details. |
epsilon |
Numeric; lower bound for p-value clamping during density estimation.
Default is |
nn |
Numeric; nearest-neighbor bandwidth for |
monotone |
Logical; if TRUE, enforces monotonicity of the estimated density with respect to p-values within bins of the surrogate variable. Default is TRUE. |
monotone_method |
Character; method for enforcing monotonicity. Options are "min" (running minimum) or "pava" (Pool Adjacent Violators Algorithm). Default is "min". We do not recommend "pava" for GWAS data sets due to LD, but it may be better for expression data. |
fp_ties |
Logical; whether to break ties in functional p-values using the original p-value ordering. Default is TRUE. |
seed |
Integer; random seed for reproducibility of rank tie-breaking, the jitter and the null downsampling in the density fit. Default is 2026. |
verbose |
Logical; print progress messages. Default is TRUE. |
... |
Additional arguments passed to |
The surrogate functional FDR (sfFDR) methodology extends the functional FDR framework to leverage multiple informative variables (e.g., functional annotations, GWAS summary statistics) for increased power while controlling false discovery rates.
Workflow:
Estimate functional pi0 (proportion of nulls) using fpi0est
Call sffdr() with p-values and estimated functional pi0
Use returned functional p-values/q-values/local FDRs for significance testing
Surrogate Variable: If not specified, the estimated functional pi0 is used as the surrogate variable.
Weights and LD:
The joint density is estimated over all tested SNPs, which is the set the
functional p-values count errors in, so sffdr() no longer reweights
it. Inverse-LD weights change the density to that of independent signals;
because associated SNPs have larger LD scores, this biases the local FDR
(conservative with inverse-LD weights, anti-conservative with weights that
are larger at associated SNPs). LD scores carry information about
association, so their place is the functional pi0 model, as a covariate
(see the vignette); as inverse weights they remain available in
pi0_model and fpi0est in place of an LD-pruned
training set.
Perfect-LD copies:
SNPs with identical p-values are copies of one signal (perfect LD). Each
set of copies contributes one point to the density fit, and the density is
then evaluated for every SNP, so copies do not narrow the effective window
in the tail (which inflated functional p-values below 1e-3 before 1.1.4).
Rounded p-values tie between unrelated SNPs too, almost only in the bulk
of the distribution where values are dense per significant digit, and
collapsing those would distort the density; so identical p-values count as
copies only when they are at least half as common below p = 0.023 as above
it, which holds for LD copies and fails for rounding. The window itself
spans about 5,000 of the rows the density is fitted on; see
kernelEstimator.
Interpretation note: Functional p-values reflect the global ranking of all SNPs by the functional local FDR. A SNP with a small functional p-value but large flfdr (> 0.5) has weak individual evidence but ranks favorably relative to null SNPs. This happens when there is strong prior evidence (fpi0) of the SNP being non-null. Thus, users should also consider flfdr (and fqvalues) alongside fpvalues when assessing individual SNP evidence. SNPs with flfdr > 0.5 should be interpreted cautiously regardless of their functional p-value.
An S3 object of class "sffdr" containing:
call |
The function call. |
pvalues |
The p-values as analysed (exact zeros clamped). |
fpvalues |
Functional p-values. |
fqvalues |
Functional q-values. |
flfdr |
Functional local false discovery rates. |
fpi0 |
Functional pi0 estimates. |
fx |
Joint density estimates at observed (p-value, surrogate) pairs. |
Andrew J. Bass
fpi0est, plot.sffdr, kernelEstimator
## Not run: # Import data data(bmi) # Separate main p-values and conditioning p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1]) # Apply pi0_model to create model # (note: use indep_snps argument to specify independent SNPs for training) fmod <- pi0_model(z) # Estimate functional pi0 # (note: use indep_snps argument to specify independent SNPs for training) fpi0_out <- fpi0est(p, fmod) fpi0 <- fpi0_out$fpi0 # Apply sffdr sffdr_out <- sffdr(p, fpi0) # Plot significance results plot(sffdr_out) # Extract functional quantities fp <- sffdr_out$fpvalues fq <- sffdr_out$fqvalues flfdr <- sffdr_out$flfdr ## End(Not run)## Not run: # Import data data(bmi) # Separate main p-values and conditioning p-values p <- sumstats$bmi z <- as.matrix(sumstats[, -1]) # Apply pi0_model to create model # (note: use indep_snps argument to specify independent SNPs for training) fmod <- pi0_model(z) # Estimate functional pi0 # (note: use indep_snps argument to specify independent SNPs for training) fpi0_out <- fpi0est(p, fmod) fpi0 <- fpi0_out$fpi0 # Apply sffdr sffdr_out <- sffdr(p, fpi0) # Plot significance results plot(sffdr_out) # Extract functional quantities fp <- sffdr_out$fpvalues fq <- sffdr_out$fqvalues flfdr <- sffdr_out$flfdr ## End(Not run)