Building QtlDataset and MultiStudyQtlDataset objects
pecotmr authors
2026-07-21
Source:vignettes/constructing-qtl-datasets.Rmd
constructing-qtl-datasets.RmdOverview
A QtlDataset is an individual-level container for a
single study: one genotype source plus a per-context
set of molecular phenotypes. MultiStudyQtlDataset bundles
several QtlDatasets (optionally alongside a
summary-statistics-only QtlSumStats) so a molecular trait
measured across multiple studies can be analyzed with downstream
pipelines together.
| Slot | Holds |
|---|---|
study |
study identifier (length-1 character) |
genotypes |
a GenotypeHandle (lazy PLINK1/PLINK2/VCF/GDS
handle) |
phenotypes |
a named list of SummarizedExperiments,
one per context; assay = traits x samples, rowRanges =
trait positions, colData = per-context covariates |
genotypeCovariates |
numeric matrix (samples x covariates), e.g. genotype PCs, applied to every context |
| QC thresholds |
mafCutoff, macCutoff,
xvarCutoff, imissCutoff,
keepIndel, keepSamples,
keepVariants, scaleResiduals — stored as
lazy filters, applied at genotype-extraction time |
There are two ways to build a QtlDataset:
-
From a manifest —
loadQtlDatasetFromManifest()reads a table (data.frameor file) that points at the phenotype / covariate / genotype files and assembles the object for you. This is the recommended route for real data. -
Manually — build the per-context
SummarizedExperiments and theGenotypeHandleyourself and call theQtlDataset()constructor. Useful when your data is already in memory.
Toy data for this vignette
We will reuse the toy_ref PLINK1 panel shipped with the
package as the genotype source, take its first few samples as our
phenotyped individuals, and write two small per-context phenotype BED
files (TensorQTL-style layout: #chr, start,
end, gene_id, then one column per sample) plus
a covariate table to a temporary directory.
genoPrefix <- sub("\\.bed$", "",
system.file("extdata", "toy_ref.bed", package = "pecotmr"))
samples <- head(getSampleIds(GenotypeHandle(plink1Prefix = genoPrefix)), 40)
workdir <- file.path(tempdir(), "qtl-vignette")
dir.create(workdir, showWarnings = FALSE)
# Two genes on chr22, positioned inside the toy_ref variant span so a cis
# window captures nearby variants.
writePhenoBed <- function(path, seed) {
set.seed(seed)
bed <- data.frame(
`#chr` = c("chr22", "chr22"),
start = c(14600000L, 14900000L),
end = c(14600001L, 14900001L),
gene_id = c("GENE1", "GENE2"),
check.names = FALSE)
for (s in samples) bed[[s]] <- round(rnorm(2), 3)
readr::write_tsv(bed, path)
path
}
muscleBed <- writePhenoBed(file.path(workdir, "Muscle.bed"), 1)
bloodBed <- writePhenoBed(file.path(workdir, "Blood.bed"), 2)
# Per-context covariate table: samples as rows, covariate columns after an ID
# column (here two genotype PCs).
covPath <- file.path(workdir, "covariates.tsv")
set.seed(3)
readr::write_tsv(
data.frame(sample = samples,
PC1 = round(rnorm(length(samples)), 3),
PC2 = round(rnorm(length(samples)), 3),
check.names = FALSE),
covPath)Constructing a QtlDataset from a manifest
The QtlDataset manifest has one row per context.
Canonical column names are camelCase; snake_case aliases such as
cond, cov_path, genotype_prefix
are also accepted and normalized on read.
| Column | Required | Meaning |
|---|---|---|
context (alias cond) |
yes | context / condition label |
phenotypePath (alias path) |
yes | bgzipped or plain BED phenotype file |
covariatePath (alias cov_path) |
no | per-context covariate table |
study |
single-valued | study identifier (may also be an argument) |
genotypePath (alias genotype_prefix) |
single-valued | genotype file or PLINK prefix (may also be an argument) |
genotypeCovariatePath |
no, single-valued | genotype-derived covariates applied to all contexts |
study and genotypePath are single-valued
for the whole dataset, so you can supply them either as a constant
manifest column or as an argument; if you give both,
they must agree.
manifest <- data.frame(
context = c("Muscle", "Blood"),
phenotypePath = c(muscleBed, bloodBed),
covariatePath = c(covPath, covPath),
study = "GTEx_toy",
genotypePath = genoPrefix,
stringsAsFactors = FALSE)
# The path columns hold absolute paths; show basenames for readability
# (the real `manifest` keeps the full paths the loader needs).
transform(manifest,
phenotypePath = basename(phenotypePath),
covariatePath = basename(covariatePath),
genotypePath = basename(genotypePath))## context phenotypePath covariatePath study genotypePath
## 1 Muscle Muscle.bed covariates.tsv GTEx_toy toy_ref
## 2 Blood Blood.bed covariates.tsv GTEx_toy toy_ref
qd <- loadQtlDatasetFromManifest(manifest)
qd## QtlDataset for study 'GTEx_toy'
## 2 context(s): Muscle, Blood
## 2 unique traits across contexts
## Genotypes: plink1 @ /tmp/Rtmpi0K25n/temp_libpatha6b6057489c/pecotmr/extdata/toy_ref
## Genotype covariates: 0 cols
## Scale residuals: TRUEThe manifest can equally be a path to a .csv (read as
CSV) or any other extension (read as TSV):
manifestPath <- file.path(workdir, "qtl_manifest.tsv")
readr::write_tsv(manifest, manifestPath)
qd <- loadQtlDatasetFromManifest(manifestPath)The QC thresholds are ordinary arguments (they are stored as lazy filters, not applied at load), for example:
qdFiltered <- loadQtlDatasetFromManifest(
manifest, mafCutoff = 0.01, imissCutoff = 0.05, keepIndel = FALSE)Constructing a QtlDataset manually
The manifest loader is a thin wrapper over the
QtlDataset() constructor. Doing it by hand shows exactly
what the object contains: a named list of
SummarizedExperiments and a
GenotypeHandle.
# One SummarizedExperiment per context. Rows are traits, columns are samples;
# rowRanges carries each trait's position, colData carries per-context
# covariates.
buildSE <- function(bedPath) {
bed <- as.data.frame(readr::read_tsv(bedPath, show_col_types = FALSE),
check.names = FALSE)
sampleCols <- setdiff(names(bed), c("#chr", "start", "end", "gene_id"))
expr <- as.matrix(bed[, sampleCols, drop = FALSE])
rownames(expr) <- bed$gene_id
rr <- GenomicRanges::GRanges(
seqnames = bed$`#chr`,
ranges = IRanges::IRanges(start = bed$start + 1L, end = bed$end))
names(rr) <- bed$gene_id
cov <- as.data.frame(readr::read_tsv(covPath, show_col_types = FALSE),
check.names = FALSE)
cd <- S4Vectors::DataFrame(cov[match(sampleCols, cov$sample),
c("PC1", "PC2"), drop = FALSE],
row.names = sampleCols)
SummarizedExperiment::SummarizedExperiment(
assays = list(expression = expr), rowRanges = rr, colData = cd)
}
phenotypes <- list(Muscle = buildSE(muscleBed), Blood = buildSE(bloodBed))
genotypes <- GenotypeHandle(plink1Prefix = genoPrefix)
qdManual <- QtlDataset(
study = "GTEx_toy",
genotypes = genotypes,
phenotypes = phenotypes)
qdManual## QtlDataset for study 'GTEx_toy'
## 2 context(s): Muscle, Blood
## 2 unique traits across contexts
## Genotypes: plink1 @ /tmp/Rtmpi0K25n/temp_libpatha6b6057489c/pecotmr/extdata/toy_ref
## Genotype covariates: 0 cols
## Scale residuals: TRUEThis approach gives users more flexibility if they need to make more
modifications to the molecular phenotype data after loading it from a
file but before making a QtlDataset
MultiStudyQtlDataset
A MultiStudyQtlDataset needs at least two studies in
total - at least two individual-level studies or at least one
individual-level study and one summary-only study. When loading from a
manifest, a qtlDatasetsManifest with
study and genotypePath must be provided where
rows are grouped by study and one QtlDataset
is built per group. A sumStatsManifest describing
summary-only studies can also be provided (see the Constructing
QtlSumStats and GwasSumStats objects vignette).
multiManifest <- data.frame(
study = c("StudyA", "StudyB"),
context = c("Muscle", "Muscle"),
phenotypePath = c(muscleBed, bloodBed),
genotypePath = genoPrefix,
stringsAsFactors = FALSE)
msd <- loadMultiStudyQtlDatasetFromManifest(multiManifest)
msd## MultiStudyQtlDataset: 2 individual-level + 0 sumstats studies
## Individual-level studies: StudyA, StudyBA MultiStudyQtlDataset can be manually constructed by building a
named list of QtlDatasets and passing it to the
constructor:
msdManual <- MultiStudyQtlDataset(
qtlDatasets = list(StudyA = qd, StudyB = qdManual))
names(msdManual@qtlDatasets)## [1] "StudyA" "StudyB"To attach a summary-only study, pass a QtlSumStats
object to the sumStats argument. manifest form, supply
sumStatsManifest, genome, and
ldSketch).
Next steps
Once built, the accessors extract data lazily, applying the stored QC filters at extraction time:
getStudy(qd)## [1] "GTEx_toy"
getContexts(qd)## [1] "Muscle" "Blood"
# Cis genotypes around a trait (+/- a window), QC applied on the fly:
dim(getGenotypes(qd, traitId = "GENE1", cisWindow = 1e6))## [1] 165 58
# Phenotype matrix for one context:
dim(SummarizedExperiment::assay(getPhenotypes(qd, contexts = "Muscle")))## [1] 2 40QtlDataset and MultiStudyQtlDataset are
direct inputs to fineMappingPipeline(),
twasWeightsPipeline() 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.4 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/runner/work/pecotmr/pecotmr/.pixi/envs/default/lib/libopenblasp-r0.3.33.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.6.7
##
## 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.0-9
## [7] fastmap_1.2.0 reshape_0.8.10
## [9] mathjaxr_2.0-0 digest_0.6.39
## [11] lifecycle_1.0.5 survival_3.8-9
## [13] magrittr_2.0.5 compiler_4.5.3
## [15] rlang_1.3.0 sass_0.4.10
## [17] tools_4.5.3 yaml_2.3.12
## [19] knitr_1.51 S4Arrays_1.10.1
## [21] htmlwidgets_1.6.4 bit_4.6.0
## [23] DelayedArray_0.36.0 plyr_1.8.9
## [25] RColorBrewer_1.1-3 abind_1.4-8
## [27] BiocParallel_1.44.0 purrr_1.2.2
## [29] numDeriv_2016.8-1.1 BiocGenerics_0.56.0
## [31] desc_1.4.3 grid_4.5.3
## [33] stats4_4.5.3 susieR_0.16.4
## [35] ggplot2_4.0.3 scales_1.4.0
## [37] MASS_7.3-66 SummarizedExperiment_1.40.0
## [39] cli_3.6.6 rmarkdown_2.31
## [41] metafor_5.0-1 crayon_1.5.3
## [43] ragg_1.5.2 generics_0.1.4
## [45] otel_0.2.0 RcppParallel_5.1.11-2
## [47] tzdb_0.5.0 cachem_1.1.0
## [49] stringr_1.6.0 splines_4.5.3
## [51] metadat_1.6-0 parallel_4.5.3
## [53] XVector_0.50.0 matrixStats_1.5.0
## [55] vctrs_0.7.3 Matrix_1.7-5
## [57] jsonlite_2.0.0 IRanges_2.44.0
## [59] hms_1.1.4 S4Vectors_0.48.0
## [61] bit64_4.8.2 mixsqp_0.3-54
## [63] irlba_2.3.7 systemfonts_1.3.2
## [65] tidyr_1.3.2 jquerylib_0.1.4
## [67] snpStats_1.60.0 glue_1.8.1
## [69] pkgdown_2.2.1 codetools_0.2-20
## [71] stringi_1.8.7 gtable_0.3.6
## [73] GenomicRanges_1.62.1 quadprog_1.5-8
## [75] tibble_3.3.1 pillar_1.11.1
## [77] htmltools_0.5.9 Seqinfo_1.0.0
## [79] R6_2.6.1 zigg_0.0.2
## [81] textshaping_1.0.5 vroom_1.7.1
## [83] evaluate_1.0.5 lattice_0.22-9
## [85] Biobase_2.70.0 readr_2.2.0
## [87] Rsamtools_2.26.0 tictoc_1.2.1
## [89] Rfast_2.1.5.2 bslib_0.11.0
## [91] Rcpp_1.1.2 SparseArray_1.10.8
## [93] nlme_3.1-170 xfun_0.60
## [95] fs_2.1.0 MatrixGenerics_1.22.0
## [97] pkgconfig_2.0.3