Genome-wide pipeline that computes per-pair enrichment
estimates by passing an outcome PIP vector and a set of annotation
credible-set posteriors to qtlEnrichment. The returned table
feeds colocPipeline via its enrichment argument.
Either side may be a QtlFineMappingResult or a
GwasFineMappingResult, so QTL-in-GWAS, QTL-in-QTL and
GWAS-in-GWAS enrichment all run through one code path. The argument names
keep the QTL / GWAS wording of the common case; what they mean generally
is that gwasFineMappingResult is the outcome whose PIPs
are scanned and qtlFineMappingResult is the annotation
whose region fits are tested for enrichment within them.
Not gene-parallelisable: the enrichment estimator runs over the full genome of outcome PIPs and the full collection of annotation fits at once.
Usage
qtlEnrichmentPipeline(
gwasFineMappingResult,
qtlFineMappingResult,
qtlEnrichmentParam = QtlEnrichmentParam(),
numThreads = 1L,
seed = NULL,
verbose = TRUE
)Arguments
- gwasFineMappingResult
The outcome side; see above.
- qtlFineMappingResult
The annotation side; see above.
- qtlEnrichmentParam
The estimator's algorithmic settings, built with
QtlEnrichmentParam:numGwas,piQtl,lambda,impN,doubleShrinkageandbesselCorrection. A bare list is refused, since it cannot be checked.- numThreads
Number of threads used by
qtlEnrichment. Default1.- seed
Integer or
NULL. Base random seed forwarded toqtlEnrichmentfor reproducible multiple imputation.NULL(default) draws a nondeterministic seed.- verbose
Logical. Print progress messages. Default
TRUE.
Value
A tibble with one row per (outcome trait, annotation unit) pair.
The identity columns are gwasStudy, gwasContext,
gwasTrait, qtlStudy, qtlContext; the axes a side does
not have are NA (gwasContext / gwasTrait for a GWAS
outcome, qtlContext for a GWAS annotation). Suitable as the
enrichment argument to colocPipeline, which joins on
those columns.
The estimates are enrichmentLogOdds, the enrichment parameter
\(a_1\) on the log-odds scale, with its standard error
enrichmentSe; enrichment is the same quantity as a
multiplicative factor, \(e^{a_1} - 1\), reported for interpretation
(so \(a_1 = 0\) means no enrichment).
enrichmentLogOddsNoShrinkage and
enrichmentSeNoShrinkage are the same estimate before shrinkage,
intercept / interceptSe are \(a_0\), and colocP1,
colocP2, colocP12 are the enrichment-informed priors the
estimator derives from \((a_0, a_1)\) exactly as fastenloc's
set_enrich_params does; those three are what
colocPipeline scores with when given this table.
effectiveMiRounds is how many multiple-imputation rounds survived
outlier filtering.
Inputs
gwasFineMappingResult: the genome-wide outcome collection. Each row's stored fit (susieFit) must carry apipvector. One PIP vector is built per outcome trait – keyed bystudyfor a GWAS collection (one trait per study) and by (study,context,trait) for a QTL one, since a variant's PIP differs between molecular traits and pooling them would collide.qtlFineMappingResult: the annotation collection. Each row's stored fit must carryalpha,pip, and prior-variance fields (V). Its region fits are pooled per (study,context) – perstudyalone for a GWAS collection, which has no context axis.
LD reference compatibility check
A GWAS outcome collection
must have a non-NULL ldReference (it should be RSS-derived). Where
both sides carry one, the two must be the same LD reference.
The two must be the same kind of reference. Two genotype panels must
share their samples and allele orientation; two LdMatrix objects
must agree on the LD they share; two GenotypeMixture objects must
average the same panels with the same weights. They need NOT carry the
same variants: QC-ing the two sides separately trims each reference to
its own surviving variants, and a partial overlap only warns. A
NULL on either side (an individual-level fit) skips the check.
Examples
data(gwasFineMappingExample)
data(qtlFineMappingExample)
qtlEnrichmentPipeline(
gwasFineMappingResult = gwasFineMappingExample,
qtlFineMappingResult = qtlFineMappingExample
)
#> 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.00141
#> 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.00071
#> 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
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 10
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 11
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 12
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 13
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 14
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 15
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 16
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 17
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 18
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 19
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 20
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 21
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 22
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 23
#> Proportion of xQTL missing from GWAS variants: 0 in MI round 24
#> EM updates completed!
#> # A tibble: 1 × 16
#> gwasStudy gwasContext gwasTrait qtlStudy qtlContext enrichment enrichmentSe
#> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl>
#> 1 study_1 NA NA study_1 context_1 -0.00000705 1.000
#> # ℹ 9 more variables: enrichmentLogOdds <dbl>,
#> # enrichmentLogOddsNoShrinkage <dbl>, enrichmentSeNoShrinkage <dbl>,
#> # intercept <dbl>, interceptSe <dbl>, colocP1 <dbl>, colocP2 <dbl>,
#> # colocP12 <dbl>, effectiveMiRounds <dbl>