---
title: "sffdr Package Vignette"
author: "Andrew J. Bass and Chris Wallace"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{sffdr: Surrogate Functional False Discovery Rates}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
```

## 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*](https://doi.org/10.1038/s43588-025-00852-3); 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:

```{r, eval = FALSE}
# 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.

```{r}
library(sffdr)
data(bmi)

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

head(sumstats)
```

### 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)$​.

```{r}
# Create model
mpi0 <- pi0_model(z)

# Estimation of functional pi0
fpi0 <- fpi0est(p, mpi0)
```

   **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.
```{r}
# Apply sfFDR
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)   

# 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 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.

```{r}
# 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)

fpi0 <- fpi0est(p, mpi0, indep_snps = indep_snps)

# Estimate quantities for all SNPs
sffdr_out <- sffdr(p, fpi0 = fpi0$fpi0)
```

   **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

### 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.

```{r, eval = FALSE}
# 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.