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.8674931. 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 : 2500 (All SNPs)
#> --------------------------------------------------
#> [Auto-Scale]: N = 10000. Bandwidth set to 100 independent SNPs per knot.
#> [bfp]: Fallback -> 1 Knot at 0.1000
#> [cholesterol]: Fallback -> 1 Knot at 0.1000
#> [triglycerides]: Fallback -> 1 Knot at 0.1000
#> --------------------------------------------------
#> ~splines::ns(bfp, knots = c(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.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.0)
#> ==================================================
#> Transforming surrogate variable...
#> Fitting kernel density...
#> High-Res Tcub (nn=0.1)... ✓
#> 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))
#> Scale for x is already present.
#> Adding another scale for x, which will replace the existing scale.
3. Incorporating linkage disequilibrium (LD)
In standard GWAS, SNPs are often correlated due to LD. To account for
this, you should provide a 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 : 2500 (All SNPs)
#> --------------------------------------------------
#> [Auto-Scale]: N = 10000. Bandwidth set to 100 independent SNPs per knot.
#> [bfp]: Fallback -> 1 Knot at 0.1000
#> [cholesterol]: Fallback -> 1 Knot at 0.1000
#> [triglycerides]: Fallback -> 1 Knot at 0.1000
#> --------------------------------------------------
#> ~splines::ns(bfp, knots = c(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.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.0)
#> ==================================================
#> Transforming surrogate variable...
#> Fitting kernel density...
#> High-Res Tcub (nn=0.1)... ✓
#> Computing functional local FDRs...
#> Computing functional p-values and q-values...
#> Done.
#> ==================================================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.