Quality Control for GWAS and QTL Summary Statistics
pecotmr authors
2026-11-02
Source:vignettes/rss-qc.Rmd
rss-qc.RmdOverview
Regression-with-summary-statistics (RSS) methods — both fine-mapping and regularized regression — all require matching variants in a summary-statistics table to an LD reference panel. Two are frequently encountered in practice:
- Allele alignment. The sumary statisticss and LD panel may flip the reference and alternate allele, encode strand-ambiguous variants inconsistently, or contain indels and duplicates that need filtering.
- LD-summary-statistic mismatch. Even when alleles align, individual z-scores can be inconsistent with the LD pattern around them because of imputation errors, ancestry differences between the GWAS and LD panel, or per-variant sample-size differences.
The summaryStatsQc() function provides a convenient
interface to performing many QC tasks to address these issues. It takes
a GwasSumStats or QtlSumStats object (see the
Building QtlSumStats and GwasSumStats objects vignette), and
returns a new collection of the same class with updated entries, the
options it ran with recorded in qcOptions(), and what it
did to each row in qcInfo(). Downstream pipelines require
qcOptions() to be filled, so QC is mandatory before
fine-mapping, TWAS-weight learning, colocalization, etc.
Running QC
summaryStatsQc() operates on a GwasSumStats
or QtlSumStats object See the Building QtlSumStats and
GwasSumStats objects vignette for details on how to make these
objects. In this example the object is GwasSumStats object
called gws, but the interface is the same for
QtlSumStats.
summaryStatsQc() provides many options for the user to
tweak:
gws_qcd <- summaryStatsQc(
gws,
sumstatsFilterParam = SumstatsFilterParam(
removeIndels = FALSE,
removeStrandAmbiguous = TRUE,
infoCutoff = 0, # 0 = off; > 0 requires INFO column
nCutoff = 5
),
panelFilterParam = PanelFilterParam(
mafCutoff = 0 # 0 = off; > 0 requires MAF column
),
skipRegion = NULL, # e.g. MHC "chr6:25000000-35000000"
ldMismatchQcMethod = "slalom", # "none", "slalom", "dentist"
alleleFlipKriging = FALSE, # opt-in pre-SLALOM LD-consistency check
impute = FALSE, # RAISS imputation against ldReference
matchMinProp = 0
)gws_qc is a new GwasSumStats of the same
class, with cleaned entries. What QC did to each row is a table,
qcInfo(), and the options it ran with are
qcOptions():
qcInfo(gws_qcd) # one row per tuple: variantsIn, matchedAgainstLdReference,
# krigingOutliersDropped, ldMismatchOutliersDropped, ..., variantsOut
qcOptions(gws_qcd)$options # the knobs you passed, kept for provenanceA non-zero mafCutoff requires every entry to carry a
MAF mcol; non-zero infoCutoff requires
INFO; non-zero nCutoff requires
N. A missing column with a non-zero cutoff will produce an
error - if you don’t want a filter to run, leave its cutoff at the
default value of 0.
The pieces QC is built from
summaryStatsQc() orchestrates the steps below. Each is
exported, because the same operations come up outside a full QC run —
preparing a variant list, auditing why something was dropped, or
applying one check without the rest.
Classifying variants
Allele strings decide how a variant is handled: indels cannot be matched to a panel the same way SNPs are, and strand-ambiguous SNPs need care.
ids <- c(
"chr22:10685239:C:T",
"chr22:10761227:G:A",
"chr22:11285239:CT:C",
"chr22:11500000:A:T"
)
classifyVariantType(ids)## [1] "SNP" "SNP" "deletion" "SNP"Matching a variant list to an LD panel
filterVariantsByLdReference() keeps only variants the
panel actually carries, and reports how many it dropped. It returns both
the surviving IDs and their positions in the input, so a parallel table
can be subset the same way.
ldMeta <- system.file(
"extdata",
"ld_reference",
"ld_meta_file.tsv",
package = "pecotmr"
)
kept <- filterVariantsByLdReference(ids, ldMeta)
kept$data## [1] "chr22:10685239:C:T" "chr22:10761227:G:A" "chr22:11285239:CT:C"
kept$idx## [1] 1 2 3keepIndel = FALSE additionally drops indels, which is
what an RSS method wanting a SNP-only design needs:
filterVariantsByLdReference(ids, ldMeta, keepIndel = FALSE)$data## [1] "chr22:10685239:C:T" "chr22:10761227:G:A"Reconciling two variant tables
mergeVariantInfo() matches two variant sets while
allowing for flipped alleles — the second table is re-oriented to the
first rather than treated as a different variant:
panel <- data.frame(
chrom = "chr22", pos = c(10685239, 10761227),
ref = c("C", "G"), alt = c("T", "A")
)
sumstats <- data.frame(
chrom = "chr22", pos = c(10685239, 10761227),
ref = c("T", "G"), alt = c("C", "A") # first variant is flipped
)
mergeVariantInfo(panel, sumstats)## # A tibble: 2 × 4
## chrom pos alt ref
## <chr> <dbl> <chr> <chr>
## 1 chr22 10685239 T C
## 2 chr22 10761227 A GThe flipped variant comes back on the panel’s orientation.
all = FALSE returns only the re-oriented second set instead
of the union.
Applying one mismatch check
ldMismatchQc() is the LD-mismatch step on its own,
taking z-scores plus either a correlation matrix (R) or
genotypes (X):
ldMismatchQc(zScore = z, R = ldMatrix, nSample = 10000, method = "slalom")and autoDecision() is the helper that reads credible-set
structure to decide whether a flagged variant should be imputed or
dropped — the choice summaryStatsQc() otherwise makes for
you.
Allele-flip QC: the SuSiE kriging prefilter
alleleFlipKriging = TRUE enables an optional prefilter
that runs before the SLALOM / DENTIST z-mismatch step.
Where SLALOM and DENTIST look for z-scores that are broadly inconsistent
with the LD pattern, the kriging prefilter specifically targets
allele-flip / LD-mismatch outliers — variants whose
observed z-score disagrees with the value the LD structure
predicts for them from their neighbours.
It uses susieR’s kriging RSS diagnostic:
-
susieR::estimate_s_rss()estimates the overall LD-mismatch scalesbetween the harmonized z-scores and the reference LD matrixR. -
susieR::kriging_rss(z, R, n, s)computes, for each variant, the leave-one-out conditional distribution ofz_igiven all the other z-scores under the LD model — i.e. the z-score the LD panel predicts for that variant. - The standardized residual between observed and predicted
(
z_std_diff) is ~N(0, 1)when the z-scores and LD panel agree. A large residual flags an allele-flip / LD-mismatch outlier; variants below a two-sided p-value threshold (default5e-8) are dropped. The number removed is recorded in theqcInfo()table askrigingOutliersDropped.
Because it needs a per-region LD matrix and the kriging diagnostic,
this is an RSS-only check: it draws R from
each row’s LD reference and requires a susieR that exports
kriging_rss() / estimate_s_rss() (otherwise it
errors — which is why it is off by default).
The same diagnostic is exposed standalone as
krigingOutlierQc() when you want to inspect the per-variant
predicted z, residual, statistic, and p-value directly:
kr <- krigingOutlierQc(zScore = z, R = ld, n = 10000)
kr$outlier # logical vector: TRUE = flagged as an outlier
head(kr$diagnostics) # variant_id, z, predicted, residual, statistic,
# p_value, outlierKriging and SLALOM/DENTIST are complementary and can be combined: the
kriging prefilter removes allele-flip outliers up front, then
ldMismatchQcMethod handles the remaining LD-mismatch
variants.
LD mismatch correction with SLALOM and DENTIST
summaryStatsQc(ldMismatchQcMethod = ...) dispatches to
one of two LD-mismatch detectors. They differ in how they decide whether
a variant’s z-score is inconsistent with the LD pattern around it.
SLALOM (Kanai et al. 2022, the recommended default once you turn ldMismatchQcMethod on):
- For each variant, compares its z-score against the lead
variant’s z-score scaled by their LD: a large
(z_i − r_{i,lead} · z_lead)flags the variant as inconsistent. - Restricted to variants in high LD with the lead
(
r2Threshold). - Stable: results do not depend on window-size parameters.
DENTIST (Chen et al. 2021):
- Iteratively imputes each z-score from all other z-scores via SVD-
truncated LD regression, then flags large
(z_obs − z_imp)^2 / (1 − r²)values as outliers. - Operates on windows of the chromosome; the window size noticeably affects how many outliers it reports.
Both methods are available directly as slalom() /
dentist() / dentistSingleWindow() if you want
to run them outside the summaryStatsQc() pipeline
(typically to tune parameters or to inspect the per-variant intermediate
output).
Dentist combines imputation and LD mismatch correction
DENTIST’s outlier rate is sensitive to its window size. On a typical ∼17k-variant test region we observe:
min_dim |
Approx outlier rate |
|---|---|
| 2,000 | ∼8 % |
| 6,000 | ∼1.5 % |
| 20,000 | ∼0.6 % |
For this reason SLALOM is the recommended default when you turn ldMismatchQcMethod on. Use DENTIST when you specifically want its iterative imputation output, or when you are reproducing a DENTIST-based pipeline.
Combining LD mismatch correction and imputation with SLALOM and RAISS
SLALOM only flags outliers; it does not impute them. To
re-impute the flagged variants from the surviving ones, pair SLALOM with
RAISS (a LD-regression imputer).
summaryStatsQc(ldMismatchQcMethod = "slalom", impute = TRUE)
does this for you. The standalone version is exposed as
raiss(); the full set of RAISS knobs lives under the
raissParam argument to summaryStatsQc(), built
with RaissParam():
gws_qcd <- summaryStatsQc(
gws,
ldMismatchQcMethod = "slalom",
impute = TRUE,
raissParam = RaissParam(
r2Threshold = 0.6,
minimumLd = 5,
lamb = 0.01,
# Bound what RAISS is asked to impute, and widen the window it reads.
mafCutoff = 0.01,
flank = 5e5
)
)QC always hands RAISS a genotype matrix, so it solves through the SVD
path and takes svdTol; raiss()’s
rcond belongs to its LD-matrix path and is deliberately not
offered here.
Next steps
QtlSumStats and GwasSumStats are direct
inputs to fineMappingPipeline(),
twasWeightsPipeline() (QtlSumStats only) and
colocboostPipeline()— see the Fine-mapping with
pecotmr, Learning TWAS weights with pecotmr and
Multi-trait colocalization with ColocBoost vignettes.
## R version 4.5.3 (2026-03-11)
## Platform: x86_64-conda-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/runner/work/pecotmr/pecotmr/.pixi/envs/default/lib/libopenblasp-r0.3.34.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pecotmr_0.9.4
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 dplyr_1.2.1
## [3] farver_2.1.2 Biostrings_2.78.0
## [5] S7_0.2.2 bitops_1.1-0
## [7] fastmap_1.2.0 reshape_0.8.10
## [9] mathjaxr_2.0-0 digest_0.6.39
## [11] lifecycle_1.0.5 magrittr_2.0.5
## [13] compiler_4.5.3 rlang_1.3.0
## [15] sass_0.4.10 tools_4.5.3
## [17] utf8_1.2.6 yaml_2.3.12
## [19] knitr_1.52 S4Arrays_1.10.1
## [21] htmlwidgets_1.6.4 bit_4.6.0
## [23] DelayedArray_0.36.1 plyr_1.8.9
## [25] RColorBrewer_1.1-3 abind_1.4-8
## [27] BiocParallel_1.44.0 withr_3.0.3
## [29] purrr_1.2.2 numDeriv_2016.8-1.1
## [31] BiocGenerics_0.56.0 desc_1.4.3
## [33] grid_4.5.3 stats4_4.5.3
## [35] susieR_0.16.6 ggplot2_4.0.3
## [37] scales_1.4.0 MASS_7.3-66
## [39] MultiAssayExperiment_1.36.2 SummarizedExperiment_1.40.0
## [41] cli_3.6.6 rmarkdown_2.32
## [43] metafor_5.0-1 crayon_1.5.3
## [45] ragg_1.5.2 generics_0.1.4
## [47] otel_0.2.0 RcppParallel_5.1.11-2
## [49] httr_1.4.9 tzdb_0.5.0
## [51] BiocBaseUtils_1.12.0 cachem_1.1.0
## [53] stringr_1.6.0 metadat_1.6-0
## [55] parallel_4.5.3 XVector_0.50.0
## [57] matrixStats_1.5.0 vctrs_0.7.3
## [59] Matrix_1.7-6 jsonlite_2.0.0
## [61] IRanges_2.44.0 hms_1.1.4
## [63] S4Vectors_0.48.1 bit64_4.8.6
## [65] mixsqp_0.3-54 irlba_2.4.1
## [67] archive_1.1.13 systemfonts_1.3.2
## [69] tidyr_1.3.2 jquerylib_0.1.4
## [71] glue_1.8.1 pkgdown_2.2.1
## [73] codetools_0.2-20 stringi_1.8.9
## [75] gtable_0.3.6 GenomeInfoDb_1.46.2
## [77] UCSC.utils_1.6.1 GenomicRanges_1.62.1
## [79] quadprog_1.5-8 tibble_3.3.1
## [81] pillar_1.11.1 htmltools_0.5.9
## [83] Seqinfo_1.0.0 R6_2.6.1
## [85] zigg_0.0.2 textshaping_1.0.5
## [87] vroom_1.7.1 evaluate_1.0.5
## [89] lattice_0.23-1 Biobase_2.70.0
## [91] readr_2.2.0 tictoc_1.2.1
## [93] backports_1.5.1 Rsamtools_2.26.0
## [95] Rfast_2.1.5.2 bslib_0.12.0
## [97] Rcpp_1.1.2 SparseArray_1.10.10
## [99] nlme_3.1-170 checkmate_2.3.3
## [101] xfun_0.61 fs_2.1.0
## [103] MatrixGenerics_1.22.0 pkgconfig_2.0.3