Building QtlSumStats and GwasSumStats objects
pecotmr authors
2026-10-22
Source:vignettes/constructing-sumstats.Rmd
constructing-sumstats.RmdOverview
GwasSumStats and QtlSumStats are
DFrame-subclass collections of summary statistics. Each
row carries one entry — a
GRanges of per-variant statistics — annotated by an
identity tuple:
| Class | Row key | One row is |
|---|---|---|
GwasSumStats |
study |
one GWAS study over a region |
QtlSumStats |
(study, context, trait) |
one molecular trait in one context |
Every entry GRanges must carry these mcols:
SNP, A1 (effect / dosage-counted allele),
A2 (other allele), Z, N;
optionally MAF, INFO, BETA,
SE, P. Two settings apply to the whole
collection and are stored in slots, not columns:
-
genome— the build label (e.g."GRCh38"). -
ldReference— the LD reference used for QC and downstream analyses: a genotype panel (aRangedSummarizedExperiment, as returned byreadGenotypes()), aGenotypeMixture, or anLdMatrixholding a correlation matrix you built (see the LD reference vignette).
Optional per-row columns: varY (phenotype variance) and,
for case/control GWAS, nCase / nControl.
QtlSumStats and GwasSumStats can be
constructed in two ways:
-
From a manifest —
loadGwasSumStatsFromManifest()/loadQtlSumStatsFromManifest()read a table pointing at the sumstats files and assemble the collection. -
Manually — build the
entryGRanges yourself and call theGwasSumStats()/QtlSumStats()constructor.
The loaders run no QC. They build the raw collection with an empty
qcInfoslot. Allele harmonization, MAF/INFO/N filtering, and LD-mismatch QC are a separate, mandatory step —summaryStatsQc()— documented in the Quality Control for GWAS and QTL Summary Statistics vignette. Downstream pipelines reject collections whoseqcInfois empty.
Toy data for this vignette
We will use the toy_ref PLINK1 panel (chr22) as the LD
reference, and a handful of its variants as our summary statistics so
the LD-overlap check succeeds.
genoPrefix <- sub(
"\\.bed$",
"",
system.file("extdata", "toy_ref.bed", package = "pecotmr")
)
ldReference <- readGenotypes(plink1Prefix = genoPrefix)
# Five real toy_ref chr22 sites (positions and alleles taken from the panel).
sumstatsDf <- data.frame(
chrom = rep("22", 5),
pos = c(14560203, 14564328, 14850625, 14870204, 14878387),
variant_id = c(
"rs11089128",
"rs7288972",
"rs11167319",
"rs8138488",
"rs2186521"
),
A1 = c("G", "C", "G", "C", "A"),
A2 = c("A", "T", "T", "T", "G"),
z = c(1.20, -0.42, 2.13, 0.05, -1.83),
n_sample = rep(10000L, 5),
stringsAsFactors = FALSE
)
workdir <- file.path(tempdir(), "sumstats-vignette")
dir.create(workdir, showWarnings = FALSE)
gwasTsv <- file.path(workdir, "study1.tsv")
readr::write_tsv(sumstatsDf, gwasTsv)GwasSumStats from a manifest
The GWAS manifest has one row per study.
genome and ldReference are collection-level,
so they are supplied as arguments (or as constant genome /
ldReferencePath columns).
| Column | Required | Meaning |
|---|---|---|
study (alias study_id) |
yes, unique | study identifier |
sumStatsPath (aliases path,
file_path) |
yes | sumstats file |
columnMapping (alias column_mapping) |
no | YAML column map (see below) |
nCase / nControl
|
no | case / control counts |
varY |
no | phenotype variance |
genome, ldReferencePath
|
single-valued | may instead be arguments |
gwasManifest <- data.frame(
study = "MyGWAS",
sumStatsPath = gwasTsv,
stringsAsFactors = FALSE
)
gws <- loadGwasSumStatsFromManifest(
gwasManifest,
genome = "GRCh38",
ldReference = ldReference
)
gws## GwasSumStats: 1 studies, genome build GRCh38
## LD reference: plink1 @ /tmp/RtmpmBLZhq/temp_libpathaae388f3498/pecotmr/extdata/toy_ref
S4Vectors::mcols(gws[[1]])## DataFrame with 5 rows and 6 columns
## SNP A1 A2 Z N z_supplied
## <character> <character> <character> <numeric> <numeric> <logical>
## 1 rs11089128 G A 1.20 10000 TRUE
## 2 rs7288972 C T -0.42 10000 TRUE
## 3 rs11167319 G T 2.13 10000 TRUE
## 4 rs8138488 C T 0.05 10000 TRUE
## 5 rs2186521 A G -1.83 10000 TRUEThe ldReference argument also accepts a path/prefix or a
chromosome-sharded genoMeta specification (see
?readGenotypes), so a genome-wide reference can be supplied
as a #chr,path meta file rather than a single panel.
GwasSumStats manually
The loader is a wrapper over the GwasSumStats()
constructor. Built by hand, the pattern is: make one
GRanges per study with the required mcols,
then call the constructor with the shared genome and
ldReference.
entry <- GenomicRanges::GRanges(
seqnames = paste0("chr", sumstatsDf$chrom),
ranges = IRanges::IRanges(start = sumstatsDf$pos, width = 1L)
)
S4Vectors::mcols(entry) <- S4Vectors::DataFrame(
SNP = sumstatsDf$variant_id,
A1 = sumstatsDf$A1,
A2 = sumstatsDf$A2,
Z = sumstatsDf$z,
N = sumstatsDf$n_sample
)
gwsManual <- GwasSumStats(
studyName = "MyGWAS",
entry = list(entry),
genome = "GRCh38",
ldReference = ldReference
)
gwsManual## GwasSumStats: 1 studies, genome build GRCh38
## LD reference: plink1 @ /tmp/RtmpmBLZhq/temp_libpathaae388f3498/pecotmr/extdata/toy_refFor a case/control study, pass nCase /
nControl; for sufficient statistics interface with
continuous traits, pass varY.
QtlSumStats from a manifest
The QTL manifest has one row per
(study, context, trait) tuple. Each row points at
a sumstats file for that trait.
| Column | Required | Meaning |
|---|---|---|
study |
yes | study identifier |
context (alias cond) |
yes | context / condition |
trait (alias gene_id) |
yes | molecular trait id |
sumStatsPath (alias path) |
yes | sumstats file for this trait |
columnMapping, varY
|
no | as for GWAS |
genome, ldReferencePath
|
single-valued | may instead be arguments |
qtlManifest <- data.frame(
study = "MyQTL",
context = "Blood",
trait = "ENSG00000001",
sumStatsPath = gwasTsv, # reuse the toy table for one gene
stringsAsFactors = FALSE
)
qts <- loadQtlSumStatsFromManifest(
qtlManifest,
genome = "GRCh38",
ldReference = ldReference
)
qts## QtlSumStats: 1 entries, genome build GRCh38
## 1 studies, 1 contexts, 1 traits
## LD reference: plink1 @ /tmp/RtmpmBLZhq/temp_libpathaae388f3498/pecotmr/extdata/toy_refQtlSumStats manually
qtsManual <- QtlSumStats(
studyName = "MyQTL",
context = "Blood",
trait = "ENSG00000001",
entry = list(entry),
genome = "GRCh38",
ldReference = ldReference
)
qtsManual## QtlSumStats: 1 entries, genome build GRCh38
## 1 studies, 1 contexts, 1 traits
## LD reference: plink1 @ /tmp/RtmpmBLZhq/temp_libpathaae388f3498/pecotmr/extdata/toy_refFor a multi-context collection (e.g. as input to mash),
pass parallel study / context /
trait vectors and a matching list of entry
GRanges — one element per tuple.
From effect matrices
Multi-context work often starts from matrices — variants by contexts — rather than per-context tables. mash in particular produces and consumes that shape, so two builders take it directly instead of making you assemble entries by hand.
qtlSumStatsFromZMatrix() takes z-scores, with variant
IDs as row names and contexts as column names:
data(mashInputExample)
z <- mashInputExample$strong.z[1:6, 1:3]
rownames(z) <- paste0("chr1:", seq(1000, by = 100, length.out = 6), ":A:G")
qtlSumStatsFromZMatrix(z, studyName = "demo")## QtlSumStats: 3 entries, genome build GRCh38
## 1 studies, 3 contexts, 1 traits
## LD reference: none (LD-free)qtlSumStatsFromBetaMatrix() is the effect-size
counterpart, taking bhat and shat:
b <- mashInputExample$strong.b[1:6, 1:3]
shat <- mashInputExample$strong.s[1:6, 1:3]
rownames(b) <- rownames(shat) <- rownames(z)
contexts(qtlSumStatsFromBetaMatrix(b, shat, studyName = "demo"))## [1] "ALL" "Ast" "End"Both default ldReference = NULL, giving an LD-free
collection. That is fine for mash,
which never needs LD, but any RSS method does — supply a panel if the
result is headed for fine-mapping or TWAS-weight learning.
The remaining arguments (genome, n,
a1, a2, trait) fill in the
metadata a matrix cannot carry. The allele defaults are placeholders, so
set them if the alleles matter downstream.
Reading from files
Column mapping (YAML)
When columns in a sumstats file don’t match the standard names,
supply a columnMapping — a YAML file of
standardName: sourceName entries (the xqtl-protocol
format), either per-row in the manifest or as a default argument.
Standard keys are chrom, pos,
variant_id, A1, A2,
z, n_sample (and optional beta,
se, p, maf,
info).
oddDf <- sumstatsDf
names(oddDf) <- c("CHR", "POS", "RSID", "EA", "OA", "ZSCORE", "SAMPLE_N")
oddTsv <- file.path(workdir, "odd_columns.tsv")
readr::write_tsv(oddDf, oddTsv)
mapPath <- file.path(workdir, "column_map.yaml")
yaml::write_yaml(
list(
chrom = "CHR",
pos = "POS",
variant_id = "RSID",
A1 = "EA",
A2 = "OA",
z = "ZSCORE",
n_sample = "SAMPLE_N"
),
mapPath
)
gwsMapped <- loadGwasSumStatsFromManifest(
data.frame(
study = "MyGWAS",
sumStatsPath = oddTsv,
columnMapping = mapPath
),
genome = "GRCh38",
ldReference = ldReference
)
length(gwsMapped[[1]])## [1] 5Region-restricted reads (tabix)
If a region is supplied and the sumstats file is
compressed with bgzip and indexed with tabix (a
<file>.tbi must exist in the same directory), only
that region is read. Plain-text files are read whole and the region is
ignored with a warning. The tabix header line must begin with
# so its column names can be recovered.
tabixDf <- sumstatsDf
names(tabixDf)[1] <- "#chrom" # tabix header convention
tabixTsv <- file.path(workdir, "study1.sorted.tsv")
readr::write_tsv(tabixDf, tabixTsv)
gz <- Rsamtools::bgzip(
tabixTsv,
dest = paste0(tabixTsv, ".gz"),
overwrite = TRUE
)
Rsamtools::indexTabix(gz, seq = 1L, start = 2L, end = 2L, comment = "#")## [1] "/tmp/Rtmp0mJTlQ/sumstats-vignette/study1.sorted.tsv.gz.tbi"
gwsRegion <- loadGwasSumStatsFromManifest(
data.frame(study = "MyGWAS", sumStatsPath = gz),
genome = "GRCh38",
ldReference = ldReference,
region = "22:14560000-14565000"
)
length(gwsRegion[[1]]) # only the two variants in the region## [1] 2GWAS-VCF
VCF summary statistics are read as GWAS-VCF: the
effect allele A1 is the ALT allele,
A2 is REF, and the statistics come from the
per-study FORMAT fields ES, SE,
LP, SS, EAF
(Z = ES / SE, N = SS, P = 10^-LP,
MAF = fold(EAF)). Reading only a specific requires the VCF
to be compressed with bgzip and indexed with tabix. BCF is not
supported — convert to a bgzipped VCF first. Support for
reading VCFs requires the VariantAnnotation package.
vcfLines <- c(
"##fileformat=VCFv4.2",
'##FORMAT=<ID=ES,Number=A,Type=Float,Description="effect size">',
'##FORMAT=<ID=SE,Number=A,Type=Float,Description="standard error">',
'##FORMAT=<ID=LP,Number=A,Type=Float,Description="-log10 p">',
'##FORMAT=<ID=SS,Number=A,Type=Float,Description="sample size">',
"##contig=<ID=22>",
paste(
"#CHROM",
"POS",
"ID",
"REF",
"ALT",
"QUAL",
"FILTER",
"INFO",
"FORMAT",
"MyGWAS",
sep = "\t"
)
)
for (i in seq_len(nrow(sumstatsDf))) {
vcfLines <- c(
vcfLines,
paste(
"22",
sumstatsDf$pos[i],
sumstatsDf$variant_id[i],
sumstatsDf$A2[i],
sumstatsDf$A1[i], # REF = A2, ALT = A1
".",
"PASS",
".",
"ES:SE:LP:SS",
sprintf("%.3f:0.100:1.000:10000", sumstatsDf$z[i] * 0.1),
sep = "\t"
)
)
}
vcfPath <- file.path(workdir, "study1.vcf")
writeLines(vcfLines, vcfPath)
gwsVcf <- loadGwasSumStatsFromManifest(
data.frame(study = "MyGWAS", sumStatsPath = vcfPath),
genome = "GRCh38",
ldReference = ldReference
)
S4Vectors::mcols(gwsVcf[[1]])[, c("SNP", "A1", "A2", "Z", "N")]## DataFrame with 5 rows and 5 columns
## SNP A1 A2 Z N
## <character> <character> <character> <numeric> <numeric>
## 1 rs11089128 G A 1.20 10000
## 2 rs7288972 C T -0.42 10000
## 3 rs11167319 G T 2.13 10000
## 4 rs8138488 C T 0.05 10000
## 5 rs2186521 A G -1.83 10000For a multi-sample GWAS-VCF, pass
sampleSelect = "<column>" to pick the study; override
the FORMAT tag names with formatMapping if they differ from
ES/SE/LP/SS/EAF.
Next steps
The collections carry no QC yet. Run summaryStatsQc()
before any downstream analysis (it populates qcInfo;
downstream pipelines require it):
gws_qcd <- summaryStatsQc(gws, ldMismatchQcMethod = "slalom")See the Quality Control for GWAS and QTL Summary Statistics
vignette for the full QC pipeline. 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.
unname(GenomeInfoDb::genome(gws))## [1] "GRCh38"
as.character(gws$study)## [1] "MyGWAS"## [1] 5## 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.3
##
## loaded via a namespace (and not attached):
## [1] DBI_1.3.0 bitops_1.1-0
## [3] httr2_1.3.0 rlang_1.3.0
## [5] magrittr_2.0.5 otel_0.2.0
## [7] matrixStats_1.5.0 susieR_0.16.6
## [9] compiler_4.5.3 RSQLite_3.53.3
## [11] GenomicFeatures_1.62.0 png_0.1-9
## [13] systemfonts_1.3.2 vctrs_0.7.3
## [15] quadprog_1.5-8 stringr_1.6.0
## [17] pkgconfig_2.0.3 crayon_1.5.3
## [19] fastmap_1.2.0 dbplyr_2.6.0
## [21] backports_1.5.1 XVector_0.50.0
## [23] Rsamtools_2.26.0 rmarkdown_2.32
## [25] tzdb_0.5.0 UCSC.utils_1.6.1
## [27] ragg_1.5.2 purrr_1.2.2
## [29] bit_4.6.0 xfun_0.61
## [31] Rfast_2.1.5.2 MultiAssayExperiment_1.36.2
## [33] cachem_1.1.0 cigarillo_1.0.0
## [35] GenomeInfoDb_1.46.2 jsonlite_2.0.0
## [37] blob_1.3.0 reshape_0.8.10
## [39] DelayedArray_0.36.1 tictoc_1.2.1
## [41] BiocParallel_1.44.0 irlba_2.4.1
## [43] parallel_4.5.3 R6_2.6.1
## [45] VariantAnnotation_1.56.0 bslib_0.12.0
## [47] stringi_1.8.9 RColorBrewer_1.1-3
## [49] rtracklayer_1.70.1 GenomicRanges_1.62.1
## [51] jquerylib_0.1.4 numDeriv_2016.8-1.1
## [53] Rcpp_1.1.2 Seqinfo_1.0.0
## [55] SummarizedExperiment_1.40.0 knitr_1.52
## [57] readr_2.2.0 IRanges_2.44.0
## [59] BiocBaseUtils_1.12.0 Matrix_1.7-6
## [61] tidyselect_1.2.1 abind_1.4-8
## [63] yaml_2.3.12 codetools_0.2-20
## [65] metafor_5.0-1 curl_8.0.0
## [67] plyr_1.8.9 lattice_0.23-1
## [69] tibble_3.3.1 withr_3.0.3
## [71] S7_0.2.2 Biobase_2.70.0
## [73] KEGGREST_1.50.0 evaluate_1.0.5
## [75] archive_1.1.13 desc_1.4.3
## [77] BiocFileCache_3.0.0 RcppParallel_5.1.11-2
## [79] Biostrings_2.78.0 filelock_1.0.3
## [81] pillar_1.11.1 MatrixGenerics_1.22.0
## [83] metadat_1.6-0 checkmate_2.3.3
## [85] stats4_4.5.3 generics_0.1.4
## [87] vroom_1.7.1 mathjaxr_2.0-0
## [89] RCurl_1.98-1.19 S4Vectors_0.48.1
## [91] hms_1.1.4 ggplot2_4.0.3
## [93] scales_1.4.0 glue_1.8.1
## [95] tools_4.5.3 BiocIO_1.20.0
## [97] BSgenome_1.78.0 GenomicAlignments_1.46.0
## [99] fs_2.1.0 XML_3.99-0.25
## [101] grid_4.5.3 PlinkMatrix_1.0.0
## [103] tidyr_1.3.2 AnnotationDbi_1.72.0
## [105] nlme_3.1-170 restfulr_0.0.17
## [107] cli_3.6.6 zigg_0.0.2
## [109] textshaping_1.0.5 mixsqp_0.3-54
## [111] S4Arrays_1.10.1 dplyr_1.2.1
## [113] gtable_0.3.6 sass_0.4.10
## [115] digest_0.6.39 BiocGenerics_0.56.0
## [117] SparseArray_1.10.10 rjson_0.2.23
## [119] htmlwidgets_1.6.4 farver_2.1.2
## [121] memoise_2.0.1 htmltools_0.5.9
## [123] pkgdown_2.2.1 lifecycle_1.0.5
## [125] httr_1.4.9 bit64_4.8.6
## [127] MASS_7.3-66