Building QtlDataset and MultiStudyQtlDataset objects
pecotmr authors
2026-09-04
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.
It is a MultiAssayExperiment,
so the whole multi-assay surface works on it directly. Each QTL context
is one experiment, and the genotypes are one more:
-
genotype—rowRangesis one range per variant (carryingA1/A2); thedosageassay is variants x samples, read lazily from the genotype files;colDataholds genotype-derived covariates such as ancestry PCs. -
one per context —
rowRangesis trait positions, the assay is traits x samples, andcolDataholds that context’s phenotype covariates.
Alongside them the object keeps a colData naming every
sample and a sampleMap recording which of them each
experiment observes — so contexts need not share a sample set — plus its
own slots:
-
study— study identifier (length-1 character). -
genotypes— the genotype panel (a lazy PLINK1/PLINK2/VCF/GDS source, read withreadGenotypes()) that the dosage assay reads through. - QC thresholds —
mafCutoff,macCutoff,xvarCutoff,imissCutoff,keepIndel,keepVariants,scaleResiduals. Stored as lazy filters, applied at genotype-extraction time.
To restrict a dataset to a set of samples, subset it:
x[, samples, ] narrows colData and
sampleMap together, and keeps the class and its slots. The
keepSamples constructor argument is shorthand for the same
thing.
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 the genotype panel yourself 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(colnames(readGenotypes(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/RtmppU8QTb/temp_libpath984aa04109/pecotmr/extdata/toy_ref
## Genotype covariates: 0 cols
## Samples: 165
## 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 genotype panel.
# 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 <- readGenotypes(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/RtmppU8QTb/temp_libpath984aa04109/pecotmr/extdata/toy_ref
## Genotype covariates: 0 cols
## Samples: 165
## 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(getQtlDatasets(msdManual))## [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.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.7.8
##
## 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-11
## [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.6
## [35] ggplot2_4.0.3 scales_1.4.0
## [37] MASS_7.3-66 MultiAssayExperiment_1.36.1
## [39] SummarizedExperiment_1.40.0 cli_3.6.6
## [41] rmarkdown_2.32 metafor_5.0-1
## [43] crayon_1.5.3 ragg_1.5.2
## [45] generics_0.1.4 otel_0.2.0
## [47] RcppParallel_5.1.11-2 tzdb_0.5.0
## [49] BiocBaseUtils_1.12.0 cachem_1.1.0
## [51] stringr_1.6.0 splines_4.5.3
## [53] metadat_1.6-0 parallel_4.5.3
## [55] XVector_0.50.0 matrixStats_1.5.0
## [57] vctrs_0.7.3 Matrix_1.7-6
## [59] jsonlite_2.0.0 IRanges_2.44.0
## [61] hms_1.1.4 S4Vectors_0.48.0
## [63] bit64_4.8.6 mixsqp_0.3-54
## [65] irlba_2.3.7 archive_1.1.13
## [67] systemfonts_1.3.2 tidyr_1.3.2
## [69] jquerylib_0.1.4 snpStats_1.60.0
## [71] glue_1.8.1 pkgdown_2.2.1
## [73] codetools_0.2-20 stringi_1.8.9
## [75] gtable_0.3.6 GenomicRanges_1.62.1
## [77] quadprog_1.5-8 tibble_3.3.1
## [79] pillar_1.11.1 htmltools_0.5.9
## [81] Seqinfo_1.0.0 R6_2.6.1
## [83] zigg_0.0.2 textshaping_1.0.5
## [85] vroom_1.7.1 evaluate_1.0.5
## [87] lattice_0.23-1 Biobase_2.70.0
## [89] readr_2.2.0 tictoc_1.2.1
## [91] Rsamtools_2.26.0 Rfast_2.1.5.2
## [93] bslib_0.12.0 Rcpp_1.1.2
## [95] SparseArray_1.10.8 nlme_3.1-170
## [97] xfun_0.60 fs_2.1.0
## [99] MatrixGenerics_1.22.0 pkgconfig_2.0.3