Implementation of enrichment analysis described in https://doi.org/10.1371/journal.pgen.1006646

Largely follows from fastenloc https://github.com/xqwen/fastenloc but uses `susieR` fitted objects as input to estimate prior for use with `coloc` package (coloc v5, aka SuSiE-coloc). The main differences are 1) now enrichment is based on all QTL variants whether or not they are inside signal clusters; 2) Causal QTL are sampled from SuSiE single effects, not signal clusters; 3) Allow a variant to be QTL for not only multiple conditions (eg cell types) but also multiple regions (eg genes). Other minor improvements include 1) Make GSL RNG thread-safe; 2) Release memory from QTL binary annotation samples immediately after they are used.

Usage

qtlEnrichment(
  gwasPip,
  susieQtlRegions,
  numGwas = NULL,
  piQtl = NULL,
  lambda = 1,
  impN = 25,
  doubleShrinkage = FALSE,
  besselCorrection = TRUE,
  numThreads = 1,
  verbose = TRUE,
  alignNames = TRUE,
  seed = NULL
)

Arguments

gwasPip

This is a vector of GWAS PIP, genome-wide.

susieQtlRegions

This is a list of SuSiE fitted objects per QTL unit analyzed

numGwas

This parameter is highly important if GWAS input does not contain all SNPs interrogated (e.g., in some cases, only fine-mapped geomic regions are included). Then users must pick a value of total_variants and estimate piGwas beforehand by: sum(gwasPip$pip)/numGwas. If numGwas is null, piGwas would be sum(gwasPip$pip)/total_variants.

piQtl

This parameter can be safely left to default if your input QTL data has enough regions to estimate it.

lambda

Similar to the shrinkage parameter used in ridge regression. It takes any non-negative value and shrinks the enrichment estimate towards 0. When it is set to 0, no shrinkage will be applied. A large value indicates strong shrinkage. The default value is set to 1.0.

impN

Rounds of multiple imputation to draw QTL from, default is 25.

doubleShrinkage

Logical. Apply the double-shrinkage correction to the enrichment estimate. Default FALSE.

besselCorrection

Logical. Apply Bessel's correction when estimating the sampling variance. Default TRUE.

numThreads

Number of Simultaneous running CPU threads for multiple imputation, default is 1.

verbose

Logical. Print progress messages. Default TRUE.

alignNames

Logical; when TRUE (default) QTL pip names are aligned to the GWAS variant-naming convention via matchVariants. Set FALSE when the caller has already aligned them (e.g. qtlEnrichmentPipeline aligns each QTL tuple once against the union GWAS panel rather than re-aligning per GWAS study); only the cheap per-study unmatched set is then recomputed, skipping the costly harmonizeAlleles pass.

seed

Integer or NULL. Base random seed for the multiple- imputation sampler; each imputation round derives its own seed from it, so a fixed seed gives reproducible results. NULL (default) draws a nondeterministic seed.

Value

A list of enrichment parameter estimates

Details

Uses output of susie from the susieR package.

See also

susie

Examples

# Simulate fake data for gwasPip
nGwasPip <- 1000
gwasPip <- runif(nGwasPip)
names(gwasPip) <- paste0("snp", 1:nGwasPip)
# Simulate fake data for a single SuSiEFit object
simulateSusiefit <- function(n, p) {
  pip <- runif(n)
  names(pip) <- paste0("snp", 1:n)
  alpha <- t(matrix(runif(n * p), nrow = n))
  alpha <- t(apply(alpha, 1, function(row) row / sum(row)))
  list(
    pip = pip,
    alpha = alpha,
    prior_variance = runif(p)
  )
}
# Simulate multiple SuSiEFit objects
nSusieFits <- 2
susieFits <- replicate(
  nSusieFits, simulateSusiefit(nGwasPip, 10), simplify = FALSE)
# Add these fits to a list, providing names to each element
names(susieFits) <- paste0("fit", seq_along(susieFits))
# Set other parameters
impN <- 10
lambda <- 1
numThreads <- 1
library(pecotmr)
en <- qtlEnrichment(
  gwasPip, susieFits, lambda = lambda, impN = impN,
  numThreads = numThreads)
#> Warning: numGwas is not provided. Estimating piGwas from the data. Note that this estimate may be biased if the input gwasPip does not contain genome-wide variants.
#> Estimated piGwas: 0.50153
#> Warning: piQtl is not provided. Estimating piQtl from the data. Note that this estimate may be biased if either 1) the input susieQtlRegions does not have enough data, or 2) the single effects only include variables inside of credible sets or signal clusters.
#> Estimated piQtl: 0.5032
#> Fine-mapped GWAS and QTL data loaded successfully for enrichment analysis!
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 0
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 1
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 2
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 3
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 4
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 5
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 6
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 7
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 8
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 9
#> EM updates completed!