Building QtlDataset and MultiStudyQtlDataset objects
pecotmr authors
2026-10-22
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/BGEN 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/RtmpmBLZhq/temp_libpathaae388f3498/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,
genotypeFilterParam = GenotypeFilterParam(
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(
studyName = "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/RtmpmBLZhq/temp_libpathaae388f3498/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(qtlDatasets(msdManual))## [1] "StudyA" "StudyB"To attach a summary-only study, pass a QtlSumStats
object to the sumStats argument. For the manifest form,
supply sumStatsManifest, genome, and
ldReference.
Excluding related samples
Both fineMappingPipeline() and
twasWeightsPipeline() fit models that assume independent
samples. Related individuals break that assumption: they inflate the
effective sample size, so cross-validated R-squared comes back
optimistic and PIPs overconfident. Relatedness is therefore a QC step
applied before the dataset is built — the pipelines do not
model it.
filterRelatedness() takes a kinship table — one row per
pair, the shape PLINK’s --genome or KING emits — and
returns the individuals to drop:
relatedness <- data.frame(
IID1 = c(samples[1], samples[3], samples[5]),
IID2 = c(samples[2], samples[4], samples[6]),
PI_HAT = c(0.50, 0.25, 0.125),
stringsAsFactors = FALSE
)
excluded <- filterRelatedness(relatedness, relatednessThreshold = 0.0625)
excluded## [1] "NA06989" "NA11843" "NA12739"## [1] 37One member of each related pair comes back, chosen so that what
remains is as large an unrelated set as possible.
relatednessThreshold defaults to 0.0625, which is
third-degree relatives; raise it to 0.125 to keep everyone but
second-degree and closer.
The column names above are the defaults. A table that spells them
differently is described with
RelatednessColumnsParam():
filterRelatedness(
relatedness,
columns = RelatednessColumnsParam(iid1 = "ID1", iid2 = "ID2", value = "KINSHIP")
)Building the description rather than passing a bare list means a misspelled field is an error at the call site instead of a default silently left in place.
The exclusion list is not applied for you — feed it to
keepSamples, the constructor argument this vignette
introduced above, or subset an existing dataset with
qd[, unrelated, ]. Both narrow colData and
sampleMap together.
analysisType changes what is being maximised. The
default, "maximizeUnrelated", keeps the largest unrelated
set. "maximizeCases" instead retains as many cases as it
can — the right objective for a case/control design, where cases are the
scarce resource — and requires phenoData naming each
sample’s status.
filterRelatedness() needs the suggested packages
igraph and plinkQC; the chunk above is
skipped when they are not installed.
Next steps
Once built, the accessors extract data lazily, applying the stored QC filters at extraction time:
studyName(qd)## [1] "GTEx_toy"
contexts(qd)## [1] "Muscle" "Blood"
# Cis genotypes around a trait (+/- a window), QC applied on the fly:
dim(genotypes(qd, traitId = "GENE1", cisWindow = 1e6))## [1] 165 58
# Phenotype matrix for one context:
dim(SummarizedExperiment::assay(molecularTraits(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.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 gridExtra_2.3.1
## [5] rlang_1.3.0 magrittr_2.0.5
## [7] otel_0.2.0 matrixStats_1.5.0
## [9] susieR_0.16.6 compiler_4.5.3
## [11] RSQLite_3.53.3 systemfonts_1.3.2
## [13] vctrs_0.7.3 quadprog_1.5-8
## [15] stringr_1.6.0 pkgconfig_2.0.3
## [17] crayon_1.5.3 fastmap_1.2.0
## [19] dbplyr_2.6.0 backports_1.5.1
## [21] XVector_0.50.0 Rsamtools_2.26.0
## [23] rmarkdown_2.32 tzdb_0.5.0
## [25] ragg_1.5.2 UpSetR_1.4.1
## [27] purrr_1.2.2 bit_4.6.0
## [29] xfun_0.61 Rfast_2.1.5.2
## [31] MultiAssayExperiment_1.36.2 randomForest_4.7-1.2
## [33] cachem_1.1.0 jsonlite_2.0.0
## [35] blob_1.3.0 DelayedArray_0.36.1
## [37] reshape_0.8.10 tictoc_1.2.1
## [39] BiocParallel_1.44.0 irlba_2.4.1
## [41] parallel_4.5.3 R6_2.6.1
## [43] bslib_0.12.0 stringi_1.8.9
## [45] RColorBrewer_1.1-3 plinkQC_1.1.0
## [47] GenomicRanges_1.62.1 jquerylib_0.1.4
## [49] numDeriv_2016.8-1.1 Rcpp_1.1.2
## [51] Seqinfo_1.0.0 SummarizedExperiment_1.40.0
## [53] knitr_1.52 readr_2.2.0
## [55] IRanges_2.44.0 BiocBaseUtils_1.12.0
## [57] Matrix_1.7-6 igraph_2.3.4
## [59] tidyselect_1.2.1 abind_1.4-8
## [61] yaml_2.3.12 codetools_0.2-20
## [63] curl_8.0.0 metafor_5.0-1
## [65] lattice_0.23-1 tibble_3.3.1
## [67] plyr_1.8.9 withr_3.0.3
## [69] Biobase_2.70.0 S7_0.2.2
## [71] evaluate_1.0.5 archive_1.1.13
## [73] desc_1.4.3 RcppParallel_5.1.11-2
## [75] BiocFileCache_3.0.0 Biostrings_2.78.0
## [77] filelock_1.0.3 pillar_1.11.1
## [79] MatrixGenerics_1.22.0 metadat_1.6-0
## [81] checkmate_2.3.3 stats4_4.5.3
## [83] generics_0.1.4 vroom_1.7.1
## [85] mathjaxr_2.0-0 S4Vectors_0.48.1
## [87] hms_1.1.4 ggplot2_4.0.3
## [89] scales_1.4.0 glue_1.8.1
## [91] tools_4.5.3 robustbase_0.99-7
## [93] fs_2.1.0 cowplot_1.2.0
## [95] grid_4.5.3 PlinkMatrix_1.0.0
## [97] optparse_1.8.2 tidyr_1.3.2
## [99] nlme_3.1-170 cli_3.6.6
## [101] zigg_0.0.2 textshaping_1.0.5
## [103] mixsqp_0.3-54 S4Arrays_1.10.1
## [105] dplyr_1.2.1 DEoptimR_1.2-1
## [107] gtable_0.3.6 sass_0.4.10
## [109] digest_0.6.39 BiocGenerics_0.56.0
## [111] SparseArray_1.10.10 htmlwidgets_1.6.4
## [113] farver_2.1.2 memoise_2.0.1
## [115] htmltools_0.5.9 pkgdown_2.2.1
## [117] lifecycle_1.0.5 bit64_4.8.6
## [119] MASS_7.3-66