diff --git a/.github/recipe/recipe.yaml b/.github/recipe/recipe.yaml index dd574a80..21666705 100644 --- a/.github/recipe/recipe.yaml +++ b/.github/recipe/recipe.yaml @@ -84,6 +84,7 @@ requirements: - r-vctrs - r-vroom - r-xgboost + - r-yaml run: - bioconductor-biocparallel - bioconductor-genomicranges @@ -145,6 +146,7 @@ requirements: - r-vctrs - r-vroom - r-xgboost + - r-yaml tests: - script: diff --git a/DESCRIPTION b/DESCRIPTION index 55e8b870..4fda7654 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -39,7 +39,8 @@ Imports: tictoc, tidyr, vctrs, - vroom + vroom, + yaml Suggests: aSPU, arrow, @@ -124,19 +125,20 @@ Collate: 'jointEngine.R' 'jointSpecification.R' 'ld.R' + 'sumstatsQc.R' + 'variantId.R' + 'qtlSumStats.R' + 'manifestLoaders.R' 'mashPipeline.R' 'mashWrapper.R' 'pvalCombine.R' 'qtlEnrichmentPipeline.R' - 'qtlSumStats.R' 'regularizedRegressionWrappers.R' 'relatednessQc.R' 'sldscPostprocessingPipeline.R' 'sldscWrapper.R' - 'sumstatsQc.R' 'twasWeights.R' 'twasWeightsPipeline.R' - 'variantId.R' 'vcfWriter.R' Depends: R (>= 3.5) diff --git a/NAMESPACE b/NAMESPACE index 995c5163..93537aa6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -82,7 +82,6 @@ export(fitMashContrast) export(fitMvsusie) export(fitMvsusieRss) export(fitSusieInfThenSusieRss) -export(fixupExampleGenotypePaths) export(formatFinemappingOutput) export(fsusieGetCs) export(fsusieWeights) @@ -185,11 +184,15 @@ export(ldMismatchQc) export(ldPruneByCorrelation) export(learnTwasWeights) export(loadGenotypeRegion) +export(loadGwasSumStatsFromManifest) export(loadLdMatrix) export(loadLdSketch) +export(loadMultiStudyQtlDatasetFromManifest) export(loadMultitaskRegionalData) export(loadMultitraitRSumstat) export(loadMultitraitTensorqtlSumstat) +export(loadQtlDatasetFromManifest) +export(loadQtlSumStatsFromManifest) export(loadRegionalAssociationData) export(loadRegionalFunctionalData) export(loadRegionalMultivariateData) diff --git a/R/GenotypeHandle.R b/R/GenotypeHandle.R index 7e882e16..8ec96674 100644 --- a/R/GenotypeHandle.R +++ b/R/GenotypeHandle.R @@ -424,6 +424,34 @@ setMethod("getFormat", "GenotypeHandle", function(x) x@format) #' @export setMethod("getPath", "GenotypeHandle", function(x) x@path) +# Resolve a portable bundled-resource reference of the form +# "pecotmr://extdata/" to a concrete filesystem path via +# `system.file()`; any ordinary path is returned unchanged. This lets the +# packaged example objects (which cannot bake an install-specific absolute +# path into their serialized GenotypeHandle) carry a portable reference that +# resolves at extraction time, on whatever machine the package is installed. +.resolveGenotypeResourcePath <- function(path) { + if (length(path) != 1L || is.na(path)) return(path) + m <- regmatches(path, regexec("^pecotmr://extdata/(.+)$", path))[[1L]] + if (length(m) != 2L) return(path) + stem <- m[[2L]] + bed <- system.file("extdata", paste0(stem, ".bed"), package = "pecotmr") + if (nzchar(bed)) return(sub("\\.bed$", "", bed)) + resolved <- system.file("extdata", stem, package = "pecotmr") + if (!nzchar(resolved)) + stop("cannot resolve bundled genotype resource '", stem, + "' under inst/extdata/ (package 'pecotmr').") + resolved +} + +# The concrete filesystem path/stem to open for a handle's data. Resolves a +# bundled-resource reference; ordinary paths pass through. Used by the block +# extractors and LD-panel keying, so `getPath()` keeps returning the stored +# (portable) value for display/provenance while file access resolves it. +.genotypeReadPath <- function(handle) { + .resolveGenotypeResourcePath(getPath(handle)) +} + #' @rdname getSampleIds #' @export setMethod("getSampleIds", "GenotypeHandle", function(x) x@sampleIds) diff --git a/R/ctwasPipeline.R b/R/ctwasPipeline.R index 2dab23ba..66f28270 100644 --- a/R/ctwasPipeline.R +++ b/R/ctwasPipeline.R @@ -1242,7 +1242,7 @@ mergeCtwasBoundaryRegions <- function(finemapResult, # @noRd .ctwasLdPanelKey <- function(handle) { fmt <- getFormat(handle) - stem <- getPath(handle) + stem <- .genotypeReadPath(handle) candidates <- switch(fmt, "plink2" = c(paste0(stem, ".pgen")), "plink1" = c(paste0(stem, ".bed")), diff --git a/R/exampleData.R b/R/exampleData.R index 6e38dc1e..bd841c9f 100644 --- a/R/exampleData.R +++ b/R/exampleData.R @@ -288,59 +288,6 @@ NULL NULL -#' @title Resolve bundled-example GenotypeHandle paths -#' -#' @description The bundled S4 example objects (\code{qtl_dataset_example}, -#' \code{qtl_sumstats_example}, \code{qtl_sumstats_multicontext_example}, -#' \code{gwas_sumstats_s4_example}, \code{multi_study_qtl_dataset_example}) -#' store a relative-style path to the bundled -#' \code{inst/extdata/toy_canonical} PLINK1 reference. -#' Call this helper once at the top of a vignette to re-point each -#' contained \code{GenotypeHandle} at the installed path resolved via -#' \code{system.file()}. -#' -#' @param x A bundled \code{QtlDataset}, \code{MultiStudyQtlDataset}, -#' \code{QtlSumStats}, or \code{GwasSumStats} example object. -#' @return The same object with \code{GenotypeHandle@@path} rewritten -#' to the resolved install path. -#' @export -fixupExampleGenotypePaths <- function(x) { - resolve <- function(handle) { - bn <- basename(handle@path) - # PLINK1 / PLINK2 are stems (no extension); resolve via the .bed - # (or .pgen) sidecar and strip the extension back off. - bedPath <- system.file("extdata", paste0(bn, ".bed"), - package = "pecotmr") - if (nzchar(bedPath)) { - handle@path <- sub("\\.bed$", "", bedPath) - return(handle) - } - new <- system.file("extdata", bn, package = "pecotmr") - if (!nzchar(new)) - stop("fixupExampleGenotypePaths: cannot resolve '", bn, - "' under inst/extdata/") - handle@path <- new - handle - } - if (methods::is(x, "QtlDataset")) { - x@genotypes <- resolve(x@genotypes) - } else if (methods::is(x, "MultiStudyQtlDataset")) { - for (nm in names(x@qtlDatasets)) - x@qtlDatasets[[nm]]@genotypes <- - resolve(x@qtlDatasets[[nm]]@genotypes) - if (!is.null(x@sumStats)) - x@sumStats@ldSketch <- resolve(x@sumStats@ldSketch) - } else if (methods::is(x, "QtlSumStats") || - methods::is(x, "GwasSumStats")) { - x@ldSketch <- resolve(x@ldSketch) - } else { - stop("fixupExampleGenotypePaths: unsupported class '", - class(x)[[1L]], "'.") - } - x -} - - #' @name multi_study_qtl_dataset_example #' #' @title Example MultiStudyQtlDataset (S4) diff --git a/R/genotypeIo.R b/R/genotypeIo.R index 8ca39f63..fdb8995c 100644 --- a/R/genotypeIo.R +++ b/R/genotypeIo.R @@ -323,7 +323,7 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { #' @keywords internal .extractBlockGds <- function(handle, snpIdx) { - .withGds(getPath(handle), allow.fork = TRUE, fn = function(gds) { + .withGds(.genotypeReadPath(handle), allow.fork = TRUE, fn = function(gds) { snpIds <- getSnpInfo(handle)$SNP[snpIdx] # Use snpgdsGetGeno for proper non-contiguous SNP selection geno <- SNPRelate::snpgdsGetGeno(gds, snp.id = snpIds, @@ -345,7 +345,7 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { ) param <- VariantAnnotation::ScanVcfParam(which = gr, fixed = NA, info = NA, geno = "GT") - vcf <- VariantAnnotation::readVcf(getPath(handle), genome = "", param = param) + vcf <- VariantAnnotation::readVcf(.genotypeReadPath(handle), genome = "", param = param) gt <- VariantAnnotation::geno(vcf)$GT # Convert GT strings to ALT dosage (A1 dosage) @@ -365,7 +365,7 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { #' @keywords internal .extractBlockPlink1 <- function(handle, snpIdx) { snpIds <- getSnpInfo(handle)$SNP[snpIdx] - pathStem <- getPath(handle) + pathStem <- .genotypeReadPath(handle) plinkData <- snpStats::read.plink( bed = paste0(pathStem, ".bed"), bim = paste0(pathStem, ".bim"), @@ -386,7 +386,7 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { # become stale), so we re-open from getPath() on the fly if the cached # pointer errors out. Opening is cheap relative to dosage extraction. ptr <- getPgenPtr(handle) - paths <- resolvePlink2Paths(getPath(handle)) + paths <- resolvePlink2Paths(.genotypeReadPath(handle)) # A sharded handle routes through a transient view with pgenPtr = NULL (one # pgen per chromosome), and a deserialized pointer is stale; open a fresh # pgen up front in those cases rather than provoking a caught read error. @@ -461,7 +461,7 @@ computeBlockLdCor <- function(handle, snpIdx, backend = "internal", #' @keywords internal .computeBlockLdGds <- function(handle, snpIdx) { - .withGds(getPath(handle), allow.fork = TRUE, fn = function(gds) { + .withGds(.genotypeReadPath(handle), allow.fork = TRUE, fn = function(gds) { snpIds <- getSnpInfo(handle)$SNP[snpIdx] ldMat <- SNPRelate::snpgdsLDMat( gds, snp.id = snpIds, method = "corr", @@ -1052,7 +1052,7 @@ loadGenotypeRegion <- function(genotype, region = NULL, keepIndel = TRUE, # --- Attach allele frequency from .afreq sidecar (plink2 only) --- if (getFormat(handle) == "plink2") { - afreq <- readAfreq(getPath(handle)) + afreq <- readAfreq(.genotypeReadPath(handle)) if (!is.null(afreq)) { afreqCols <- intersect(c("id", "alt_freq", "obs_ct"), colnames(afreq)) variantInfo <- merge(variantInfo, afreq[, afreqCols, drop = FALSE], diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R new file mode 100644 index 00000000..a6bcca5f --- /dev/null +++ b/R/manifestLoaders.R @@ -0,0 +1,966 @@ +# ============================================================================= +# Manifest loaders +# ----------------------------------------------------------------------------- +# Reader functions that build the individual-level / summary-statistic S4 +# containers from an on-disk (or in-memory) manifest: +# +# loadQtlDatasetFromManifest -> QtlDataset +# loadQtlSumStatsFromManifest -> QtlSumStats +# loadGwasSumStatsFromManifest -> GwasSumStats +# loadMultiStudyQtlDatasetFromManifest -> MultiStudyQtlDataset +# +# These are pure readers: they never run summaryStatsQc() (the loaded +# sumstats objects carry qcInfo = list()) and never filter genotypes (the +# QtlDataset QC knobs are stored as lazy slots, applied at extraction time). +# +# A manifest is either a data.frame or a path to a delimited file. A ".csv" +# extension is read with readr::read_csv; anything else is read with +# readr::read_tsv. Both snake_case and camelCase column names are accepted at +# the parsing boundary and canonicalised to camelCase; snake_case is never +# propagated internally. +# ============================================================================= + +#' @include GenotypeHandle.R QtlDataset.R gwasSumStats.R qtlSumStats.R +#' @include MultiStudyQtlDataset.R variantId.R sumstatsQc.R +NULL + +# ============================================================================= +# Manifest ingestion + column canonicalisation +# ============================================================================= + +# Read a manifest into a base data.frame. A data.frame is passed through; a +# single path is read by extension (.csv -> read_csv, else read_tsv). +.readManifest <- function(manifest) { + if (is.data.frame(manifest)) { + return(as.data.frame(manifest, stringsAsFactors = FALSE, check.names = FALSE)) + } + if (!is.character(manifest) || length(manifest) != 1L) { + stop("`manifest` must be a data.frame or a single file path.") + } + if (!file.exists(manifest)) { + stop("manifest file not found: ", manifest) + } + parsed <- if (grepl("\\.csv$", manifest, ignore.case = TRUE)) { + readr::read_csv(manifest, show_col_types = FALSE, progress = FALSE) + } else { + readr::read_tsv(manifest, show_col_types = FALSE, progress = FALSE) + } + as.data.frame(parsed, stringsAsFactors = FALSE, check.names = FALSE) +} + +# Directory a manifest's relative paths resolve against: the file's own +# directory when the manifest was a path, otherwise NULL (paths used as-is). +.manifestBase <- function(manifest) { + if (is.character(manifest) && length(manifest) == 1L && file.exists(manifest)) { + dirname(normalizePath(manifest)) + } else { + NULL + } +} + +# Resolve a possibly-relative path against a manifest base directory. +.resolveRel <- function(p, base) { + if (is.null(p) || length(p) != 1L || is.na(p) || !nzchar(p)) return(p) + if (grepl("^(/|[A-Za-z]:|~)", p) || is.null(base) || file.exists(p)) return(p) + file.path(base, p) +} + +# Rename accepted aliases to the canonical camelCase name and validate that all +# `required` canonical columns are present. `aliases` is a named list keyed by +# canonical name; each value is the character vector of accepted source names +# (including the canonical name itself). The first alias present wins. +.canonManifestCols <- function(df, aliases, required, label) { + for (canon in names(aliases)) { + if (canon %in% names(df)) next + hit <- intersect(aliases[[canon]], names(df)) + if (length(hit) >= 1L) { + names(df)[match(hit[[1L]], names(df))] <- canon + } + } + missingCols <- setdiff(required, names(df)) + if (length(missingCols) > 0L) { + stop(label, " manifest is missing required column(s): ", + paste(missingCols, collapse = ", ")) + } + df +} + +# Resolve a single collection-level scalar (e.g. genome) that may be supplied +# both as a constant manifest column and as a function argument. The column (if +# present) must be constant; a supplied argument must agree with it. +.reconcileScalar <- function(colValues, arg, name, required = TRUE) { + colVal <- NULL + if (!is.null(colValues)) { + v <- unique(as.character(colValues[!is.na(colValues) & + nzchar(as.character(colValues))])) + if (length(v) > 1L) { + stop("`", name, "` column must be constant across the manifest; got: ", + paste(v, collapse = ", ")) + } + if (length(v) == 1L) colVal <- v + } + if (!is.null(arg) && !is.null(colVal)) { + if (!identical(as.character(arg), colVal)) { + stop("`", name, "` argument ('", arg, "') disagrees with the manifest ", + "column value ('", colVal, "').") + } + return(arg) + } + if (!is.null(arg)) return(arg) + if (!is.null(colVal)) return(colVal) + if (required) { + stop("`", name, "` must be provided as an argument or a manifest column.") + } + NULL +} + +# ============================================================================= +# Genotype-format detection + LD-sketch resolution +# ============================================================================= + +# Build a single-file / prefix GenotypeHandle, detecting the backend from the +# path. Recognised: .vcf[.gz|.bgz]/.bcf/.gds single files, .bed/.pgen single +# files, and bare PLINK1/PLINK2 prefixes (probed by sidecar existence). +.detectGenotypeFormat <- function(path) { + lower <- tolower(path) + if (grepl("\\.vcf(\\.b?gz)?$", lower) || endsWith(lower, ".bcf") || + endsWith(lower, ".gds")) { + return(GenotypeHandle(path = path)) + } + if (endsWith(lower, ".bed")) { + return(GenotypeHandle(plink1Prefix = sub("\\.bed$", "", path, ignore.case = TRUE))) + } + if (endsWith(lower, ".pgen")) { + return(GenotypeHandle(plink2Prefix = sub("\\.pgen$", "", path, ignore.case = TRUE))) + } + if (file.exists(paste0(path, ".bed"))) return(GenotypeHandle(plink1Prefix = path)) + if (file.exists(paste0(path, ".pgen"))) return(GenotypeHandle(plink2Prefix = path)) + stop("Could not determine genotype format for '", path, + "' (expected .bed/.pgen/.vcf[.gz]/.bcf/.gds or a PLINK prefix).") +} + +# Turn an ldSketch specification into a GenotypeHandle. Accepts a prebuilt +# handle, a named chrom -> path vector (genoMeta), or a single string that is +# either a genotype file/prefix or a chrom-sharded genoMeta meta file. +.resolveLdSketch <- function(spec) { + if (is.null(spec)) return(NULL) + if (methods::is(spec, "GenotypeHandle")) return(spec) + if (is.character(spec) && !is.null(names(spec))) { + return(GenotypeHandle(genoMeta = spec)) + } + if (is.character(spec) && length(spec) == 1L) { + lower <- tolower(spec) + looksLikeGenotype <- grepl("\\.vcf(\\.b?gz)?$", lower) || + endsWith(lower, ".bcf") || endsWith(lower, ".gds") || + endsWith(lower, ".bed") || endsWith(lower, ".pgen") || + file.exists(paste0(spec, ".bed")) || file.exists(paste0(spec, ".pgen")) + if (looksLikeGenotype) return(.detectGenotypeFormat(spec)) + # Otherwise treat it as a chrom-sharded genoMeta meta file. + return(GenotypeHandle(genoMeta = spec)) + } + stop("`ldSketch` must be a GenotypeHandle, a genotype path/prefix, or a ", + "genoMeta spec (named chrom->path vector or meta-file path).") +} + +# Resolve the LD sketch from the (argument, ldSketchPath column) pair. +.resolveLdSketchInput <- function(df, ldSketch, base) { + if (methods::is(ldSketch, "GenotypeHandle")) return(ldSketch) + colSpec <- NULL + if (!is.null(df$ldSketchPath)) { + v <- unique(as.character(df$ldSketchPath[!is.na(df$ldSketchPath) & + nzchar(as.character(df$ldSketchPath))])) + if (length(v) > 1L) stop("`ldSketchPath` column must be constant across the manifest.") + if (length(v) == 1L) colSpec <- .resolveRel(v, base) + } + spec <- ldSketch + if (!is.null(spec) && !is.null(colSpec) && + !identical(as.character(spec), as.character(colSpec))) { + stop("`ldSketch` argument disagrees with the manifest `ldSketchPath` column.") + } + spec <- spec %||% colSpec + if (is.null(spec)) { + stop("`ldSketch` must be provided as an argument or an `ldSketchPath` column.") + } + .resolveLdSketch(spec) +} + +# ============================================================================= +# LD-sketch containment check (validation only; no filtering) +# ============================================================================= + +# Canonical chromosomes present in an LD sketch: for a sharded handle these are +# the chromPaths names; for a single-file handle the unique snpInfo chromosomes. +.ldSketchChroms <- function(ldSketch) { + if (length(ldSketch@chromPaths) > 0L) { + canonChrom(names(ldSketch@chromPaths)) + } else { + unique(canonChrom(as.character(ldSketch@snpInfo$CHR))) + } +} + +# Two-tier containment check for one entry GRanges against the LD sketch: +# * chromosome tier (hard error): any entry chromosome absent from the sketch. +# * variant tier: match entry variants against the sketch's variants on the +# shared chromosomes. Zero overlap -> error; overlap fraction below +# `minLdOverlapWarn` -> warning. Never filters. +.checkLdContainment <- function(ldSketch, entryGr, minLdOverlapWarn, label) { + if (length(entryGr) == 0L) return(invisible(NULL)) + reqChrom <- unique(canonChrom(as.character(GenomicRanges::seqnames(entryGr)))) + availChrom <- .ldSketchChroms(ldSketch) + missingChrom <- setdiff(reqChrom, availChrom) + if (length(missingChrom) > 0L) { + stop(label, ": chromosome(s) ", paste(missingChrom, collapse = ", "), + " are absent from the LD sketch (available: ", + paste(availChrom, collapse = ", "), ").") + } + snpInfo <- ldSketch@snpInfo + if (is.null(snpInfo) || nrow(snpInfo) == 0L) { + warning(label, ": LD sketch carries no variant metadata; ", + "skipping the variant-overlap check.") + return(invisible(NULL)) + } + keepSketch <- canonChrom(as.character(snpInfo$CHR)) %in% reqChrom + snpInfo <- snpInfo[keepSketch, , drop = FALSE] + entryVid <- formatVariantId( + chrom = canonChrom(as.character(GenomicRanges::seqnames(entryGr))), + pos = GenomicRanges::start(entryGr), + A2 = as.character(entryGr$A2), + A1 = as.character(entryGr$A1)) + sketchVid <- formatVariantId( + chrom = canonChrom(as.character(snpInfo$CHR)), + pos = as.integer(snpInfo$BP), + A2 = as.character(snpInfo$A2), + A1 = as.character(snpInfo$A1)) + km <- matchVariants(entryVid, sketchVid) + nOverlap <- length(km$idxA) + if (nOverlap == 0L) { + stop(label, ": none of the ", length(entryGr), " summary-statistic ", + "variants are present in the LD sketch.") + } + frac <- nOverlap / length(entryGr) + if (frac < minLdOverlapWarn) { + warning(sprintf( + "%s: only %.1f%% of summary-statistic variants (%d/%d) are present in the LD sketch.", + label, 100 * frac, nOverlap, length(entryGr))) + } + invisible(NULL) +} + +# ============================================================================= +# Summary-statistic file readers +# ============================================================================= + +# Coerce a region argument (chr:start-end string, GRanges, or one-row +# data.frame with chrom/start/end) into a GRanges. +.asGRegion <- function(region) { + if (methods::is(region, "GRanges")) return(region) + if (is.data.frame(region)) { + return(GenomicRanges::GRanges( + seqnames = as.character(region$chrom), + ranges = IRanges::IRanges(start = as.integer(region$start), + end = as.integer(region$end)))) + } + if (is.character(region) && length(region) == 1L) { + m <- regmatches(region, regexec("^([^:]+):([0-9,]+)-([0-9,]+)$", region))[[1L]] + if (length(m) != 4L) { + stop("`region` must be a 'chr:start-end' string, a GRanges, or a ", + "one-row data.frame with chrom/start/end.") + } + return(GenomicRanges::GRanges( + seqnames = m[[2L]], + ranges = IRanges::IRanges(start = as.integer(gsub(",", "", m[[3L]])), + end = as.integer(gsub(",", "", m[[4L]]))))) + } + stop("`region` must be a 'chr:start-end' string, a GRanges, or a ", + "one-row data.frame with chrom/start/end.") +} + +# Remap a region's seqnames to the seqnames actually present in a tabix file +# (tolerant of chr-prefix differences). Returns NULL when no seqname matches. +.matchRegionSeqnames <- function(gr, available) { + want <- as.character(GenomicRanges::seqnames(gr)) + mapped <- vapply(want, function(w) { + hit <- available[canonChrom(available) == canonChrom(w)] + if (length(hit) >= 1L) hit[[1L]] else NA_character_ + }, character(1)) + if (anyNA(mapped)) return(NULL) + GenomicRanges::GRanges( + seqnames = mapped, + ranges = IRanges::IRanges(start = GenomicRanges::start(gr), + end = GenomicRanges::end(gr))) +} + +# Read the region of a bgzipped + tabix-indexed delimited file into a +# data.frame, using the file's last header (#...) line for column names. +.readTabixRegion <- function(path, region) { + tf <- Rsamtools::TabixFile(path) + hdr <- Rsamtools::headerTabix(tf)$header + colLine <- if (length(hdr) > 0L) { + sub("^#", "", hdr[[length(hdr)]]) + } else { + NULL + } + gr <- .matchRegionSeqnames(.asGRegion(region), Rsamtools::seqnamesTabix(tf)) + emptyDf <- if (!is.null(colLine)) { + as.data.frame(readr::read_tsv(I(colLine), show_col_types = FALSE, + progress = FALSE), check.names = FALSE)[0, , drop = FALSE] + } else { + data.frame() + } + if (is.null(gr)) return(emptyDf) + lines <- unlist(Rsamtools::scanTabix(tf, param = gr), use.names = FALSE) + if (length(lines) == 0L) return(emptyDf) + txt <- paste(c(colLine, lines), collapse = "\n") + as.data.frame(readr::read_tsv(I(txt), show_col_types = FALSE, progress = FALSE), + check.names = FALSE) +} + +# Default source-column aliases used when no explicit columnMapping is given +# (mirrors gwas_sumstats_construct.R). +.sumstatColumnAliases <- list( + chrom = c("chrom", "#chrom", "chr", "#chr"), + pos = c("pos", "position", "BP", "bp"), + variant_id = c("variant_id", "SNP", "snp", "rsid"), + A1 = c("A1", "a1"), + A2 = c("A2", "a2"), + Z = c("z", "Z"), + N = c("n_sample", "N", "n"), + BETA = c("beta", "BETA"), + SE = c("se", "SE"), + P = c("p", "P", "pvalue", "pval"), + MAF = c("maf", "MAF", "effect_allele_frequency"), + INFO = c("info", "INFO")) + +# Accepted columnMapping *key* spellings per canonical field. A mapping file +# uses the xqtl-protocol standard-key names (e.g. `z`, `n_sample`), which +# differ from the canonical mcol names (`Z`, `N`); this bridges the two. +.sumstatMappingKeys <- list( + chrom = c("chrom", "chr"), + pos = c("pos", "position"), + variant_id = c("variant_id", "snp", "SNP"), + A1 = c("A1", "a1"), + A2 = c("A2", "a2"), + Z = c("z", "Z"), + N = c("n_sample", "N", "n"), + BETA = c("beta", "BETA"), + SE = c("se", "SE"), + P = c("p", "P", "pvalue"), + MAF = c("maf", "MAF", "effect_allele_frequency"), + INFO = c("info", "INFO")) + +# Read a columnMapping specification into a named character vector mapping a +# standard key (chrom/pos/variant_id/...) to the source column name. Accepts a +# named list/vector, or a path to a YAML file of `standardName: sourceName` +# entries (the xqtl-protocol column-mapping format). +.readColumnMapping <- function(columnMapping) { + if (is.null(columnMapping)) return(NULL) + if (is.list(columnMapping) || (is.character(columnMapping) && + !is.null(names(columnMapping)))) { + return(vapply(columnMapping, as.character, character(1))) + } + if (is.character(columnMapping) && length(columnMapping) == 1L) { + if (!file.exists(columnMapping)) { + stop("columnMapping file not found: ", columnMapping) + } + mapping <- yaml::read_yaml(columnMapping) + if (!is.list(mapping) || is.null(names(mapping)) || + any(!nzchar(names(mapping)))) { + stop("columnMapping file '", columnMapping, + "' must be a YAML mapping of standardName: sourceName.") + } + vapply(mapping, as.character, character(1)) + } else { + stop("`columnMapping` must be a named list/vector or a path to a YAML ", + "mapping file.") + } +} + +# Standardise the columns of a raw sumstats data.frame into the canonical +# schema expected by .dfToEntryGranges: chrom, pos, SNP, A1, A2, Z, N (+ +# optional BETA, SE, P, MAF, INFO). `mapping` (named std -> source) overrides +# the default aliases per key. +.resolveSumstatCols <- function(df, columnMapping, label) { + # Strip a leading '#' from the first column name so a '#CHR'-style header + # (common in GWAS TSVs) resolves the same way on the plain-text and tabix + # read paths, and so explicit columnMapping values match the bare name. + if (ncol(df) > 0L) names(df)[1L] <- sub("^#", "", names(df)[1L]) + mapping <- .readColumnMapping(columnMapping) + resolveKey <- function(key) { + # Look the field up in the mapping under any accepted standard-key spelling + # (`z`/`Z`, `n_sample`/`N`, ...). `intersect` avoids the "subscript out of + # bounds" a named character vector throws for an absent `[[` key. + if (!is.null(mapping)) { + cand <- intersect(.sumstatMappingKeys[[key]], names(mapping)) + if (length(cand) >= 1L) { + src <- mapping[[cand[[1L]]]] + if (!(src %in% names(df))) { + stop(label, ": columnMapping['", cand[[1L]], "'] = '", src, + "' is not a column in the sumstats file.") + } + return(src) + } + } + intersect(.sumstatColumnAliases[[key]], names(df))[1L] + } + required <- c("chrom", "pos", "variant_id", "A1", "A2", "Z", "N") + resolved <- setNames(lapply(required, resolveKey), required) + missingKeys <- required[vapply(resolved, function(x) is.na(x) || is.null(x), logical(1))] + if (length(missingKeys) > 0L) { + stop(label, ": sumstats file is missing required field(s): ", + paste(missingKeys, collapse = ", "), + " (supply a columnMapping if the source names differ).") + } + out <- data.frame( + chrom = as.character(df[[resolved$chrom]]), + pos = as.integer(df[[resolved$pos]]), + SNP = as.character(df[[resolved$variant_id]]), + A1 = as.character(df[[resolved$A1]]), + A2 = as.character(df[[resolved$A2]]), + Z = as.numeric(df[[resolved$Z]]), + N = as.numeric(df[[resolved$N]]), + stringsAsFactors = FALSE) + for (key in c("BETA", "SE", "P", "MAF", "INFO")) { + src <- resolveKey(key) + if (!is.na(src) && !is.null(src)) out[[key]] <- as.numeric(df[[src]]) + } + out +} + +# Default GWAS-VCF FORMAT tag mapping (canonical stat -> FORMAT field). +.gwasVcfFormatDefaults <- list(BETA = "ES", SE = "SE", P = "LP", + N = "SS", MAF = "EAF") + +# Pick the FORMAT sample column of a GWAS-VCF to read. Defaults to the single +# sample when there is exactly one; otherwise `sampleSelect` must name it. +.pickVcfSample <- function(samples, sampleSelect, label) { + if (!is.null(sampleSelect)) { + if (!(sampleSelect %in% samples)) { + stop(label, ": sampleSelect '", sampleSelect, + "' is not a sample in the VCF (have: ", + paste(samples, collapse = ", "), ").") + } + return(sampleSelect) + } + if (length(samples) == 1L) return(samples[[1L]]) + stop(label, ": the VCF has ", length(samples), " samples (", + paste(samples, collapse = ", "), + "); pass `sampleSelect` to choose the study column.") +} + +# Convert a VCF (GWAS-VCF format) into the canonical sumstats data.frame. The +# effect allele (A1) is ALT; the other allele (A2) is REF. Stats come from the +# per-study FORMAT fields ES/SE/LP/SS/EAF, with Z = ES / SE. +.vcfToSumstatDf <- function(vcf, sampleSelect, formatMapping, label) { + fmap <- modifyList(.gwasVcfFormatDefaults, + if (is.null(formatMapping)) list() else as.list(formatMapping)) + rr <- SummarizedExperiment::rowRanges(vcf) + altList <- VariantAnnotation::alt(vcf) + alt <- vapply(as.list(altList), function(a) as.character(a)[[1L]], character(1)) + ref <- as.character(VariantAnnotation::ref(vcf)) + geno <- VariantAnnotation::geno(vcf) + sample <- .pickVcfSample(colnames(vcf), sampleSelect, label) + getField <- function(tag) { + if (!is.null(tag) && tag %in% names(geno)) { + as.numeric(geno[[tag]][, sample]) + } else { + NULL + } + } + es <- getField(fmap$BETA) + se <- getField(fmap$SE) + lp <- getField(fmap$P) + ss <- getField(fmap$N) + eaf <- getField(fmap$MAF) + if (is.null(es) || is.null(se)) { + stop(label, ": GWAS-VCF must carry ES and SE FORMAT fields to derive Z ", + "(configure via `formatMapping`).") + } + if (is.null(ss)) { + stop(label, ": GWAS-VCF must carry an SS (sample size) FORMAT field.") + } + out <- data.frame( + chrom = as.character(GenomicRanges::seqnames(rr)), + pos = as.integer(GenomicRanges::start(rr)), + SNP = names(rr), + A1 = alt, + A2 = ref, + Z = es / se, + N = ss, + BETA = es, + SE = se, + stringsAsFactors = FALSE) + if (!is.null(lp)) out$P <- 10^(-lp) + if (!is.null(eaf)) out$MAF <- pmin(eaf, 1 - eaf) + out +} + +# Read one GWAS-VCF sumstats file (region-restricted only when bgzipped + +# tabix-indexed) into the canonical sumstats data.frame. +.readSumStatsVcf <- function(path, region, sampleSelect, formatMapping, label) { + if (!requireNamespace("VariantAnnotation", quietly = TRUE)) { + stop(label, ": reading VCF sumstats requires the 'VariantAnnotation' ", + "package; please install it.") + } + hasTbi <- file.exists(paste0(path, ".tbi")) + if (!is.null(region) && hasTbi) { + tf <- Rsamtools::TabixFile(path) + gr <- .matchRegionSeqnames(.asGRegion(region), Rsamtools::seqnamesTabix(tf)) + vcf <- if (is.null(gr)) { + VariantAnnotation::readVcf(path)[0, ] + } else { + VariantAnnotation::readVcf(tf, param = gr) + } + } else { + if (!is.null(region) && !hasTbi) { + warning(label, ": VCF '", path, "' is not bgzipped + tabix-indexed; ", + "region ignored, reading the whole file.") + } + vcf <- VariantAnnotation::readVcf(path) + } + .vcfToSumstatDf(vcf, sampleSelect, formatMapping, label) +} + +# Read one delimited-text sumstats file. Region reads only when a .tbi +# sidecar exists; otherwise the whole file is read and any region is ignored. +.readSumStatsText <- function(path, region, columnMapping, label) { + hasTbi <- file.exists(paste0(path, ".tbi")) + raw <- if (!is.null(region) && hasTbi) { + .readTabixRegion(path, region) + } else { + if (!is.null(region) && !hasTbi) { + warning(label, ": '", path, "' is not tabix-indexed; ", + "region ignored, reading the whole file.") + } + as.data.frame(readr::read_tsv(path, show_col_types = FALSE, progress = FALSE), + check.names = FALSE) + } + .resolveSumstatCols(raw, columnMapping, label) +} + +# Dispatch a sumstats file to the text or GWAS-VCF reader. BCF is rejected. +.readSumStatsFile <- function(path, region, columnMapping, sampleSelect, + formatMapping, label) { + lower <- tolower(path) + if (endsWith(lower, ".bcf")) { + stop(label, ": BCF sumstats are not supported; convert '", path, + "' to a bgzipped, tabix-indexed VCF.") + } + if (grepl("\\.vcf(\\.b?gz)?$", lower)) { + .readSumStatsVcf(path, region, sampleSelect, formatMapping, label) + } else { + .readSumStatsText(path, region, columnMapping, label) + } +} + +# Read one sumstats source into an entry GRanges (canonical mcols). +.loadSumStatsEntry <- function(path, region, columnMapping, sampleSelect, + formatMapping, label) { + df <- .readSumStatsFile(path, region, columnMapping, sampleSelect, + formatMapping, label) + .dfToEntryGranges(df) +} + +# ============================================================================= +# Phenotype SummarizedExperiment builder (individual-level QTL) +# ============================================================================= + +# Read a covariate TSV into a samples-as-rows numeric matrix. The first column +# holds row identifiers; `transpose = TRUE` handles QTLtools-format inputs +# (covariates as rows, samples as columns). +.readCovariateMatrix <- function(path, transpose = FALSE) { + raw <- as.data.frame(readr::read_tsv(path, show_col_types = FALSE, + progress = FALSE), check.names = FALSE) + rn <- as.character(raw[[1L]]) + m <- as.matrix(raw[, -1L, drop = FALSE]) + rownames(m) <- rn + storage.mode(m) <- "double" + if (transpose) m <- t(m) + m +} + +# Build one per-context SummarizedExperiment from a bgzipped BED phenotype +# file (QTLtools layout: #chr/start/end/ID + sample columns) and an optional +# per-context covariate TSV (attached as colData). +.buildContextSe <- function(bedPath, covPath, transposeCov) { + bed <- as.data.frame(readr::read_tsv(bedPath, show_col_types = FALSE, + progress = FALSE, comment = ""), + check.names = FALSE) + chrCol <- intersect(c("#chr", "#chrom", "chrom", "chr"), names(bed))[1L] + startCol <- intersect(c("start", "Start"), names(bed))[1L] + endCol <- intersect(c("end", "End"), names(bed))[1L] + geneCol <- intersect(c("gene_id", "ID", "phenotype_id"), names(bed))[1L] + if (any(is.na(c(chrCol, startCol, endCol, geneCol)))) { + stop("phenotype BED is missing one of chrom/start/end/gene_id columns: ", + bedPath) + } + meta <- bed[, c(chrCol, startCol, endCol, geneCol)] + names(meta) <- c("chrom", "start", "end", "gene_id") + sampleCols <- setdiff(names(bed), + c("#chr", "#chrom", "chrom", "chr", "start", "Start", "end", "End", + "gene_id", "ID", "phenotype_id", "strand", "Strand")) + expr <- as.matrix(bed[, sampleCols, drop = FALSE]) + storage.mode(expr) <- "double" + rownames(expr) <- meta$gene_id + rr <- GenomicRanges::GRanges( + seqnames = meta$chrom, + ranges = IRanges::IRanges(start = meta$start + 1L, end = meta$end)) + names(rr) <- meta$gene_id + cd <- if (!is.null(covPath) && nzchar(covPath)) { + pcov <- .readCovariateMatrix(covPath, transposeCov) + common <- intersect(rownames(pcov), colnames(expr)) + if (length(common) == 0L) { + stop("No shared samples between phenotype and covariate file: ", bedPath) + } + expr <- expr[, common, drop = FALSE] + S4Vectors::DataFrame(pcov[common, , drop = FALSE], row.names = common) + } else { + S4Vectors::DataFrame(row.names = colnames(expr)) + } + SummarizedExperiment::SummarizedExperiment( + assays = list(expression = expr), rowRanges = rr, colData = cd) +} + +# Aliases for the QtlDataset (phenotype) manifest. +.qtlDatasetManifestAliases <- list( + context = c("context", "cond", "condition"), + phenotypePath = c("phenotypePath", "phenotype_path", "path", + "phenotypeFile", "phenotype_file"), + covariatePath = c("covariatePath", "covariate_path", "cov_path", + "covPath"), + study = c("study", "study_id", "studyId"), + genotypePath = c("genotypePath", "genotype_path", "genotypePrefix", + "genotype_prefix", "genotypeFile"), + genotypeCovariatePath = c("genotypeCovariatePath", "genotype_covariate_path", + "genotype_covariates", "genotypeCovariates")) + +# Build one QtlDataset from an already-canonicalised set of manifest rows for a +# single study. `genotypesOverride` / `genoCovOverride` short-circuit the +# genotypePath / genotypeCovariatePath columns when supplied as arguments. +.buildQtlDatasetFromRows <- function(rows, study, base, transposeCov, qc, + genotypesOverride = NULL, + genoCovOverride = NULL) { + contexts <- unique(as.character(rows$context)) + phenotypes <- setNames(lapply(contexts, function(cx) { + sub <- rows[as.character(rows$context) == cx, , drop = FALSE] + pths <- unique(sub$phenotypePath[!is.na(sub$phenotypePath)]) + if (length(pths) != 1L) { + stop("Context '", cx, "' (study '", study, + "') must reference exactly one phenotypePath; got: ", + paste(pths, collapse = ", ")) + } + covPath <- NULL + if ("covariatePath" %in% names(sub)) { + covs <- unique(sub$covariatePath[!is.na(sub$covariatePath) & + nzchar(as.character(sub$covariatePath))]) + if (length(covs) > 1L) { + stop("Context '", cx, "' (study '", study, + "') references multiple covariatePath values: ", + paste(covs, collapse = ", ")) + } + if (length(covs) == 1L) covPath <- .resolveRel(covs[[1L]], base) + } + .buildContextSe(.resolveRel(pths[[1L]], base), covPath, transposeCov) + }), contexts) + + genoHandle <- if (methods::is(genotypesOverride, "GenotypeHandle")) { + genotypesOverride + } else { + spec <- if (is.character(genotypesOverride)) { + genotypesOverride + } else { + v <- unique(as.character(rows$genotypePath[!is.na(rows$genotypePath)])) + if (length(v) != 1L) { + stop("study '", study, "' must reference exactly one genotypePath; got: ", + paste(v, collapse = ", ")) + } + v[[1L]] + } + .detectGenotypeFormat(.resolveRel(spec, base)) + } + + genoCov <- if (is.matrix(genoCovOverride)) { + genoCovOverride + } else { + spec <- if (is.character(genoCovOverride)) { + genoCovOverride + } else if ("genotypeCovariatePath" %in% names(rows)) { + v <- unique(as.character(rows$genotypeCovariatePath[ + !is.na(rows$genotypeCovariatePath) & + nzchar(as.character(rows$genotypeCovariatePath))])) + if (length(v) > 1L) { + stop("study '", study, "' references multiple genotypeCovariatePath values.") + } + if (length(v) == 1L) v[[1L]] else NULL + } else { + NULL + } + if (!is.null(spec)) { + .readCovariateMatrix(.resolveRel(spec, base), transposeCov) + } else { + matrix(numeric(0), nrow = 0, ncol = 0) + } + } + + do.call(QtlDataset, c(list( + study = study, genotypes = genoHandle, phenotypes = phenotypes, + genotypeCovariates = genoCov), qc)) +} + +# ============================================================================= +# Public loaders +# ============================================================================= + +#' @title Load a QtlDataset from a manifest +#' @description Build a single-study \code{\link{QtlDataset}} from a manifest +#' describing one row per QTL context (or per trait x context). The manifest +#' is a data.frame or a path to a delimited file (\code{.csv} -> CSV, else +#' TSV). Column names may be snake_case or camelCase. +#' @param manifest A data.frame or a path to a manifest file. Recognised +#' columns (canonical camelCase; snake_case aliases accepted): +#' \code{context} (required), \code{phenotypePath} (required, bgzipped BED), +#' \code{covariatePath} (optional, per-context covariates), and the +#' single-valued \code{study} / \code{genotypePath} / +#' \code{genotypeCovariatePath}. +#' @param study Study identifier; reconciled with a \code{study} column. +#' @param genotypes A \code{\link{GenotypeHandle}} or a genotype path/prefix; +#' reconciled with a \code{genotypePath} column. +#' @param genotypeCovariates A numeric matrix (samples x covariates) or a path +#' to a covariate TSV; reconciled with a \code{genotypeCovariatePath} column. +#' @param scaleResiduals,mafCutoff,macCutoff,xvarCutoff,imissCutoff,keepSamples,keepVariants,keepIndel +#' Pass-through \code{\link{QtlDataset}} arguments (stored as lazy QC slots). +#' @param transposeCovariates Transpose covariate TSVs (QTLtools layout) before +#' treating them as samples-as-rows. +#' @return A \code{QtlDataset} object. +#' @export +loadQtlDatasetFromManifest <- function(manifest, study = NULL, + genotypes = NULL, + genotypeCovariates = NULL, + scaleResiduals = TRUE, + mafCutoff = 0, macCutoff = 0, + xvarCutoff = 0, imissCutoff = 0, + keepSamples = character(0), + keepVariants = character(0), + keepIndel = TRUE, + transposeCovariates = FALSE) { + base <- .manifestBase(manifest) + df <- .canonManifestCols(.readManifest(manifest), + .qtlDatasetManifestAliases, + required = c("context", "phenotypePath"), + label = "QtlDataset") + study <- .reconcileScalar(df$study, study, "study") + qc <- list(scaleResiduals = scaleResiduals, mafCutoff = mafCutoff, + macCutoff = macCutoff, xvarCutoff = xvarCutoff, + imissCutoff = imissCutoff, keepSamples = keepSamples, + keepVariants = keepVariants, keepIndel = keepIndel) + .buildQtlDatasetFromRows(df, study, base, transposeCovariates, qc, + genotypesOverride = genotypes, + genoCovOverride = genotypeCovariates) +} + +# Aliases shared by the sumstats manifests. +.gwasSumStatsManifestAliases <- list( + study = c("study", "study_id", "studyId"), + sumStatsPath = c("sumStatsPath", "sumstats_path", "path", "file_path", + "gwas_tsv", "filePath"), + columnMapping = c("columnMapping", "column_mapping", "column_mapping_file", + "columnMappingFile"), + nCase = c("nCase", "n_case"), + nControl = c("nControl", "n_control"), + varY = c("varY", "var_y"), + genome = c("genome"), + ldSketchPath = c("ldSketchPath", "ld_sketch_path", "ld_sketch")) + +.qtlSumStatsManifestAliases <- list( + study = c("study", "study_id", "studyId"), + context = c("context", "cond", "condition"), + trait = c("trait", "gene_id", "geneId", "molecular_trait_id"), + sumStatsPath = c("sumStatsPath", "sumstats_path", "path", "file_path", + "filePath"), + columnMapping = c("columnMapping", "column_mapping", "column_mapping_file", + "columnMappingFile"), + varY = c("varY", "var_y"), + genome = c("genome"), + ldSketchPath = c("ldSketchPath", "ld_sketch_path", "ld_sketch")) + +#' @title Load a GwasSumStats collection from a manifest +#' @description Build a \code{\link{GwasSumStats}} from a manifest with one row +#' per study. No QC is run (the result carries \code{qcInfo = list()}). +#' @param manifest A data.frame or path. Columns (snake_case aliases accepted): +#' \code{study} (required, unique), \code{sumStatsPath} (required), +#' \code{columnMapping} (optional), \code{nCase} / \code{nControl} +#' (optional), \code{varY} (optional), and the single-valued \code{genome} / +#' \code{ldSketchPath}. +#' @param genome Genome build; reconciled with a \code{genome} column. +#' @param ldSketch A \code{\link{GenotypeHandle}} or a spec +#' (path/prefix/genoMeta); reconciled with an \code{ldSketchPath} column. +#' @param region Optional \code{chr:start-end} string, GRanges, or one-row +#' data.frame restricting the variants read (honoured only for tabix-indexed +#' text and bgzipped+tabixed VCF; ignored with a warning otherwise). +#' @param minLdOverlapWarn Warn when the fraction of sumstats variants present +#' in the LD sketch falls below this (default 0.5). +#' @param columnMapping Optional default column mapping (a named list/vector, +#' or a path to a YAML file of \code{standardName: sourceName} entries) +#' applied when a row has no \code{columnMapping}. +#' @param sampleSelect Optional GWAS-VCF FORMAT sample (study) column to read. +#' @param formatMapping Optional GWAS-VCF FORMAT tag mapping (canonical stat -> +#' FORMAT field), overriding ES/SE/LP/SS/EAF defaults. +#' @return A \code{GwasSumStats} object. +#' @export +loadGwasSumStatsFromManifest <- function(manifest, genome = NULL, + ldSketch = NULL, region = NULL, + minLdOverlapWarn = 0.5, + columnMapping = NULL, + sampleSelect = NULL, + formatMapping = NULL) { + base <- .manifestBase(manifest) + df <- .canonManifestCols(.readManifest(manifest), + .gwasSumStatsManifestAliases, + required = c("study", "sumStatsPath"), + label = "GwasSumStats") + if (anyDuplicated(as.character(df$study))) { + stop("GwasSumStats manifest `study` values must be unique.") + } + genome <- .reconcileScalar(df$genome, genome, "genome") + ldSketch <- .resolveLdSketchInput(df, ldSketch, base) + + entries <- lapply(seq_len(nrow(df)), function(i) { + label <- paste0("GwasSumStats[study=", df$study[[i]], "]") + mapping <- if ("columnMapping" %in% names(df) && + !is.na(df$columnMapping[[i]]) && + nzchar(as.character(df$columnMapping[[i]]))) { + .resolveRel(as.character(df$columnMapping[[i]]), base) + } else { + columnMapping + } + gr <- .loadSumStatsEntry(.resolveRel(as.character(df$sumStatsPath[[i]]), base), + region, mapping, sampleSelect, formatMapping, label) + .checkLdContainment(ldSketch, gr, minLdOverlapWarn, label) + gr + }) + + args <- list(study = as.character(df$study), entry = entries, + genome = genome, ldSketch = ldSketch) + if ("nCase" %in% names(df)) args$nCase <- as.numeric(df$nCase) + if ("nControl" %in% names(df)) args$nControl <- as.numeric(df$nControl) + if ("varY" %in% names(df)) args$varY <- as.numeric(df$varY) + do.call(GwasSumStats, args) +} + +# Build the entry list + tuple vectors for a QtlSumStats manifest. +.loadQtlSumStatsEntries <- function(df, base, region, ldSketch, + minLdOverlapWarn, columnMapping, + sampleSelect, formatMapping) { + lapply(seq_len(nrow(df)), function(i) { + label <- paste0("QtlSumStats[", df$study[[i]], "/", df$context[[i]], "/", + df$trait[[i]], "]") + mapping <- if ("columnMapping" %in% names(df) && + !is.na(df$columnMapping[[i]]) && + nzchar(as.character(df$columnMapping[[i]]))) { + .resolveRel(as.character(df$columnMapping[[i]]), base) + } else { + columnMapping + } + gr <- .loadSumStatsEntry(.resolveRel(as.character(df$sumStatsPath[[i]]), base), + region, mapping, sampleSelect, formatMapping, label) + if (!is.null(ldSketch)) { + .checkLdContainment(ldSketch, gr, minLdOverlapWarn, label) + } + gr + }) +} + +#' @title Load a QtlSumStats collection from a manifest +#' @description Build a \code{\link{QtlSumStats}} from a manifest with one row +#' per \code{(study, context, trait)} tuple. No QC is run. +#' @param manifest A data.frame or path. Columns (snake_case aliases accepted): +#' \code{study}, \code{context}, \code{trait} (required), \code{sumStatsPath} +#' (required), \code{columnMapping} (optional), \code{varY} (optional), and +#' the single-valued \code{genome} / \code{ldSketchPath}. +#' @param genome Genome build; reconciled with a \code{genome} column. +#' @param ldSketch A \code{\link{GenotypeHandle}} or spec; reconciled with an +#' \code{ldSketchPath} column. +#' @param region,minLdOverlapWarn,columnMapping,sampleSelect,formatMapping As +#' for \code{\link{loadGwasSumStatsFromManifest}}. +#' @return A \code{QtlSumStats} object. +#' @export +loadQtlSumStatsFromManifest <- function(manifest, genome = NULL, + ldSketch = NULL, region = NULL, + minLdOverlapWarn = 0.5, + columnMapping = NULL, + sampleSelect = NULL, + formatMapping = NULL) { + base <- .manifestBase(manifest) + df <- .canonManifestCols(.readManifest(manifest), + .qtlSumStatsManifestAliases, + required = c("study", "context", "trait", + "sumStatsPath"), + label = "QtlSumStats") + genome <- .reconcileScalar(df$genome, genome, "genome") + ldSketch <- .resolveLdSketchInput(df, ldSketch, base) + entries <- .loadQtlSumStatsEntries(df, base, region, ldSketch, + minLdOverlapWarn, columnMapping, + sampleSelect, formatMapping) + args <- list(study = as.character(df$study), + context = as.character(df$context), + trait = as.character(df$trait), + entry = entries, genome = genome, ldSketch = ldSketch) + if ("varY" %in% names(df)) args$varY <- as.numeric(df$varY) + do.call(QtlSumStats, args) +} + +#' @title Load a MultiStudyQtlDataset from manifests +#' @description Build a \code{\link{MultiStudyQtlDataset}} from a QtlDataset +#' manifest (individual-level studies) and, optionally, a QtlSumStats +#' manifest (summary-only studies). +#' @param qtlDatasetsManifest A data.frame or path. The QtlDataset schema with +#' \code{study} and \code{genotypePath} \strong{mandatory} (grouping key = +#' \code{study}); one \code{QtlDataset} is built per study group. +#' @param sumStatsManifest Optional QtlSumStats manifest (see +#' \code{\link{loadQtlSumStatsFromManifest}}). +#' @param genome,ldSketch,region,minLdOverlapWarn,columnMapping,sampleSelect,formatMapping +#' Passed to the QtlSumStats loader for \code{sumStatsManifest}. +#' @param transposeCovariates Transpose covariate TSVs (QTLtools layout). +#' @param scaleResiduals,mafCutoff,macCutoff,xvarCutoff,imissCutoff,keepSamples,keepVariants,keepIndel +#' Pass-through \code{\link{QtlDataset}} arguments applied to every study. +#' @return A \code{MultiStudyQtlDataset} object. +#' @export +loadMultiStudyQtlDatasetFromManifest <- function(qtlDatasetsManifest, + sumStatsManifest = NULL, + genome = NULL, ldSketch = NULL, + region = NULL, + minLdOverlapWarn = 0.5, + columnMapping = NULL, + sampleSelect = NULL, + formatMapping = NULL, + transposeCovariates = FALSE, + scaleResiduals = TRUE, + mafCutoff = 0, macCutoff = 0, + xvarCutoff = 0, imissCutoff = 0, + keepSamples = character(0), + keepVariants = character(0), + keepIndel = TRUE) { + base <- .manifestBase(qtlDatasetsManifest) + df <- .canonManifestCols(.readManifest(qtlDatasetsManifest), + .qtlDatasetManifestAliases, + required = c("study", "genotypePath", "context", + "phenotypePath"), + label = "MultiStudyQtlDataset (qtlDatasetsManifest)") + qc <- list(scaleResiduals = scaleResiduals, mafCutoff = mafCutoff, + macCutoff = macCutoff, xvarCutoff = xvarCutoff, + imissCutoff = imissCutoff, keepSamples = keepSamples, + keepVariants = keepVariants, keepIndel = keepIndel) + studies <- unique(as.character(df$study)) + qtlDatasets <- setNames(lapply(studies, function(s) { + .buildQtlDatasetFromRows(df[as.character(df$study) == s, , drop = FALSE], + s, base, transposeCovariates, qc) + }), studies) + + sumStats <- NULL + if (!is.null(sumStatsManifest)) { + sumStats <- loadQtlSumStatsFromManifest( + sumStatsManifest, genome = genome, ldSketch = ldSketch, region = region, + minLdOverlapWarn = minLdOverlapWarn, columnMapping = columnMapping, + sampleSelect = sampleSelect, formatMapping = formatMapping) + } + MultiStudyQtlDataset(qtlDatasets = qtlDatasets, sumStats = sumStats) +} diff --git a/R/twasWeightsPipeline.R b/R/twasWeightsPipeline.R index dc057bb8..5c9c9c01 100644 --- a/R/twasWeightsPipeline.R +++ b/R/twasWeightsPipeline.R @@ -424,6 +424,24 @@ names(.fineMappingMethodCapabilities) } +# Canonical fine-mapping tokens actually present as methods in a +# FineMappingResult, matched tolerantly across canonical / camelCase / +# snake_case spellings (mirrors the candidate logic in .twasFineMappingFits). +# @noRd +.twasFineMappingMethodsPresent <- function(fineMappingResult) { + if (is.null(fineMappingResult)) return(character(0)) + methods <- tolower(as.character(fineMappingResult$method)) + present <- character(0) + for (canonical in .twasFineMappingTokens()) { + candidates <- tolower(c( + canonical, + paste0(tolower(substring(canonical, 1L, 1L)), substring(canonical, 2L)), + gsub("([A-Z])", "_\\1", canonical))) + if (any(methods %in% candidates)) present <- c(present, canonical) + } + present +} + # Reject fine-mapping methods (susie / susieInf / susieAsh / mvsusie / # fsusie) when no FineMappingResult is supplied. twasWeightsPipeline is # not allowed to re-fit fine-mapping models from scratch; users must run @@ -459,6 +477,16 @@ if (!is(fineMappingResult, "FineMappingResultBase")) { stop("`fineMappingResult` must be a FineMappingResult or NULL.") } + # A supplied fineMappingResult must actually contain each requested + # fine-mapping method; twasWeightsPipeline never re-fits them from scratch. + missingMethods <- setdiff(fmTokens, + .twasFineMappingMethodsPresent(fineMappingResult)) + if (length(missingMethods) > 0L) { + stop(sprintf( + "twasWeightsPipeline: method(s) %s were requested but the supplied fineMappingResult contains no such fine-mapping fit. Run fineMappingPipeline() with method(s) %s first and pass the result via `fineMappingResult = `.", + paste(unique(missingMethods), collapse = ", "), + paste(unique(missingMethods), collapse = ", "))) + } invisible(NULL) } diff --git a/_pkgdown.yml b/_pkgdown.yml index 14d4ba82..383dd0c3 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -45,6 +45,13 @@ navbar: aria-label: View source on GitHub articles: + - title: "Building input objects" + desc: > + Constructing the individual-level and summary-statistic S4 containers + from manifests or by hand. + contents: + - constructing-qtl-datasets + - constructing-sumstats - title: "Quality Control" desc: > Allele alignment and LD-mismatch QC for GWAS and QTL summary @@ -122,6 +129,17 @@ reference: - TwasWeights - TwasWeightsEntry + - title: "Manifest loaders" + desc: > + Build the individual-level and summary-statistic containers from a + manifest (a data.frame or a delimited file) that points at the + phenotype / genotype / summary-statistic sources. + contents: + - loadQtlDatasetFromManifest + - loadMultiStudyQtlDatasetFromManifest + - loadQtlSumStatsFromManifest + - loadGwasSumStatsFromManifest + - title: "Class methods" desc: > Accessor and behaviour methods grouped by the class they're @@ -499,14 +517,12 @@ reference: - metaSldscRandom - qtlEnrichment - - title: "Bundled datasets and example helpers" + - title: "Bundled example datasets" desc: > Synthetic data shipped with the package for vignettes and - tests. `fixupExampleGenotypePaths()` re-points the bundled - GenotypeHandle paths at the install-time `inst/extdata/` - location. + tests. Their bundled GenotypeHandle paths resolve to the installed + `inst/extdata/` location automatically at extraction time. contents: - - fixupExampleGenotypePaths - qtl_dataset_example - qtl_sumstats_example - qtl_sumstats_multicontext_example diff --git a/data-raw/build_examples.R b/data-raw/build_examples.R index e8d2ca06..40fd19a2 100644 --- a/data-raw/build_examples.R +++ b/data-raw/build_examples.R @@ -245,21 +245,23 @@ qtl_sumstats_multicontext_example <- QtlSumStats( qcInfo = list(prebuilt = "synthetic multi-context example; QC bypassed")) # ----------------------------------------------------------------------------- -# 7. Strip the GenotypeHandle paths down to bare basenames so the bundled -# .rda objects don't carry source-tree-absolute paths. Vignettes / -# users resolve them at use-time via fixupExampleGenotypePaths(). +# 7. Rewrite the GenotypeHandle paths to a portable bundled-resource reference +# ("pecotmr://extdata/") so the .rda objects carry no source-tree +# absolute path. The genotype readers resolve this via system.file() at +# extraction time, so no user-facing fixup step is required. # ----------------------------------------------------------------------------- -qtl_dataset_example@genotypes@path <- basename( +asResource <- function(p) paste0("pecotmr://extdata/", basename(p)) +qtl_dataset_example@genotypes@path <- asResource( qtl_dataset_example@genotypes@path) -qtl_sumstats_example@ldSketch@path <- basename( +qtl_sumstats_example@ldSketch@path <- asResource( qtl_sumstats_example@ldSketch@path) -qtl_sumstats_multicontext_example@ldSketch@path <- basename( +qtl_sumstats_multicontext_example@ldSketch@path <- asResource( qtl_sumstats_multicontext_example@ldSketch@path) -gwas_sumstats_s4_example@ldSketch@path <- basename( +gwas_sumstats_s4_example@ldSketch@path <- asResource( gwas_sumstats_s4_example@ldSketch@path) for (nm in names(multi_study_qtl_dataset_example@qtlDatasets)) multi_study_qtl_dataset_example@qtlDatasets[[nm]]@genotypes@path <- - basename(multi_study_qtl_dataset_example@qtlDatasets[[nm]]@genotypes@path) + asResource(multi_study_qtl_dataset_example@qtlDatasets[[nm]]@genotypes@path) # ----------------------------------------------------------------------------- # 8. Save. diff --git a/data/gwas_sumstats_s4_example.rda b/data/gwas_sumstats_s4_example.rda index d59338e9..8232c763 100644 Binary files a/data/gwas_sumstats_s4_example.rda and b/data/gwas_sumstats_s4_example.rda differ diff --git a/data/multi_study_qtl_dataset_example.rda b/data/multi_study_qtl_dataset_example.rda index c17c6f3e..73f52791 100644 Binary files a/data/multi_study_qtl_dataset_example.rda and b/data/multi_study_qtl_dataset_example.rda differ diff --git a/data/qtl_dataset_example.rda b/data/qtl_dataset_example.rda index 8b7b2d31..b6ed7994 100644 Binary files a/data/qtl_dataset_example.rda and b/data/qtl_dataset_example.rda differ diff --git a/data/qtl_sumstats_example.rda b/data/qtl_sumstats_example.rda index f1bcaf7d..63078f0e 100644 Binary files a/data/qtl_sumstats_example.rda and b/data/qtl_sumstats_example.rda differ diff --git a/data/qtl_sumstats_multicontext_example.rda b/data/qtl_sumstats_multicontext_example.rda index e352ed3b..5018e503 100644 Binary files a/data/qtl_sumstats_multicontext_example.rda and b/data/qtl_sumstats_multicontext_example.rda differ diff --git a/man/fixupExampleGenotypePaths.Rd b/man/fixupExampleGenotypePaths.Rd deleted file mode 100644 index 8819d929..00000000 --- a/man/fixupExampleGenotypePaths.Rd +++ /dev/null @@ -1,26 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/exampleData.R -\name{fixupExampleGenotypePaths} -\alias{fixupExampleGenotypePaths} -\title{Resolve bundled-example GenotypeHandle paths} -\usage{ -fixupExampleGenotypePaths(x) -} -\arguments{ -\item{x}{A bundled \code{QtlDataset}, \code{MultiStudyQtlDataset}, -\code{QtlSumStats}, or \code{GwasSumStats} example object.} -} -\value{ -The same object with \code{GenotypeHandle@path} rewritten - to the resolved install path. -} -\description{ -The bundled S4 example objects (\code{qtl_dataset_example}, -\code{qtl_sumstats_example}, \code{qtl_sumstats_multicontext_example}, -\code{gwas_sumstats_s4_example}, \code{multi_study_qtl_dataset_example}) -store a relative-style path to the bundled -\code{inst/extdata/toy_canonical} PLINK1 reference. -Call this helper once at the top of a vignette to re-point each -contained \code{GenotypeHandle} at the installed path resolved via -\code{system.file()}. -} diff --git a/man/loadGwasSumStatsFromManifest.Rd b/man/loadGwasSumStatsFromManifest.Rd new file mode 100644 index 00000000..48c03c65 --- /dev/null +++ b/man/loadGwasSumStatsFromManifest.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/manifestLoaders.R +\name{loadGwasSumStatsFromManifest} +\alias{loadGwasSumStatsFromManifest} +\title{Load a GwasSumStats collection from a manifest} +\usage{ +loadGwasSumStatsFromManifest( + manifest, + genome = NULL, + ldSketch = NULL, + region = NULL, + minLdOverlapWarn = 0.5, + columnMapping = NULL, + sampleSelect = NULL, + formatMapping = NULL +) +} +\arguments{ +\item{manifest}{A data.frame or path. Columns (snake_case aliases accepted): +\code{study} (required, unique), \code{sumStatsPath} (required), +\code{columnMapping} (optional), \code{nCase} / \code{nControl} +(optional), \code{varY} (optional), and the single-valued \code{genome} / +\code{ldSketchPath}.} + +\item{genome}{Genome build; reconciled with a \code{genome} column.} + +\item{ldSketch}{A \code{\link{GenotypeHandle}} or a spec +(path/prefix/genoMeta); reconciled with an \code{ldSketchPath} column.} + +\item{region}{Optional \code{chr:start-end} string, GRanges, or one-row +data.frame restricting the variants read (honoured only for tabix-indexed +text and bgzipped+tabixed VCF; ignored with a warning otherwise).} + +\item{minLdOverlapWarn}{Warn when the fraction of sumstats variants present +in the LD sketch falls below this (default 0.5).} + +\item{columnMapping}{Optional default column mapping (a named list/vector, +or a path to a YAML file of \code{standardName: sourceName} entries) +applied when a row has no \code{columnMapping}.} + +\item{sampleSelect}{Optional GWAS-VCF FORMAT sample (study) column to read.} + +\item{formatMapping}{Optional GWAS-VCF FORMAT tag mapping (canonical stat -> +FORMAT field), overriding ES/SE/LP/SS/EAF defaults.} +} +\value{ +A \code{GwasSumStats} object. +} +\description{ +Build a \code{\link{GwasSumStats}} from a manifest with one row + per study. No QC is run (the result carries \code{qcInfo = list()}). +} diff --git a/man/loadMultiStudyQtlDatasetFromManifest.Rd b/man/loadMultiStudyQtlDatasetFromManifest.Rd new file mode 100644 index 00000000..ac0b5a5a --- /dev/null +++ b/man/loadMultiStudyQtlDatasetFromManifest.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/manifestLoaders.R +\name{loadMultiStudyQtlDatasetFromManifest} +\alias{loadMultiStudyQtlDatasetFromManifest} +\title{Load a MultiStudyQtlDataset from manifests} +\usage{ +loadMultiStudyQtlDatasetFromManifest( + qtlDatasetsManifest, + sumStatsManifest = NULL, + genome = NULL, + ldSketch = NULL, + region = NULL, + minLdOverlapWarn = 0.5, + columnMapping = NULL, + sampleSelect = NULL, + formatMapping = NULL, + transposeCovariates = FALSE, + scaleResiduals = TRUE, + mafCutoff = 0, + macCutoff = 0, + xvarCutoff = 0, + imissCutoff = 0, + keepSamples = character(0), + keepVariants = character(0), + keepIndel = TRUE +) +} +\arguments{ +\item{qtlDatasetsManifest}{A data.frame or path. The QtlDataset schema with +\code{study} and \code{genotypePath} \strong{mandatory} (grouping key = +\code{study}); one \code{QtlDataset} is built per study group.} + +\item{sumStatsManifest}{Optional QtlSumStats manifest (see +\code{\link{loadQtlSumStatsFromManifest}}).} + +\item{genome, ldSketch, region, minLdOverlapWarn, columnMapping, sampleSelect, formatMapping}{Passed to the QtlSumStats loader for \code{sumStatsManifest}.} + +\item{transposeCovariates}{Transpose covariate TSVs (QTLtools layout).} + +\item{scaleResiduals, mafCutoff, macCutoff, xvarCutoff, imissCutoff, keepSamples, keepVariants, keepIndel}{Pass-through \code{\link{QtlDataset}} arguments applied to every study.} +} +\value{ +A \code{MultiStudyQtlDataset} object. +} +\description{ +Build a \code{\link{MultiStudyQtlDataset}} from a QtlDataset + manifest (individual-level studies) and, optionally, a QtlSumStats + manifest (summary-only studies). +} diff --git a/man/loadQtlDatasetFromManifest.Rd b/man/loadQtlDatasetFromManifest.Rd new file mode 100644 index 00000000..d3cd56c9 --- /dev/null +++ b/man/loadQtlDatasetFromManifest.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/manifestLoaders.R +\name{loadQtlDatasetFromManifest} +\alias{loadQtlDatasetFromManifest} +\title{Load a QtlDataset from a manifest} +\usage{ +loadQtlDatasetFromManifest( + manifest, + study = NULL, + genotypes = NULL, + genotypeCovariates = NULL, + scaleResiduals = TRUE, + mafCutoff = 0, + macCutoff = 0, + xvarCutoff = 0, + imissCutoff = 0, + keepSamples = character(0), + keepVariants = character(0), + keepIndel = TRUE, + transposeCovariates = FALSE +) +} +\arguments{ +\item{manifest}{A data.frame or a path to a manifest file. Recognised +columns (canonical camelCase; snake_case aliases accepted): +\code{context} (required), \code{phenotypePath} (required, bgzipped BED), +\code{covariatePath} (optional, per-context covariates), and the +single-valued \code{study} / \code{genotypePath} / +\code{genotypeCovariatePath}.} + +\item{study}{Study identifier; reconciled with a \code{study} column.} + +\item{genotypes}{A \code{\link{GenotypeHandle}} or a genotype path/prefix; +reconciled with a \code{genotypePath} column.} + +\item{genotypeCovariates}{A numeric matrix (samples x covariates) or a path +to a covariate TSV; reconciled with a \code{genotypeCovariatePath} column.} + +\item{scaleResiduals, mafCutoff, macCutoff, xvarCutoff, imissCutoff, keepSamples, keepVariants, keepIndel}{Pass-through \code{\link{QtlDataset}} arguments (stored as lazy QC slots).} + +\item{transposeCovariates}{Transpose covariate TSVs (QTLtools layout) before +treating them as samples-as-rows.} +} +\value{ +A \code{QtlDataset} object. +} +\description{ +Build a single-study \code{\link{QtlDataset}} from a manifest + describing one row per QTL context (or per trait x context). The manifest + is a data.frame or a path to a delimited file (\code{.csv} -> CSV, else + TSV). Column names may be snake_case or camelCase. +} diff --git a/man/loadQtlSumStatsFromManifest.Rd b/man/loadQtlSumStatsFromManifest.Rd new file mode 100644 index 00000000..259b8a7b --- /dev/null +++ b/man/loadQtlSumStatsFromManifest.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/manifestLoaders.R +\name{loadQtlSumStatsFromManifest} +\alias{loadQtlSumStatsFromManifest} +\title{Load a QtlSumStats collection from a manifest} +\usage{ +loadQtlSumStatsFromManifest( + manifest, + genome = NULL, + ldSketch = NULL, + region = NULL, + minLdOverlapWarn = 0.5, + columnMapping = NULL, + sampleSelect = NULL, + formatMapping = NULL +) +} +\arguments{ +\item{manifest}{A data.frame or path. Columns (snake_case aliases accepted): +\code{study}, \code{context}, \code{trait} (required), \code{sumStatsPath} +(required), \code{columnMapping} (optional), \code{varY} (optional), and +the single-valued \code{genome} / \code{ldSketchPath}.} + +\item{genome}{Genome build; reconciled with a \code{genome} column.} + +\item{ldSketch}{A \code{\link{GenotypeHandle}} or spec; reconciled with an +\code{ldSketchPath} column.} + +\item{region, minLdOverlapWarn, columnMapping, sampleSelect, formatMapping}{As +for \code{\link{loadGwasSumStatsFromManifest}}.} +} +\value{ +A \code{QtlSumStats} object. +} +\description{ +Build a \code{\link{QtlSumStats}} from a manifest with one row + per \code{(study, context, trait)} tuple. No QC is run. +} diff --git a/pixi.toml b/pixi.toml index 7b95def4..419b6bee 100644 --- a/pixi.toml +++ b/pixi.toml @@ -103,3 +103,4 @@ r45 = {features = ["r45"]} "r-vctrs" = "*" "r-vroom" = "*" "r-xgboost" = "*" +"r-yaml" = "*" diff --git a/tests/testthat/test_ctwasPipeline.R b/tests/testthat/test_ctwasPipeline.R index 343fa3c2..46fb9ce0 100644 --- a/tests/testthat/test_ctwasPipeline.R +++ b/tests/testthat/test_ctwasPipeline.R @@ -899,8 +899,8 @@ test_that("ctwasPipeline: real-engine end-to-end on the bundled example panel", skip_if_not_installed("ctwas") data(gwas_sumstats_s4_example) data(qtl_dataset_example) - gss <- fixupExampleGenotypePaths(gwas_sumstats_s4_example) - qd <- fixupExampleGenotypePaths(qtl_dataset_example) + gss <- gwas_sumstats_s4_example + qd <- qtl_dataset_example gh <- qd@genotypes # 5-variant synthetic gene from the bundled panel. @@ -1015,7 +1015,7 @@ test_that(".ctwasFilterVariants: returns NULL when no variants survive", { test_that(".ctwasBuildWeights: maxNumVariants caps the per-gene weight matrix", { data(qtl_dataset_example) - qd <- fixupExampleGenotypePaths(qtl_dataset_example) + qd <- qtl_dataset_example gh <- qd@genotypes vids <- gh@snpInfo$SNP[1:5] ent <- TwasWeightsEntry( @@ -1035,7 +1035,7 @@ test_that(".ctwasBuildWeights: maxNumVariants caps the per-gene weight matrix", test_that(".ctwasBuildWeights: twasWeightCutoff drops low-magnitude variants", { data(qtl_dataset_example) - qd <- fixupExampleGenotypePaths(qtl_dataset_example) + qd <- qtl_dataset_example gh <- qd@genotypes vids <- gh@snpInfo$SNP[1:5] ent <- TwasWeightsEntry( diff --git a/tests/testthat/test_manifestLoaders.R b/tests/testthat/test_manifestLoaders.R new file mode 100644 index 00000000..65647861 --- /dev/null +++ b/tests/testthat/test_manifestLoaders.R @@ -0,0 +1,807 @@ +# Tests for R/manifestLoaders.R + +# --------------------------------------------------------------------------- +# Fixtures +# --------------------------------------------------------------------------- + +# The toy_ref PLINK1 panel (chr22) serves as the LD sketch. These variants +# are drawn verbatim from toy_ref.bim so LD-containment succeeds. +.toyRefPrefix <- function() { + sub("\\.bed$", "", system.file("extdata", "toy_ref.bed", package = "pecotmr")) +} + +.toyLdSketch <- function() { + GenotypeHandle(plink1Prefix = .toyRefPrefix()) +} + +# A small GWAS sumstats data.frame over real toy_ref chr22 variants. +.toyGwasDf <- function(n = 5) { + base <- 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.2, -0.4, 2.1, 0.05, -1.8), + n_sample = rep(10000L, 5), + stringsAsFactors = FALSE) + base[seq_len(n), , drop = FALSE] +} + +# Write a data.frame to a TSV whose header line starts with "#" (so a tabix +# index can treat it as a comment/header line). +.writeSumstatsTsv <- function(df, path) { + hdr <- names(df); hdr[[1L]] <- paste0("#", hdr[[1L]]) + con <- file(path, "w") + writeLines(paste(hdr, collapse = "\t"), con) + utils::write.table(df, con, sep = "\t", quote = FALSE, row.names = FALSE, + col.names = FALSE) + close(con) + path +} + +# --------------------------------------------------------------------------- +# Column canonicalisation + scalar reconciliation +# --------------------------------------------------------------------------- + +test_that(".canonManifestCols renames snake_case aliases to camelCase", { + df <- data.frame(study_id = "s1", sumstats_path = "x.tsv", + stringsAsFactors = FALSE) + out <- pecotmr:::.canonManifestCols(df, pecotmr:::.gwasSumStatsManifestAliases, + required = c("study", "sumStatsPath"), + label = "GwasSumStats") + expect_true(all(c("study", "sumStatsPath") %in% names(out))) + expect_false(any(c("study_id", "sumstats_path") %in% names(out))) +}) + +test_that(".canonManifestCols errors on missing required columns", { + df <- data.frame(study = "s1", stringsAsFactors = FALSE) + expect_error( + pecotmr:::.canonManifestCols(df, pecotmr:::.gwasSumStatsManifestAliases, + required = c("study", "sumStatsPath"), + label = "GwasSumStats"), + "missing required column") +}) + +test_that(".reconcileScalar resolves arg/column and flags conflicts", { + expect_equal(pecotmr:::.reconcileScalar(NULL, "hg38", "genome"), "hg38") + expect_equal(pecotmr:::.reconcileScalar(c("hg38", "hg38"), NULL, "genome"), "hg38") + expect_equal(pecotmr:::.reconcileScalar(c("hg38", "hg38"), "hg38", "genome"), "hg38") + expect_error(pecotmr:::.reconcileScalar(c("hg19", "hg38"), NULL, "genome"), + "must be constant") + expect_error(pecotmr:::.reconcileScalar("hg19", "hg38", "genome"), "disagrees") + expect_error(pecotmr:::.reconcileScalar(NULL, NULL, "genome"), "must be provided") +}) + +# --------------------------------------------------------------------------- +# Genotype-format detection +# --------------------------------------------------------------------------- + +test_that(".detectGenotypeFormat builds a PLINK1 handle from a prefix", { + h <- pecotmr:::.detectGenotypeFormat(.toyRefPrefix()) + expect_s4_class(h, "GenotypeHandle") + expect_equal(h@format, "plink1") +}) + +test_that("BCF sumstats are rejected", { + expect_error( + pecotmr:::.readSumStatsFile("x.bcf", NULL, NULL, NULL, NULL, "lbl"), + "BCF sumstats are not supported") +}) + +# --------------------------------------------------------------------------- +# LD-sketch containment +# --------------------------------------------------------------------------- + +test_that(".checkLdContainment passes for overlapping variants", { + gr <- pecotmr:::.dfToEntryGranges(.toyGwasDf(5)) + expect_silent(pecotmr:::.checkLdContainment(.toyLdSketch(), gr, 0.5, "lbl")) +}) + +test_that(".checkLdContainment errors on a chromosome absent from the sketch", { + df <- .toyGwasDf(3); df$chrom <- "1" + gr <- pecotmr:::.dfToEntryGranges(df) + expect_error(pecotmr:::.checkLdContainment(.toyLdSketch(), gr, 0.5, "lbl"), + "absent from the LD sketch") +}) + +test_that(".checkLdContainment errors on zero overlap", { + df <- .toyGwasDf(2); df$pos <- c(999999001, 999999002) + gr <- pecotmr:::.dfToEntryGranges(df) + expect_error(pecotmr:::.checkLdContainment(.toyLdSketch(), gr, 0.5, "lbl"), + "none of the") +}) + +test_that(".checkLdContainment warns on low overlap", { + df <- .toyGwasDf(4) + df$pos[3:4] <- c(999999001, 999999002) # 2/4 present -> 50% + gr <- pecotmr:::.dfToEntryGranges(df) + expect_warning(pecotmr:::.checkLdContainment(.toyLdSketch(), gr, 0.9, "lbl"), + "are present in the LD sketch") +}) + +# --------------------------------------------------------------------------- +# GwasSumStats loader +# --------------------------------------------------------------------------- + +test_that("loadGwasSumStatsFromManifest builds from a data.frame manifest", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "study1.tsv")) + manifest <- data.frame(study = "study1", sumStatsPath = ssPath, + stringsAsFactors = FALSE) + obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_s4_class(obj, "GwasSumStats") + expect_true(methods::validObject(obj)) + expect_equal(as.character(obj$study), "study1") + expect_equal(getGenome(obj), "hg38") + expect_equal(length(obj$entry[[1L]]), 5L) + expect_true(all(c("SNP", "A1", "A2", "Z", "N") %in% + colnames(S4Vectors::mcols(obj$entry[[1L]])))) + expect_equal(length(obj@qcInfo), 0L) # loaders run no QC +}) + +test_that("loadGwasSumStatsFromManifest reads a manifest file and reconciles genome", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "study1.tsv")) + mfPath <- file.path(tmp, "gwas_manifest.tsv") + readr::write_tsv(data.frame(study = "study1", sumStatsPath = ssPath, + genome = "hg38", stringsAsFactors = FALSE), mfPath) + # genome column present; matching arg is fine. + obj <- loadGwasSumStatsFromManifest(mfPath, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_equal(getGenome(obj), "hg38") + # conflicting arg errors. + expect_error(loadGwasSumStatsFromManifest(mfPath, genome = "hg19", + ldSketch = .toyLdSketch()), + "disagrees") +}) + +test_that("loadGwasSumStatsFromManifest rejects duplicate studies", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(3), file.path(tmp, "s.tsv")) + manifest <- data.frame(study = c("a", "a"), + sumStatsPath = c(ssPath, ssPath), + stringsAsFactors = FALSE) + expect_error(loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()), + "must be unique") +}) + +# --------------------------------------------------------------------------- +# Region reading (tabix vs plain text) +# --------------------------------------------------------------------------- + +test_that("tabix-indexed text honours the region", { + skip_if_not_installed("Rsamtools") + tmp <- withr::local_tempdir() + tsv <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "s.tsv")) + gz <- Rsamtools::bgzip(tsv, dest = paste0(tsv, ".gz"), overwrite = TRUE) + Rsamtools::indexTabix(gz, seq = 1L, start = 2L, end = 2L, comment = "#") + df <- pecotmr:::.readSumStatsText(gz, "22:14560000-14565000", NULL, "lbl") + expect_equal(nrow(df), 2L) # first two variants only + expect_true(all(df$pos >= 14560000 & df$pos <= 14565000)) +}) + +test_that("plain-text region is ignored with a warning", { + tmp <- withr::local_tempdir() + tsv <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "s.tsv")) + expect_warning( + df <- pecotmr:::.readSumStatsText(tsv, "22:14560000-14565000", NULL, "lbl"), + "region ignored") + expect_equal(nrow(df), 5L) +}) + +# --------------------------------------------------------------------------- +# YAML column mapping +# --------------------------------------------------------------------------- + +test_that("loadGwasSumStatsFromManifest resolves a YAML column mapping", { + tmp <- withr::local_tempdir() + # Sumstats with non-standard source column names. + df <- .toyGwasDf(5) + names(df) <- c("CHR", "POS", "RSID", "EA", "OA", "ZSCORE", "SAMPLE_N") + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "study1.tsv")) + mapPath <- file.path(tmp, "mapping.yaml") + yaml::write_yaml(list(chrom = "CHR", pos = "POS", variant_id = "RSID", + A1 = "EA", A2 = "OA", z = "ZSCORE", + n_sample = "SAMPLE_N"), mapPath) + manifest <- data.frame(study = "study1", sumStatsPath = ssPath, + columnMapping = mapPath, stringsAsFactors = FALSE) + obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_s4_class(obj, "GwasSumStats") + expect_equal(length(obj$entry[[1L]]), 5L) + expect_equal(as.character(S4Vectors::mcols(obj$entry[[1L]])$A1), + c("G", "C", "G", "C", "A")) +}) + +test_that(".readColumnMapping accepts a named list and a YAML file alike", { + tmp <- withr::local_tempdir() + mp <- file.path(tmp, "m.yaml") + yaml::write_yaml(list(chrom = "CHR", z = "ZSCORE"), mp) + fromFile <- pecotmr:::.readColumnMapping(mp) + fromList <- pecotmr:::.readColumnMapping(list(chrom = "CHR", z = "ZSCORE")) + expect_equal(fromFile[["chrom"]], "CHR") + expect_equal(fromFile[["z"]], "ZSCORE") + expect_equal(fromList, fromFile) +}) + +# --------------------------------------------------------------------------- +# QtlSumStats loader +# --------------------------------------------------------------------------- + +test_that("loadQtlSumStatsFromManifest builds a per-tuple collection", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "geneA.tsv")) + manifest <- data.frame(study = "eqtl", context = "Whole_Blood", + trait = "geneA", sumStatsPath = ssPath, + stringsAsFactors = FALSE) + obj <- loadQtlSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_s4_class(obj, "QtlSumStats") + expect_true(methods::validObject(obj)) + expect_equal(as.character(obj$trait), "geneA") + expect_equal(as.character(obj$context), "Whole_Blood") +}) + +# --------------------------------------------------------------------------- +# QtlDataset loader +# --------------------------------------------------------------------------- + +# Build a tiny 2-trait phenotype BED for one context. +.writePhenoBed <- function(path, samples = paste0("S", 1:6)) { + df <- data.frame( + `#chr` = c("chr22", "chr22"), + start = c(14560000L, 14870000L), + end = c(14560001L, 14870001L), + gene_id = c("geneA", "geneB"), + check.names = FALSE, stringsAsFactors = FALSE) + for (s in samples) df[[s]] <- stats::rnorm(2) + readr::write_tsv(df, path) + path +} + +test_that("loadQtlDatasetFromManifest builds a single-study dataset", { + tmp <- withr::local_tempdir() + bed <- .writePhenoBed(file.path(tmp, "ctxA.bed")) + manifest <- data.frame(context = "ctxA", phenotypePath = bed, + study = "studyX", genotypePath = .toyRefPrefix(), + stringsAsFactors = FALSE) + qd <- loadQtlDatasetFromManifest(manifest) + expect_s4_class(qd, "QtlDataset") + expect_true(methods::validObject(qd)) + expect_equal(getStudy(qd), "studyX") + expect_equal(getContexts(qd), "ctxA") +}) + +test_that("loadQtlDatasetFromManifest accepts study/genotypes as arguments", { + tmp <- withr::local_tempdir() + bed <- .writePhenoBed(file.path(tmp, "ctxA.bed")) + manifest <- data.frame(context = "ctxA", phenotypePath = bed, + stringsAsFactors = FALSE) + qd <- loadQtlDatasetFromManifest(manifest, study = "studyY", + genotypes = .toyLdSketch()) + expect_equal(getStudy(qd), "studyY") +}) + +# --------------------------------------------------------------------------- +# MultiStudyQtlDataset loader +# --------------------------------------------------------------------------- + +test_that("loadMultiStudyQtlDatasetFromManifest builds from >=2 studies", { + tmp <- withr::local_tempdir() + bedA <- .writePhenoBed(file.path(tmp, "a.bed")) + bedB <- .writePhenoBed(file.path(tmp, "b.bed")) + manifest <- data.frame( + study = c("study1", "study2"), + context = c("ctxA", "ctxA"), + phenotypePath = c(bedA, bedB), + genotypePath = rep(.toyRefPrefix(), 2), + stringsAsFactors = FALSE) + msd <- loadMultiStudyQtlDatasetFromManifest(manifest) + expect_s4_class(msd, "MultiStudyQtlDataset") + expect_true(methods::validObject(msd)) + expect_equal(sort(names(msd@qtlDatasets)), c("study1", "study2")) +}) + +test_that("loadMultiStudyQtlDatasetFromManifest attaches a summary-only study", { + tmp <- withr::local_tempdir() + bedA <- .writePhenoBed(file.path(tmp, "a.bed")) + qdMan <- data.frame(study = "study1", context = "ctxA", phenotypePath = bedA, + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "geneA.tsv")) + ssMan <- data.frame(study = "eqtl2", context = "Blood", trait = "geneA", + sumStatsPath = ssPath, stringsAsFactors = FALSE) + msd <- loadMultiStudyQtlDatasetFromManifest( + qdMan, sumStatsManifest = ssMan, genome = "hg38", ldSketch = .toyLdSketch()) + expect_true(methods::validObject(msd)) + expect_s4_class(msd@sumStats, "QtlSumStats") + expect_equal(names(msd@qtlDatasets), "study1") +}) + +# =========================================================================== +# Coverage: manifest ingestion, path resolution, scalar reconciliation +# =========================================================================== + +test_that(".readManifest handles CSV, missing files, and bad input", { + tmp <- withr::local_tempdir() + csv <- file.path(tmp, "m.csv") + ssPath <- .writeSumstatsTsv(.toyGwasDf(3), file.path(tmp, "s.tsv")) + readr::write_csv(data.frame(study = "s1", sumStatsPath = ssPath), csv) + m <- pecotmr:::.readManifest(csv) + expect_true(is.data.frame(m) && m$study == "s1") + expect_error(pecotmr:::.readManifest("/no/such/manifest.tsv"), "not found") + expect_error(pecotmr:::.readManifest(42L), "data.frame or a single file path") +}) + +test_that(".resolveRel resolves relative paths against a base only", { + expect_equal(pecotmr:::.resolveRel("a.tsv", "/base"), "/base/a.tsv") + expect_equal(pecotmr:::.resolveRel("/abs/a.tsv", "/base"), "/abs/a.tsv") + expect_equal(pecotmr:::.resolveRel("a.tsv", NULL), "a.tsv") + expect_equal(pecotmr:::.resolveRel("", "/base"), "") +}) + +test_that("loadGwasSumStatsFromManifest resolves manifest-relative paths", { + tmp <- withr::local_tempdir() + .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "study1.tsv")) + mfPath <- file.path(tmp, "manifest.tsv") + # Relative sumStatsPath resolves against the manifest's own directory. + readr::write_tsv(data.frame(study = "study1", sumStatsPath = "study1.tsv"), + mfPath) + obj <- loadGwasSumStatsFromManifest(mfPath, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_equal(length(obj$entry[[1L]]), 5L) +}) + +test_that(".reconcileScalar returns NULL for an optional unset scalar", { + expect_null(pecotmr:::.reconcileScalar(NULL, NULL, "x", required = FALSE)) +}) + +# =========================================================================== +# Coverage: genotype-format detection and LD-sketch resolution +# =========================================================================== + +test_that(".detectGenotypeFormat dispatches by extension", { + expect_error(pecotmr:::.detectGenotypeFormat("x.vcf.gz")) # vcf branch + expect_error(pecotmr:::.detectGenotypeFormat("x.pgen")) # plink2 branch + expect_error(pecotmr:::.detectGenotypeFormat("x.unknown"), + "Could not determine genotype format") + h <- pecotmr:::.detectGenotypeFormat(paste0(.toyRefPrefix(), ".bed")) + expect_equal(h@format, "plink1") +}) + +test_that(".resolveLdSketch accepts a genoMeta vector, a path, and rejects bad input", { + sharded <- pecotmr:::.resolveLdSketch(c("22" = .toyRefPrefix())) + expect_s4_class(sharded, "GenotypeHandle") + expect_true("22" %in% names(sharded@chromPaths)) + expect_s4_class(pecotmr:::.resolveLdSketch(.toyRefPrefix()), "GenotypeHandle") + expect_error(pecotmr:::.resolveLdSketch("nonexistent_meta.tsv")) + expect_error(pecotmr:::.resolveLdSketch(42L), "must be a GenotypeHandle") +}) + +test_that("ldSketch is resolved from an ldSketchPath column and conflicts error", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "s.tsv")) + manifest <- data.frame(study = "s1", sumStatsPath = ssPath, + ldSketchPath = .toyRefPrefix(), stringsAsFactors = FALSE) + obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38") + expect_s4_class(obj, "GwasSumStats") + # Conflicting arg vs column. + expect_error( + loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = "/some/other/prefix"), + "disagrees") + # Neither arg nor column. + bare <- data.frame(study = "s1", sumStatsPath = ssPath, stringsAsFactors = FALSE) + expect_error(loadGwasSumStatsFromManifest(bare, genome = "hg38"), + "must be provided as an argument or an .ldSketchPath") +}) + +test_that("LD containment uses a sharded sketch's chromPaths", { + sharded <- GenotypeHandle(genoMeta = c("22" = .toyRefPrefix())) + gr <- pecotmr:::.dfToEntryGranges(.toyGwasDf(5)) + expect_silent(pecotmr:::.checkLdContainment(sharded, gr, 0.5, "lbl")) +}) + +test_that(".checkLdContainment short-circuits on an empty entry", { + expect_null(pecotmr:::.checkLdContainment(.toyLdSketch(), + GenomicRanges::GRanges(), 0.5, "lbl")) +}) + +test_that(".checkLdContainment warns when the sketch has no variant metadata", { + h <- methods::new("GenotypeHandle", path = "meta", format = "plink1", + snpInfo = data.frame(), nSamples = 0L, + sampleIds = character(), pgenPtr = NULL, + chromPaths = c("22" = "/p")) + gr <- pecotmr:::.dfToEntryGranges(.toyGwasDf(3)) + expect_warning(pecotmr:::.checkLdContainment(h, gr, 0.5, "lbl"), + "no variant metadata") +}) + +# =========================================================================== +# Coverage: region coercion and column-mapping errors +# =========================================================================== + +test_that(".asGRegion coerces strings, GRanges, and data.frames", { + gr <- GenomicRanges::GRanges("22", IRanges::IRanges(1, 100)) + expect_identical(pecotmr:::.asGRegion(gr), gr) + fromStr <- pecotmr:::.asGRegion("22:1-100") + expect_equal(GenomicRanges::start(fromStr), 1L) + fromDf <- pecotmr:::.asGRegion(data.frame(chrom = "22", start = 1, end = 100)) + expect_equal(as.character(GenomicRanges::seqnames(fromDf)), "22") + expect_error(pecotmr:::.asGRegion("not-a-region"), "chr:start-end") + expect_error(pecotmr:::.asGRegion(42L), "chr:start-end") +}) + +test_that(".readColumnMapping errors on bad inputs", { + tmp <- withr::local_tempdir() + expect_error(pecotmr:::.readColumnMapping("/no/map.yaml"), "not found") + bad <- file.path(tmp, "bad.yaml") + yaml::write_yaml(list("a", "b"), bad) # unnamed sequence + expect_error(pecotmr:::.readColumnMapping(bad), "standardName: sourceName") + expect_error(pecotmr:::.readColumnMapping(42L), "named list/vector or a path") +}) + +test_that(".resolveSumstatCols errors on bad mapping and missing fields", { + df <- .toyGwasDf(3) + names(df) <- c("chrom", "pos", "variant_id", "A1", "A2", "z", "n_sample") + expect_error( + pecotmr:::.resolveSumstatCols(df, list(chrom = "NOPE"), "lbl"), + "is not a column") + expect_error( + pecotmr:::.resolveSumstatCols(df[, c("chrom", "pos")], NULL, "lbl"), + "missing required field") +}) + +# =========================================================================== +# Coverage: nCase/nControl/varY columns +# =========================================================================== + +test_that("loadGwasSumStatsFromManifest attaches nCase/nControl/varY", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "s.tsv")) + manifest <- data.frame(study = "cc", sumStatsPath = ssPath, + nCase = 1200, nControl = 3400, varY = 0.19, + stringsAsFactors = FALSE) + obj <- loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_equal(as.numeric(obj$nCase), 1200) + expect_equal(as.numeric(obj$nControl), 3400) + expect_equal(as.numeric(obj$varY), 0.19) +}) + +test_that("loadQtlSumStatsFromManifest attaches varY and honours a region arg", { + tmp <- withr::local_tempdir() + ssPath <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "geneA.tsv")) + manifest <- data.frame(study = "eqtl", context = "Blood", trait = "geneA", + sumStatsPath = ssPath, varY = 1.7, + stringsAsFactors = FALSE) + obj <- suppressWarnings(loadQtlSumStatsFromManifest( + manifest, genome = "hg38", ldSketch = .toyLdSketch(), + region = data.frame(chrom = "22", start = 1, end = 5e7))) + expect_equal(as.numeric(obj$varY), 1.7) +}) + +# =========================================================================== +# Coverage: QtlDataset covariates and error branches +# =========================================================================== + +.writeCovTsv <- function(path, samples, npc = 3, transpose = FALSE) { + set.seed(1) + m <- matrix(stats::rnorm(length(samples) * npc), nrow = length(samples), + dimnames = list(samples, paste0("PC", seq_len(npc)))) + if (transpose) { + out <- data.frame(covariate = colnames(m), t(m), check.names = FALSE) + } else { + out <- data.frame(sample = rownames(m), m, check.names = FALSE) + } + readr::write_tsv(out, path) + path +} + +test_that("loadQtlDatasetFromManifest attaches per-context and genotype covariates", { + tmp <- withr::local_tempdir() + samples <- paste0("S", 1:6) + bed <- .writePhenoBed(file.path(tmp, "ctx.bed"), samples) + pcov <- .writeCovTsv(file.path(tmp, "pcov.tsv"), samples, npc = 3) + gcov <- .writeCovTsv(file.path(tmp, "gcov.tsv"), samples, npc = 2) + manifest <- data.frame(context = "ctx", phenotypePath = bed, + covariatePath = pcov, study = "s", + genotypePath = .toyRefPrefix(), + genotypeCovariatePath = gcov, stringsAsFactors = FALSE) + qd <- loadQtlDatasetFromManifest(manifest) + expect_equal(ncol(getPhenotypeCovariates(qd, "ctx")[["ctx"]]), 3L) + expect_equal(ncol(getGenotypeCovariates(qd)), 2L) +}) + +test_that("loadQtlDatasetFromManifest reads transposed (QTLtools) covariates", { + tmp <- withr::local_tempdir() + samples <- paste0("S", 1:6) + bed <- .writePhenoBed(file.path(tmp, "ctx.bed"), samples) + pcov <- .writeCovTsv(file.path(tmp, "pcovT.tsv"), samples, npc = 4, + transpose = TRUE) + manifest <- data.frame(context = "ctx", phenotypePath = bed, + covariatePath = pcov, study = "s", + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + qd <- loadQtlDatasetFromManifest(manifest, transposeCovariates = TRUE) + expect_equal(ncol(getPhenotypeCovariates(qd, "ctx")[["ctx"]]), 4L) +}) + +test_that("loadQtlDatasetFromManifest accepts a genotype prefix and matrix covariates", { + tmp <- withr::local_tempdir() + samples <- paste0("S", 1:6) + bed <- .writePhenoBed(file.path(tmp, "ctx.bed"), samples) + manifest <- data.frame(context = "ctx", phenotypePath = bed, study = "s", + stringsAsFactors = FALSE) + gmat <- matrix(0, nrow = 6, ncol = 2, dimnames = list(samples, c("A", "B"))) + qd <- loadQtlDatasetFromManifest(manifest, genotypes = .toyRefPrefix(), + genotypeCovariates = gmat) + expect_equal(ncol(getGenotypeCovariates(qd)), 2L) +}) + +test_that("QtlDataset builder errors on inconsistent per-context paths", { + tmp <- withr::local_tempdir() + bed1 <- .writePhenoBed(file.path(tmp, "c1.bed")) + bed2 <- .writePhenoBed(file.path(tmp, "c2.bed")) + # Same context, two different phenotype paths. + manifest <- data.frame(context = c("ctx", "ctx"), + phenotypePath = c(bed1, bed2), study = "s", + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + expect_error(loadQtlDatasetFromManifest(manifest), + "exactly one phenotypePath") +}) + +# =========================================================================== +# Coverage: GWAS-VCF reading +# =========================================================================== + +# Write a minimal GWAS-VCF (ES/SE/LP/SS/EAF FORMAT) over toy_ref chr22 sites. +# A1 (effect) = ALT, A2 = REF, matching the toy_ref alleles. +.writeGwasVcf <- function(path, samples = "study1", n = 5, dropSs = FALSE, + dropEs = FALSE) { + pos <- c(14560203, 14564328, 14850625, 14870204, 14878387)[seq_len(n)] + ids <- c("rs11089128", "rs7288972", "rs11167319", "rs8138488", "rs2186521")[seq_len(n)] + ref <- c("A", "T", "T", "T", "G")[seq_len(n)] + alt <- c("G", "C", "G", "C", "A")[seq_len(n)] + es <- c(0.12, -0.04, 0.21, 0.005, -0.18)[seq_len(n)] + se <- rep(0.1, n) + lp <- c(1.3, 0.2, 2.1, 0.05, 1.0)[seq_len(n)] + eaf <- c(0.3, 0.4, 0.2, 0.5, 0.1)[seq_len(n)] + # Build the FORMAT tag list and matching per-record cells dynamically. + tags <- c(if (!dropEs) c("ES", "SE"), "LP", if (!dropSs) "SS", "EAF") + hdrFor <- c( + ES = '##FORMAT=', + SE = '##FORMAT=', + LP = '##FORMAT=', + SS = '##FORMAT=', + EAF = '##FORMAT=') + meta <- c("##fileformat=VCFv4.2", unname(hdrFor[tags]), "##contig=", + paste(c("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", + "INFO", "FORMAT", samples), collapse = "\t")) + rows <- vapply(seq_len(n), function(i) { + vals <- c(ES = sprintf("%.4f", es[i]), SE = sprintf("%.4f", se[i]), + LP = sprintf("%.4f", lp[i]), SS = "10000", + EAF = sprintf("%.4f", eaf[i])) + cell <- paste(vals[tags], collapse = ":") + geno <- paste(rep(cell, length(samples)), collapse = "\t") + paste(c("22", pos[i], ids[i], ref[i], alt[i], ".", "PASS", ".", + paste(tags, collapse = ":"), geno), collapse = "\t") + }, character(1)) + writeLines(c(meta, rows), path) + path +} + +test_that("loadGwasSumStatsFromManifest reads a GWAS-VCF (A1=ALT, Z=ES/SE)", { + skip_if_not_installed("VariantAnnotation") + tmp <- withr::local_tempdir() + vcf <- .writeGwasVcf(file.path(tmp, "study1.vcf")) + manifest <- data.frame(study = "study1", sumStatsPath = vcf, + stringsAsFactors = FALSE) + obj <- suppressWarnings(loadGwasSumStatsFromManifest( + manifest, genome = "hg38", ldSketch = .toyLdSketch())) + mc <- S4Vectors::mcols(obj$entry[[1L]]) + expect_equal(as.character(mc$A1), c("G", "C", "G", "C", "A")) # ALT + expect_equal(as.character(mc$A2), c("A", "T", "T", "T", "G")) # REF + expect_equal(mc$Z, c(0.12, -0.04, 0.21, 0.005, -0.18) / 0.1, tolerance = 1e-3) + expect_equal(mc$N, rep(10000, 5)) +}) + +test_that("GWAS-VCF sample selection is enforced for multi-sample files", { + skip_if_not_installed("VariantAnnotation") + tmp <- withr::local_tempdir() + vcf <- .writeGwasVcf(file.path(tmp, "multi.vcf"), + samples = c("study1", "study2")) + manifest <- data.frame(study = "study1", sumStatsPath = vcf, + stringsAsFactors = FALSE) + expect_error( + suppressWarnings(loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch())), + "pass .sampleSelect") + obj <- suppressWarnings(loadGwasSumStatsFromManifest( + manifest, genome = "hg38", ldSketch = .toyLdSketch(), + sampleSelect = "study2")) + expect_equal(length(obj$entry[[1L]]), 5L) + expect_error( + suppressWarnings(loadGwasSumStatsFromManifest( + manifest, genome = "hg38", ldSketch = .toyLdSketch(), + sampleSelect = "nope")), + "not a sample") +}) + +test_that("GWAS-VCF without SS errors", { + skip_if_not_installed("VariantAnnotation") + tmp <- withr::local_tempdir() + vcf <- .writeGwasVcf(file.path(tmp, "noss.vcf"), dropSs = TRUE) + manifest <- data.frame(study = "study1", sumStatsPath = vcf, + stringsAsFactors = FALSE) + expect_error( + suppressWarnings(loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch())), + "SS") +}) + +test_that("GWAS-VCF region read: tabixed honours it, plain warns", { + skip_if_not_installed("VariantAnnotation") + tmp <- withr::local_tempdir() + vcf <- .writeGwasVcf(file.path(tmp, "study1.vcf")) + gz <- Rsamtools::bgzip(vcf, dest = paste0(vcf, ".gz"), overwrite = TRUE) + Rsamtools::indexTabix(gz, format = "vcf") + df <- suppressWarnings(pecotmr:::.readSumStatsVcf( + gz, "22:14560000-14565000", NULL, NULL, "lbl")) + expect_equal(nrow(df), 2L) + # Plain (non-indexed) VCF ignores the region with a warning. + expect_warning( + pecotmr:::.readSumStatsVcf(vcf, "22:14560000-14565000", NULL, NULL, "lbl"), + "region ignored") +}) + +test_that("GWAS-VCF without ES/SE errors", { + skip_if_not_installed("VariantAnnotation") + tmp <- withr::local_tempdir() + vcf <- .writeGwasVcf(file.path(tmp, "noes.vcf"), dropEs = TRUE) + manifest <- data.frame(study = "study1", sumStatsPath = vcf, + stringsAsFactors = FALSE) + expect_error( + suppressWarnings(loadGwasSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch())), + "ES and SE") +}) + +test_that("GWAS-VCF region on an absent chromosome yields no rows", { + skip_if_not_installed("VariantAnnotation") + tmp <- withr::local_tempdir() + vcf <- .writeGwasVcf(file.path(tmp, "s.vcf")) + gz <- Rsamtools::bgzip(vcf, dest = paste0(vcf, ".gz"), overwrite = TRUE) + Rsamtools::indexTabix(gz, format = "vcf") + df <- suppressWarnings( + pecotmr:::.readSumStatsVcf(gz, "99:1-100", NULL, NULL, "lbl")) + expect_equal(nrow(df), 0L) +}) + +# =========================================================================== +# Coverage: remaining error/edge branches +# =========================================================================== + +test_that("optional sumstats columns (BETA/SE/P/MAF) are attached", { + tmp <- withr::local_tempdir() + df <- .toyGwasDf(5) + df$beta <- df$z * 0.1; df$se <- rep(0.1, 5) + df$p <- rep(0.05, 5); df$maf <- rep(0.2, 5) + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "s.tsv")) + obj <- loadGwasSumStatsFromManifest( + data.frame(study = "s1", sumStatsPath = ssPath, stringsAsFactors = FALSE), + genome = "hg38", ldSketch = .toyLdSketch()) + mc <- S4Vectors::mcols(obj$entry[[1L]]) + expect_true(all(c("BETA", "SE", "P", "MAF") %in% colnames(mc))) +}) + +test_that(".detectGenotypeFormat probes a PLINK2 prefix", { + tmp <- withr::local_tempdir() + file.create(file.path(tmp, "g.pgen")) + expect_error(pecotmr:::.detectGenotypeFormat(file.path(tmp, "g"))) +}) + +test_that(".resolveLdSketch passes through NULL and a handle", { + expect_null(pecotmr:::.resolveLdSketch(NULL)) + h <- .toyLdSketch() + expect_identical(pecotmr:::.resolveLdSketch(h), h) +}) + +test_that("non-constant ldSketchPath column errors", { + tmp <- withr::local_tempdir() + ss1 <- .writeSumstatsTsv(.toyGwasDf(3), file.path(tmp, "a.tsv")) + ss2 <- .writeSumstatsTsv(.toyGwasDf(3), file.path(tmp, "b.tsv")) + manifest <- data.frame(study = c("a", "b"), sumStatsPath = c(ss1, ss2), + ldSketchPath = c(.toyRefPrefix(), "other"), + stringsAsFactors = FALSE) + expect_error(loadGwasSumStatsFromManifest(manifest, genome = "hg38"), + "must be constant") +}) + +test_that("tabix region with absent chrom or empty range returns no rows", { + tmp <- withr::local_tempdir() + tsv <- .writeSumstatsTsv(.toyGwasDf(5), file.path(tmp, "s.tsv")) + gz <- Rsamtools::bgzip(tsv, dest = paste0(tsv, ".gz"), overwrite = TRUE) + Rsamtools::indexTabix(gz, seq = 1L, start = 2L, end = 2L, comment = "#") + expect_equal(nrow(pecotmr:::.readSumStatsText(gz, "99:1-100", NULL, "l")), 0L) + expect_equal(nrow(pecotmr:::.readSumStatsText(gz, "22:1-2", NULL, "l")), 0L) +}) + +test_that("phenotype BED missing required columns errors", { + tmp <- withr::local_tempdir() + bad <- file.path(tmp, "bad.bed") + readr::write_tsv( + data.frame(`#chr` = "chr22", start = 1L, end = 2L, check.names = FALSE), bad) + manifest <- data.frame(context = "ctx", phenotypePath = bad, study = "s", + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + expect_error(loadQtlDatasetFromManifest(manifest), "missing one of") +}) + +test_that("covariate file with no shared samples errors", { + tmp <- withr::local_tempdir() + bed <- .writePhenoBed(file.path(tmp, "c.bed"), paste0("S", 1:6)) + cov <- .writeCovTsv(file.path(tmp, "cov.tsv"), paste0("Z", 1:6)) + manifest <- data.frame(context = "ctx", phenotypePath = bed, + covariatePath = cov, study = "s", + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + expect_error(loadQtlDatasetFromManifest(manifest), "No shared samples") +}) + +test_that("genotypeCovariates accepts a path string", { + tmp <- withr::local_tempdir() + samples <- paste0("S", 1:6) + bed <- .writePhenoBed(file.path(tmp, "c.bed"), samples) + gcov <- .writeCovTsv(file.path(tmp, "g.tsv"), samples, npc = 2) + manifest <- data.frame(context = "ctx", phenotypePath = bed, study = "s", + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + qd <- loadQtlDatasetFromManifest(manifest, genotypeCovariates = gcov) + expect_equal(ncol(getGenotypeCovariates(qd)), 2L) +}) + +test_that("QtlSumStats manifest resolves a per-row columnMapping", { + tmp <- withr::local_tempdir() + df <- .toyGwasDf(5) + names(df) <- c("CHR", "POS", "RSID", "EA", "OA", "ZSCORE", "SAMPLE_N") + ssPath <- .writeSumstatsTsv(df, file.path(tmp, "g.tsv")) + mapPath <- file.path(tmp, "m.yaml") + yaml::write_yaml(list(chrom = "CHR", pos = "POS", variant_id = "RSID", + A1 = "EA", A2 = "OA", z = "ZSCORE", + n_sample = "SAMPLE_N"), mapPath) + manifest <- data.frame(study = "eqtl", context = "Blood", trait = "geneA", + sumStatsPath = ssPath, columnMapping = mapPath, + stringsAsFactors = FALSE) + obj <- loadQtlSumStatsFromManifest(manifest, genome = "hg38", + ldSketch = .toyLdSketch()) + expect_equal(length(obj$entry[[1L]]), 5L) +}) + +test_that("QtlDataset builder errors on conflicting per-context/study paths", { + tmp <- withr::local_tempdir() + samples <- paste0("S", 1:6) + bed <- .writePhenoBed(file.path(tmp, "c.bed"), samples) + cov1 <- .writeCovTsv(file.path(tmp, "cov1.tsv"), samples) + cov2 <- .writeCovTsv(file.path(tmp, "cov2.tsv"), samples) + # Two rows, same context + phenotype, conflicting covariatePath. + m1 <- data.frame(context = c("ctx", "ctx"), phenotypePath = c(bed, bed), + covariatePath = c(cov1, cov2), study = "s", + genotypePath = .toyRefPrefix(), stringsAsFactors = FALSE) + expect_error(loadQtlDatasetFromManifest(m1), "multiple covariatePath") + # Conflicting genotypeCovariatePath. + m2 <- data.frame(context = c("ctx", "ctx"), phenotypePath = c(bed, bed), + study = "s", genotypePath = .toyRefPrefix(), + genotypeCovariatePath = c(cov1, cov2), stringsAsFactors = FALSE) + expect_error(loadQtlDatasetFromManifest(m2), "multiple genotypeCovariatePath") +}) + +test_that("MultiStudy builder errors on conflicting per-study genotypePath", { + tmp <- withr::local_tempdir() + bedA <- .writePhenoBed(file.path(tmp, "a.bed")) + bedB <- .writePhenoBed(file.path(tmp, "b.bed")) + man <- data.frame(study = c("study1", "study1"), context = c("c1", "c2"), + phenotypePath = c(bedA, bedB), + genotypePath = c(.toyRefPrefix(), "other"), + stringsAsFactors = FALSE) + expect_error(loadMultiStudyQtlDatasetFromManifest(man), + "exactly one genotypePath") +}) diff --git a/tests/testthat/test_twasWeightsPipeline.R b/tests/testthat/test_twasWeightsPipeline.R index d0cd5a78..04765e5f 100644 --- a/tests/testthat/test_twasWeightsPipeline.R +++ b/tests/testthat/test_twasWeightsPipeline.R @@ -2315,19 +2315,17 @@ test_that("twasWeightsPipeline(QtlSumStats): multivariate mr.mash returns a colu expect_equal(nrow(res), 2L) # one row per context }) -test_that("twasWeightsPipeline(QtlSumStats): mvsusie with no matching FM fit warns and skips", { +test_that("twasWeightsPipeline(QtlSumStats): fine-mapping method absent from fineMappingResult errors", { ss <- .tp_makeQtlSumStats(n_entries = 2L) - # FineMappingResult that does NOT contain an mvsusie fit for (s1, t1). + # FineMappingResult that does NOT contain an mvsusie fit (only susie). fmr <- .tp_makeStubFineMappingResult(study = "s1", contexts = "c1", traits = "t1", method = "susie") - do.call(local_mocked_bindings, - c(list(extractBlockGenotypes = .tp_mockExtractor()), - .tp_mockSumstatWeights(), list(.package = "pecotmr"))) + # The gate rejects the request before any fitting: a requested fine-mapping + # method must be present in the supplied fineMappingResult. expect_error( - suppressWarnings(suppressMessages( - twasWeightsPipeline(ss, methods = "mvsusie", fineMappingResult = fmr, - verbose = 0))), - "no entries produced weights") + twasWeightsPipeline(ss, methods = "mvsusie", fineMappingResult = fmr, + verbose = 0), + "contains no such fine-mapping fit") }) # ----------------------------------------------------------------------------- diff --git a/vignettes/coloc-pipeline.Rmd b/vignettes/coloc-pipeline.Rmd index c1e2666f..6e78d579 100644 --- a/vignettes/coloc-pipeline.Rmd +++ b/vignettes/coloc-pipeline.Rmd @@ -24,10 +24,10 @@ across many QTL outcomes in one model, use Inputs: -- **QTL side** — a `QtlFineMappingResult` produced by +- **QTL input** — a `QtlFineMappingResult` produced by [`fineMappingPipeline`](fine-mapping.html). -- **GWAS side** — either a `GwasFineMappingResult` (used directly) or a - `GwasSumStats` (the pipeline fine-maps it inline first). +- **GWAS input** — either a `GwasSumStats` (fine-mapping run internally) or a +`GwasFineMappingResult` (used directly). The output is a long-format data frame with one row per `(QTL tuple, GWAS tuple, credible-set pair)` combination and the @@ -41,8 +41,6 @@ library(pecotmr) ```{r bundled} data(qtl_dataset_example, gwas_sumstats_s4_example) -qtl_dataset_example <- fixupExampleGenotypePaths(qtl_dataset_example) -gwas_sumstats_s4_example <- fixupExampleGenotypePaths(gwas_sumstats_s4_example) # Build a QtlFineMappingResult from the bundled QtlDataset. qtlFmr <- fineMappingPipeline(qtl_dataset_example, methods = "susie", @@ -51,9 +49,9 @@ qtlFmr <- fineMappingPipeline(qtl_dataset_example, methods = "susie", ## QTL FMR + GWAS sumstats (inline GWAS fine-mapping) -The simplest call: hand `colocPipeline()` the QTL fine-mapping result -and a GWAS `SumStats`. The pipeline runs `fineMappingPipeline()` -internally on the GWAS sumstats (with `finemappingMethods`), then runs +The recommend approach is to provide `colocPipeline()` a `QtlFineMappingResult` +and a `GwasSumStats`. The pipeline runs `fineMappingPipeline()` +internally on the GWAS sumstats (with `finemappingMethods`) over the exact window used in the `QtlFineMappingResult`, then runs `coloc.bf_bf` per (QTL, GWAS) pair. ```{r coloc-quickstart} @@ -75,11 +73,11 @@ res2 <- colocPipeline( attr(res2, "gwasFineMapping") ``` +If the QTL `FineMappingResult` has an LD sketch because it was run on summary statistics it must match the LD sketch for the `GWASSumStats`. + ## QTL FMR + pre-computed GWAS FMR -If you already have a `GwasFineMappingResult` (e.g. computed once and -reused), pass it directly — the pipeline skips the inline fine-mapping -step: +If you already have a `GwasFineMappingResult` (typically computed over LD blocks), you can use that fine-mapping result directly: ```{r coloc-precomputed-fmr, eval=FALSE} gwasFmr <- fineMappingPipeline(gwas_sumstats_s4_example, methods = "susie") @@ -87,28 +85,22 @@ colocPipeline(qtlFineMappingResult = qtlFmr, gwasInput = gwasFmr) ``` +The disadvantage of this approach is that LD blocks will not align with cis-QTL windows. Ultimately this means that the variable +selection in the QTL and GWAS fine-mapping will have been done over different variant sets, which may reduce the sensitivity for +colocalizations, particuarly in areas with complex LD. + ## Filtering effects fed into coloc -`coloc::coloc.bf_bf` accepts one LBF matrix per side. By default, every +`coloc::coloc.bf_bf` accepts one LBF matrix for the QTL and one for hte GWAS. By default, every effect with non-trivial prior variance (`V > priorTol`) gets included. -Two opt-in filters tighten this: - -```{r coloc-filters-table, echo=FALSE} -df <- data.frame( - Knob = c("`filterLbfCs`", "`filterLbfCsSecondary`", - "`filterLbfCsConcentration`", "`priorTol`"), - When = c("`filterLbfCsSecondary` is NULL", - "`filterLbfCsSecondary` is non-NULL (overrides `filterLbfCs`)", - "only with `filterLbfCsSecondary`", - "default filter (no `filterLbfCs*`)"), - What = c("TRUE keeps only effects that produced a credible set (`trimmedFit$sets$cs_index`).", - "Numeric coverage in (0,1). Runs a CS-concentration filter at that coverage.", - "Numeric in (0,1). Concentration factor: a CS at coverage X is kept only if it spans fewer than `nVariants * X * concentration` variants.", - "Drop effects whose estimated prior variance is at or below this threshold."), - Default = c("FALSE", "NULL", "0.5", "1e-9"), - stringsAsFactors = FALSE) -knitr::kable(df) -``` +Two opt-in filters can adjust this: + +| Filter | Description | Default | +|---|---|---| +| `filterLbfCs` | TRUE keeps only effects that produced a credible set (`trimmedFit$sets$cs_index`). | FALSE | +| `filterLbfCsSecondary` | Numeric coverage in (0,1). Runs a CS-concentration filter at that coverage. | NULL | +| `filterLbfCsConcentration` | Numeric in (0,1). Concentration factor: a CS at coverage X is kept only if it spans fewer than `nVariants * X * concentration` variants. | 0.5 | +| `priorTol` | Drop effects whose estimated prior variance is at or below this threshold. | 1e-9 | ```{r coloc-filter-example, eval=FALSE} colocPipeline( @@ -150,7 +142,7 @@ colocPipeline( When the QTL and GWAS fine-mapping results were fit on different variant sets (common when the user declines RAISS imputation, or when the GWAS FMR carries variants the QTL FMR doesn't), `adjustPips = TRUE` -(default) renormalizes each side's PIPs to the cross-FMR variant +(default) renormalizes the PIPs for each objectl to the cross-FMR variant intersection before any per-pair inference. Pass `FALSE` to use the FMRs as supplied. @@ -161,13 +153,6 @@ colocPipeline( adjustPips = TRUE) # default ``` -## LD-sketch consistency - -The QTL `FineMappingResult` and the GWAS `SumStats`/FMR must share an -LD sketch (or the QTL side's `ldSketch` must be `NULL` — indicating an -individual-level fit, in which case the check is skipped on the QTL -side). Mismatch is a hard error. - ## Common parameters ```{r coloc-params, eval=FALSE} @@ -199,11 +184,3 @@ Each row is one CS-pair combination; in particular the `idx1` / `idx2` columns identify the QTL / GWAS effect indices within the respective LBF matrices, and `PP.H4.abf` is the standard coloc shared-signal posterior. - -## Next steps - -- For multi-trait colocalization across many QTL outcomes - simultaneously, see - [`colocboostPipeline`](colocboost-pipeline.html). -- Couple coloc evidence with per-gene TWAS Z + MR via - [`causalInferencePipeline`](twas-zscore.html). diff --git a/vignettes/colocboost-pipeline.Rmd b/vignettes/colocboost-pipeline.Rmd index fa68869a..f62d70aa 100644 --- a/vignettes/colocboost-pipeline.Rmd +++ b/vignettes/colocboost-pipeline.Rmd @@ -18,8 +18,8 @@ knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, `colocboostPipeline()` runs the [colocboost](https://github.com/StatFunGen/colocboost) multi-trait variable-selection model on QTL data — optionally with one -or more GWAS studies layered on top. It dispatches on the QTL input -class: +or more GWAS studies layered on top. It accepts inputs from the QTL input +classes: | QTL input | What it is | |---|---| @@ -43,30 +43,21 @@ data(qtl_dataset_example, qtl_sumstats_example, gwas_sumstats_s4_example, multi_study_qtl_dataset_example) -qtl_dataset_example <- fixupExampleGenotypePaths(qtl_dataset_example) -qtl_sumstats_example <- fixupExampleGenotypePaths(qtl_sumstats_example) -gwas_sumstats_s4_example <- fixupExampleGenotypePaths(gwas_sumstats_s4_example) -multi_study_qtl_dataset_example <- fixupExampleGenotypePaths(multi_study_qtl_dataset_example) ``` ## Three analysis variants ColocBoost can run three different model flavours in one call. -```{r cb-flags-table, echo=FALSE} -df <- data.frame( - Flag = c("`xqtlColoc`", "`jointGwas`", "`separateGwas`"), - Default = c("TRUE", "FALSE", "FALSE"), - What = c("QTL-only colocboost across the QTL contexts.", - "Single non-focal colocboost model that stacks all QTL contexts/studies with every supplied GWAS study.", - "One focal colocboost model per GWAS study (the GWAS is the focal outcome)."), - stringsAsFactors = FALSE) -knitr::kable(df) -``` +| Flag | Default | What | +|---|---|---| +| `xqtlColoc` | TRUE | QTL-only colocboost across the QTL contexts. | +| `jointGwas` | FALSE | Single non-focal colocboost model that stacks all QTL contexts/studies with every supplied GWAS study. | +| `separateGwas` | FALSE | One focal colocboost model per GWAS study (the GWAS is the focal outcome). | ## xQTL-only colocalization -The default — run colocboost over the contexts in the QTL data with no +The default call runs colocboost over the contexts in the QTL data with no GWAS layer. ```{r xqtl-only} @@ -167,7 +158,7 @@ the `QtlDataset` (`keepSamples`, `keepVariants`, the per-context `SummarizedExperiment` row selection); the pipeline parameters above are the per-call refinements. -## QC contract +## QC requirements - **Individual-level (`QtlDataset` / `MultiStudyQtlDataset`)**: QC (MAF / MAC / X-variance / per-sample missingness, sample / variant @@ -178,11 +169,3 @@ are the per-call refinements. variant filters, panel harmonization, LD-mismatch detection, and RAISS imputation live in [`summaryStatsQc()`](rss-qc.html). The pipeline rejects inputs whose `getQcInfo()` is empty. - -## Next steps - -- For pairwise QTL ↔ GWAS coloc with the classic - `coloc.bf_bf` model (instead of multi-trait ColocBoost), see - [`colocPipeline`](coloc-pipeline.html). -- To compose colocboost results with a per-gene TWAS Z + MR analysis, - see [`causalInferencePipeline`](twas-zscore.html). diff --git a/vignettes/constructing-qtl-datasets.Rmd b/vignettes/constructing-qtl-datasets.Rmd new file mode 100644 index 00000000..060cbe93 --- /dev/null +++ b/vignettes/constructing-qtl-datasets.Rmd @@ -0,0 +1,237 @@ +--- +title: "Building QtlDataset and MultiStudyQtlDataset objects" +author: "pecotmr authors" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Building QtlDataset and MultiStudyQtlDataset objects} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, + fig.width = 7, fig.height = 4.5) +``` + +## Overview + +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 `QtlDataset`s (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 `SummarizedExperiment`s**, 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.frame` + or 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 `SummarizedExperiment`s and the + `GenotypeHandle` yourself and call the `QtlDataset()` constructor. Useful when + your data is already in memory. + +```{r load-packages} +library(pecotmr) +``` + +## 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. + +```{r toy-data} +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. + +```{r qtl-from-manifest} +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)) + +qd <- loadQtlDatasetFromManifest(manifest) +qd +``` + +The manifest can equally be a path to a `.csv` (read as CSV) or any other +extension (read as TSV): + +```{r qtl-manifest-file} +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: + +```{r qtl-manifest-qc} +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 +`SummarizedExperiment`s and a `GenotypeHandle`. + +```{r qtl-manual} +# 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 +``` + +This 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). + +```{r multistudy-manifest} +multiManifest <- data.frame( + study = c("StudyA", "StudyB"), + context = c("Muscle", "Muscle"), + phenotypePath = c(muscleBed, bloodBed), + genotypePath = genoPrefix, + stringsAsFactors = FALSE) + +msd <- loadMultiStudyQtlDatasetFromManifest(multiManifest) +msd +``` + +A MultiStudyQtlDataset can be manually constructed by building a named list of `QtlDataset`s and passing it to the constructor: + +```{r multistudy-manual} +msdManual <- MultiStudyQtlDataset( + qtlDatasets = list(StudyA = qd, StudyB = qdManual)) +names(msdManual@qtlDatasets) +``` + +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: + +```{r accessors} +getStudy(qd) +getContexts(qd) +# Cis genotypes around a trait (+/- a window), QC applied on the fly: +dim(getGenotypes(qd, traitId = "GENE1", cisWindow = 1e6)) +# Phenotype matrix for one context: +dim(SummarizedExperiment::assay(getPhenotypes(qd, contexts = "Muscle"))) +``` + +`QtlDataset` 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 session-info} +sessionInfo() +``` diff --git a/vignettes/constructing-sumstats.Rmd b/vignettes/constructing-sumstats.Rmd new file mode 100644 index 00000000..b61cbdec --- /dev/null +++ b/vignettes/constructing-sumstats.Rmd @@ -0,0 +1,294 @@ +--- +title: "Building QtlSumStats and GwasSumStats objects" +author: "pecotmr authors" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Building QtlSumStats and GwasSumStats objects} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, + fig.width = 7, fig.height = 4.5) +hasVA <- requireNamespace("VariantAnnotation", quietly = TRUE) +``` + +## Overview + +`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"`). +* `ldSketch` — a `GenotypeHandle` giving the LD reference used for QC and downstream analyses. + +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 `entry` GRanges yourself and call the + `GwasSumStats()` / `QtlSumStats()` constructor. + +> **The loaders run no QC.** They build the raw collection with an empty +> `qcInfo` slot. 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 whose `qcInfo` is empty. + +```{r load-packages} +library(pecotmr) +``` + +## 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. + +```{r toy-data} +genoPrefix <- sub("\\.bed$", "", + system.file("extdata", "toy_ref.bed", package = "pecotmr")) +ldSketch <- GenotypeHandle(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 `ldSketch` are +collection-level, so they are supplied as arguments (or as constant +`genome` / `ldSketchPath` 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`, `ldSketchPath` | single-valued | may instead be arguments | + +```{r gwas-from-manifest} +gwasManifest <- data.frame( + study = "MyGWAS", + sumStatsPath = gwasTsv, + stringsAsFactors = FALSE) + +gws <- loadGwasSumStatsFromManifest( + gwasManifest, genome = "GRCh38", ldSketch = ldSketch) +gws +S4Vectors::mcols(gws$entry[[1]]) +``` + +The `ldSketch` argument also accepts a path/prefix or a chromosome-sharded +`genoMeta` specification (see `?GenotypeHandle`), so a genome-wide reference can +be supplied as a `#chr,path` meta file rather than a single handle. + +## 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 `ldSketch`. + +```{r gwas-manual} +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( + study = "MyGWAS", + entry = list(entry), + genome = "GRCh38", + ldSketch = ldSketch) +gwsManual +``` + +For 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`, `ldSketchPath` | single-valued | may instead be arguments | + +```{r qtl-from-manifest} +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", ldSketch = ldSketch) +qts +``` + +## QtlSumStats manually + +```{r qtl-manual} +qtsManual <- QtlSumStats( + study = "MyQTL", + context = "Blood", + trait = "ENSG00000001", + entry = list(entry), + genome = "GRCh38", + ldSketch = ldSketch) +qtsManual +``` + +For 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. + +## 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`). + +```{r column-mapping} +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", ldSketch = ldSketch) +length(gwsMapped$entry[[1]]) +``` + +### Region-restricted reads (tabix) + +If a `region` is supplied and the sumstats file is **compressed with bgzip and +indexed with tabix** (a `.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. + +```{r tabix-region} +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 = "#") + +gwsRegion <- loadGwasSumStatsFromManifest( + data.frame(study = "MyGWAS", sumStatsPath = gz), + genome = "GRCh38", ldSketch = ldSketch, + region = "22:14560000-14565000") +length(gwsRegion$entry[[1]]) # only the two variants in the region +``` + +### GWAS-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. + +```{r gwas-vcf, eval=hasVA} +vcfLines <- c( + "##fileformat=VCFv4.2", + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + "##contig=", + 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", ldSketch = ldSketch) +S4Vectors::mcols(gwsVcf$entry[[1]])[, c("SNP", "A1", "A2", "Z", "N")] +``` + +For a multi-sample GWAS-VCF, pass `sampleSelect = ""` 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): + +```{r qc, eval=FALSE} +gws_qcd <- summaryStatsQc(gws, zMismatchQc = "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. + +```{r accessors} +getGenome(gws) +as.character(gws$study) +length(getSumStats(gws)) # variants in the first entry +``` + +```{r session-info} +sessionInfo() +``` diff --git a/vignettes/fine-mapping.Rmd b/vignettes/fine-mapping.Rmd index 26c882f6..048e4a3d 100644 --- a/vignettes/fine-mapping.Rmd +++ b/vignettes/fine-mapping.Rmd @@ -1,10 +1,10 @@ --- -title: "Fine-mapping with the pecotmr S4 pipeline" +title: "Fine-mapping with pecotmr" author: "pecotmr authors" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{Fine-mapping with the pecotmr S4 pipeline} + %\VignetteIndexEntry{Fine-mapping with pecotmr} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -16,17 +16,16 @@ knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, ## Overview -`fineMappingPipeline()` is a single S4 generic that dispatches on the -input class: +`fineMappingPipeline()` is a single pipeline that can accept any of the following inputs: | Input class | What it represents | Methods supported | |---|---|---| -| `QtlDataset` | Single-study individual-level data: genotypes + per-context phenotypes | susie, susieInf, susieAsh, mvsusie, fsusie | -| `MultiStudyQtlDataset` | Multiple `QtlDataset`s (optionally with an embedded `QtlSumStats`) | Same as above; runs per study | -| `QtlSumStats` | QTL summary statistics keyed by `(study, context, trait)` | susie, susieInf, susieAsh (RSS implementation auto-dispatched) | -| `GwasSumStats` | GWAS summary statistics keyed by `study` | susie, susieInf, susieAsh (RSS implementation auto-dispatched) | +| `QtlDataset` | Single-study individual-level data | susie, susieInf, susieAsh, mvsusie, fsusie, susieSer | +| `MultiStudyQtlDataset` | Multiple `QtlDataset` objects (optionally with an embedded `QtlSumStats`) | susie, susieInf, susieAsh, mvsusie, fsusie (individual-level only), susieSer (individual-level only)| +| `QtlSumStats` | QTL summary statistics annotated by `(study, context, trait)` | susie, susieInf, susieAsh, mvsusie | +| `GwasSumStats` | GWAS summary statistics annotated by `study` | susie, susieInf, susieAsh, mvsusie | -Method tokens are the same across input classes; `fineMappingPipeline()` +Method arguments are the same across input classes; `fineMappingPipeline()` routes `methods = "susie"` to `susieR::susie` for individual-level data and to `susieR::susie_rss` for summary-statistics inputs automatically. @@ -44,12 +43,6 @@ data(qtl_dataset_example, qtl_sumstats_example, gwas_sumstats_s4_example, multi_study_qtl_dataset_example) -# The bundled objects carry relative GenotypeHandle paths; this helper -# resolves them to the installed inst/extdata location. -qtl_dataset_example <- fixupExampleGenotypePaths(qtl_dataset_example) -qtl_sumstats_example <- fixupExampleGenotypePaths(qtl_sumstats_example) -gwas_sumstats_s4_example <- fixupExampleGenotypePaths(gwas_sumstats_s4_example) -multi_study_qtl_dataset_example <- fixupExampleGenotypePaths(multi_study_qtl_dataset_example) qtl_dataset_example ``` @@ -99,7 +92,7 @@ head(getTopLoci(fmr, study = "study1", context = "brain", ## Multi-study fine-mapping -A `MultiStudyQtlDataset` runs the per-study pipeline in turn and +A `MultiStudyQtlDataset` runs the pipeline per-study and concatenates the results into a single `QtlFineMappingResult`: ```{r multistudy-fm} @@ -117,7 +110,7 @@ fmrSs ``` The `QtlSumStats` collection carries an `ldSketch` `GenotypeHandle`; the -pipeline pulls the per-block LD from there automatically. +pipeline pulls the per-region LD from there automatically. ## GWAS summary-statistics fine-mapping @@ -127,8 +120,7 @@ fmrGwas <- fineMappingPipeline(gwas_sumstats_s4_example, fmrGwas ``` -The GWAS path is essentially the QTL-sumstats path with the 1-tuple -`(study)` key instead of the 4-tuple `(study, context, trait, method)`. +The GWAS path is similar to the QTL-sumstats path but only annotates by study. ## Common parameters @@ -163,7 +155,7 @@ fineMappingPipeline( ## Joint multi-context / multi-trait fits `jointSpecification` triggers the dispatch for cross-axis joint -methods (`mvsusie` for cross-context, `fsusie` for ordered-trait): +methods (`mvsusie` for cross-context, `mvsusie` or `fsusie` for multi-trait): ```{r fm-joint, eval=FALSE} fineMappingPipeline( @@ -173,7 +165,7 @@ fineMappingPipeline( contexts = c("brain", "liver"))) ``` -See `?jointSpecification` for the spec grammar. +See `?jointSpecification` for the specification grammar. ## Reading the result @@ -199,11 +191,7 @@ Each entry carries: ## Next steps -- Hand a `QtlFineMappingResult` to - [`colocPipeline`](coloc-pipeline.html) for cross-study colocalization. -- Pass `fineMappingPipeline()` output to - [`twasWeightsPipeline`](twas-weights.html) so the susie weights are - contributed to the TWAS ensemble alongside other weight methods. -- Run the full - [`causalInferencePipeline`](twas-zscore.html) (TWAS + MR) on top of a - fine-mapping result and GWAS sumstats. +- `QtlFineMappingResult` is used by [`colocPipeline`](coloc-pipeline.html) for QTL-GWAS colocalization + and by [`twasWeightsPipeline`](twas-weights.html) so that fine-mapping weights are + contributed to the TWAS ensemble weights model alongside regularized regression methods. +- Run the full [`causalInferencePipeline`](twas-zscore.html) (TWAS + MR) with a `FineMappingResult` and `GwasSumStats`. diff --git a/vignettes/rss-qc.Rmd b/vignettes/rss-qc.Rmd index f52b020e..c589c10e 100644 --- a/vignettes/rss-qc.Rmd +++ b/vignettes/rss-qc.Rmd @@ -16,12 +16,9 @@ knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, ## Overview -Regression-with-summary-statistics (RSS) methods — SuSiE-RSS, ColocBoost, -DENTIST, lassosum, PRS-CS, SDPR — all require a summary-statistics table -(`Z` or `BETA`/`SE`) and an LD reference panel that agree about which -variant is which. Two things go wrong in practice: +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: -1. **Allele alignment.** Sumstats and LD panel may flip A1/A2, encode +1. **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. 2. **LD-summary-statistic mismatch.** Even when alleles align, individual @@ -29,116 +26,22 @@ variant is which. Two things go wrong in practice: of imputation errors, ancestry differences between the GWAS and LD panel, or per-variant sample-size differences. -In the current S4 API both stages are unified behind a single pipeline -entry point: - -| Stage | Mechanism | -|---|---| -| Variant-content QC (column standardization, indels, strand-ambiguous, MAF/INFO/N filters, p-value & effect-size sanity checks, optional dbSNP / liftover) | Delegated to `MungeSumstats::format_sumstats()` from inside `summaryStatsQc()`. | -| Pecotmr-specific QC (`skipRegion`, panel harmonization against the `ldSketch`, optional PIP screen on the harmonized variants, optional SLALOM/DENTIST LD-mismatch QC, optional RAISS imputation) | Implemented inside `summaryStatsQc()`. | - -`summaryStatsQc()` is the **only** QC entry point. It takes a -`GwasSumStats` or `QtlSumStats` (a DFrame collection that already carries -its LD reference via the `ldSketch` slot), runs the full pipeline above, -and returns a new collection of the same class with cleaned entries and -a populated `qcInfo` slot. Downstream pipelines refuse SumStats -collections where `length(getQcInfo(x)) == 0L`, so QC is mandatory before +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 and +a populated `qcInfo` slot. Downstream pipelines require the `qcInfo` to be populated, so QC is mandatory before fine-mapping, TWAS-weight learning, colocalization, etc. ```{r load-packages} library(pecotmr) ``` -## Reading summary statistics into pecotmr - -There is no `readSumstats()` or `loadRssData()` in the new API. The S4 -design is constructor-only: build an in-memory `GRanges` (one per -study/context/trait tuple) carrying the variant metadata in its `mcols`, -then pass a list of those `GRanges` to `GwasSumStats()` or -`QtlSumStats()` together with the LD reference (`ldSketch`). - -The minimum required mcols on each entry are `SNP`, `A1`, `A2`, `Z`, -`N`; optional columns are `MAF`, `INFO`, `BETA`, `SE`, `P`. The -`ldSketch` is a `GenotypeHandle` covering the same variant universe -(can be built from a VCF, GDS, PLINK1/2 file, or — new in this release -— from an LD-meta TSV row via `GenotypeHandle(ldMeta = ..., region = ...)`). - -### Worked example: GWAS sumstats from a TSV - -```{r build-gwas-sumstats, eval=FALSE} -# TODO(data): replace paths with your own; this block is illustrative. - -# Read your sumstats with MungeSumstats first to harmonize column names -# and genome build (summaryStatsQc() will re-run it with the curated knobs, -# but it's convenient to standardize before building the GRanges). -df <- MungeSumstats::format_sumstats( - path = "your_gwas.tsv.gz", - ref_genome = "GRCh38", - return_data = TRUE, - return_format = "data.frame") - -gr <- GenomicRanges::GRanges( - seqnames = paste0("chr", df$CHR), - ranges = IRanges::IRanges(start = df$BP, width = 1L)) -S4Vectors::mcols(gr) <- S4Vectors::DataFrame( - SNP = df$SNP, A1 = df$A1, A2 = df$A2, - Z = df$BETA / df$SE, - N = df$N) - -ldSketch <- GenotypeHandle(ldMeta = "ld_meta.tsv", - region = "chr1:1-1000000") - -gws <- GwasSumStats( - study = "MyStudy", - entry = list(gr), - genome = "GRCh38", - ldSketch = ldSketch) -``` - -### Worked example: QTL sumstats (tabix-indexed per-gene tables) - -```{r build-qtl-sumstats, eval=FALSE} -# TODO(data): point at your tabix-indexed per-gene QTL tables and column -# map. The pattern below works for tensorQTL / APEX nominal output once -# columns are renamed to SNP / A1 / A2 / Z / N. - -readOneGene <- function(tabix_path, gene_id, region) { - rows <- Rsamtools::scanTabix(Rsamtools::TabixFile(tabix_path), - param = GenomicRanges::GRanges(region))[[1L]] - df <- read.table(text = rows, sep = "\t", - col.names = c("chrom","pos","SNP","A1","A2","BETA","SE","N")) - df <- df[df$gene_id == gene_id, , drop = FALSE] - gr <- GenomicRanges::GRanges( - seqnames = paste0("chr", df$chrom), - ranges = IRanges::IRanges(start = df$pos, width = 1L)) - S4Vectors::mcols(gr) <- S4Vectors::DataFrame( - SNP = df$SNP, A1 = df$A1, A2 = df$A2, - Z = df$BETA / df$SE, - N = df$N) - gr -} - -ldSketch <- GenotypeHandle(ldMeta = "ld_meta.tsv", - region = "chr1:1-1000000") - -qts <- QtlSumStats( - study = "MyQTL", - context = "brain", - trait = c("ENSG00000001", "ENSG00000002"), - entry = list(readOneGene("brain.tsv.gz", "ENSG00000001", "chr1:1-1000000"), - readOneGene("brain.tsv.gz", "ENSG00000002", "chr1:1-1000000")), - genome = "GRCh38", - ldSketch = ldSketch) -``` - -The binary-trait variance formula `varY = n / (n - 1) * phi * (1 - phi)` -(with `phi = nCase / n`) is computed in user code and passed via the -constructor's `varY = c()` argument. +## Running QC -## Running the 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`. -Once the SumStats collection is built, every QC step lives in -`summaryStatsQc()`: +`summaryStatsQc()` provides many options for the user to tweak: ```{r run-qc, eval=FALSE} gws_qcd <- summaryStatsQc( @@ -156,7 +59,7 @@ gws_qcd <- summaryStatsQc( matchMinProp = 0) ``` -`gws_qcd` is a new `GwasSumStats` of the same class, with cleaned +`gws_qc` is a new `GwasSumStats` of the same class, with cleaned entries and a populated `qcInfo` slot. The per-entry audit record is accessible via `getQcInfo()`: @@ -168,15 +71,55 @@ str(qc$entryAudit[[1]]) # per-entry counts: variantsIn, mungeSumstatsDropped, qc$options # the curated knobs you passed (frozen for provenance) ``` -## Column-availability contract - A non-zero `mafCutoff` requires every entry to carry a `MAF` mcol; non-zero `infoCutoff` requires `INFO`; non-zero `nCutoff` requires `N`. -Missing column with a non-zero cutoff is a hard error — there is no -silent fallback. If you don't want a given filter to run, leave its -cutoff at 0 (its default). +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. + +## 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: + +1. `susieR::estimate_s_rss()` estimates the overall LD-mismatch scale `s` + between the harmonized z-scores and the reference LD matrix `R`. +2. `susieR::kriging_rss(z, R, n, s)` computes, for each variant, the + leave-one-out **conditional distribution** of `z_i` given all the other + z-scores under the LD model — i.e. the z-score the LD panel predicts for + that variant. +3. 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 (default `5e-8`) are dropped. The number removed is recorded in + the `qcInfo` audit as `krigingOutliersDropped`. + +Because it needs a per-region LD matrix and the kriging diagnostic, this is an +**RSS-only** check: it draws `R` from the collection's `ldSketch` 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: + +```{r kriging, eval=FALSE} +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, outlier +``` + +Kriging and SLALOM/DENTIST are complementary and can be combined: the kriging +prefilter removes allele-flip outliers up front, then `zMismatchQc` handles the +remaining LD-mismatch variants. -## Method internals: SLALOM vs DENTIST +## LD mismatch correction with SLALOM and DENTIST `summaryStatsQc(zMismatchQc = ...)` dispatches to one of two LD-mismatch detectors. They differ in how they decide whether a @@ -204,7 +147,7 @@ Both methods are available directly as `slalom()` / `dentist()` / `summaryStatsQc()` pipeline (typically to tune parameters or to inspect the per-variant intermediate output). -### When DENTIST window size matters +### 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: @@ -220,7 +163,7 @@ zMismatchQc on. Use DENTIST when you specifically want its iterative imputation output, or when you are reproducing a DENTIST-based pipeline. -### SLALOM + RAISS as the imputation alternative +### 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 @@ -238,43 +181,11 @@ gws_qcd <- summaryStatsQc( minimumLd = 5, lamb = 0.01)) ``` -## Where else RAISS is used - -QC re-imputation (above) is one of three places RAISS is wired into the -pipelines. The other two impute summary statistics *outside* the QC -context, and are documented in the TWAS vignettes: - -- **Filling GWAS variants missing from the LD reference**, so that TWAS - weight variants with LD neighbors but no observed GWAS hit are not - silently dropped at the weight-vs-GWAS intersection. See the - *Inferring TWAS Z-scores from GWAS summary statistics* vignette. -- **Filling QTL variants missing from the LD reference** before - summary-statistics TWAS weight training, so the weight learners - (lassosum-RSS, PRS-CS, SDPR, SuSiE-RSS, etc.) see a richer variant - set. See the *Learning TWAS weights* vignette. +## Next steps -In all three cases RAISS imputes **summary statistics only**, never -weights. Imputed variants with `R^2 < r2Threshold` are dropped by -RAISS's internal filter. Weight variants with no LD-reference -counterpart at all are warned about and dropped (they have no neighbors -to impute from). - -## Summary - -| Function | Use it when | -|---|---| -| `summaryStatsQc()` | The single QC entry point. Takes a `GwasSumStats` / `QtlSumStats`, runs MungeSumstats variant-content QC + pecotmr-specific QC (allele harmonization, SLALOM/DENTIST, RAISS), returns a new SumStats with `qcInfo` populated. Mandatory before fine-mapping / TWAS / colocalization. | -| `getQcInfo()` | Inspect the per-entry audit record left behind by `summaryStatsQc()`. | -| `slalom()` / `dentist()` / `dentistSingleWindow()` / `raiss()` | Direct access to the underlying methods for inspection or parameter tuning. | -| `ldMismatchQc()` | One-shot SLALOM/DENTIST flagging when you already have aligned `z` and `R` and want just the per-variant outlier flag. | - -Downstream pipelines (`fineMappingPipeline()`, `twasWeightsPipeline()`, -`colocboostPipeline()`, etc.) enforce that `length(getQcInfo(x)) > 0L` -before they accept a SumStats input, so QC is never silently skipped. -The two TWAS pipelines additionally accept `imputeMissing = TRUE` to -fill missing variants from the LD reference at the harmonization stage -(see the TWAS vignettes). +`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 session-info} sessionInfo() diff --git a/vignettes/twas-weights.Rmd b/vignettes/twas-weights.Rmd index a51c4d95..9ddd04f7 100644 --- a/vignettes/twas-weights.Rmd +++ b/vignettes/twas-weights.Rmd @@ -1,10 +1,10 @@ --- -title: "Learning TWAS weights with the pecotmr S4 pipeline" +title: "Learning TWAS weights with pecotmr" author: "pecotmr authors" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{Learning TWAS weights with the pecotmr S4 pipeline} + %\VignetteIndexEntry{Learning TWAS weights with pecotmr} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -18,20 +18,19 @@ knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, TWAS weights are per-variant coefficients that predict a molecular phenotype (gene expression, splicing, protein abundance) from genotype. -Once learned, the weights are applied to GWAS summary statistics to +Once learned, the weights are combined with GWAS summary statistics to compute per-gene TWAS Z-scores (covered in the [TWAS Z-score + causal inference vignette](twas-zscore.html)). -`twasWeightsPipeline()` is a single S4 generic that dispatches on the -input class. The output is always a `TwasWeights` collection (DFrame +`twasWeightsPipeline()` is a single that accepts multiple input types. The output is always a `TwasWeights` collection (DFrame subclass keyed by `(study, context, trait, method)` with per-tuple `TwasWeightsEntry` payloads). | Input class | Methods supported | |---|---| -| `QtlDataset` | Univariate: susie, susieInf, glmnet (lasso/enet/ridge), mr.ash, MCP, SCAD, L0Learn, Bayesian alphabet methods. Multivariate (multi-context): mr.mash, mvsusie | -| `MultiStudyQtlDataset` | Iterates the embedded studies; same methods as `QtlDataset` | -| `QtlSumStats` | RSS-only: susieRss, lassosum-RSS, PRS-CS, SDPR, mr.ash-RSS, P+T, MCP-RSS, ... | +| `QtlDataset` | Univariate: LASSO, elastic net, mr.ash, MCP, SCAD, L0Learn, Bayesian alphabet, DPR. Multivariate: mr.mash, mvsusie | +| `MultiStudyQtlDataset` | Iterates over embedded studies; same methods as `QtlDataset` | +| `QtlSumStats` | RSS-only: LASSO, mr.ash, MCP, SCAD, L0Learn, DPR, PRS-CS, P+T. Multivariate: mr.mash, mvsusie | ```{r load-packages} library(pecotmr) @@ -43,25 +42,18 @@ library(pecotmr) data(qtl_dataset_example, qtl_sumstats_example, multi_study_qtl_dataset_example) -qtl_dataset_example <- fixupExampleGenotypePaths(qtl_dataset_example) -qtl_sumstats_example <- fixupExampleGenotypePaths(qtl_sumstats_example) -multi_study_qtl_dataset_example <- fixupExampleGenotypePaths(multi_study_qtl_dataset_example) ``` ## Individual-level weights -The simplest call learns a single regression-based weight set from a -`QtlDataset`. Cross-validation is on by default (`cvFolds = 5`). Like +The simplest individual-level call learns a single regression-based weight set from a +`QtlDataset`. Cross-validation is enabled by default (`cvFolds = 5`). Like `fineMappingPipeline()`, you need to tell it which genotype window to extract — `cisWindow` (per trait) is the most common choice. Fine-mapping methods (`susie`, `susieInf`, `susieAsh`, `mvsusie`, -`fsusie`) are intentionally NOT runnable here on their own: -`twasWeightsPipeline()` will not re-fit fine-mapping models from -scratch. Run `fineMappingPipeline()` first and pass the result via -`fineMappingResult = …` (see the next section). The quickstart below -uses `lasso`, a TWAS-regression method that has no fine-mapping -dependency. +`fsusie`) are not re-fit by `twasWeightsPipeline()` - instead run fine-mapping models with `fineMappingPipeline()` first and pass the result via +`fineMappingResult = …` (see the next section). ```{r individual-quickstart} tw <- twasWeightsPipeline(qtl_dataset_example, @@ -93,13 +85,11 @@ length(getVariantIds(entry)) ## Including a fine-mapping susie fit in the TWAS ensemble -A susie/susieInf/susieAsh fine-mapping fit is itself a set of per-variant +The posterior effect estimates of SuSiE fine-mapping models are themselves a set of per-variant weight estimates. When you pass an existing `FineMappingResult` to `twasWeightsPipeline()`, the trimmed susie fit for each matching -`(study, context, trait)` tuple is added to the ensemble as a `"susie"` -weights row instead of being re-fit from scratch — so the susie-derived -predictions can be compared and combined against the other weight -methods (glmnet, mr.ash, …) without paying the susie fit cost twice. +`(study, context, trait)` tuple is added to the ensemble as an additional +weights row instead of being re-fit. ```{r individual-with-fmr} fmr <- fineMappingPipeline(qtl_dataset_example, methods = "susie", @@ -110,10 +100,6 @@ tw2 <- twasWeightsPipeline(qtl_dataset_example, fineMappingResult = fmr) ``` -The matching is on the canonical method token (`susie`, `susieInf`, -`susieAsh`) and the `(study, context, trait)` tuple; tuples without a -matching fine-mapping row fit the susie weights from scratch as normal. - ## Multi-study weights ```{r multistudy} @@ -123,13 +109,14 @@ msTw <- twasWeightsPipeline(multi_study_qtl_dataset_example, table(msTw$study) ``` +The call for iterating over the studies of a MultiStudyQtlDataset is identical to the call for a single QtlDataset. + ## Multiple weight methods at once Pass a character vector to fit several methods on the same fold split and return them as separate `(method)` rows in the collection. Fine-mapping methods (`susie`, `susieInf`, `susieAsh`, `mvsusie`) can -ride along when a `fineMappingResult` is supplied; without one, restrict -the list to the TWAS-regression methods: +only be included if their results are present in the object passed to the `fineMappingResult` argument: ```{r multimethod, eval=FALSE} # TWAS-regression only — no fineMappingResult needed. @@ -154,17 +141,13 @@ twSs <- twasWeightsPipeline(qtl_sumstats_example, methods = "lasso") twSs ``` -The RSS path uses the `ldSketch` `GenotypeHandle` carried on the -`QtlSumStats` collection to compute / fetch LD as needed. The fine- -mapping methods (`susie`, `susieInf`, `susieAsh`, `mvsusie`) follow the -same `fineMappingResult` gate as the individual-level path — call -`fineMappingPipeline()` on the `QtlSumStats` first, then pass the -result into `twasWeightsPipeline(..., fineMappingResult = ...)`. +The RSS path uses the `ldSketch` `GenotypeHandle` inside the +`QtlSumStats` collection to compute LD as needed. The same ## Multivariate (multi-context) weights When you have multiple contexts for the same trait (e.g. multi-tissue -expression), pass `methods = "mrmash"` or `"mvsusie"`: +expression), pass `methods = "mrmash"` to do a joint multivariate fit: ```{r multivariate, eval=FALSE} twasWeightsPipeline( @@ -209,8 +192,8 @@ twasWeightsPipeline( ``` - `ensemble = TRUE` adds an additional `ensembleWeights` row to the - collection — a quadprog-fitted weighted combination of the individual - methods' CV predictions, gated by `ensembleR2Threshold`. + collection — a weighted combination of the individual + methods' CV predictions filterd by `ensembleR2Threshold`. - `naAction` is forwarded to `getResidualizedPhenotypes()` (see the [fine-mapping vignette](fine-mapping.html#common-parameters) for semantics). @@ -244,8 +227,4 @@ Each entry exposes: - Combine the learned weights with a GWAS via [`causalInferencePipeline`](twas-zscore.html) to compute per-gene - TWAS Z-scores and MR estimates. -- Skip the standalone TWAS step if you want - [`colocPipeline`](coloc-pipeline.html) or - [`colocboostPipeline`](colocboost-pipeline.html) directly on the - fine-mapping result. + TWAS Z-scores. diff --git a/vignettes/twas-zscore.Rmd b/vignettes/twas-zscore.Rmd index 0214ca62..82f10e1c 100644 --- a/vignettes/twas-zscore.Rmd +++ b/vignettes/twas-zscore.Rmd @@ -1,10 +1,10 @@ --- -title: "Causal inference: TWAS Z-scores + Wald-ratio MR" +title: "Causal inference with pecotmr: TWAS Z-scores and MR" author: "pecotmr authors" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{Causal inference: TWAS Z-scores + Wald-ratio MR} + %\VignetteIndexEntry{Causal inference with pecotmr: TWAS Z-scores and MR} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -16,13 +16,12 @@ knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, ## Overview -`causalInferencePipeline()` pairs QTL-side weights and/or a QTL -fine-mapping result with one or more GWAS studies to produce, per +`causalInferencePipeline()` pairs QTL fine-mapping and/or regularized regression weights one or more GWAS studies to produce, per `(qtlStudy, context, trait, method, gwasStudy)` tuple: - a **TWAS Z-score** (and p-value), and -- when `fineMappingResult` is supplied, a **Wald-ratio Mendelian - Randomization** estimate over the QTL credible sets. +- when `fineMappingResult` is supplied, a **Mendelian + Randomization Wald-ratio** estimate over the QTL credible sets. Three valid input combinations: @@ -40,8 +39,6 @@ library(pecotmr) ```{r bundled} data(qtl_dataset_example, gwas_sumstats_s4_example) -qtl_dataset_example <- fixupExampleGenotypePaths(qtl_dataset_example) -gwas_sumstats_s4_example <- fixupExampleGenotypePaths(gwas_sumstats_s4_example) # Pre-compute the QTL-side products from the bundled QtlDataset. fmr <- fineMappingPipeline(qtl_dataset_example, methods = "susie", @@ -61,7 +58,7 @@ res <- causalInferencePipeline( head(S4Vectors::mcols(res)) ``` -The result is a `GRanges` positioned at the variant span of each QTL +The result is a `GRanges` object positioned at the variant span of each QTL weight set, with identity columns (`qtlStudy`, `context`, `trait`, `method`, `gwasStudy`) and the TWAS columns (`twasZ`, `twasPval`). MR columns are `NA` in this mode. @@ -85,7 +82,7 @@ S4Vectors::mcols(res2)[, c("twasZ", "twasPval", "waldRatio", ## MR variants: per-variant vs CS-aware IVW -| `mrMethod` | What it does | Extra knobs | +| `mrMethod` | What it does | Extra filters | |---|---|---| | `"ivwPerVariant"` (default) | Treats every `topLoci` variant with `pip > mrPipCutoff` as an IV; pools their Wald ratios via inverse-variance weighting. | `mrPipCutoff` | | `"csAware"` | Groups variants by credible set (column `cs` in `topLoci`), forms a PIP-weighted composite Wald ratio per CS (subject to per-CS cumulative-PIP cutoff `mrCpipCutoff`), then IVW-pools across CSs. Reports Cochran's Q and I-squared. | `mrCpipCutoff` | @@ -115,14 +112,14 @@ causalInferencePipeline( combineMethods = c("acat", "hmp")) ``` -Accepted method tokens for `combineMethods`: `acat`, `hmp`, +Accepted values for `combineMethods`: `acat`, `hmp`, `bonferroni`, `fisher`, `stouffer`, `invchisq`, `gbj`, `bj`, `hc`, `ghc`, `minp`, `gbj_omni`, `aspu`, `gates` (see `?combinePValues`). ## LD-sketch consistency Both `twasWeights` and `fineMappingResult` carry an `ldSketch` slot. -When they're non-`NULL`, they must reference the same LD panel as +When provided, they must reference the same LD panel as `gwasSumStats`. The pipeline errors on mismatch. A QTL input with `ldSketch = NULL` (e.g. weights learned from @@ -154,12 +151,3 @@ mc <- as.data.frame(S4Vectors::mcols(res2)) # Significant TWAS hits at Bonferroni alpha = 0.05 across rows mc[mc$twasPval < 0.05 / nrow(mc), ] ``` - -## Next steps - -- Pair TWAS hits with coloc evidence via - [`colocPipeline`](coloc-pipeline.html) or - [`colocboostPipeline`](colocboost-pipeline.html). -- Refine the input weight ensemble by re-running - [`twasWeightsPipeline`](twas-weights.html) with a richer `methods` - vector and the same `fineMappingResult`.