diff --git a/DESCRIPTION b/DESCRIPTION index 09da490c..1d69cc1d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -93,8 +93,10 @@ Collate: 'GenotypeHandle.R' 'AllClasses.R' 'AnnotationMatrix.R' - 'FineMappingEntry.R' 'tupleSelectors.R' + 'CtwasResult.R' + 'CtwasResultEntry.R' + 'FineMappingEntry.R' 'GwasFineMappingResult.R' 'H2Estimate.R' 'JointGroup.R' @@ -113,6 +115,7 @@ Collate: 'colocPipeline.R' 'colocboostPipeline.R' 'cpp11.R' + 'crossValidation.R' 'ctwasPipeline.R' 'deprecated.R' 'exampleData.R' diff --git a/NAMESPACE b/NAMESPACE index 93537aa6..feea5479 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,8 @@ S3method(postprocessFinemappingFit,susie) S3method(postprocessFinemappingFit,susieInf) S3method(postprocessFinemappingFit,susieRss) export(AnnotationMatrix) +export(CtwasResult) +export(CtwasResultEntry) export(FineMappingEntry) export(GenotypeHandle) export(GwasFineMappingResult) @@ -23,6 +25,8 @@ export(TwasWeightsEntry) export(adjustPips) export(alignVariantNames) export(alleleQc) +export(asCtwasResult) +export(asGranges) export(assembleCtwasInputs) export(autoDecision) export(bLassoWeights) @@ -43,7 +47,9 @@ export(colocPipeline) export(colocPostProcessor) export(colocWrapper) export(colocboostPipeline) +export(combineFineMappingResults) export(combinePValues) +export(combineTwasWeights) export(computeBlockLdCor) export(computeCoefficientsGlasso) export(computeCovDiag) @@ -99,12 +105,14 @@ export(getContexts) export(getCorrelation) export(getCs) export(getCtwasMetaData) +export(getCtwasParam) export(getCvFits) export(getCvResult) export(getDataType) export(getEigenList) export(getEnrichment) export(getFineMappingResult) +export(getFinemap) export(getFits) export(getFormat) export(getFrqData) @@ -138,6 +146,7 @@ export(getQcInfo) export(getQtlDatasets) export(getRefPanel) export(getRefVariantInfo) +export(getRegion) export(getResidualizedGenotypes) export(getResidualizedPhenotypes) export(getSampleIds) @@ -150,6 +159,7 @@ export(getStandardized) export(getStudy) export(getSumStats) export(getSumstatDf) +export(getSusieAlpha) export(getSusieFit) export(getSusieResult) export(getTauBlocks) @@ -233,6 +243,7 @@ export(prsCs) export(prsCsWeights) export(qtlEnrichment) export(qtlEnrichmentPipeline) +export(qtlSumStatsFromZMatrix) export(raiss) export(readAfreq) export(readAnnotations) @@ -244,6 +255,7 @@ export(regionDataToIndInput) export(regionDataToRssInput) export(regionToDf) export(regionsOverlap) +export(resolveWeights) export(rssAnalysisPipeline) export(sanitizeMashData) export(scadRssWeights) @@ -253,6 +265,7 @@ export(sdpr) export(sdprWeights) export(slalom) export(sldscPostprocessingPipeline) +export(sldscSubsetMeta) export(sliceMashData) export(standardiseSumstatsColumns) export(standardizeSldscTrait) @@ -304,12 +317,14 @@ exportMethods(getBlocks) exportMethods(getContexts) exportMethods(getCorrelation) exportMethods(getCs) +exportMethods(getCtwasParam) exportMethods(getCvFits) exportMethods(getCvResult) exportMethods(getDataType) exportMethods(getEigenList) exportMethods(getEnrichment) exportMethods(getFineMappingResult) +exportMethods(getFinemap) exportMethods(getFits) exportMethods(getFormat) exportMethods(getFrqData) @@ -342,6 +357,7 @@ exportMethods(getQcDiagnostics) exportMethods(getQcInfo) exportMethods(getQtlDatasets) exportMethods(getRefPanel) +exportMethods(getRegion) exportMethods(getResidualizedGenotypes) exportMethods(getResidualizedPhenotypes) exportMethods(getSampleIds) @@ -354,6 +370,7 @@ exportMethods(getStandardized) exportMethods(getStudy) exportMethods(getSumStats) exportMethods(getSumstatDf) +exportMethods(getSusieAlpha) exportMethods(getSusieFit) exportMethods(getTauBlocks) exportMethods(getTopLoci) @@ -371,6 +388,7 @@ exportMethods(hasGenotypes) exportMethods(nSnps) exportMethods(readAnnotations) exportMethods(readGenotypes) +exportMethods(resolveWeights) exportMethods(show) exportMethods(subsetChr) exportMethods(twasWeightsPipeline) @@ -444,6 +462,7 @@ importFrom(stats,complete.cases) importFrom(stats,cor) importFrom(stats,cutree) importFrom(stats,hclust) +importFrom(stats,lm) importFrom(stats,median) importFrom(stats,na.omit) importFrom(stats,optim) diff --git a/R/AllClasses.R b/R/AllClasses.R index d756316d..37e2d2ee 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -178,9 +178,12 @@ setGeneric(".fmrSelectEntry", function(x, ...) standardGeneric(".fmrSelectEntry" setMethod("getCs", "FineMappingResultBase", function(x, study = NULL, context = NULL, trait = NULL, method = NULL, region = NULL, coverage = 0.95, ...) { - getCs(.fmrSelectEntry(x, study = study, context = context, trait = trait, - method = method, region = region), - coverage = coverage) + # Selectors pinning one entry -> that entry's bare credible-set table; no / + # partial selectors -> aggregate every matching entry's credible sets, + # prefixed with the row identity (study/context/trait/region_id/method). + .fmrAggregateView(x, study = study, context = context, trait = trait, + method = method, region = region, + perEntry = function(e) getCs(e, coverage = coverage)) }) #' @rdname getTopLoci @@ -189,9 +192,27 @@ setMethod("getTopLoci", "FineMappingResultBase", function(x, type = c("data.frame", "GRanges"), signalCutoff = 0.025, study = NULL, context = NULL, trait = NULL, method = NULL, region = NULL, ...) { - getTopLoci(.fmrSelectEntry(x, study = study, context = context, - trait = trait, method = method, region = region), - type = match.arg(type), signalCutoff = signalCutoff) + type <- match.arg(type) + # type = "GRanges" can only be honored for a single pinned entry; the + # aggregate (identity-prefixed) form is data.frame-only. + if (type == "GRanges") { + sel <- tryCatch( + .fmrSelectEntry(x, study = study, context = context, trait = trait, + method = method, region = region), + error = function(e) e) + if (inherits(sel, "error")) { + stop("getTopLoci: aggregating across multiple entries requires ", + "type = 'data.frame'.") + } + return(getTopLoci(sel, type = "GRanges", signalCutoff = signalCutoff)) + } + # data.frame: bare per-variant table when selectors pin one entry, else the + # per-variant tables of every matching row stacked with row identity + # (study / context / trait / region_id / method) columns. + .fmrAggregateView(x, study = study, context = context, trait = trait, + method = method, region = region, + perEntry = function(e) getTopLoci(e, type = "data.frame", + signalCutoff = signalCutoff)) }) #' @rdname getMarginalEffects @@ -200,11 +221,18 @@ setMethod("getMarginalEffects", "FineMappingResultBase", function(x, maxPval = NULL, study = NULL, context = NULL, trait = NULL, method = NULL, region = NULL, ...) { - getMarginalEffects(.fmrSelectEntry(x, study = study, context = context, - trait = trait, method = method, - region = region), maxPval = maxPval) + # Selectors pinning one entry -> that entry's bare marginal table; no / + # partial selectors -> aggregate every matching entry's marginals, prefixed + # with the row identity (study/context/trait/region_id/method). + .fmrAggregateView(x, study = study, context = context, trait = trait, + method = method, region = region, + perEntry = function(e) getMarginalEffects(e, maxPval = maxPval)) }) +#' @rdname getRegion +#' @export +setMethod("getRegion", "FineMappingResultBase", function(x, ...) .getRegionColumn(x)) + #' @rdname getSusieFit #' @export setMethod("getSusieFit", "FineMappingResultBase", diff --git a/R/AllGenerics.R b/R/AllGenerics.R index e435c458..55c6091f 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -486,6 +486,50 @@ setGeneric("getCs", function(x, ...) standardGeneric("getCs")) #' @export setGeneric("getWeights", function(x, ...) standardGeneric("getWeights")) +#' @title Resolve Per-Variant Weights From a Weight Source +#' @description Return an aligned \code{(variantIds, weights)} pair from a single +#' weight-source entry, so cTWAS and \code{\link{causalInferencePipeline}} +#' extract weights identically whether the source is a +#' \code{TwasWeightsEntry} (its learned weight vector) or a +#' \code{FineMappingEntry} (its topLoci posterior effect). +#' @param x A \code{TwasWeightsEntry} or \code{FineMappingEntry}. +#' @param ... Reserved for future use. +#' @return A \code{list} with \code{variantIds} (character) and \code{weights} +#' (numeric) of equal length; both empty when no usable weights are present. +#' @export +setGeneric("resolveWeights", function(x, ...) standardGeneric("resolveWeights")) + +#' @title Get cTWAS Fine-mapping Posteriors +#' @description Return the per-gene (and per-SNP) fine-mapping posterior table +#' from a \code{\link{CtwasResultEntry}} or a \code{\link{CtwasResult}} +#' collection (aggregated across rows, tagged with run identity). +#' @param x A \code{CtwasResultEntry} or \code{CtwasResult}. +#' @param ... Class-specific selection arguments. +#' @return A \code{data.frame} of posteriors (or \code{NULL} when absent). +#' @export +setGeneric("getFinemap", function(x, ...) standardGeneric("getFinemap")) + +#' @title Get cTWAS Per-effect Susie Alpha Table +#' @description Return the full per-effect susie alpha table +#' (\code{ctwas::finemap_regions} \code{susie_alpha_res} shape) from a +#' \code{\link{CtwasResultEntry}} or a \code{\link{CtwasResult}} collection +#' (aggregated across rows, tagged with run identity). This is the fuller +#' cTWAS output retained so the raw run is reconstructable. +#' @param x A \code{CtwasResultEntry} or \code{CtwasResult}. +#' @param ... Class-specific selection arguments. +#' @return A \code{data.frame} of per-effect alphas (or \code{NULL} when absent). +#' @export +setGeneric("getSusieAlpha", function(x, ...) standardGeneric("getSusieAlpha")) + +#' @title Get cTWAS Group Prior Parameters +#' @description Return the estimated \code{group_prior} / \code{group_prior_var} +#' of a \code{\link{CtwasResultEntry}}. +#' @param x A \code{CtwasResultEntry}. +#' @param ... Reserved for future use. +#' @return The stored parameter object (typically a list), or \code{NULL}. +#' @export +setGeneric("getCtwasParam", function(x, ...) standardGeneric("getCtwasParam")) + #' @title Get Standardized Flag #' @description Check whether weights are on the standardized scale. #' @param x A \code{TwasWeightsEntry} or \code{TwasWeights}. @@ -591,6 +635,18 @@ setGeneric("getContexts", function(x) standardGeneric("getContexts")) #' @export setGeneric("getTraits", function(x) standardGeneric("getTraits")) +#' @title Get Per-Row Genomic Regions +#' @description Return the genomic anchor of each row of a per-tuple collection +#' as a \code{GRanges} with one range per row -- the trait's own region for a +#' \code{TwasWeights}, or the fine-mapping window for a \code{FineMappingResult}. +#' This is location provenance (e.g. for cTWAS LD-block placement). +#' @param x The object. +#' @param ... Reserved for future use. +#' @return A \code{GRanges} with one range per row of \code{x}, or an empty +#' \code{GRanges} when the collection carries no region provenance. +#' @export +setGeneric("getRegion", function(x, ...) standardGeneric("getRegion")) + #' @title Get Residualized Genotypes #' @description Residualize the genotype matrix against the per-context #' phenotype covariates and the genotype covariates, optionally diff --git a/R/CtwasResult.R b/R/CtwasResult.R new file mode 100644 index 00000000..e6cd2384 --- /dev/null +++ b/R/CtwasResult.R @@ -0,0 +1,154 @@ +# ============================================================================= +# CtwasResult S4 class +# ----------------------------------------------------------------------------- +# DFrame-subclass collection of cTWAS runs, keyed by the identity tuple +# (gwasStudy, study, context, method). Each row holds a CtwasResultEntry +# payload (fine-mapping posteriors + group priors + region metadata) for one +# run. Unlike the QTL family, `trait` is NOT a key -- cTWAS is multi-gene per +# run, so genes live inside the payload. Optional joint columns (jointStudies, +# jointContexts) tag rows born from a multi-study / multi-context cTWAS run and +# participate in the uniqueness key, exactly as in the TwasWeights / +# FineMappingResult families. +# ============================================================================= + +#' @include AllGenerics.R tupleSelectors.R +NULL + +setClass("CtwasResult", + contains = "DFrame", + validity = function(object) { + errors <- character() + required <- c("gwasStudy", "study", "context", "method", "entry") + missingCols <- setdiff(required, names(object)) + if (length(missingCols) > 0L) + errors <- c(errors, paste("missing columns:", + paste(missingCols, collapse = ", "))) + if (length(errors) == 0L) { + if (length(object$entry) != nrow(object)) + errors <- c(errors, "length(entry) must equal nrow(.) for CtwasResult") + entryTypes <- vapply(object$entry, + function(e) methods::is(e, "CtwasResultEntry"), + logical(1)) + if (!all(entryTypes)) + errors <- c(errors, + "every element of the `entry` column must be a CtwasResultEntry") + jointCols <- intersect(c("jointStudies", "jointContexts"), names(object)) + for (jc in jointCols) { + if (!is.character(object[[jc]])) + errors <- c(errors, sprintf("'%s' column must be character (got %s)", + jc, class(object[[jc]])[[1L]])) + } + keyCols <- c("gwasStudy", "study", "context", "method", jointCols) + keyDf <- as.data.frame( + lapply(keyCols, function(cn) object[[cn]]), + col.names = keyCols, stringsAsFactors = FALSE) + if (anyDuplicated(keyDf)) + errors <- c(errors, + "(gwasStudy, study, context, method[, joint*]) tuple uniqueness violated") + } + if (length(errors) == 0L) TRUE else errors + } +) + +#' @title Create a CtwasResult Collection +#' @description Construct a \code{CtwasResult} DFrame-subclass collection of +#' cTWAS runs from per-tuple vectors and a list of +#' \code{\link{CtwasResultEntry}} payloads (one per run). A single-context run +#' is one row; a multi-context run emits per-context rows sharing the same +#' \code{jointContexts} set and jointly-estimated group priors. +#' @param gwasStudy Character vector of GWAS (disease) study identifiers. +#' @param study Character vector of QTL study identifiers. +#' @param context Character vector of context labels (scalar per row). +#' @param method Character vector of weight-method names. +#' @param entry List / \code{SimpleList} of \code{CtwasResultEntry} objects. +#' @param jointStudies Optional character vector (length \code{length(gwasStudy)}) +#' listing the joined studies for a multi-study run, \code{NA_character_} +#' otherwise. \code{NULL} (default) omits the column. +#' @param jointContexts Optional character vector for multi-context runs. Same +#' shape as \code{jointStudies}. +#' @return A \code{CtwasResult} object. +#' @export +CtwasResult <- function(gwasStudy, study, context, method, entry, + jointStudies = NULL, jointContexts = NULL) { + n <- length(gwasStudy) + if (length(study) != n || length(context) != n || length(method) != n || + length(entry) != n) { + stop("`gwasStudy`, `study`, `context`, `method`, and `entry` must all ", + "have the same length.") + } + cols <- list( + gwasStudy = as.character(gwasStudy), + study = as.character(study), + context = as.character(context), + method = as.character(method), + entry = S4Vectors::SimpleList(entry) + ) + for (nm in c("jointStudies", "jointContexts")) { + val <- get(nm) + if (is.null(val)) next + if (length(val) != n) + stop("`", nm, "` must have the same length as `gwasStudy`.") + cols[[nm]] <- as.character(val) + } + df <- do.call(S4Vectors::DataFrame, c(cols, list(check.names = FALSE))) + obj <- new("CtwasResult", df) + validObject(obj) + obj +} + +#' @rdname getMethodNames +#' @export +setMethod("getMethodNames", "CtwasResult", + function(x) unique(as.character(x$method))) + +#' @rdname getStudy +#' @export +setMethod("getStudy", "CtwasResult", + function(x) unique(as.character(x$study))) + +#' @rdname getContexts +#' @export +setMethod("getContexts", "CtwasResult", + function(x) unique(as.character(x$context))) + +# Aggregate a per-entry table (finemap or susieAlpha) across all rows, prefixing +# each with the row's (gwasStudy, study, context, method) run identity. +# @noRd +.ctwasAggregateRows <- function(x, getter) { + n <- nrow(x) + if (n == 0L) return(NULL) + parts <- lapply(seq_len(n), function(i) { + tb <- getter(x$entry[[i]]) + if (is.null(tb) || nrow(tb) == 0L) return(NULL) + tb <- as.data.frame(tb) + id <- data.frame(gwasStudy = as.character(x$gwasStudy)[[i]], + study = as.character(x$study)[[i]], + context = as.character(x$context)[[i]], + method = as.character(x$method)[[i]], + stringsAsFactors = FALSE) + cbind(id[rep(1L, nrow(tb)), , drop = FALSE], tb, row.names = NULL) + }) + parts <- parts[!vapply(parts, is.null, logical(1L))] + if (length(parts) == 0L) return(NULL) + .rbindAligned(parts) +} + +#' @rdname getFinemap +#' @export +setMethod("getFinemap", "CtwasResult", + function(x, ...) .ctwasAggregateRows(x, getFinemap)) + +#' @rdname getSusieAlpha +#' @export +setMethod("getSusieAlpha", "CtwasResult", + function(x, ...) .ctwasAggregateRows(x, getSusieAlpha)) + +#' @export +setMethod("show", "CtwasResult", function(object) { + cat(sprintf("CtwasResult: %d run(s)\n", nrow(object))) + if (nrow(object) > 0L) + cat(sprintf(" %d GWAS stud(y/ies), %d method(s): %s\n", + length(unique(object$gwasStudy)), + length(unique(object$method)), + paste(unique(as.character(object$method)), collapse = ", "))) +}) diff --git a/R/CtwasResultEntry.R b/R/CtwasResultEntry.R new file mode 100644 index 00000000..4167759b --- /dev/null +++ b/R/CtwasResultEntry.R @@ -0,0 +1,55 @@ +# ============================================================================= +# CtwasResultEntry S4 class +# ----------------------------------------------------------------------------- +# Per-run cTWAS payload: the per-gene (+ per-SNP) posterior table from +# ctwas::ctwas_sumstats, the jointly-estimated group prior + prior variance, +# and per-region metadata. One entry sits in every row of a CtwasResult +# collection. +# ============================================================================= + +#' @include AllGenerics.R +NULL + +setClass("CtwasResultEntry", + representation( + finemap = "ANY", # per-gene/SNP posterior summary (ctwas finemap_res) + susieAlpha = "ANY", # per-effect susie alpha table (ctwas susie_alpha_res) + param = "ANY", # group_prior / group_prior_var for this run + regionInfo = "ANY" # per-region metadata (optional) + ), + prototype = prototype(finemap = NULL, susieAlpha = NULL, param = NULL, + regionInfo = NULL)) + +#' @title Create a CtwasResultEntry +#' @description Per-run cTWAS payload wrapping the fine-mapping posterior table, +#' the full per-effect susie alpha table, the estimated group prior(s), and +#' region metadata. Held in every row of a \code{\link{CtwasResult}} +#' collection. +#' @param finemap The per-gene (and, when SNPs are retained, per-SNP) posterior +#' summary table (\code{ctwas::finemap_regions} \code{finemap_res} shape), or +#' \code{NULL}. +#' @param susieAlpha The per-effect susie alpha table +#' (\code{ctwas::finemap_regions} \code{susie_alpha_res} shape) — the fuller +#' cTWAS output retained so the raw run is reconstructable, or \code{NULL}. +#' @param param The estimated \code{group_prior} / \code{group_prior_var} for +#' this run, or \code{NULL}. +#' @param regionInfo Per-region metadata, or \code{NULL}. +#' @return A \code{CtwasResultEntry} object. +#' @export +CtwasResultEntry <- function(finemap = NULL, susieAlpha = NULL, param = NULL, + regionInfo = NULL) { + new("CtwasResultEntry", finemap = finemap, susieAlpha = susieAlpha, + param = param, regionInfo = regionInfo) +} + +#' @rdname getFinemap +#' @export +setMethod("getFinemap", "CtwasResultEntry", function(x, ...) x@finemap) + +#' @rdname getSusieAlpha +#' @export +setMethod("getSusieAlpha", "CtwasResultEntry", function(x, ...) x@susieAlpha) + +#' @rdname getCtwasParam +#' @export +setMethod("getCtwasParam", "CtwasResultEntry", function(x, ...) x@param) diff --git a/R/FineMappingEntry.R b/R/FineMappingEntry.R index 3a3f892b..1d41d89d 100644 --- a/R/FineMappingEntry.R +++ b/R/FineMappingEntry.R @@ -113,6 +113,44 @@ FineMappingEntry <- function(variantIds, susieFit, topLoci, cvResult = NULL) { setMethod("getVariantIds", "FineMappingEntry", function(x, ...) x@variantIds) +# cTWAS weight-source accessors for a FineMappingEntry. resolveWeights returns +# the topLoci posterior effect colSums(alpha * mu) (buildTopLoci), which is on +# the STANDARDIZED scale -- unlike susieWeights / the cTWAS renorm, it does NOT +# divide by X_column_scale_factors. So: +# getStandardized -> TRUE skip cTWAS's per-variant variance scaling, which +# would double-standardize the effect (RSS and +# individual fits alike -- posterior_mean is always +# colSums(alpha*mu)). +# getFits -> NULL skip cTWAS's alpha renormalisation, which recomputes +# an UNstandardized weight (÷ scale factors) that is +# inconsistent with the standardized posterior effect. +# This lands on the same scale as the susie-TwasWeights path: an unstandardized +# susie weight (mu / scale) x sqrt(var) = mu = the standardized posterior effect. +# (getSusieFit still exposes the fit for non-weight uses.) +#' @rdname getFits +#' @export +setMethod("getFits", "FineMappingEntry", function(x, ...) NULL) + +#' @rdname getStandardized +#' @export +setMethod("getStandardized", "FineMappingEntry", function(x, ...) TRUE) + +#' @rdname resolveWeights +#' @export +setMethod("resolveWeights", "FineMappingEntry", function(x, ...) { + empty <- list(variantIds = character(0), weights = numeric(0)) + # getTopLoci projects the posterior effect to `beta` (= posterior_mean); + # use it as the per-variant weight, aligned with variant_id. + tl <- getTopLoci(x) + if (is.null(tl) || nrow(tl) == 0L || + !all(c("variant_id", "beta") %in% names(tl))) return(empty) + vids <- as.character(tl$variant_id) + w <- as.numeric(tl$beta) + ok <- !is.na(vids) & !is.na(w) + if (!any(ok)) return(empty) + list(variantIds = vids[ok], weights = w[ok]) +}) + #' @rdname getSusieFit #' @export setMethod("getSusieFit", "FineMappingEntry", diff --git a/R/GenotypeHandle.R b/R/GenotypeHandle.R index 8ec96674..43baff5f 100644 --- a/R/GenotypeHandle.R +++ b/R/GenotypeHandle.R @@ -298,10 +298,36 @@ GenotypeHandle <- function(path = NULL, tryCatch(handle@chromPaths, error = function(e) character(0)) } +# Case-insensitive match of the first of `aliases` present in `cols`; falls +# back to the `default` column position when none of the aliases are present. +.metaMatchCol <- function(cols, aliases, default) { + hit <- which(tolower(cols) %in% tolower(aliases)) + if (length(hit) > 0L) hit[[1L]] else default +} + +# Enforce one payload per chromosome key: LD sketches / sharded genotypes are +# genome-wide or per-chromosome, never sub-chromosomal, so a chromosome mapping +# to two distinct payloads is an error. Returns the data.frame unchanged. +.metaCheckUniqueChrom <- function(df) { + byChr <- split(df$path, df$chrom) + multi <- names(byChr)[vapply(byChr, function(p) length(unique(p)) > 1L, + logical(1L))] + if (length(multi) > 0L) + stop("GenotypeHandle(genoMeta): chromosome(s) ", + paste(multi, collapse = ", "), + " map to multiple genotype payloads; each chromosome must map to ", + "exactly one payload (no sub-chromosomal blocks).") + df +} + # Parse the genoMeta input into a data.frame(chrom, path). Accepts either a -# path to a #chr/path meta file (whitespace- or tab-delimited, with header) -# or a named character vector (names = chromosomes). Relative payload paths -# in a meta file are resolved against the meta file's own directory. +# path to a meta file (whitespace- or tab-delimited, with header) or a named +# character vector (names = chromosomes). In a meta file the chromosome and +# path columns are matched BY NAME -- chromosome from #chr/#chrom/chr/chrom and +# path from path/payload/prefix/genotype -- so extra columns (e.g. legacy +# start/end) and non-standard column order are tolerated; the first two columns +# are used only when the header carries none of those aliases. Relative payload +# paths are resolved against the meta file's own directory. #' @keywords internal .parseChromMeta <- function(genoMeta) { isMetaFile <- is.character(genoMeta) && length(genoMeta) == 1L && @@ -313,18 +339,25 @@ GenotypeHandle <- function(path = NULL, if (ncol(meta) < 2L) stop("GenotypeHandle(genoMeta): meta file '", genoMeta, "' must have at least 2 columns (chromosome, path).") - chrom <- as.character(meta[[1L]]) - pth <- as.character(meta[[2L]]) + chromCol <- .metaMatchCol(names(meta), + c("#chr", "#chrom", "chr", "chrom"), default = 1L) + pathCol <- .metaMatchCol(names(meta), + c("path", "payload", "prefix", "genotype"), + default = 2L) + chrom <- as.character(meta[[chromCol]]) + pth <- as.character(meta[[pathCol]]) base <- dirname(normalizePath(genoMeta)) pth <- vapply(pth, function(p) { if (grepl("^(/|[A-Za-z]:)", p) || file.exists(p)) p else file.path(base, p) }, character(1), USE.NAMES = FALSE) - return(data.frame(chrom = chrom, path = pth, stringsAsFactors = FALSE)) + return(.metaCheckUniqueChrom( + data.frame(chrom = chrom, path = pth, stringsAsFactors = FALSE))) } if (is.character(genoMeta) && length(genoMeta) >= 1L && !is.null(names(genoMeta))) { - return(data.frame(chrom = names(genoMeta), path = unname(genoMeta), - stringsAsFactors = FALSE)) + return(.metaCheckUniqueChrom( + data.frame(chrom = names(genoMeta), path = unname(genoMeta), + stringsAsFactors = FALSE))) } stop("GenotypeHandle(genoMeta): expected a path to a `#chr,path` meta file ", "or a named character vector (names = chromosomes, values = paths).") diff --git a/R/GwasFineMappingResult.R b/R/GwasFineMappingResult.R index 8ba4f587..99bc5fdb 100644 --- a/R/GwasFineMappingResult.R +++ b/R/GwasFineMappingResult.R @@ -36,6 +36,7 @@ setClass("GwasFineMappingResult", if (!all(entryTypes)) errors <- c(errors, "every element of the `entry` column must be a FineMappingEntry") + errors <- c(errors, .validateRegionColumn(object)) # Extract key columns directly rather than via `object[, keyCols]`: # column-subsetting preserves the GwasFineMappingResult class while # dropping the required `entry` column, and older S4Vectors revalidates @@ -102,11 +103,18 @@ setClass("GwasFineMappingResult", #' (\code{study}, \code{method}, \code{region_id}) triple is unique. #' Supplying meaningful labels (e.g. \code{"chr22_10516173_17414263"}) #' is preferred for downstream consumers that join on region. +#' @param region Optional \code{GRanges} (length \code{length(study)}) genomic +#' anchor per row. When \code{NULL} (default) it is derived from +#' \code{region_id} (parsing \code{chr_start_end} / \code{chr:start-end}; a +#' 0-width \code{chrUn} sentinel for ids that do not encode coordinates). +#' Carried forward as provenance (e.g. for cTWAS LD-block placement); not +#' part of the identity key. #' @param ldSketch An optional \code{GenotypeHandle}. #' @return A \code{GwasFineMappingResult} object. #' @export GwasFineMappingResult <- function(study, method, entry, region_id = NULL, + region = NULL, ldSketch = NULL) { n <- length(study) if (length(method) != n || length(entry) != n) { @@ -124,6 +132,8 @@ GwasFineMappingResult <- function(study, method, entry, region_id = as.character(region_id), entry = S4Vectors::SimpleList(entry) ) + if (is.null(region)) region <- .regionFromIds(region_id) + cols <- .appendRegionCol(cols, region, n) df <- do.call(S4Vectors::DataFrame, c(cols, list(check.names = FALSE))) obj <- new("GwasFineMappingResult", df, ldSketch = ldSketch) @@ -131,6 +141,28 @@ GwasFineMappingResult <- function(study, method, entry, obj } +# Internal: parse region_id strings (chr_start_end / chr:start-end) into a +# per-row GRanges, using the canonical `asGranges` parser; a 0-width chrUn +# sentinel for ids (e.g. synthetic "region_1") that carry no coordinates. +# Built from vectors so mixed seqlevels do not emit a seqinfo-merge warning. +.regionFromIds <- function(ids) { + n <- length(ids) + chrom <- character(n); start <- integer(n); end <- integer(n) + for (i in seq_len(n)) { + g <- tryCatch( + asGranges(sub("_([0-9]+)_([0-9]+)$", ":\\1-\\2", as.character(ids[[i]]))), + error = function(e) NULL) + if (!is.null(g) && length(g) >= 1L) { + chrom[[i]] <- as.character(GenomicRanges::seqnames(g))[[1L]] + start[[i]] <- GenomicRanges::start(g)[[1L]] + end[[i]] <- GenomicRanges::end(g)[[1L]] + } else { + chrom[[i]] <- "chrUn"; start[[i]] <- 1L; end[[i]] <- 0L + } + } + GenomicRanges::GRanges(chrom, IRanges::IRanges(start = start, end = end)) +} + # Internal: return integer row indices of `x` where every (column, value) # pair in `keys` matches as.character(x[[column]]) == value. Shared # building block for all of pecotmr's tuple-keyed row selectors and diff --git a/R/QtlDataset.R b/R/QtlDataset.R index 497174e6..4d660845 100644 --- a/R/QtlDataset.R +++ b/R/QtlDataset.R @@ -285,6 +285,27 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals) isTRUE(tryCatch(x@keepIndel, error = function(e) TRUE)) } +# Internal: return a copy of a QtlDataset with the supplied filter cutoffs / +# keep-lists REPLACING the stored slot values (NULL = leave the stored value +# untouched). This lets a pipeline accept per-call filter overrides as ordinary +# arguments instead of forcing callers to mutate @slots directly (which bypasses +# the class's validity checks). Applied against a validated copy. +.qtlApplyFilterOverrides <- function(data, + mafCutoff = NULL, macCutoff = NULL, + xvarCutoff = NULL, imissCutoff = NULL, + keepIndel = NULL, + keepSamples = NULL, keepVariants = NULL) { + if (!is.null(mafCutoff)) data@mafCutoff <- as.numeric(mafCutoff) + if (!is.null(macCutoff)) data@macCutoff <- as.numeric(macCutoff) + if (!is.null(xvarCutoff)) data@xvarCutoff <- as.numeric(xvarCutoff) + if (!is.null(imissCutoff)) data@imissCutoff <- as.numeric(imissCutoff) + if (!is.null(keepIndel)) data@keepIndel <- as.logical(keepIndel) + if (!is.null(keepSamples)) data@keepSamples <- as.character(keepSamples) + if (!is.null(keepVariants)) data@keepVariants <- as.character(keepVariants) + methods::validObject(data) + data +} + # Internal: extract the panel dosage block (samples x variants) for the # requested region, narrow to the requested sample set, and apply lazy QC # (per-sample imiss filter, then per-variant max(mafCutoff, diff --git a/R/QtlFineMappingResult.R b/R/QtlFineMappingResult.R index 17761df8..c8d2ec13 100644 --- a/R/QtlFineMappingResult.R +++ b/R/QtlFineMappingResult.R @@ -33,6 +33,7 @@ setClass("QtlFineMappingResult", if (!all(entryTypes)) errors <- c(errors, "every element of the `entry` column must be a FineMappingEntry") + errors <- c(errors, .validateRegionColumn(object)) jointCols <- intersect( c("jointStudies", "jointContexts", "jointTraits"), names(object)) for (jc in jointCols) { @@ -96,6 +97,10 @@ setClass("QtlFineMappingResult", #' Same shape as \code{jointStudies}. #' @param jointTraits Optional character vector for cross-trait joints. #' Same shape as \code{jointStudies}. +#' @param region Optional \code{GRanges} (length \code{length(study)}) giving +#' the genomic anchor of each row's trait (its own coordinates). Carried +#' forward as provenance (e.g. for cTWAS LD-block placement); not part of the +#' identity key. \code{NULL} (default) omits the column. #' @param ldSketch An optional \code{GenotypeHandle} (the LD reference for #' RSS-derived fits), or \code{NULL} for individual-level fits. #' @return A \code{QtlFineMappingResult} object. @@ -104,6 +109,7 @@ QtlFineMappingResult <- function(study, context, trait, method, entry, jointStudies = NULL, jointContexts = NULL, jointTraits = NULL, + region = NULL, ldSketch = NULL) { n <- length(study) if (length(context) != n || length(trait) != n || length(method) != n || @@ -127,6 +133,7 @@ QtlFineMappingResult <- function(study, context, trait, method, entry, stop("`", pair[[1L]], "` must have the same length as `study`.") cols[[pair[[2L]]]] <- as.character(val) } + cols <- .appendRegionCol(cols, region, n) df <- do.call(S4Vectors::DataFrame, c(cols, list(check.names = FALSE))) obj <- new("QtlFineMappingResult", df, ldSketch = ldSketch) diff --git a/R/TwasWeightsEntry.R b/R/TwasWeightsEntry.R index d70c7fbd..94f99263 100644 --- a/R/TwasWeightsEntry.R +++ b/R/TwasWeightsEntry.R @@ -116,6 +116,16 @@ setMethod("getWeights", "TwasWeightsEntry", setMethod("getVariantIds", "TwasWeightsEntry", function(x, ...) x@variantIds) +#' @rdname resolveWeights +#' @export +setMethod("resolveWeights", "TwasWeightsEntry", function(x, ...) { + vids <- as.character(getVariantIds(x)) + w <- as.numeric(getWeights(x)) + if (length(vids) == 0L || length(vids) != length(w)) + return(list(variantIds = character(0), weights = numeric(0))) + list(variantIds = vids, weights = w) +}) + #' @rdname getFits #' @export setMethod("getFits", "TwasWeightsEntry", diff --git a/R/causalInferencePipeline.R b/R/causalInferencePipeline.R index cf008851..6b3734f9 100644 --- a/R/causalInferencePipeline.R +++ b/R/causalInferencePipeline.R @@ -404,9 +404,11 @@ causalInferencePipeline <- function(gwasSumStats, trait = trait, method = method))) > 0L } -# Extract the per-tuple weights vector. From TwasWeights: read the -# TwasWeightsEntry. From FineMappingResult: extract the SuSiE-style -# coefficient (betahat) from topLoci. +# Fetch the per-tuple weight-source entry (a TwasWeightsEntry from +# `twasWeights`, or a FineMappingEntry from `fineMappingResult`) and resolve it +# to an aligned (variantIds, weights) pair via the shared `resolveWeights` -- +# the SuSiE-style posterior effect for the FMR path. Returns NULL when the tuple +# is absent or has no usable weights. .cipExtractWeights <- function(twasWeights, fineMappingResult, study, context, trait, method, useFmr) { if (!useFmr) { @@ -414,28 +416,18 @@ causalInferencePipeline <- function(gwasSumStats, list(study = study, context = context, trait = trait, method = method))) == 0L) return(NULL) - twEntry <- getTwasWeights(twasWeights, study = study, context = context, - trait = trait, method = method) - vids <- getVariantIds(twEntry) - w <- as.numeric(getWeights(twEntry)) - if (length(vids) != length(w) || length(vids) == 0L) return(NULL) - return(list(variantIds = vids, weights = w)) + ent <- getTwasWeights(twasWeights, study = study, context = context, + trait = trait, method = method) + } else { + if (!.cipFmrHasTuple(fineMappingResult, study, context, trait, method)) + return(NULL) + ent <- getFineMappingResult(fineMappingResult, study = study, + context = context, trait = trait, + method = method) } - # FMR-based weights: pull from the entry's topLoci$betahat column. - if (!.cipFmrHasTuple(fineMappingResult, study, context, trait, method)) - return(NULL) - ent <- getFineMappingResult(fineMappingResult, study = study, - context = context, trait = trait, - method = method) - tl <- getTopLoci(ent) - if (is.null(tl) || nrow(tl) == 0L) return(NULL) - betaCol <- .cipTlCols(tl)$beta - if (length(betaCol) == 0L) return(NULL) - vids <- as.character(tl$variant_id) - w <- as.numeric(tl[[betaCol[[1L]]]]) - ok <- !is.na(vids) & !is.na(w) - if (sum(ok) == 0L) return(NULL) - list(variantIds = vids[ok], weights = w[ok]) + wr <- resolveWeights(ent) + if (length(wr$variantIds) == 0L) return(NULL) + wr } # Compute the per-tuple TWAS Z from a single GwasSumStats tuple's diff --git a/R/crossValidation.R b/R/crossValidation.R new file mode 100644 index 00000000..8d9c3320 --- /dev/null +++ b/R/crossValidation.R @@ -0,0 +1,203 @@ +# ============================================================================= +# Cross-validation engine (shared by twasWeightsCv + fine-mapping CV) +# ----------------------------------------------------------------------------- +# A single K-fold CV harness backs both the TWAS predictive-weight CV +# (twasWeightsCv, twasWeights.R) and the fine-mapping CV (.fmWeightsCv, +# fineMappingPipeline.R). Callers supply only a per-fold fit; partitioning, +# the parallel fold loop, prediction, aggregation, the metric block, and the +# output key format live here so the two paths cannot drift. +# ============================================================================= + +# ============================================================================= +# Shared cross-validation engine +# ============================================================================= +# One CV harness backs both twasWeightsCv() (predictive weight methods) and the +# fine-mapping CV (.fmWeightsCv, susie/mvsusie/fsusie). Callers differ only in +# the per-fold fit; partitioning, the parallel fold loop, prediction, fold +# aggregation, the metric block, and the output key format live here. + +# Sample/Fold partition: shuffle samples, then cut into `fold` contiguous blocks. +# @noRd +.cvSamplePartition <- function(sampleNames, fold) { + idx <- sample(length(sampleNames)) + folds <- cut(seq_along(sampleNames), breaks = fold, labels = FALSE) + data.frame(Sample = sampleNames[idx], Fold = folds, stringsAsFactors = FALSE) +} + +# Canonical CV output key: "_predicted" / "_performance". +# @noRd +.cvOutputKey <- function(methodKey, suffix) paste0(methodKey, "_", suffix) + +# One CV metric row (corr, rsq, adj_rsq, pval, RMSE, MAE) for a single outcome. +# Drops NA in either vector; needs >= 3 valid points and non-constant predictions. +# @noRd +.cvMetricRow <- function(pred, actual) { + out <- setNames(rep(NA_real_, 6L), + c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE")) + ok <- !is.na(pred) & !is.na(actual) + pred <- pred[ok]; actual <- actual[ok] + if (length(pred) < 3L || stats::sd(pred) == 0) return(out) + lmFit <- stats::lm(actual ~ pred); s <- summary(lmFit) + out["corr"] <- stats::cor(actual, pred) + out["rsq"] <- s$r.squared + out["adj_rsq"] <- s$adj.r.squared + out["pval"] <- if (nrow(s$coefficients) >= 2L) s$coefficients[2L, 4L] else NA_real_ + res <- actual - pred + out["RMSE"] <- sqrt(mean(res^2)) + out["MAE"] <- mean(abs(res)) + out +} + +# Shared K-fold cross-validation engine. +# +# `fitFold(Xtrain, Ytrain, foldIndex)` must return +# list(weights = (variants x outcomes) weight matrix, +# rownames indexing colnames(Xtrain)>, +# fits = fitted model or NULL>) +# The engine drops zero-variance training columns before calling fitFold, +# predicts held-out samples (Xtest[, common] %*% W[common, ]), aggregates the +# out-of-fold predictions into full-Y matrices, scores them with .cvMetricRow, +# and keys the output via .cvOutputKey(). `maxNumVariants` (optional) randomly +# subsamples variants up front to bound compute; `numThreads` parallelises the +# fold loop (-1 = all cores, 0/1 = serial). +#' @importFrom BiocParallel bplapply bpworkers MulticoreParam +#' @importFrom stats sd lm cor +#' @noRd +.crossValidateWeights <- function(X, Y, fold = NULL, samplePartitions = NULL, + fitFold, numThreads = 1, maxNumVariants = NULL, + variantsToKeep = NULL, retainFits = FALSE, + verbose = 1) { + if (!is.null(fold) && (!is.numeric(fold) || fold <= 0)) { + stop("Invalid value for 'fold'. It must be a positive integer.") + } + if (!is.matrix(X) || (!is.matrix(Y) && !is.vector(Y))) { + stop("X must be a matrix and Y must be a matrix or a vector.") + } + if (is.vector(Y)) { + Y <- matrix(Y, ncol = 1) + if (verbose >= 1) message(sprintf( + "Y converted to matrix of %d rows and %d columns.", nrow(Y), ncol(Y))) + } + if (nrow(X) != nrow(Y)) { + stop("The number of rows in X and Y must be the same.") + } + sampleNames <- if (!is.null(rownames(Y))) rownames(Y) + else if (!is.null(rownames(X))) rownames(X) + else paste0("sample_", seq_len(nrow(X))) + if (is.null(rownames(X))) rownames(X) <- sampleNames + if (is.null(rownames(Y))) rownames(Y) <- sampleNames + if (is.null(colnames(X))) colnames(X) <- paste0("variable_", seq_len(ncol(X))) + if (is.null(colnames(Y))) colnames(Y) <- paste0("context_", seq_len(ncol(Y))) + + # Optional variant subsample (compute saver; e.g. TWAS weight CV). + if (!is.null(maxNumVariants) && ncol(X) > maxNumVariants) { + if (!is.null(variantsToKeep) && length(variantsToKeep) > 0) { + variantsToKeep <- intersect(variantsToKeep, colnames(X)) + remaining <- setdiff(colnames(X), variantsToKeep) + if (length(variantsToKeep) < maxNumVariants) { + additional <- sample(remaining, maxNumVariants - length(variantsToKeep)) + selected <- union(variantsToKeep, additional) + if (verbose >= 1) message(sprintf( + "Including %d specified variants and randomly selecting %d additional variants, for a total of %d variants out of %d for cross-validation purpose.", + length(variantsToKeep), length(additional), length(selected), ncol(X))) + } else { + selected <- sample(variantsToKeep, maxNumVariants) + if (verbose >= 1) message(sprintf( + "Randomly selecting %d out of %d input variants for cross validation purpose.", + length(selected), length(variantsToKeep))) + } + } else { + selected <- sort(sample(ncol(X), maxNumVariants)) + if (verbose >= 1) message(sprintf( + "Randomly selecting %d out of %d variants for cross validation purpose.", + length(selected), ncol(X))) + } + X <- X[, selected, drop = FALSE] + } + + # Fold partition: reuse a provided one, else build it. + if (!is.null(samplePartitions)) { + if (!all(samplePartitions$Sample %in% sampleNames)) { + stop("Some samples in 'samplePartitions' do not match 'X' and 'Y'.") + } + if (!is.null(fold) && verbose >= 1 && + fold != length(unique(samplePartitions$Fold))) { + message(sprintf(paste0( + "fold number provided does not match with sample partition, performing ", + "%d fold cross validation based on provided sample partition. "), + length(unique(samplePartitions$Fold)))) + } + samplePartition <- samplePartitions + } else if (!is.null(fold)) { + samplePartition <- .cvSamplePartition(sampleNames, fold) + } else { + stop("Either 'fold' or 'samplePartitions' must be provided.") + } + foldIds <- sort(unique(samplePartition$Fold)) + + st <- proc.time() + runFold <- function(j) { + if (verbose >= 1) message(sprintf(" CV fold %s/%s ...", j, length(foldIds))) + testIds <- samplePartition$Sample[samplePartition$Fold == j] + isTest <- rownames(X) %in% testIds + if (all(isTest) || !any(isTest)) return(list(preds = list(), fits = list())) + Xtr <- X[!isTest, , drop = FALSE] + Xte <- X[isTest, , drop = FALSE] + Ytr <- Y[!isTest, , drop = FALSE] + keep <- .nonzeroVarColumns(Xtr) + Xtr <- Xtr[, keep, drop = FALSE] + ff <- fitFold(Xtr, Ytr, j) + preds <- lapply(ff$weights, function(W) { + if (is.null(W)) return(NULL) + W[is.na(W)] <- 0 + common <- intersect(colnames(Xte), rownames(W)) + if (length(common) == 0L) return(NULL) + yhat <- Xte[, common, drop = FALSE] %*% W[common, , drop = FALSE] + rownames(yhat) <- rownames(Xte) + yhat + }) + list(preds = preds, fits = if (isTRUE(retainFits)) ff$fits else list()) + } + + numCores <- if (numThreads == -1) bpworkers(MulticoreParam()) else numThreads + numCores <- min(numCores, bpworkers(MulticoreParam())) + foldResults <- if (numCores >= 2) { + bplapply(foldIds, runFold, + BPPARAM = MulticoreParam(workers = numCores, RNGseed = 1L)) + } else { + lapply(foldIds, runFold) + } + + metricNames <- c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE") + methodKeys <- unique(unlist(lapply(foldResults, function(fr) names(fr$preds)))) + prediction <- list() + performance <- list() + for (mk in methodKeys) { + predMat <- matrix(NA_real_, nrow(Y), ncol(Y), dimnames = dimnames(Y)) + for (fr in foldResults) { + yh <- fr$preds[[mk]] + if (!is.null(yh)) predMat[rownames(yh), ] <- yh + } + prediction[[.cvOutputKey(mk, "predicted")]] <- predMat + perf <- t(vapply(seq_len(ncol(Y)), function(r) { + if (verbose >= 1) { + s <- suppressWarnings(stats::sd(predMat[, r], na.rm = TRUE)) + if (is.na(s) || s == 0) message(sprintf(paste0( + "Predicted values for condition %d using %s have zero variance. ", + "Filling performance metric with NAs"), r, mk)) + } + .cvMetricRow(predMat[, r], Y[, r]) + }, setNames(numeric(6L), metricNames))) + rownames(perf) <- colnames(Y) + performance[[.cvOutputKey(mk, "performance")]] <- perf + } + + foldFits <- setNames(lapply(foldResults, function(fr) { + ff <- fr$fits + ff[!vapply(ff, is.null, logical(1))] + }), paste0("fold_", foldIds)) + + list(samplePartition = samplePartition, prediction = prediction, + performance = performance, foldFits = foldFits, + timeElapsed = proc.time() - st) +} diff --git a/R/ctwasPipeline.R b/R/ctwasPipeline.R index 66f28270..1d846bda 100644 --- a/R/ctwasPipeline.R +++ b/R/ctwasPipeline.R @@ -24,12 +24,19 @@ #' @param gwasSumStats NAMED LIST of \code{\link{GwasSumStats}} keyed #' by \code{region_id} (at least two entries). Each must have #' \code{getQcInfo()} non-empty. -#' @param twasWeights NAMED LIST of \code{\link{TwasWeights}} keyed by -#' \code{region_id}. Keys must be a SUBSET of \code{gwasSumStats}'s -#' keys: blocks without any TWAS weights still contribute their -#' SNP-level signal to ctwas's joint group prior estimate (matches -#' the legacy whole-chromosome pattern where only a few of many LD -#' blocks carry gene weights). +#' @param twasWeights The per-gene weight source. Either (a) a FLAT +#' \code{\link{TwasWeights}} / \code{QtlFineMappingResult} (or a +#' homogeneous list of them) carrying \code{region} provenance — each +#' gene is placed into its home LD block internally by +#' \code{start(region)} (matching cTWAS's \code{p0} assignment rule); +#' or (b) a pre-bucketed NAMED LIST keyed by \code{region_id} (keys a +#' SUBSET of \code{gwasSumStats}'s), used as-is. Blocks without any +#' TWAS weights still contribute their SNP-level signal to ctwas's +#' joint group prior estimate (the legacy whole-chromosome pattern +#' where only a few of many LD blocks carry gene weights). A gene whose +#' cis span straddles a block boundary is homed by its single anchor; +#' the cross-block signal is cTWAS's boundary-gene concern +#' (\code{\link{mergeCtwasBoundaryRegions}}), not placement. #' @param twasZ Optional \code{GRanges} of TWAS Z-scores (output of #' \code{\link{causalInferencePipeline}}). When supplied, the #' per-(trait, context) Z is used as the \code{z_gene} input to @@ -44,10 +51,11 @@ #' @param method Optional character (length 1). Picks which TWAS #' method's weights to feed into ctwas for each (study, context, #' trait) gene. When \code{NULL} (default): use \code{"ensemble"} if -#' that method is present across the TwasWeights; otherwise use the -#' sole method when only one is present; otherwise error. Passing -#' the name explicitly (e.g. \code{"mrash"}) overrides the default -#' resolution. +#' that method is present across the weight sources; otherwise use the +#' sole method when only one is present; otherwise run \strong{every} +#' method as an independent cTWAS run (one \code{CtwasResult} row-set +#' per method). Passing the name explicitly (e.g. \code{"mrash"}) +#' restricts the run to that single method. #' @param thin,niterPrefit,niter,L Pass-throughs to #' \code{ctwas::ctwas_sumstats}. #' @param groupPriorVarStructure Pass-through (defaults @@ -76,10 +84,42 @@ #' NaN failure is recovered by falling back to the prefit estimates #' (mirrors the legacy ctwas_2 workaround on underpowered data). #' Default \code{FALSE}. +#' @param keepSnps Logical (length 1). When \code{TRUE}, retain the +#' context-agnostic SNP background of each run as one extra +#' \code{CtwasResult} row (\code{study = context = "SNP"}, mirroring +#' cTWAS's own \code{"SNP"} group) so the full ctwas output is +#' reconstructable from \code{\link{getFinemap}} / \code{getSusieAlpha}. +#' Default \code{FALSE} — the SNP rows are the null background and are +#' dropped from the structured gene-level result. +#' @param mergeBoundary Logical (length 1). When \code{TRUE}, run +#' \code{\link{mergeCtwasBoundaryRegions}} after fine-mapping each +#' run: a high-PIP gene whose cis window straddles an LD-block boundary +#' has its adjacent regions merged and re-fine-mapped (the legacy +#' default-off \code{ctwas_3} post-processing). Default \code{FALSE}. +#' @param mergePipThresh Numeric (length 1). PIP threshold for selecting +#' which boundary genes to merge (\code{\link{mergeCtwasBoundaryRegions}} +#' \code{pipThresh}). Default \code{0.5}. Ignored unless +#' \code{mergeBoundary = TRUE}. +#' @param mergeFilterCs Logical (length 1). Require the boundary gene to +#' be in a credible set to be selected. Default \code{FALSE}. Ignored +#' unless \code{mergeBoundary = TRUE}. +#' @param mergeMaxSNP Numeric (length 1). Per-merged-region SNP cap. +#' Default \code{Inf}. Ignored unless \code{mergeBoundary = TRUE}. #' @param ... Additional arguments forwarded to #' \code{ctwas::ctwas_sumstats}. -#' @return Whatever \code{ctwas::ctwas_sumstats} returns (a list with -#' \code{susie_alpha_res}, \code{param}, and other diagnostics). +#' @return A \code{\link{CtwasResult}} collection: one row per +#' \code{(gwasStudy, study, context, method)}. A single-context run is +#' one row per method; a multi-context (joint) run emits per-context +#' rows sharing the same \code{jointContexts} set and the +#' jointly-estimated group priors. Each row's +#' \code{\link{CtwasResultEntry}} payload carries that context's +#' per-gene fine-mapping posteriors (\code{finemap}), the run's +#' \code{param}, and its \code{regionInfo}. For the raw +#' \code{ctwas::finemap_regions} list (e.g. to feed +#' \code{\link{mergeCtwasBoundaryRegions}}), call the granular +#' \code{\link{assembleCtwasInputs}} \eqn{\to} \code{\link{estCtwasParam}} +#' \eqn{\to} \code{\link{screenCtwasRegions}} \eqn{\to} +#' \code{\link{finemapCtwasRegions}} path instead. #' @export ctwasPipeline <- function(gwasSumStats, twasWeights, @@ -101,37 +141,58 @@ ctwasPipeline <- function(gwasSumStats, minPipCutoff = 0, maxNumVariants = Inf, fallbackToPrefit = FALSE, + keepSnps = FALSE, + mergeBoundary = FALSE, + mergePipThresh = 0.5, + mergeFilterCs = FALSE, + mergeMaxSNP = Inf, ...) { groupPriorVarStructure <- match.arg(groupPriorVarStructure) - inputs <- assembleCtwasInputs( - gwasSumStats = gwasSumStats, - twasWeights = twasWeights, - twasZ = twasZ, - fineMappingResult = fineMappingResult, - method = method, - twasWeightCutoff = twasWeightCutoff, - csMinCor = csMinCor, - minPipCutoff = minPipCutoff, - maxNumVariants = maxNumVariants) - est <- estCtwasParam( - inputs, - thin = thin, - niterPrefit = niterPrefit, - niter = niter, - groupPriorVarStructure = groupPriorVarStructure, - ncore = ncore, - fallbackToPrefit = fallbackToPrefit, - ...) - screened <- screenCtwasRegions( - est, - L = L, - ncore = ncore, - ...) - finemapCtwasRegions( - screened, - L = L, - ncore = ncore, - ...) + .ctwasRequireNamedLists(gwasSumStats, twasWeights) + methods <- .ctwasResolveMethods(twasWeights, method) + gwasStudy <- .ctwasGwasStudy(gwasSumStats) + + # One independent cTWAS run per method (weights are homogeneous within a run), + # each decomposed into per-context CtwasResult row-specs. + rows <- list() + for (m in methods) { + inputs <- assembleCtwasInputs( + gwasSumStats = gwasSumStats, + twasWeights = twasWeights, + twasZ = twasZ, + fineMappingResult = fineMappingResult, + method = m, + twasWeightCutoff = twasWeightCutoff, + csMinCor = csMinCor, + minPipCutoff = minPipCutoff, + maxNumVariants = maxNumVariants) + est <- estCtwasParam( + inputs, + thin = thin, + niterPrefit = niterPrefit, + niter = niter, + groupPriorVarStructure = groupPriorVarStructure, + ncore = ncore, + fallbackToPrefit = fallbackToPrefit, + ...) + screened <- screenCtwasRegions(est, L = L, ncore = ncore, ...) + finemap <- finemapCtwasRegions(screened, L = L, ncore = ncore, ...) + # Optional boundary-gene region merging: a gene whose cis window straddles an + # LD-block boundary is split across two regions in the first-pass fine-map; + # this merges each high-PIP straddler's adjacent regions and re-fine-maps. + # The decomposition below is merge-transparent (it keys by gene id). + if (mergeBoundary) + finemap <- mergeCtwasBoundaryRegions( + finemap, pipThresh = mergePipThresh, filterCs = mergeFilterCs, + maxSNP = mergeMaxSNP, L = L, ncore = ncore, ...) + rows <- c(rows, .ctwasRunToRows(finemap, gwasStudy = gwasStudy, method = m, + keepSnps = keepSnps)) + } + if (length(rows) == 0L) + stop("ctwasPipeline: no genes were modeled (the weight sources produced no ", + "usable gene weights for method(s): ", + paste(methods, collapse = ", "), ").") + .ctwasRowsToResult(rows) } #' Assemble cTWAS inputs from S4 GwasSumStats / TwasWeights @@ -174,15 +235,17 @@ assembleCtwasInputs <- function(gwasSumStats, twasWeights, "region_id (got ", class(gwasSumStats)[[1L]], "). cTWAS's EM ", "requires multi-block context to converge; single-block calls ", "are no longer supported.") - if (missing(twasWeights) || !is.list(twasWeights) || - methods::is(twasWeights, "TwasWeights")) - stop("`twasWeights` must be a NAMED LIST of TwasWeights keyed by ", - "region_id.") if (is.null(names(gwasSumStats)) || any(!nzchar(names(gwasSumStats)))) stop("`gwasSumStats` must be a named list keyed by region_id (got an ", "unnamed or empty-named list).") + if (missing(twasWeights) || is.null(twasWeights)) + stop("`twasWeights` is required (a TwasWeights / QtlFineMappingResult ", + "weight source, or a per-region named list keyed by region_id).") + # Place a FLAT weight source into its home LD blocks by start(region) + # (cTWAS's p0 rule); a pre-bucketed per-region named list passes through. + twasWeights <- .ctwasResolveWeightBuckets(twasWeights, gwasSumStats) if (is.null(names(twasWeights)) || any(!nzchar(names(twasWeights)))) - stop("`twasWeights` must be a named list keyed by region_id (got an ", + stop("`twasWeights` must resolve to a named list keyed by region_id (got an ", "unnamed or empty-named list).") extra_tw_keys <- setdiff(names(twasWeights), names(gwasSumStats)) if (length(extra_tw_keys) > 0L) @@ -200,8 +263,10 @@ assembleCtwasInputs <- function(gwasSumStats, twasWeights, "']] has no QC record. Call summaryStatsQc() first.") } for (rid in names(twasWeights)) { - if (!methods::is(twasWeights[[rid]], "TwasWeights")) - stop("twasWeights[['", rid, "']] is not a TwasWeights.") + if (!methods::is(twasWeights[[rid]], "TwasWeights") && + !methods::is(twasWeights[[rid]], "QtlFineMappingResult")) + stop("twasWeights[['", rid, "']] must be a TwasWeights or ", + "QtlFineMappingResult (the per-gene weight source).") } if (!is.null(twasZ) && !methods::is(twasZ, "GRanges")) stop("`twasZ` must be a GRanges (output of causalInferencePipeline) ", @@ -743,6 +808,359 @@ mergeCtwasBoundaryRegions <- function(finemapResult, "pick one (e.g. method = \"mrash\").") } +# Fail fast on the two cTWAS inputs. `gwasSumStats` defines the LD-block grid, so +# it must be a NAMED LIST keyed by region_id (a bare S4 collection is the common +# mistake). `twasWeights` may be a FLAT weight source (a single TwasWeights / +# QtlFineMappingResult, or a list of them) — placed into blocks internally by +# `assembleCtwasInputs` — or a pre-bucketed per-region named list. +# @noRd +.ctwasRequireNamedLists <- function(gwasSumStats, twasWeights) { + if (!is.list(gwasSumStats) || methods::is(gwasSumStats, "GwasSumStats")) + stop("`gwasSumStats` must be a NAMED LIST of GwasSumStats keyed by ", + "region_id (got ", class(gwasSumStats)[[1L]], "). cTWAS's EM ", + "requires multi-block context to converge; single-block calls ", + "are no longer supported.") + okTw <- methods::is(twasWeights, "TwasWeights") || + methods::is(twasWeights, "QtlFineMappingResult") || is.list(twasWeights) + if (!okTw) + stop("`twasWeights` must be a TwasWeights / QtlFineMappingResult weight ", + "source (placed into LD blocks internally by region), or a per-region ", + "named list keyed by region_id (got ", class(twasWeights)[[1L]], ").") +} + +# Strip a leading "chr" and case so region_id-derived seqnames ("chr22") and +# phenotype rowRanges seqnames ("22" / "chr22") compare equal for placement. +# @noRd +.ctwasChrKey <- function(x) sub("^chr", "", tolower(as.character(x))) + +# TRUE when `tw` is ALREADY a per-region named list (the pre-bucketed contract): +# a NAMED plain list (not an S4 collection) of weight collections. Its keys are +# validated against the block grid downstream (the `extra_tw_keys` check). A +# flat collection or an UNNAMED list is instead treated as a flat weight source +# to place internally by region. +# @noRd +.ctwasIsPreBucketed <- function(tw, gwasSumStats) { + is.list(tw) && !methods::is(tw, "DFrame") && length(tw) > 0L && + !is.null(names(tw)) && all(nzchar(names(tw))) && + all(vapply(tw, function(x) methods::is(x, "TwasWeights") || + methods::is(x, "QtlFineMappingResult"), logical(1))) +} + +# Combine a flat weight source into ONE collection. Accepts a single +# TwasWeights / QtlFineMappingResult, or a homogeneous list of one kind. +# @noRd +.ctwasCombineWeightSources <- function(weights) { + if (methods::is(weights, "TwasWeights") || + methods::is(weights, "QtlFineMappingResult")) + return(weights) + if (is.list(weights)) { + weights <- weights[!vapply(weights, is.null, logical(1))] + if (length(weights) == 0L) + stop("assembleCtwasInputs: `twasWeights` is an empty weight source.") + if (all(vapply(weights, methods::is, logical(1), "TwasWeights"))) + return(do.call(combineTwasWeights, unname(weights))) + if (all(vapply(weights, methods::is, logical(1), "QtlFineMappingResult"))) + return(do.call(combineFineMappingResults, unname(weights))) + } + stop("assembleCtwasInputs: `twasWeights` must be a TwasWeights or ", + "QtlFineMappingResult (or a homogeneous list of one kind), or a ", + "per-region named list keyed by region_id.") +} + +# Place each gene (row) of a flat weight source into its home LD block. The +# anchor is start(region) -- matching cTWAS's own `assign_region_data` rule, +# which homes a gene by its p0 (single point) into the block where +# p0 in [start, stop). Returns a region_id per row (NA when the anchor falls in +# no block). A gene whose cis SPAN straddles a boundary is still homed by its +# single anchor here; the cross-block signal is cTWAS's boundary-gene concern +# (get_boundary_genes / postprocess_region_merging), not placement. +# @noRd +.ctwasPlaceByAnchor <- function(region, gwasSumStats) { + ids <- names(gwasSumStats) + blockGr <- .regionFromIds(ids) + aChr <- .ctwasChrKey(as.character(GenomicRanges::seqnames(region))) + aPos <- GenomicRanges::start(region) + bChr <- .ctwasChrKey(as.character(GenomicRanges::seqnames(blockGr))) + bS <- GenomicRanges::start(blockGr) + bE <- GenomicRanges::end(blockGr) + vapply(seq_along(aPos), function(i) { + if (is.na(aPos[[i]])) return(NA_character_) + hit <- which(bChr == aChr[[i]] & aPos[[i]] >= bS & aPos[[i]] < bE) + if (length(hit) > 0L) ids[[hit[[1L]]]] else NA_character_ + }, character(1)) +} + +# Bucket a flat weight source into a per-region named list keyed to the block +# grid, homing each gene by start(region). Each per-block sub-collection carries +# that block's GWAS LD sketch (the panel its weights are harmonized against, and +# what the downstream match-check expects). +# @noRd +.ctwasBucketWeights <- function(weights, gwasSumStats) { + combined <- .ctwasCombineWeightSources(weights) + if (!("region" %in% names(combined))) + stop("assembleCtwasInputs: the weight source carries no `region` ", + "provenance, which internal LD-block placement requires (produced by ", + "twasWeightsPipeline / fineMappingPipeline). Supply a pre-bucketed ", + "per-region named list if placement was done upstream.") + home <- .ctwasPlaceByAnchor(getRegion(combined), gwasSumStats) + unplaced <- sum(is.na(home)) + if (unplaced > 0L) + warning(sprintf(paste0("assembleCtwasInputs: %d gene(s) whose region anchor ", + "fell in no LD block were dropped."), unplaced)) + out <- list() + for (rid in names(gwasSumStats)) { + idx <- which(home == rid) + if (length(idx) == 0L) next + sub <- combined[idx, ] + sub@ldSketch <- getLdSketch(gwasSumStats[[rid]]) + out[[rid]] <- sub + } + if (length(out) == 0L) + stop("assembleCtwasInputs: no gene placed into any LD block. Check that the ", + "weight `region`s and the gwasSumStats region_id keys share a ", + "coordinate system (e.g. 'chr22_1_1000000').") + out +} + +# Resolve `twasWeights` to a per-region named list: pass a pre-bucketed list +# through, otherwise place a flat weight source by region. +# @noRd +.ctwasResolveWeightBuckets <- function(twasWeights, gwasSumStats) { + if (.ctwasIsPreBucketed(twasWeights, gwasSumStats)) return(twasWeights) + .ctwasBucketWeights(twasWeights, gwasSumStats) +} + +# Resolve the LIST of TWAS methods a `ctwasPipeline` run should iterate over +# (one independent cTWAS run per method — weights are homogeneous within a run). +# - explicit `method`: exactly that one (validated present). +# - NULL + an "ensemble" method present: just "ensemble" (the pre-combined +# weight, the historical default). +# - NULL + a single method present: that one. +# - NULL + MULTIPLE methods, no "ensemble": ALL of them (the singular +# `.ctwasResolveMethod` errors here; the pipeline instead fans out). +# @noRd +.ctwasResolveMethods <- function(twasWeightsList, method = NULL) { + methodsOf <- function(tw) if (is.null(tw)) NULL else as.character(tw$method) + available <- unique( + if (methods::is(twasWeightsList, "TwasWeights") || + methods::is(twasWeightsList, "QtlFineMappingResult")) + methodsOf(twasWeightsList) # a flat weight source + else unlist(lapply(twasWeightsList, methodsOf))) # a list of them + if (length(available) == 0L) + stop("ctwasPipeline: weight sources carry no method entries.") + if (!is.null(method) && nzchar(method)) { + if (!method %in% available) + stop("ctwasPipeline: method '", method, "' not present in the weight ", + "sources (available: ", paste(available, collapse = ", "), ").") + return(method) + } + if ("ensemble" %in% available) return("ensemble") + available # single -> length-1 (one run); multiple -> iterate over all +} + +# The single GWAS (disease) study a cTWAS run models. cTWAS solves one disease +# per run (the z_snp carries a single z per SNP), so the input blocks must all +# reference the same GWAS study. +# @noRd +.ctwasGwasStudy <- function(gwasSumStats) { + studies <- unique(unlist(lapply(gwasSumStats, function(g) + as.character(g$study)))) + studies <- studies[!is.na(studies) & nzchar(studies)] + if (length(studies) == 0L) return(NA_character_) + if (length(studies) > 1L) + stop("ctwasPipeline: the input blocks reference multiple GWAS studies (", + paste(studies, collapse = ", "), "); cTWAS models one disease per run.") + studies +} + +# Parse the cTWAS gene ids (`region|study|context|trait|method`) that name the +# assembled weights list into their identity components. `method` is the LAST +# field and `trait` everything between context and method, so a trait that +# itself contains "|" is preserved. +# @noRd +.ctwasParseGeneIds <- function(ids) { + parts <- strsplit(ids, "|", fixed = TRUE) + n <- lengths(parts) + if (any(n < 5L)) + stop("ctwasPipeline: malformed cTWAS gene id(s): ", + paste(ids[n < 5L], collapse = ", "), + " (expected 'region|study|context|trait|method').") + pick <- function(i) vapply(parts, function(p) p[[i]], character(1)) + data.frame( + id = ids, + rid = pick(1L), + study = pick(2L), + context = pick(3L), + trait = mapply(function(p, k) paste(p[4:(k - 1L)], collapse = "|"), + parts, n, USE.NAMES = FALSE), + method = mapply(function(p, k) p[[k]], parts, n, USE.NAMES = FALSE), + stringsAsFactors = FALSE) +} + +# Enforce the multi-context joint-model invariant: every context in a run must +# carry the SAME set of genes (traits). A cTWAS joint fit couples the contexts +# through shared group priors; contexts with disjoint gene sets make the joint +# model meaningless. No-op for a single-context run. +# @noRd +.ctwasAssertSharedGenes <- function(parsed) { + byCtx <- split(parsed$trait, parsed$context) + if (length(byCtx) < 2L) return(invisible()) + geneSets <- lapply(byCtx, function(g) sort(unique(g))) + ref <- geneSets[[1L]] + mismatch <- names(geneSets)[ + !vapply(geneSets, function(s) identical(s, ref), logical(1))] + if (length(mismatch) > 0L) + stop("ctwasPipeline: multi-context cTWAS requires the SAME gene set in ", + "every context (the joint model is only meaningful when genes are ", + "shared across contexts). Context(s) differing from the reference '", + names(geneSets)[[1L]], "' (", length(ref), " gene(s)): ", + paste(mismatch, collapse = ", "), ".") + invisible() +} + +# Subset a ctwas result frame (finemap_res / susie_alpha_res) to the rows whose +# `id` is in `ids`. Returns NULL when the frame is absent or nothing matches. +# @noRd +.ctwasSubsetById <- function(df, ids) { + if (is.null(df)) return(NULL) + sub <- df[as.character(df$id) %in% ids, , drop = FALSE] + if (nrow(sub) == 0L) NULL else `rownames<-`(sub, NULL) +} + +# Subset a ctwas result frame to its SNP rows (type == "SNP"; anno_susie tags +# the non-gene background this way). Returns NULL when absent or none present. +# @noRd +.ctwasSubsetSnp <- function(df) { + if (is.null(df) || is.null(df$type)) return(NULL) + sub <- df[as.character(df$type) == "SNP", , drop = FALSE] + if (nrow(sub) == 0L) NULL else `rownames<-`(sub, NULL) +} + +# Decompose one cTWAS run (a `finemapCtwasRegions` output) into per-context +# row-specs for a CtwasResult. The row skeleton comes from the ASSEMBLED weights +# (so every modeled (study, context) appears even if no gene reached +# fine-mapping); each row's finemap / susieAlpha payloads are the subsets whose +# gene id belongs to that context. Multi-context runs are annotated with the +# shared `jointContexts` set and share the jointly-estimated `param`. +# +# `keepSnps` (default FALSE) additionally retains the context-agnostic SNP +# background as ONE extra row (study = context = "SNP"), mirroring cTWAS's own +# "SNP" group in `group_prior`. Kept off by default because the SNP rows are the +# null background and bloat the structured gene-level result; when on, the full +# ctwas run is reconstructable from `getFinemap()` / `getSusieAlpha()`. +# @noRd +.ctwasRunToRows <- function(runResult, gwasStudy, method, keepSnps = FALSE) { + geneIds <- names(runResult$weights) + if (is.null(geneIds) || length(geneIds) == 0L) return(list()) + parsed <- .ctwasParseGeneIds(geneIds) + contexts <- unique(parsed$context) + .ctwasAssertSharedGenes(parsed) + jointStr <- if (length(contexts) > 1L) + paste(sort(unique(contexts)), collapse = ",") else NA_character_ + fmDf <- if (is.null(runResult$finemap_res)) NULL else + as.data.frame(runResult$finemap_res) + saDf <- if (is.null(runResult$susie_alpha_res)) NULL else + as.data.frame(runResult$susie_alpha_res) + mkEntry <- function(fm, sa) CtwasResultEntry( + finemap = fm, susieAlpha = sa, param = runResult$param, + regionInfo = runResult$region_info) + + rows <- lapply(contexts, function(cx) { + inCx <- parsed$context == cx + studyCx <- unique(parsed$study[inCx]) + if (length(studyCx) != 1L) + stop("ctwasPipeline: context '", cx, "' mixes multiple QTL studies (", + paste(studyCx, collapse = ", "), "); one study per context.") + idsCx <- parsed$id[inCx] + list(gwasStudy = gwasStudy, study = studyCx, context = cx, method = method, + jointContexts = jointStr, + entry = mkEntry(.ctwasSubsetById(fmDf, idsCx), + .ctwasSubsetById(saDf, idsCx))) + }) + + if (keepSnps) { + snpFm <- .ctwasSubsetSnp(fmDf) + snpSa <- .ctwasSubsetSnp(saDf) + if (!is.null(snpFm) || !is.null(snpSa)) + rows <- c(rows, list(list( + gwasStudy = gwasStudy, study = "SNP", context = "SNP", method = method, + jointContexts = jointStr, entry = mkEntry(snpFm, snpSa)))) + } + rows +} + +# Assemble accumulated per-run row-specs (from .ctwasRunToRows) into a single +# CtwasResult; the jointContexts column is omitted when every row is +# single-context. Shared by ctwasPipeline (across methods) and asCtwasResult. +# @noRd +.ctwasRowsToResult <- function(rows) { + if (length(rows) == 0L) + stop("cTWAS: no genes were modeled (the weight source produced no usable ", + "gene weights).") + jointContexts <- vapply(rows, function(r) r$jointContexts, character(1)) + CtwasResult( + gwasStudy = vapply(rows, function(r) r$gwasStudy, character(1)), + study = vapply(rows, function(r) r$study, character(1)), + context = vapply(rows, function(r) r$context, character(1)), + method = vapply(rows, function(r) r$method, character(1)), + entry = lapply(rows, function(r) r$entry), + jointContexts = if (any(!is.na(jointContexts))) jointContexts else NULL) +} + +# The single GWAS study a finemap result models, read from z_snp$study (which +# `.ctwasBuildZSnp` stamps per row). Errors on multiple; NA when absent. +# @noRd +.ctwasGwasStudyFromZSnp <- function(zSnp) { + if (is.null(zSnp) || is.null(zSnp$study)) return(NA_character_) + s <- unique(as.character(zSnp$study)) + s <- s[!is.na(s) & nzchar(s)] + if (length(s) == 0L) return(NA_character_) + if (length(s) > 1L) + stop("asCtwasResult: z_snp references multiple GWAS studies (", + paste(s, collapse = ", "), "); cTWAS models one disease per run.") + s +} + +# The single weight method a finemap result was built for, read from the +# assembled weight ids. Errors on a mix (the granular path is one method/run). +# @noRd +.ctwasMethodFromWeights <- function(weights) { + if (is.null(weights) || length(weights) == 0L) + stop("asCtwasResult: the finemap result carries no weights to derive a ", + "method from.") + m <- unique(.ctwasParseGeneIds(names(weights))$method) + if (length(m) != 1L) + stop("asCtwasResult: the finemap result mixes weight methods (", + paste(m, collapse = ", "), "); expected one per run.") + m +} + +#' @title Structure a granular cTWAS finemap result as a CtwasResult +#' @description Decompose the raw list returned by +#' \code{\link{finemapCtwasRegions}} (optionally after +#' \code{\link{mergeCtwasBoundaryRegions}}) into the structured, per-(study, +#' context) \code{\link{CtwasResult}} — the same decomposition +#' \code{\link{ctwasPipeline}} applies to its one-shot output, exposed for the +#' granular \code{assembleCtwasInputs} \eqn{\to} \code{estCtwasParam} \eqn{\to} +#' \code{screenCtwasRegions} \eqn{\to} \code{finemapCtwasRegions} path. The +#' GWAS study is read from \code{z_snp$study} and the (single) weight method +#' from the gene ids. +#' @param finemapResult A list from \code{\link{finemapCtwasRegions}} or +#' \code{\link{mergeCtwasBoundaryRegions}}. +#' @param keepSnps Logical (length 1). Retain the context-agnostic SNP +#' background as one extra \code{study = context = "SNP"} row. Default +#' \code{FALSE}. See \code{\link{ctwasPipeline}}. +#' @return A \code{\link{CtwasResult}}. +#' @seealso \code{\link{ctwasPipeline}}, \code{\link{finemapCtwasRegions}} +#' @export +asCtwasResult <- function(finemapResult, keepSnps = FALSE) { + gwasStudy <- .ctwasGwasStudyFromZSnp(finemapResult$z_snp) + method <- .ctwasMethodFromWeights(finemapResult$weights) + rows <- .ctwasRunToRows(finemapResult, gwasStudy = gwasStudy, method = method, + keepSnps = keepSnps) + .ctwasRowsToResult(rows) +} + # Subset a TwasWeights collection to rows whose `method` matches the # resolved method. Used to enforce the "one ctwas gene per (study, # context, trait)" semantics — the legacy pipeline fed a single @@ -752,13 +1170,11 @@ mergeCtwasBoundaryRegions <- function(finemapResult, .ctwasFilterMethod <- function(tw, method) { keep <- which(as.character(tw$method) == method) if (length(keep) == 0L) return(NULL) - TwasWeights( - study = as.character(tw$study)[keep], - context = as.character(tw$context)[keep], - trait = as.character(tw$trait)[keep], - method = as.character(tw$method)[keep], - entry = as.list(tw$entry[keep]), - ldSketch = getLdSketch(tw)) + # Row subset carries every column forward (joint* / region / ...); the old + # hand-listed rebuild silently dropped them. + out <- tw[keep, ] + out@ldSketch <- getLdSketch(tw) + out } # Build the per-variant Z data.frame ctwas expects from a GwasSumStats. @@ -979,9 +1395,12 @@ mergeCtwasBoundaryRegions <- function(finemapResult, out <- list() for (i in seq_len(nrow(twasWeights))) { entry <- twasWeights$entry[[i]] - origVids <- getVariantIds(entry) - origW <- as.numeric(getWeights(entry)) - if (length(origVids) == 0L || length(origVids) != length(origW)) next + # Shared resolver: aligned (variantIds, weights) from either a + # TwasWeightsEntry or a FineMappingEntry (topLoci posterior effect). + wr <- resolveWeights(entry) + origVids <- wr$variantIds + origW <- wr$weights + if (length(origVids) == 0L) next # --- Step 1: allele-harmonize against the LD panel ------------- # Parses chr:pos:A2:A1 IDs into the data.frame `harmonizeAlleles` diff --git a/R/fineMappingPipeline.R b/R/fineMappingPipeline.R index d2af774b..e03e0ccc 100644 --- a/R/fineMappingPipeline.R +++ b/R/fineMappingPipeline.R @@ -131,8 +131,10 @@ #' \code{NULL} (all contexts). #' @param traitId Optional character vector of trait names to restrict #' processing to. -#' @param region Optional \code{GRanges} for QtlDataset trait -#' selection. Mutually exclusive with \code{traitId}. +#' @param region Optional variant window for QtlDataset trait selection: a +#' \code{GRanges}, a \code{"chr:start-end"} string, or a one-row data.frame +#' with \code{chrom}/\code{start}/\code{end}. Mutually exclusive with +#' \code{traitId}. #' @param cisWindow For QtlDataset: cis-window (bp) around each trait's #' genomic position when extracting variants. Required when #' \code{traitId} is supplied. Mutually exclusive with \code{region}. @@ -170,6 +172,10 @@ #' this partition and feeds these predictions into the SR-TWAS ensemble. #' Individual-level (\code{QtlDataset} / \code{MultiStudyQtlDataset}) #' input only; ignored for sumstat inputs. +#' @param cvThreads Integer. Number of parallel workers for the +#' cross-validation fold refits (passed to the shared CV engine's +#' \code{numThreads}). Default \code{1} (serial); \code{-1} uses all cores. +#' Only consulted when \code{cvFolds > 1}. #' @param samplePartition Optional pre-defined CV partition #' \code{data.frame} with columns \code{Sample} and \code{Fold}. When #' supplied (and \code{cvFolds > 1}), every method reuses this exact @@ -565,26 +571,8 @@ setGeneric("fineMappingPipeline", ldSketch = ldSketch) } -# Combine an optional joint column across two collections. Returns NULL -# when neither input carries the column (so the rebuilt collection -# omits it too); otherwise pads the missing side with NA_character_ so -# both halves contribute a same-length character vector. -# @noRd -.combineJointCol <- function(a, b, colName) { - hasA <- colName %in% names(a) - hasB <- colName %in% names(b) - if (!hasA && !hasB) return(NULL) - aVals <- if (hasA) as.character(a[[colName]]) - else rep(NA_character_, nrow(a)) - bVals <- if (hasB) as.character(b[[colName]]) - else rep(NA_character_, nrow(b)) - c(aVals, bVals) -} - -# Concatenate two FineMappingResult collections row-wise. Routes to the -# right constructor based on input class. rbind() on a DFrame subclass -# does not reliably preserve the ldSketch slot, so this helper rebuilds -# the collection via the constructor. +# Concatenate two same-class FineMappingResult collections row-wise, carrying +# forward every column (delegates to the generic `.rbindCollections`). # @noRd .rbindFineMappingResult <- function(a, b, ldSketch = NULL) { if (!is(a, "FineMappingResultBase") || !is(b, "FineMappingResultBase")) { @@ -594,25 +582,32 @@ setGeneric("fineMappingPipeline", stop(".rbindFineMappingResult: inputs must be the same concrete class ", "(got '", class(a)[[1L]], "' and '", class(b)[[1L]], "').") } - if (is(a, "QtlFineMappingResult")) { - QtlFineMappingResult( - study = c(as.character(a$study), as.character(b$study)), - context = c(as.character(a$context), as.character(b$context)), - trait = c(as.character(a$trait), as.character(b$trait)), - method = c(as.character(a$method), as.character(b$method)), - entry = c(as.list(a$entry), as.list(b$entry)), - jointStudies = .combineJointCol(a, b, "jointStudies"), - jointContexts = .combineJointCol(a, b, "jointContexts"), - jointTraits = .combineJointCol(a, b, "jointTraits"), - ldSketch = ldSketch) - } else { - GwasFineMappingResult( - study = c(as.character(a$study), as.character(b$study)), - method = c(as.character(a$method), as.character(b$method)), - region_id = c(as.character(a$region_id), as.character(b$region_id)), - entry = c(as.list(a$entry), as.list(b$entry)), - ldSketch = ldSketch) - } + # Carry forward every column (region_id / joint* / ...) via the generic + # combine; the concrete class (QTL vs GWAS) is preserved automatically. + .rbindCollections(list(a, b), ldSketch = ldSketch) +} + +#' Combine FineMappingResult collections +#' +#' Row-bind two or more fine-mapping result collections of the SAME concrete +#' class (all \code{\link{QtlFineMappingResult}} or all +#' \code{\link{GwasFineMappingResult}}) into one -- e.g. per-block GWAS results +#' into a genome-wide collection for cTWAS. Mixing the two concrete classes is +#' an error. +#' +#' @param ... Two or more \code{FineMappingResultBase} objects, or a single +#' \code{list} of them. +#' @param ldSketch Optional \code{\link{GenotypeHandle}} to attach to the +#' combined collection. Default \code{NULL}. Applied when combining two or +#' more inputs; a single input is returned unchanged. +#' @return A single combined fine-mapping result of the shared concrete class. +#' @seealso \code{\link{combineTwasWeights}} +#' @export +combineFineMappingResults <- function(..., ldSketch = NULL) { + parts <- .asCombineList(list(...), "FineMappingResultBase", + "combineFineMappingResults") + Reduce(function(a, b) .rbindFineMappingResult(a, b, ldSketch = ldSketch), + parts) } @@ -743,6 +738,9 @@ setGeneric("fineMappingPipeline", # fineMapping & twas methods. #' @keywords internal .makeXRegions <- function(region, jointRegions) { + # Accept a "chr:start-end" string / one-row data.frame as well as a GRanges + # (a GRanges passes through unchanged), so pipeline callers need not pre-parse. + if (!is.null(region)) region <- .asGRegion(region) if (is.null(region)) { list(NULL) } else if (isTRUE(jointRegions)) { @@ -1071,15 +1069,6 @@ setGeneric("fineMappingPipeline", # the TWAS snake method name (adapter methodKey) for a drop-in merge. # ============================================================================= -# Generate a Sample/Fold partition over the rows of X, matching the scheme in -# twasWeightsCv() (shuffle samples, then cut into `fold` contiguous blocks). -# @noRd -.fmMakeSamplePartition <- function(sampleNames, fold) { - idx <- sample(length(sampleNames)) - folds <- cut(seq_along(sampleNames), breaks = fold, labels = FALSE) - data.frame(Sample = sampleNames[idx], Fold = folds, stringsAsFactors = FALSE) -} - # Snake method key (e.g. "susie_inf") for a fine-mapping token, taken from the # shared adapter registry so fineMapping CV keys match the TwasWeights `method` # column and twasWeightsCv()'s prediction keys. @@ -1090,26 +1079,6 @@ setGeneric("fineMappingPipeline", sub("_weights$", "", adapter$methodKey) } -# Compact CV metric row (corr, rsq, adj_rsq, pval, RMSE, MAE) for one outcome, -# mirroring the metric block of twasWeightsCv(). -# @noRd -.fmCvMetricRow <- function(pred, actual) { - out <- setNames(rep(NA_real_, 6L), - c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE")) - ok <- !is.na(pred) & !is.na(actual) - pred <- pred[ok]; actual <- actual[ok] - if (length(pred) < 3L || stats::sd(pred) == 0) return(out) - lmFit <- stats::lm(actual ~ pred); s <- summary(lmFit) - out["corr"] <- stats::cor(actual, pred) - out["rsq"] <- s$r.squared - out["adj_rsq"] <- s$adj.r.squared - out["pval"] <- if (nrow(s$coefficients) >= 2L) s$coefficients[2L, 4L] else NA_real_ - res <- actual - pred - out["RMSE"] <- sqrt(mean(res^2)) - out["MAE"] <- mean(abs(res)) - out -} - # Fit one fine-mapping method on (Xtr, Ytr) for a CV fold and return a # variants x outcomes weight matrix (rownames = colnames(Xtr)). susie-family # tokens are fit independently (no chained init) per fold, matching @@ -1160,51 +1129,33 @@ setGeneric("fineMappingPipeline", NULL } -# Cross-validate a homogeneous set of fine-mapping `tokens` over (X, Y). For -# univariate tokens Y is a single column; for mvsusie/fsusie Y carries one -# column per condition/feature (and fsusie additionally needs `pos`). Returns -# a list(samplePartition, prediction, performance) shaped like twasWeightsCv(). +# Cross-validate a homogeneous set of fine-mapping `tokens` over (X, Y) via the +# shared .crossValidateWeights() engine. For univariate tokens Y is a single +# column; for mvsusie/fsusie Y carries one column per condition/feature (and +# fsusie additionally needs `pos`). Each token's per-fold fit is refit here; the +# engine owns partitioning, the (optionally parallel) fold loop, prediction, and +# the metric block. Returns list(samplePartition, prediction, performance), +# keyed identically to twasWeightsCv(). # @noRd -.fmCrossValidate <- function(X, Y, tokens, methodArgs, fold, - samplePartition = NULL, coverage = 0.95, - pos = NULL, verbose = 1, mvPrior = NULL, - mvPriorCv = NULL) { +.fmWeightsCv <- function(X, Y, tokens, methodArgs, fold, + samplePartition = NULL, coverage = 0.95, + pos = NULL, verbose = 1, mvPrior = NULL, + mvPriorCv = NULL, numThreads = 1) { if (length(tokens) == 0L) return(NULL) - if (!is.matrix(Y)) { - Y <- matrix(Y, ncol = 1L, - dimnames = list(rownames(X), NULL)) - } - if (is.null(rownames(Y))) rownames(Y) <- rownames(X) - sampleNames <- rownames(X) - if (is.null(samplePartition)) { - samplePartition <- .fmMakeSamplePartition(sampleNames, fold) - } - foldIds <- sort(unique(samplePartition$Fold)) - - preds <- setNames( - lapply(tokens, function(tk) { - matrix(NA_real_, nrow(Y), ncol(Y), dimnames = dimnames(Y)) - }), tokens) - - for (j in foldIds) { - testIds <- samplePartition$Sample[samplePartition$Fold == j] - isTest <- rownames(X) %in% testIds - if (all(isTest) || !any(isTest)) next - Xtr <- X[!isTest, , drop = FALSE] - Xte <- X[isTest, , drop = FALSE] - Ytr <- Y[!isTest, , drop = FALSE] + # Per-fold fit: the engine has already dropped zero-variance training columns. + # Weights are keyed by the canonical method key so _predicted / + # _performance line up with the TwasWeights method column. + fitFold <- function(Xtr, Ytr, j) { # Honest per-fold mvSuSiE prior when supplied (the fold's own mr.mash-derived # prior); otherwise the single full-data prior is reused on every fold. mvPriorThisFold <- if (!is.null(mvPriorCv)) { p <- mvPriorCv[[as.character(j)]] if (is.null(p)) mvPrior else p } else mvPrior - # Drop columns with zero variance in this training fold. - keepCol <- .nonzeroVarColumns(Xtr) - XtrK <- Xtr[, keepCol, drop = FALSE] + weights <- list() for (tk in tokens) { - W <- tryCatch( - .fmFoldWeights(tk, XtrK, Ytr, coverage, methodArgs[[tk]], pos, + weights[[.fmTwasMethodKey(tk)]] <- tryCatch( + .fmFoldWeights(tk, Xtr, Ytr, coverage, methodArgs[[tk]], pos, mvPriorThisFold), error = function(e) { if (verbose >= 1) @@ -1212,30 +1163,17 @@ setGeneric("fineMappingPipeline", j, tk, conditionMessage(e))) NULL }) - if (is.null(W)) next - common <- intersect(colnames(Xte), rownames(W)) - if (length(common) == 0L) next - yhat <- Xte[, common, drop = FALSE] %*% W[common, , drop = FALSE] - preds[[tk]][rownames(Xte), ] <- yhat } + list(weights = weights, fits = list()) } - - prediction <- list() - performance <- list() - for (tk in tokens) { - key <- .fmTwasMethodKey(tk) - prediction[[paste0(key, "_predicted")]] <- preds[[tk]] - perf <- t(vapply(seq_len(ncol(Y)), function(r) { - .fmCvMetricRow(preds[[tk]][, r], Y[, r]) - }, numeric(6L))) - rownames(perf) <- colnames(Y) - performance[[paste0(key, "_performance")]] <- perf - } - list(samplePartition = samplePartition, - prediction = prediction, performance = performance) + res <- .crossValidateWeights( + X, Y, fold = fold, samplePartitions = samplePartition, + fitFold = fitFold, numThreads = numThreads, verbose = verbose) + list(samplePartition = res$samplePartition, + prediction = res$prediction, performance = res$performance) } -# Slice a full .fmCrossValidate() result down to one method's payload, keeping +# Slice a full .fmWeightsCv() result down to one method's payload, keeping # the shared samplePartition. Stored on that method's FineMappingEntry. # @noRd .fmSliceCv <- function(cv, token) { @@ -1273,6 +1211,15 @@ setMethod("fineMappingPipeline", "QtlDataset", traitId = NULL, region = NULL, cisWindow = NULL, + # Per-call genotype-filter overrides; NULL = use the QtlDataset's + # construct-time slot value (applied lazily at extraction). + mafCutoff = NULL, + macCutoff = NULL, + xvarCutoff = NULL, + imissCutoff = NULL, + keepIndel = NULL, + keepSamples = NULL, + keepVariants = NULL, jointRegions = FALSE, jointSpecification = NULL, addSusieInf = TRUE, @@ -1285,6 +1232,7 @@ setMethod("fineMappingPipeline", "QtlDataset", medianAbsCorr = NULL, fineMappingResult = NULL, cvFolds = 0, + cvThreads = 1, samplePartition = NULL, pipCutoffToSkip = 0, usePCA = FALSE, @@ -1302,6 +1250,11 @@ setMethod("fineMappingPipeline", "QtlDataset", ...) { naAction <- match.arg(naAction) if (!is.null(seed)) set.seed(as.integer(seed)) + # Apply any per-call filter overrides to a validated copy of the dataset + # (replaces the construct-time slot values for this call only). + data <- .qtlApplyFilterOverrides(data, mafCutoff, macCutoff, xvarCutoff, + imissCutoff, keepIndel, keepSamples, + keepVariants) # `cisWindow` expands a trait's own coordinates; `region` is taken # literally. Supplying both signals a misunderstanding -> reject. if (!is.null(region) && !is.null(cisWindow)) { @@ -1330,7 +1283,7 @@ setMethod("fineMappingPipeline", "QtlDataset", methodArgs = methodArgs, xRegions = xRegions, twasWeights = twasWeights, dataDrivenPriorWeightsCutoff = dataDrivenPriorWeightsCutoff, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, pipCutoffToSkip = pipCutoffToSkip, fineMappingResult = fineMappingResult) tokens <- setdiff(tokens, c("mvsusie", "fsusie")) @@ -1471,7 +1424,7 @@ setMethod("fineMappingPipeline", "QtlDataset", .fmFitXBlock(X, y, toRun, addSusieInf, coverage, secondaryCoverage, signalCutoff, minAbsCorr, methodArgs, verbose, ctx, tid, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, af = afVec) }) @@ -1522,7 +1475,7 @@ setMethod("fineMappingPipeline", "QtlDataset", .fmFitXBlock(Xb, pcY[common], "susie", FALSE, coverage, secondaryCoverage, signalCutoff, minAbsCorr, methodArgs, verbose, ctx, pcName, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, af = afVec) }) ents <- lapply(blockEntries, function(be) be[["susie"]]) @@ -1550,7 +1503,7 @@ setMethod("fineMappingPipeline", "QtlDataset", methodArgs = methodArgs, xRegions = xRegions, twasWeights = twasWeights, dataDrivenPriorWeightsCutoff = dataDrivenPriorWeightsCutoff, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, pipCutoffToSkip = pipCutoffToSkip, fineMappingResult = fineMappingResult) jointResult <- if (is.null(jointResult)) autoJoint @@ -1600,6 +1553,7 @@ setMethod("fineMappingPipeline", "MultiStudyQtlDataset", twasWeights = NULL, dataDrivenPriorWeightsCutoff = 1e-10, cvFolds = 0, + cvThreads = 1, samplePartition = NULL, pipCutoffToSkip = 0, seed = NULL, @@ -1662,7 +1616,7 @@ setMethod("fineMappingPipeline", "MultiStudyQtlDataset", jointSpecification = NULL, addSusieInf = addSusieInf, coverage = coverage, secondaryCoverage = secondaryCoverage, signalCutoff = signalCutoff, minAbsCorr = minAbsCorr, fineMappingResult = fineMappingResult, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, pipCutoffToSkip = pipCutoffToSkip, seed = seed, naAction = naAction, verbose = verbose), dotArgs)) } diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index 28bf8fb7..7beabfa8 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -1862,7 +1862,8 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) { .fmFitXBlock <- function(X, y, toRun, addSusieInf, coverage, secondaryCoverage, signalCutoff, minAbsCorr, methodArgs, verbose, ctx, tid, - cvFolds = 0, samplePartition = NULL, af = NULL) { + cvFolds = 0, cvThreads = 1, samplePartition = NULL, + af = NULL) { chainLocal <- .fmResolveSusieChain(toRun, addSusieInf) infFit <- NULL if (chainLocal$runInf) { @@ -1898,9 +1899,10 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) { if (verbose >= 1) message(sprintf("Cross-validating (%d folds) for (context='%s', trait='%s') ...", cvFolds, ctx, tid)) - cv <- .fmCrossValidate(X, y, names(out), methodArgs, cvFolds, + cv <- .fmWeightsCv(X, y, names(out), methodArgs, cvFolds, samplePartition = samplePartition, - coverage = coverage, verbose = verbose) + coverage = coverage, verbose = verbose, + numThreads = cvThreads) for (tk in names(out)) { out[[tk]] <- .fmAttachCv(out[[tk]], .fmSliceCv(cv, tk)) } diff --git a/R/jointEngine.R b/R/jointEngine.R index 9ec0e747..26cd37b6 100644 --- a/R/jointEngine.R +++ b/R/jointEngine.R @@ -69,31 +69,46 @@ NULL foldFits = cvRes$foldFits) } -# Mutable accumulator for the joint rows the engine assembles; each add() -# appends one fitted group as one result row, deriving its identity + joint* -# members from the group's conditions. +# Mutable accumulator for the joint rows the engine assembles. Each add(...) +# stores one fitted group as a per-row RECORD (a named list of column values +# whose names match the target collection constructor's parameters). Nothing +# here enumerates columns, so adding a column needs only that it be passed to +# add(); construct() folds the records into a collection via .buildJointResult. .jointRows <- function() { e <- new.env(parent = emptyenv()) - e$study <- character(0); e$context <- character(0); e$trait <- character(0) - e$method <- character(0); e$entries <- list() - e$jointStudies <- character(0); e$jointContexts <- character(0) - e$jointTraits <- character(0) - e$add <- function(study, context, trait, method, entry, - jointStudies = NA_character_, - jointContexts = NA_character_, - jointTraits = NA_character_) { - e$study <- c(e$study, study) - e$context <- c(e$context, context) - e$trait <- c(e$trait, trait) - e$method <- c(e$method, method) - e$entries[[length(e$entries) + 1L]] <- entry - e$jointStudies <- c(e$jointStudies, jointStudies) - e$jointContexts <- c(e$jointContexts, jointContexts) - e$jointTraits <- c(e$jointTraits, jointTraits) - } + e$records <- list() + e$add <- function(...) e$records[[length(e$records) + 1L]] <- list(...) e } +# Genomic anchor of one trait for the `region` provenance column: the trait's +# own coordinates for a QtlDataset (phenotype rowRanges), or the variant-span +# window of the matching entry for a QtlSumStats (an approximation -- summary +# stats carry no gene coordinate). Returns a length-1 GRanges, or NULL when the +# anchor cannot be determined (the accumulator then records a chrUn sentinel). +.jointTraitRegion <- function(data, context, trait) { + if (methods::is(data, "QtlDataset")) { + se <- tryCatch(getPhenotypes(data, contexts = context), + error = function(e) NULL) + if (is.null(se)) return(NULL) + rr <- SummarizedExperiment::rowRanges(se) + if (!(trait %in% names(rr))) return(NULL) + return(GenomicRanges::granges(rr[trait])[1L]) + } + if (methods::is(data, "QtlSumStats")) { + idx <- which(as.character(data$context) == context & + as.character(data$trait) == trait) + if (length(idx) == 0L) return(NULL) + gr <- data$entry[[idx[[1L]]]] + if (length(gr) == 0L) return(NULL) + return(GenomicRanges::GRanges( + as.character(GenomicRanges::seqnames(gr))[[1L]], + IRanges::IRanges(min(GenomicRanges::start(gr)), + max(GenomicRanges::end(gr))))) + } + NULL +} + # ---- fitters (fitJointGroup) ------------------------------------------------ # (individual, fine-mapping) -> mvSuSiE joint fit + honest per-fold CV prior. @@ -122,10 +137,11 @@ setMethod("fitJointGroup", signature("IndividualJointGroup", "FmJointPipeline"), cvM <- NULL cvFolds <- if (is.null(cfg$cvFolds)) 0L else cfg$cvFolds if (cvFolds > 1L) { - cv <- .fmCrossValidate(Xc, Yc, "fsusie", args$methodArgs, cvFolds, + cv <- .fmWeightsCv(Xc, Yc, "fsusie", args$methodArgs, cvFolds, samplePartition = cfg$samplePartition, coverage = cfg$coverage, pos = group@pos, - verbose = verbose) + verbose = verbose, + numThreads = if (is.null(cfg$cvThreads)) 1L else cfg$cvThreads) cvM <- .fmSliceCv(cv, "fsusie") } return(lapply(seq_len(nCond), function(r) { @@ -187,10 +203,11 @@ setMethod("fitJointGroup", signature("IndividualJointGroup", "FmJointPipeline"), sp <- cfg$samplePartition if (is.null(sp) && !is.null(mvCv)) sp <- mvCv$samplePartition mvPriorCv <- .fmBuildMvsusiePriorCv(mvCv, mvFitParts, colnames(Ys), ddCut) - cv <- .fmCrossValidate(Xc, Ys, "mvsusie", args$methodArgs, cvFolds, + cv <- .fmWeightsCv(Xc, Ys, "mvsusie", args$methodArgs, cvFolds, samplePartition = sp, coverage = cfg$coverage, verbose = verbose, mvPrior = mvPrior, - mvPriorCv = mvPriorCv) + mvPriorCv = mvPriorCv, + numThreads = if (is.null(cfg$cvThreads)) 1L else cfg$cvThreads) cvM <- .fmSliceCv(cv, "mvsusie") } # One per-context entry per ORIGINAL condition: NULL for screened-out columns @@ -420,27 +437,33 @@ setMethod("fitJointGroup", signature("SumStatsJointGroup", "TwasJointPipeline"), # Both pipelines assemble identically-shaped joint rows; only the result # collection differs (the axis-3 divergence the markers encode). Only the joint* # columns for axes that actually vary are attached. -.constructJointArgs <- function(pipeline, rows) { - a <- list(study = rows$study, context = rows$context, trait = rows$trait, - method = rows$method, entry = rows$entries, - ldSketch = pipeline@config$ldSketch) - if (any(!is.na(rows$jointStudies))) a$jointStudies <- rows$jointStudies - if (any(!is.na(rows$jointContexts))) a$jointContexts <- rows$jointContexts - if (any(!is.na(rows$jointTraits))) a$jointTraits <- rows$jointTraits - a +# Fold the accumulator's per-row records into one collection: build each record +# into a 1-row collection via `constructor`, then union them with the generic +# .rbindCollections (which aligns optional columns and sets the ldSketch slot). +# The only per-row transforms are: wrap `entry` in a list, and drop an NA joint* +# value so the 1-row collection omits that column (the union re-adds it, padded, +# only when some row is joint). Every other field -- study/context/trait/method, +# region, and any future column -- flows through untouched. +.buildJointResult <- function(constructor, records, ldSketch = NULL) { + if (length(records) == 0L) return(NULL) + parts <- lapply(records, function(rec) { + rec$entry <- list(rec$entry) + for (jc in c("jointStudies", "jointContexts", "jointTraits")) + if (!is.null(rec[[jc]]) && length(rec[[jc]]) == 1L && is.na(rec[[jc]])) + rec[[jc]] <- NULL + do.call(constructor, rec) + }) + .rbindCollections(parts, ldSketch = ldSketch) } setMethod("construct", "FmJointPipeline", - function(pipeline, rows, ...) { - if (length(rows$entries) == 0L) return(NULL) - do.call(QtlFineMappingResult, .constructJointArgs(pipeline, rows)) - }) + function(pipeline, rows, ...) + .buildJointResult(QtlFineMappingResult, rows$records, + pipeline@config$ldSketch)) setMethod("construct", "TwasJointPipeline", - function(pipeline, rows, ...) { - if (length(rows$entries) == 0L) return(NULL) - do.call(TwasWeights, .constructJointArgs(pipeline, rows)) - }) + function(pipeline, rows, ...) + .buildJointResult(TwasWeights, rows$records, pipeline@config$ldSketch)) # ---- enumerators (pattern x dataForm -> list) -------------------- @@ -721,11 +744,14 @@ setMethod("construct", "TwasJointPipeline", for (i in seq_len(min(length(entries), nrow(cond)))) { e <- entries[[i]] if (is.null(e)) next + ctx <- as.character(cond$context[[i]]) + tid <- as.character(cond$trait[[i]]) + reg <- .jointTraitRegion(data, ctx, tid) rows$add(study = as.character(cond$study[[i]]), - context = as.character(cond$context[[i]]), - trait = as.character(cond$trait[[i]]), + context = ctx, trait = tid, method = method, entry = e, - jointStudies = js, jointContexts = jc, jointTraits = jt) + jointStudies = js, jointContexts = jc, jointTraits = jt, + region = reg) } } # Twas: resolve this group's fine-mapping fits + CV (keyed on its first diff --git a/R/jointSpecification.R b/R/jointSpecification.R index b46d5575..1850f54d 100644 --- a/R/jointSpecification.R +++ b/R/jointSpecification.R @@ -948,7 +948,8 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { xRegions = list(NULL), twasWeights = NULL, dataDrivenPriorWeightsCutoff = 1e-10, - cvFolds = 0, samplePartition = NULL, + cvFolds = 0, cvThreads = 1, + samplePartition = NULL, pipCutoffToSkip = 0, fineMappingResult = NULL) { # Run the joint dispatch once per region block, then merge per @@ -961,7 +962,7 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { methodArgs = methodArgs, region = rg, twasWeights = twasWeights, dataDrivenPriorWeightsCutoff = dataDrivenPriorWeightsCutoff, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, pipCutoffToSkip = pipCutoffToSkip, fineMappingResult = fineMappingResult) }) @@ -981,7 +982,8 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { region = NULL, twasWeights = NULL, dataDrivenPriorWeightsCutoff = 1e-10, - cvFolds = 0, samplePartition = NULL, + cvFolds = 0, cvThreads = 1, + samplePartition = NULL, pipCutoffToSkip = 0, fineMappingResult = NULL) { # Engine routing (jointEngine.R); one region block (the caller loops regions). @@ -990,7 +992,7 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { coverage = coverage, secondaryCoverage = secondaryCoverage, signalCutoff = signalCutoff, minAbsCorr = minAbsCorr, dataDrivenPriorWeightsCutoff = dataDrivenPriorWeightsCutoff, - cvFolds = cvFolds, samplePartition = samplePartition, + cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, verbose = verbose, ldSketch = NULL)) .runJointSpecs(parsedJointSpec, data, dataForm = "individual", pipeline = pipeline, jointMethods = intersect(methods, c("mvsusie", "fsusie")), diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R index a6bcca5f..80e83c93 100644 --- a/R/manifestLoaders.R +++ b/R/manifestLoaders.R @@ -311,7 +311,12 @@ NULL 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), + # Read every column as character; .resolveSumstatCols() coerces each field to + # its target type. This prevents readr from type-guessing an allele/ID column + # that is uniformly T/F within the region as logical (the allele "T" -> TRUE), + # which would corrupt the variant id and break LD-containment matching. + as.data.frame(readr::read_tsv(I(txt), show_col_types = FALSE, progress = FALSE, + col_types = readr::cols(.default = readr::col_character())), check.names = FALSE) } @@ -531,7 +536,11 @@ NULL 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), + # Read every column as character; .resolveSumstatCols() coerces each field + # to its target type. This prevents readr from type-guessing an allele/ID + # column that is uniformly T/F as logical (the allele "T" -> TRUE). + as.data.frame(readr::read_tsv(path, show_col_types = FALSE, progress = FALSE, + col_types = readr::cols(.default = readr::col_character())), check.names = FALSE) } .resolveSumstatCols(raw, columnMapping, label) diff --git a/R/mashWrapper.R b/R/mashWrapper.R index 4456c57b..52bae2a2 100644 --- a/R/mashWrapper.R +++ b/R/mashWrapper.R @@ -232,6 +232,96 @@ mergeMashData <- function(resData, oneData) { return(combinedData) } +#' @title Build a QtlSumStats from a Z-score matrix +#' @description Assemble a per-condition \code{\link{QtlSumStats}} from a +#' \code{variants x conditions} Z-score matrix -- the input shape the mash +#' pipeline uses when only Z is available. Each column is one condition, and +#' conditions are distinguished by \code{context}, \code{trait}, or both: the +#' columns may be different cell types / tissues (contexts), different +#' molecular phenotypes (traits), or arbitrary context x trait pairs. +#' Chromosome / position are decoded from the row (variant) identifiers via +#' \code{\link{parseVariantId}} (with a synthetic-position fallback for ids +#' that do not encode coordinates); \code{A1} / \code{A2} / \code{N} are +#' placeholders because a Z-only input carries no alleles or sample sizes +#' (mash reads only Z). A pass-through \code{qcInfo} record is set so the +#' result clears the mash QC gate. +#' @param z Numeric matrix (variants x conditions). \code{rownames(z)} are +#' variant ids (ideally \code{chr:pos:A2:A1}); \code{colnames(z)} label the +#' conditions. +#' @param study Study identifier (recycled across conditions). +#' @param ldSketch A \code{\link{GenotypeHandle}} embedded in the collection. +#' @param context Condition context label(s): a single value recycled across +#' every column, or a length-\code{ncol(z)} vector (one per condition). +#' Defaults to \code{colnames(z)} -- one context per column. +#' @param trait Condition trait label(s): a single value recycled across every +#' column, or a length-\code{ncol(z)} vector. Default \code{"mash"}. Pass +#' \code{colnames(z)} here (with a constant \code{context}) when the columns +#' are traits rather than contexts. +#' @param genome Genome build. Default \code{"GRCh38"}. +#' @param n Placeholder per-variant sample size. Default \code{1000}. +#' @param a1,a2 Placeholder alleles. Defaults \code{"A"} / \code{"G"}. +#' @param role Tag stored in the \code{qcInfo} record. Default \code{"mash"}. +#' @return A \code{\link{QtlSumStats}} with one entry per condition (column). +#' @importFrom GenomicRanges GRanges +#' @importFrom IRanges IRanges +#' @export +qtlSumStatsFromZMatrix <- function(z, study, ldSketch, + context = colnames(z), trait = "mash", + genome = "GRCh38", n = 1000L, + a1 = "A", a2 = "G", role = "mash") { + if (!is.matrix(z) || !is.numeric(z)) { + stop("qtlSumStatsFromZMatrix: `z` must be a numeric variants x conditions matrix.") + } + nCond <- ncol(z) + context <- .qszmRecycle(context, nCond, "context") + trait <- .qszmRecycle(trait, nCond, "trait") + vids <- rownames(z) + if (is.null(vids)) vids <- paste0("var", seq_len(nrow(z))) + # Decode chrom/pos from the variant ids; synthesise where they do not parse. + parsed <- tryCatch(suppressWarnings(parseVariantId(vids)), + error = function(e) NULL) + chrom <- if (!is.null(parsed)) as.character(parsed$chrom) + else rep(NA_character_, length(vids)) + pos <- if (!is.null(parsed)) suppressWarnings(as.integer(parsed$pos)) + else rep(NA_integer_, length(vids)) + chrom[is.na(chrom) | !nzchar(chrom)] <- "chr1" + pos[is.na(pos)] <- seq_along(pos)[is.na(pos)] + entries <- lapply(seq_len(nCond), function(j) { + gr <- GenomicRanges::GRanges( + seqnames = chrom, ranges = IRanges::IRanges(start = pos, width = 1L)) + S4Vectors::mcols(gr) <- S4Vectors::DataFrame( + SNP = vids, A1 = rep(a1, length(vids)), A2 = rep(a2, length(vids)), + Z = as.numeric(z[, j]), N = rep(as.integer(n), length(vids))) + gr + }) + QtlSumStats( + study = rep(as.character(study), nCond), + context = context, + trait = trait, + entry = entries, + genome = genome, + ldSketch = ldSketch, + qcInfo = list(role = role, entryAudit = vector("list", nCond))) +} + +# Internal: recycle a condition-label argument (context / trait) to one value +# per column of the Z matrix. Accepts a single value (recycled to every column) +# or a vector of length ncol(z); errors otherwise -- including NULL, which is +# how `context = colnames(z)` arrives when the matrix has no column names. +.qszmRecycle <- function(v, n, what) { + if (is.null(v)) { + stop(sprintf(paste0("qtlSumStatsFromZMatrix: `%s` is NULL; pass a length-1 ", + "or length-%d value (or give `z` column names)."), + what, n)) + } + v <- as.character(v) + if (length(v) == 1L) return(rep(v, n)) + if (length(v) == n) return(v) + stop(sprintf( + "qtlSumStatsFromZMatrix: `%s` must be length 1 or ncol(z)=%d, got %d.", + what, n, length(v))) +} + # Internal: convert a single SumStats object (post-QC) into a (Bhat, Shat) # pair of matrices keyed by context. Each row of the matrix corresponds to # one (variantId × (study, trait)) cell from the SumStats entries; each diff --git a/R/sldscWrapper.R b/R/sldscWrapper.R index c8e7812e..51ae6e74 100644 --- a/R/sldscWrapper.R +++ b/R/sldscWrapper.R @@ -516,3 +516,53 @@ metaSldscRandom <- function(perTraitEstimates, category, list(summary = newDf) }) } + +#' Random-effects meta-analysis over a subset of sLDSC traits +#' +#' Re-run the DerSimonian-Laird random-effects meta-analysis on a chosen subset +#' of the per-trait standardised tables produced by +#' \code{\link{sldscPostprocessingPipeline}} -- no regression is re-run, only the +#' already-standardised per-trait estimates are re-meta'd. Powers a "meta on a +#' subset of traits" workflow: pick traits, pick target annotation categories, +#' and get the per-category tau* / enrichment / enrichstat meta results back. +#' +#' @param postprocessResult The list returned by +#' \code{\link{sldscPostprocessingPipeline}}. Must carry a \code{$per_trait} +#' element; \code{$params$target_categories} is used when +#' \code{targetCategories} is \code{NULL}. +#' @param subsetTraits Character vector of trait ids to meta over; every id must +#' be present in \code{postprocessResult$per_trait}. +#' @param targetCategories Optional character vector of target annotation names. +#' Defaults to \code{postprocessResult$params$target_categories}. +#' @return A list with \code{tau_star_single}, \code{tau_star_joint}, +#' \code{enrichment}, and \code{enrichstat}; each is a per-category named list +#' of \code{\link{metaSldscRandom}} results. +#' @seealso \code{\link{sldscPostprocessingPipeline}}, \code{\link{metaSldscRandom}} +#' @export +sldscSubsetMeta <- function(postprocessResult, subsetTraits, + targetCategories = NULL) { + perTrait <- postprocessResult$per_trait + if (is.null(perTrait)) + stop("sldscSubsetMeta: `postprocessResult` has no `per_trait` element.") + if (is.null(targetCategories)) + targetCategories <- postprocessResult$params$target_categories + if (is.null(targetCategories) || length(targetCategories) == 0L) + stop("sldscSubsetMeta: no `targetCategories` supplied and none found in ", + "`postprocessResult$params$target_categories`.") + missingTraits <- setdiff(subsetTraits, names(perTrait)) + if (length(missingTraits) > 0L) + stop("sldscSubsetMeta: trait(s) absent from `per_trait`: ", + paste(missingTraits, collapse = ", ")) + sub <- perTrait[subsetTraits] + viewSingle <- .sldscViewForMeta(sub, "single") + viewJoint <- .sldscViewForMeta(sub, "joint") + perCategory <- function(view, quantity) + setNames(lapply(targetCategories, + function(cat) metaSldscRandom(view, cat, quantity)), + targetCategories) + list( + tau_star_single = perCategory(viewSingle, "tauStar"), + tau_star_joint = perCategory(viewJoint, "tauStar"), + enrichment = perCategory(viewSingle, "enrichment"), + enrichstat = perCategory(viewSingle, "enrichstat")) +} diff --git a/R/tupleSelectors.R b/R/tupleSelectors.R index 73186b46..f528084b 100644 --- a/R/tupleSelectors.R +++ b/R/tupleSelectors.R @@ -86,3 +86,172 @@ } idx[[1L]] } + +# Internal: row indices of a FineMappingResultBase collection matching the +# given selectors. NULL selectors, and selectors naming a column the +# collection lacks, are ignored -- QTL rows carry study/context/trait/method, +# GWAS rows carry study/method/region_id -- so the same call works on either. +# `region` matches the `region_id` column. Returns all rows when nothing +# constrains. Unlike the single-row selectors above this never errors on +# ambiguity; it is the aggregate counterpart used by getTopLoci. +.fmrRowsMatching <- function(x, study = NULL, context = NULL, trait = NULL, + method = NULL, region = NULL) { + keys <- list(study = study, context = context, trait = trait, + method = method, region_id = region) + keys <- keys[!vapply(keys, is.null, logical(1L))] + keys <- keys[names(keys) %in% names(x)] + if (length(keys) == 0L) return(seq_len(nrow(x))) + ok <- rep(TRUE, nrow(x)) + for (k in names(keys)) { + ok <- ok & as.character(x[[k]]) %in% as.character(keys[[k]]) + } + which(ok) +} + +# Internal: per-row identity metadata for a FineMappingResultBase collection, +# as a data.frame with a stable column set (study, context, trait, region_id, +# method). Columns the collection lacks are NA-filled so QTL and GWAS results +# yield the same metadata shape. Safe for zero-row collections. +.fmrRowMetadata <- function(x) { + cols <- c("study", "context", "trait", "region_id", "method") + n <- nrow(x) + vals <- lapply(cols, function(cc) { + if (cc %in% names(x)) as.character(x[[cc]]) else rep(NA_character_, n) + }) + names(vals) <- cols + as.data.frame(vals, stringsAsFactors = FALSE) +} + +# Internal: rbind data.frames whose column sets may differ, aligning on the +# union of columns (missing cells become NA). When every input already shares +# the same columns -- the common case -- column types are preserved. +.rbindAligned <- function(parts) { + if (length(parts) == 1L) return(parts[[1L]]) + cols <- Reduce(union, lapply(parts, names)) + aligned <- lapply(parts, function(p) { + for (m in setdiff(cols, names(p))) p[[m]] <- NA + p[cols] + }) + do.call(rbind, aligned) +} + +# Internal: read the per-row `region` GRanges column of a collection, or an +# empty GRanges when the collection carries none. Shared by getRegion on +# TwasWeights and FineMappingResultBase (region is a uniform column across the +# family; no derivation from topLoci / region_id). +.getRegionColumn <- function(x) { + if ("region" %in% names(x)) x[["region"]] else GenomicRanges::GRanges() +} + +# Internal: append a validated `region` GRanges to a constructor's column list +# (no-op when region is NULL). Shared by the TwasWeights / FineMappingResult +# constructors so the provenance column is added identically everywhere. +.appendRegionCol <- function(cols, region, n) { + if (is.null(region)) return(cols) + if (!methods::is(region, "GRanges")) + stop("`region` must be a GRanges (one range per row) or NULL.") + if (length(region) != n) + stop("`region` must have the same length as `study`.") + cols[["region"]] <- region + cols +} + +# Internal: validate the optional `region` GRanges column (one range per row), +# shared by the TwasWeights / FineMappingResult validity methods. Returns a +# character vector of errors (empty when valid or the column is absent). +.validateRegionColumn <- function(object) { + if (!("region" %in% names(object))) return(character(0)) + if (!methods::is(object[["region"]], "GRanges")) + return("'region' column must be a GRanges") + if (length(object[["region"]]) != nrow(object)) + return("'region' column must have one range per row") + character(0) +} + +# Internal: a length-n "unknown" column matching the type of `exemplar`, used +# to pad a collection that lacks an optional column before a union row-bind. +# Character -> NA; GRanges -> a 0-width chrUn sentinel range; else NA. +.naLikeColumn <- function(exemplar, n) { + if (methods::is(exemplar, "GRanges")) + return(GenomicRanges::GRanges(rep("chrUn", n), + IRanges::IRanges(start = 1L, width = 0L))) + if (is.character(exemplar)) return(rep(NA_character_, n)) + rep(NA, n) +} + +# Internal: row-bind two or more per-tuple collection objects (TwasWeights / +# FineMappingResult subclasses), carrying forward EVERY column. Columns are +# unioned, and a collection lacking an optional column (e.g. jointContexts or +# region) is padded with `.naLikeColumn` -- so adding a column to a class flows +# through combines automatically rather than being hand-listed at each site. +# Preserves the concrete class and sets the `ldSketch` slot explicitly (rbind +# does not reliably carry it). Returns NULL when given no inputs. +.rbindCollections <- function(parts, ldSketch = NULL) { + parts <- parts[!vapply(parts, is.null, logical(1L))] + if (length(parts) == 0L) return(NULL) + cls <- class(parts[[1L]])[[1L]] + allCols <- Reduce(union, lapply(parts, names)) + exemplar <- function(cn) { + for (p in parts) if (cn %in% names(p)) return(p[[cn]]) + NULL + } + combined <- lapply(allCols, function(cn) { + pieces <- lapply(parts, function(p) + if (cn %in% names(p)) p[[cn]] else .naLikeColumn(exemplar(cn), nrow(p))) + suppressWarnings(do.call(c, pieces)) # GRanges concat may warn on seqinfo + }) + names(combined) <- allCols + df <- do.call(S4Vectors::DataFrame, c(combined, list(check.names = FALSE))) + out <- new(cls, df, ldSketch = ldSketch) + validObject(out) + out +} + +# Internal: aggregate a per-entry accessor across every row of a +# FineMappingResultBase collection that matches the given selectors, prefixing +# each entry's rows with the row identity (study/context/trait/region_id/method) +# so rows stay attributable to their source entry. Shared by getTopLoci / getCs +# / getMarginalEffects on FineMappingResultBase. +# +# When at least one selector is supplied AND resolves to exactly one entry, this +# short-circuits to `onSingle(entry)` -- the bare per-entry view with no identity +# columns -- preserving the single-entry contract. A no-selector call always +# aggregates (even a one-row collection) so the output shape is uniform. +# +# perEntry(entry) -> data.frame : per-entry view used when aggregating. Any +# column whose name collides with an identity column is dropped before the +# metadata prefix is attached (so a stamped `method` never duplicates). +# onSingle(entry) -> any : returned verbatim when selectors pin one +# entry (defaults to perEntry). +.fmrAggregateView <- function(x, study = NULL, context = NULL, trait = NULL, + method = NULL, region = NULL, perEntry, + onSingle = perEntry) { + anySelector <- !is.null(study) || !is.null(context) || !is.null(trait) || + !is.null(method) || !is.null(region) + sel <- NULL + if (anySelector) { + sel <- tryCatch( + .fmrSelectEntry(x, study = study, context = context, trait = trait, + method = method, region = region), + error = function(e) e) + if (!inherits(sel, "error")) return(onSingle(sel)) + } + idx <- .fmrRowsMatching(x, study = study, context = context, trait = trait, + method = method, region = region) + meta <- .fmrRowMetadata(x) + if (length(idx) == 0L) { + # An explicit selector that matched nothing re-raises the informative + # selection error; otherwise return an empty identity-only frame. + if (inherits(sel, "error")) stop(conditionMessage(sel)) + return(meta[integer(0), , drop = FALSE]) + } + parts <- lapply(idx, function(i) { + v <- perEntry(x$entry[[i]]) + if (is.null(v) || nrow(v) == 0L) return(NULL) + v <- v[, setdiff(names(v), names(meta)), drop = FALSE] + cbind(meta[rep(i, nrow(v)), , drop = FALSE], v, row.names = NULL) + }) + parts <- parts[!vapply(parts, is.null, logical(1L))] + if (length(parts) == 0L) return(meta[integer(0), , drop = FALSE]) + .rbindAligned(parts) +} diff --git a/R/twasWeights.R b/R/twasWeights.R index c67b312f..b92274ff 100644 --- a/R/twasWeights.R +++ b/R/twasWeights.R @@ -33,6 +33,8 @@ setClass("TwasWeights", if (!all(entryTypes)) errors <- c(errors, "every element of the `entry` column must be a TwasWeightsEntry") + # Optional `region` provenance (one genomic anchor per row; non-key). + errors <- c(errors, .validateRegionColumn(object)) jointCols <- intersect( c("jointStudies", "jointContexts", "jointTraits"), names(object)) for (jc in jointCols) { @@ -137,6 +139,12 @@ setClass("TwasWeights", #' Same shape as \code{jointStudies}. #' @param jointTraits Optional character vector for cross-trait joints. #' Same shape as \code{jointStudies}. +#' @param region Optional \code{GRanges} (length \code{length(study)}) giving +#' the genomic anchor of each row's trait -- the trait's own coordinates when +#' built from a \code{QtlDataset}, or the summary-stat variant-span window +#' when built from a \code{QtlSumStats}. Carried forward as provenance (e.g. +#' for cTWAS LD-block placement); not part of the identity key. \code{NULL} +#' (default) omits the column. #' @param ldSketch An optional \code{GenotypeHandle}, or \code{NULL} for #' individual-level fits. #' @return A \code{TwasWeights} object. @@ -145,6 +153,7 @@ TwasWeights <- function(study, context, trait, method, entry, jointStudies = NULL, jointContexts = NULL, jointTraits = NULL, + region = NULL, ldSketch = NULL) { n <- length(study) if (length(context) != n || length(trait) != n || length(method) != n || @@ -168,6 +177,7 @@ TwasWeights <- function(study, context, trait, method, entry, stop("`", pair[[1L]], "` must have the same length as `study`.") cols[[pair[[2L]]]] <- as.character(val) } + cols <- .appendRegionCol(cols, region, n) df <- do.call(S4Vectors::DataFrame, c(cols, list(check.names = FALSE))) obj <- new("TwasWeights", df, ldSketch = ldSketch) @@ -175,6 +185,10 @@ TwasWeights <- function(study, context, trait, method, entry, obj } +#' @rdname getRegion +#' @export +setMethod("getRegion", "TwasWeights", function(x, ...) .getRegionColumn(x)) + #' @title Get a Single TWAS Weights Entry #' @description Return the \code{TwasWeightsEntry} for one #' \code{(study, context, trait, method)} row of a \code{TwasWeights} @@ -617,298 +631,86 @@ setMethod("show", "TwasWeights", function(object) { #' @importFrom quadprog solve.QP #' @export twasWeightsCv <- function(X, Y, fold = NULL, samplePartitions = NULL, weightMethods = NULL, maxNumVariants = NULL, variantsToKeep = NULL, numThreads = 1, verbose = 1, retainFits = FALSE, ...) { - splitData <- function(X, Y, samplePartition, fold) { - testIds <- samplePartition[which(samplePartition$Fold == fold), "Sample"] - Xtrain <- X[!(rownames(X) %in% testIds), , drop = FALSE] - Ytrain <- Y[!(rownames(Y) %in% testIds), , drop = FALSE] - Xtest <- X[rownames(X) %in% testIds, , drop = FALSE] - Ytest <- Y[rownames(Y) %in% testIds, , drop = FALSE] - if (nrow(Xtrain) == 0 || nrow(Ytrain) == 0 || nrow(Xtest) == 0 || nrow(Ytest) == 0) { - stop("Error: One of the datasets (train or test) has zero rows.") - } - return(list(Xtrain = Xtrain, Ytrain = Ytrain, Xtest = Xtest, Ytest = Ytest)) - } - - # Validation checks - if (!is.null(fold) && (!is.numeric(fold) || fold <= 0)) { - stop("Invalid value for 'fold'. It must be a positive integer.") - } - - if (!is.matrix(X) || (!is.matrix(Y) && !is.vector(Y))) { - stop("X must be a matrix and Y must be a matrix or a vector.") - } - - if (is.vector(Y)) { - Y <- matrix(Y, ncol = 1) - if (verbose >= 1) message(paste("Y converted to matrix of", nrow(Y), "rows and", ncol(Y), "columns.")) - } - - if (nrow(X) != nrow(Y)) { - stop("The number of rows in X and Y must be the same.") - } - if (!is.null(rownames(X)) && !is.null(rownames(Y))) { - if (!identical(rownames(X), rownames(Y))) { - rownames(X) <- rownames(Y) - } - sampleNames <- rownames(Y) - } else if (!is.null(rownames(Y))) { - sampleNames <- rownames(Y) - } else if (!is.null(rownames(X))) { - sampleNames <- rownames(X) - } else { - sampleNames <- paste0("sample_", 1:nrow(X)) - } - if (is.null(rownames(X))) { - rownames(X) <- sampleNames - } - if (is.null(rownames(Y))) { - rownames(Y) <- sampleNames - } - - if (is.null(colnames(X))) { - colnames(X) <- paste0("variable_", 1:ncol(X)) - } - if (is.null(colnames(Y))) { - colnames(Y) <- paste0("context_", 1:ncol(Y)) - } - if (is.character(weightMethods)) { weightMethods <- .twasMethodLookup(weightMethods) } - - if (!exists(".Random.seed")) { - if (verbose >= 1) message("! No seed has been set. Please set seed for reproducable result. ") - } - - # Select variants if necessary - if (!is.null(maxNumVariants) && ncol(X) > maxNumVariants) { - if (!is.null(variantsToKeep) && length(variantsToKeep) > 0) { - variantsToKeep <- intersect(variantsToKeep, colnames(X)) - remainingColumns <- setdiff(colnames(X), variantsToKeep) - if (length(variantsToKeep) < maxNumVariants) { - additionalColumns <- sample(remainingColumns, maxNumVariants - length(variantsToKeep), replace = FALSE) - selectedColumns <- union(variantsToKeep, additionalColumns) - if (verbose >= 1) message(sprintf( - "Including %d specified variants and randomly selecting %d additional variants, for a total of %d variants out of %d for cross-validation purpose.", - length(variantsToKeep), length(additionalColumns), length(selectedColumns), ncol(X) - )) - } else { - selectedColumns <- sample(variantsToKeep, maxNumVariants, replace = FALSE) - if (verbose >= 1) message(paste("Randomly selecting", length(selectedColumns), "out of", length(variantsToKeep), "input variants for cross validation purpose.")) - } - } else { - selectedColumns <- sort(sample(ncol(X), maxNumVariants, replace = FALSE)) - if (verbose >= 1) message(paste("Randomly selecting", length(selectedColumns), "out of", ncol(X), "variants for cross validation purpose.")) - } - X <- X[, selectedColumns, drop = FALSE] - } - - # Create or use provided folds - if (!is.null(fold)) { - if (!is.null(samplePartitions)) { - if (fold != length(unique(samplePartitions$Fold))) { - if (verbose >= 1) message(paste0( - "fold number provided does not match with sample partition, performing ", length(unique(samplePartitions$Fold)), - " fold cross validation based on provided sample partition. " - )) - } - - folds <- samplePartitions$Fold - samplePartition <- samplePartitions - } else { - sampleIndices <- sample(nrow(X)) - folds <- cut(seq(1, nrow(X)), breaks = fold, labels = FALSE) - samplePartition <- data.frame(Sample = sampleNames[sampleIndices], Fold = folds, stringsAsFactors = FALSE) - } - } else if (!is.null(samplePartitions)) { - if (!all(samplePartitions$Sample %in% sampleNames)) { - stop("Some samples in 'samplePartitions' do not match the samples in 'X' and 'Y'.") - } - samplePartition <- samplePartitions - fold <- length(unique(samplePartition$Fold)) - } else { - stop("Either 'fold' or 'samplePartitions' must be provided.") + if (!exists(".Random.seed") && verbose >= 1) { + message("! No seed has been set. Please set seed for reproducable result. ") } - - st <- proc.time() - if (is.null(weightMethods)) { - return(list(samplePartition = samplePartition)) - } else { - # Hardcoded vector of multivariate weightMethods (accept both snake and camel). - # fSuSiE is excluded from the per-fold CV refit path: it is functional and - # cannot be refit from a bare (X, y) fold split, so its cross-validated - # predictions are supplied by fineMappingPipeline (FineMappingResult - # cvResult) rather than recomputed here. - multivariateWeightMethods <- c("mrmash_weights", "mvsusie_weights", - "mrmashWeights", "mvsusieWeights") - - # Determine the number of cores to use - numCores <- ifelse(numThreads == -1, - bpworkers(MulticoreParam()), - numThreads) - numCores <- min(numCores, - bpworkers(MulticoreParam())) - - cvArgs <- list(...) - - # Perform CV with parallel processing - computeMethodPredictions <- function(j) { - if (verbose >= 1) { - message(sprintf(" CV fold %d/%d ...", j, fold)) - tic() - } - datSplit <- splitData(X, Y, samplePartition = samplePartition, fold = j) - Xtrain <- datSplit$Xtrain - Ytrain <- datSplit$Ytrain - Xtest <- datSplit$Xtest - Ytest <- datSplit$Ytest - - # Remove columns with zero variance in the training fold. - # NOTE: Y was already NA-handled at the pipeline layer - # (getResidualizedPhenotypes naAction = "drop"/"impute"), so there - # is no need to mask X by Y NA rows here. - validColumns <- .nonzeroVarColumns(Xtrain) - Xtrain <- Xtrain[, validColumns, drop = FALSE] - validColumns <- colnames(Xtrain) - # Xtest <- Xtest[, validColumns, drop=FALSE] - foldWeightMethods <- .prepareSusieWeightMethods(Xtrain, Ytrain, weightMethods) - - foldOut <- setNames(lapply(names(foldWeightMethods), function(method) { - args <- foldWeightMethods[[method]] - fnName <- .resolveMethodFunction(method, args) - capturedFit <- NULL - - if (method %in% multivariateWeightMethods) { - # Apply multivariate method to entire Y for this fold. Per-fold priors - # bind to the fitter's camelCase args: the fold's mr.mash data-driven - # prior, and the fold's reweighted mvSuSiE prior. - if (!is.null(cvArgs$data_driven_prior_matrices_cv) && - method %in% c("mrmash_weights", "mrmashWeights")) { - args$dataDrivenPriorMatrices <- cvArgs$data_driven_prior_matrices_cv[[j]] - } - if (!is.null(cvArgs$reweightedMixturePriorCv) && - method %in% c("mvsusie_weights", "mvsusieWeights")) { - args$prior_variance <- cvArgs$reweightedMixturePriorCv[[j]] - } - # Retain the per-fold fitted model (e.g. the fold's mr.mash fit) when - # asked and supported, so a caller can reuse it as that fold's prior - # (full-CV mvSuSiE). No extra fitting -- it is the model fit here. - if (isTRUE(retainFits) && "retainFit" %in% names(formals(fnName))) { - args$retainFit <- TRUE - } - weightsMatrix <- if (verbose < 2) { - .quietEval(do.call(fnName, c(list(X = Xtrain, Y = Ytrain), args))) - } else { - do.call(fnName, c(list(X = Xtrain, Y = Ytrain), args)) - } - capturedFit <- attr(weightsMatrix, "fit") - attr(weightsMatrix, "fit") <- NULL - rownames(weightsMatrix) <- colnames(Xtrain) - fullWeightsMatrix <- .embedWeights(weightsMatrix[validColumns, , drop = FALSE], validColumns, ncol(X), ncol(Y), colnames(X), colnames(Y)) - Ypred <- Xtest %*% fullWeightsMatrix - rownames(Ypred) <- rownames(Xtest) - } else { - Ypred <- sapply(1:ncol(Ytrain), function(k) { - weights <- if (verbose < 2) { - .quietEval(do.call(fnName, c(list(X = Xtrain, y = Ytrain[, k]), args))) - } else { - do.call(fnName, c(list(X = Xtrain, y = Ytrain[, k]), args)) - } - fullWeights <- rep(0, ncol(X)) - names(fullWeights) <- colnames(X) - fullWeights[validColumns] <- weights - # Handle NAs in weights - fullWeights[is.na(fullWeights)] <- 0 - Xtest %*% fullWeights - }) - rownames(Ypred) <- rownames(Xtest) + cvArgs <- list(...) + # Multivariate weight methods (snake + camel) are fit on the whole Y for a + # fold; univariate methods are fit per Y column. fSuSiE is intentionally + # absent -- it is functional and cannot be refit from a bare (X, y) fold + # split, so its cross-validated predictions are supplied by fineMappingPipeline. + multivariateWeightMethods <- c("mrmash_weights", "mvsusie_weights", + "mrmashWeights", "mvsusieWeights") + + # Per-fold fit passed to the shared engine. Weights are keyed by the canonical + # method key (drives the _predicted / _performance output); captured + # fits keep the full method name (foldFits back-compat). + weightFitFold <- function(Xtr, Ytr, j) { + foldWeightMethods <- .prepareSusieWeightMethods(Xtr, Ytr, weightMethods) + weights <- list() + fits <- list() + for (method in names(foldWeightMethods)) { + args <- foldWeightMethods[[method]] + fnName <- .resolveMethodFunction(method, args) + mk <- sub("_weights$|Weights$", "", method) + capturedFit <- NULL + if (method %in% multivariateWeightMethods) { + # Per-fold priors bind to the fitter's camelCase args. + if (!is.null(cvArgs$data_driven_prior_matrices_cv) && + method %in% c("mrmash_weights", "mrmashWeights")) { + args$dataDrivenPriorMatrices <- cvArgs$data_driven_prior_matrices_cv[[j]] } - list(pred = Ypred, fit = capturedFit) - }), names(foldWeightMethods)) - if (verbose >= 1) { - elapsed <- toc(quiet = TRUE) - message(sprintf(" CV fold %d/%d done in %.1fs", j, fold, elapsed$toc - elapsed$tic)) - } - list(preds = lapply(foldOut, `[[`, "pred"), - fits = lapply(foldOut, `[[`, "fit")) - } - - if (numCores >= 2) { - bpParam <- MulticoreParam(workers = numCores, - RNGseed = 1L) - foldResults <- bplapply(1:fold, - computeMethodPredictions, BPPARAM = bpParam) - } else { - foldResults <- map(1:fold, computeMethodPredictions) - } - - # Reorganize into Ypred - # After cross validation, each sample should have been in - # test set at some point, and therefore has predicted value. - # The prediction matrix is therefore exactly the same dimension as input Y - Ypred <- setNames(lapply(weightMethods, function(x) `dimnames<-`(matrix(NA, nrow(Y), ncol(Y)), dimnames(Y))), names(weightMethods)) - for (j in seq_along(foldResults)) { - for (method in names(weightMethods)) { - Ypred[[method]][rownames(foldResults[[j]]$preds[[method]]), ] <- foldResults[[j]]$preds[[method]] - } - } - - names(Ypred) <- .renameSuffix(names(Ypred), "predicted") - - # Compute rsq, adj rsq, p-value, RMSE, and MAE for each method - metricsTable <- list() - - for (m in names(weightMethods)) { - metricsTable[[m]] <- matrix(NA, nrow = ncol(Y), ncol = 6) - colnames(metricsTable[[m]]) <- c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE") - rownames(metricsTable[[m]]) <- colnames(Y) - - for (r in 1:ncol(Y)) { - methodPredictions <- Ypred[[.renameSuffix(m, "predicted")]][, r] - actualValues <- Y[, r] - # Remove missing values in the first place - naIndx <- which(is.na(actualValues)) - if (length(naIndx) != 0) { - methodPredictions <- methodPredictions[-naIndx] - actualValues <- actualValues[-naIndx] + if (!is.null(cvArgs$reweightedMixturePriorCv) && + method %in% c("mvsusie_weights", "mvsusieWeights")) { + args$prior_variance <- cvArgs$reweightedMixturePriorCv[[j]] } - if (sd(methodPredictions) != 0) { - lmFit <- lm(actualValues ~ methodPredictions) - - # Calculate raw correlation and and adjusted R-squared - metricsTable[[m]][r, "corr"] <- cor(actualValues, methodPredictions) - - metricsTable[[m]][r, "rsq"] <- summary(lmFit)$r.squared - metricsTable[[m]][r, "adj_rsq"] <- summary(lmFit)$adj.r.squared - - # Calculate p-value - metricsTable[[m]][r, "pval"] <- summary(lmFit)$coefficients[2, 4] - - # Calculate RMSE - residuals <- actualValues - methodPredictions - metricsTable[[m]][r, "RMSE"] <- sqrt(mean(residuals^2)) - - # Calculate MAE - metricsTable[[m]][r, "MAE"] <- mean(abs(residuals)) + if (isTRUE(retainFits) && "retainFit" %in% names(formals(fnName))) { + args$retainFit <- TRUE + } + W <- if (verbose < 2) { + .quietEval(do.call(fnName, c(list(X = Xtr, Y = Ytr), args))) } else { - metricsTable[[m]][r, ] <- NA - if (verbose >= 1) message(paste0( - "Predicted values for condition ", r, " using ", m, - " have zero variance. Filling performance metric with NAs" - )) + do.call(fnName, c(list(X = Xtr, Y = Ytr), args)) } + capturedFit <- attr(W, "fit") + attr(W, "fit") <- NULL + rownames(W) <- colnames(Xtr) + } else { + Wcols <- lapply(seq_len(ncol(Ytr)), function(k) { + w <- if (verbose < 2) { + .quietEval(do.call(fnName, c(list(X = Xtr, y = Ytr[, k]), args))) + } else { + do.call(fnName, c(list(X = Xtr, y = Ytr[, k]), args)) + } + as.numeric(w) + }) + W <- do.call(cbind, Wcols) + rownames(W) <- colnames(Xtr) } + weights[[mk]] <- W + fits[[method]] <- capturedFit } - names(metricsTable) <- .renameSuffix(names(metricsTable), "performance") - # Per-fold retained fits (e.g. mr.mash), keyed [[fold]][[method]] with NULL - # entries dropped; empty per fold when retainFits = FALSE. Lets full-CV - # callers reuse each fold's fit as that fold's prior. - foldFits <- setNames(lapply(seq_along(foldResults), function(j) { - ff <- foldResults[[j]]$fits - ff[!vapply(ff, is.null, logical(1))] - }), paste0("fold_", seq_along(foldResults))) - return(list(samplePartition = samplePartition, prediction = Ypred, performance = metricsTable, foldFits = foldFits, timeElapsed = proc.time() - st)) + list(weights = weights, fits = fits) } + + # No weight methods: the caller only wants the fold partition. + if (is.null(weightMethods)) { + res <- .crossValidateWeights( + X, Y, fold = fold, samplePartitions = samplePartitions, + fitFold = function(Xtr, Ytr, j) list(weights = list(), fits = list()), + numThreads = numThreads, maxNumVariants = maxNumVariants, + variantsToKeep = variantsToKeep, retainFits = retainFits, verbose = verbose) + return(list(samplePartition = res$samplePartition)) + } + + .crossValidateWeights( + X, Y, fold = fold, samplePartitions = samplePartitions, + fitFold = weightFitFold, numThreads = numThreads, + maxNumVariants = maxNumVariants, variantsToKeep = variantsToKeep, + retainFits = retainFits, verbose = verbose) } #' Run multiple TWAS weight methods diff --git a/R/twasWeightsPipeline.R b/R/twasWeightsPipeline.R index 5c9c9c01..d338cf27 100644 --- a/R/twasWeightsPipeline.R +++ b/R/twasWeightsPipeline.R @@ -2,27 +2,50 @@ # Helpers + S4 dispatch surface for twasWeightsPipeline # ============================================================================= -# Concatenate two TwasWeights collections row-wise. `rbind` on DFrame -# subclasses does not reliably preserve the `ldSketch` slot, so this -# helper rebuilds via the constructor. Optional joint columns -# (`jointStudies`, `jointContexts`, `jointTraits`) are carried through -# via .combineJointCol() so a mixed rbind of joint + non-joint rows -# pads the non-joint side with NA_character_. +# Concatenate two TwasWeights collections row-wise, carrying forward every +# column (delegates to the generic `.rbindCollections`, which unions columns +# and pads a side lacking an optional column such as joint* / region). # @noRd .rbindTwasWeights <- function(a, b, ldSketch = NULL) { if (!is(a, "TwasWeights") || !is(b, "TwasWeights")) { stop(".rbindTwasWeights expects two TwasWeights inputs.") } - TwasWeights( - study = c(as.character(a$study), as.character(b$study)), - context = c(as.character(a$context), as.character(b$context)), - trait = c(as.character(a$trait), as.character(b$trait)), - method = c(as.character(a$method), as.character(b$method)), - entry = c(as.list(a$entry), as.list(b$entry)), - jointStudies = .combineJointCol(a, b, "jointStudies"), - jointContexts = .combineJointCol(a, b, "jointContexts"), - jointTraits = .combineJointCol(a, b, "jointTraits"), - ldSketch = ldSketch) + # Carry forward every column (joint*, region, ...) via the generic combine. + .rbindCollections(list(a, b), ldSketch = ldSketch) +} + +# Normalize combine() varargs: accept either N objects or a single list of +# them; drop NULLs; require at least one input of the expected class `cls`. +.asCombineList <- function(parts, cls, fn) { + if (length(parts) == 1L && is.list(parts[[1L]]) && + !methods::is(parts[[1L]], cls)) { + parts <- parts[[1L]] + } + parts <- parts[!vapply(parts, is.null, logical(1L))] + if (length(parts) == 0L) + stop(fn, ": nothing to combine (need at least one ", cls, ").") + if (!all(vapply(parts, function(p) methods::is(p, cls), logical(1L)))) + stop(fn, ": every input must be a ", cls, ".") + parts +} + +#' Combine TwasWeights collections +#' +#' Row-bind two or more \code{\link{TwasWeights}} collections into one -- e.g. +#' assembling per-gene weight sets into a single per-region collection for +#' cTWAS. Joint-specification metadata columns are carried through. +#' +#' @param ... Two or more \code{TwasWeights} objects, or a single \code{list} +#' of them. +#' @param ldSketch Optional \code{\link{GenotypeHandle}} to attach to the +#' combined collection. Default \code{NULL}. Applied when combining two or +#' more inputs; a single input is returned unchanged. +#' @return A single combined \code{TwasWeights}. +#' @seealso \code{\link{combineFineMappingResults}} +#' @export +combineTwasWeights <- function(..., ldSketch = NULL) { + parts <- .asCombineList(list(...), "TwasWeights", "combineTwasWeights") + Reduce(function(a, b) .rbindTwasWeights(a, b, ldSketch = ldSketch), parts) } # --- Multi-region (jointRegions) helpers for the QtlDataset method ---------- @@ -716,20 +739,20 @@ #' @param traitId Optional character vector of trait identifiers to #' restrict processing to (QtlDataset / QtlSumStats inputs). Default #' \code{NULL}. -#' @param region Optional \code{GRanges} for QtlDataset trait selection. -#' Mutually exclusive with \code{traitId}. +#' @param region Optional variant window for QtlDataset trait selection: a +#' \code{GRanges}, a \code{"chr:start-end"} string, or a one-row data.frame +#' with \code{chrom}/\code{start}/\code{end}. Mutually exclusive with +#' \code{traitId}. #' @param cisWindow For QtlDataset: cis-window (bp) around each trait's #' genomic position when extracting variants. Required when #' \code{traitId} is supplied. Mutually exclusive with \code{region}. -#' @param minTwasMaf For QtlDataset: optional minimum minor-allele frequency -#' applied to the variant set used for TWAS weight learning, on top of the -#' dataset's construct-time \code{mafCutoff} (the effective cutoff is the -#' larger of the two). Lets the TWAS pass use a stricter MAF threshold than -#' fine mapping. \code{NULL} (default) leaves the construct-time cutoff in -#' place. -#' @param minTwasXvar As \code{minTwasMaf} but for the per-variant genotype -#' variance cutoff (\code{xvarCutoff}). \code{NULL} (default) leaves the -#' construct-time cutoff in place. +#' @param mafCutoff,macCutoff,xvarCutoff,imissCutoff,keepIndel,keepSamples,keepVariants +#' For QtlDataset: optional per-call genotype-filter overrides. Each replaces +#' the corresponding construct-time \code{\link{QtlDataset}} slot for this +#' call only (applied to a validated copy); \code{NULL} (default) leaves the +#' stored value in place. Variant QC is a property of the data, so these are +#' applied identically here and in \code{\link{fineMappingPipeline}} -- there +#' is deliberately no TWAS-specific variant filter. #' @param jointRegions For QtlDataset with a multi-range \code{region}: #' \code{FALSE} (default) learns weights for each range independently and #' concatenates them into one entry per (study, context, trait, method); @@ -815,8 +838,16 @@ setMethod("twasWeightsPipeline", "QtlDataset", traitId = NULL, region = NULL, cisWindow = NULL, - minTwasMaf = NULL, - minTwasXvar = NULL, + # Per-call genotype-filter overrides; NULL = use the QtlDataset's + # construct-time slot value. Variant QC is a data property, so these + # match fineMappingPipeline exactly -- there is no TWAS-specific filter. + mafCutoff = NULL, + macCutoff = NULL, + xvarCutoff = NULL, + imissCutoff = NULL, + keepIndel = NULL, + keepSamples = NULL, + keepVariants = NULL, jointRegions = FALSE, jointSpecification = NULL, fineMappingResult = NULL, @@ -853,14 +884,12 @@ setMethod("twasWeightsPipeline", "QtlDataset", "coordinates, whereas `region` is the literal variant window.") } xRegions <- .makeXRegions(region, jointRegions) - # TWAS-specific variant filters: tighten the QtlDataset maf/xvar cutoffs for - # weight learning (distinct from the construct-time / fine-mapping cutoffs). - # Modifying the local `data` copy elevates them everywhere downstream - # (runOne, runMultivariate, and the jointSpec dispatcher all extract from it). - if (!is.null(minTwasMaf)) - data@mafCutoff <- max(data@mafCutoff, as.numeric(minTwasMaf)) - if (!is.null(minTwasXvar)) - data@xvarCutoff <- max(data@xvarCutoff, as.numeric(minTwasXvar)) + # Per-call filter overrides replace the construct-time slot values on a + # validated copy. Variant QC is a data property applied identically to + # fine-mapping and TWAS -- there is no TWAS-specific variant filter. + data <- .qtlApplyFilterOverrides(data, mafCutoff, macCutoff, xvarCutoff, + imissCutoff, keepIndel, keepSamples, + keepVariants) parsedJointSpec <- parseJointSpecification(jointSpecification, data) norm <- .twasNormalizeMethods(methods) .twasCheckMethodCapabilities(norm$tokens, "QtlDataset") diff --git a/R/variantId.R b/R/variantId.R index de220d75..74f3c4a5 100644 --- a/R/variantId.R +++ b/R/variantId.R @@ -675,7 +675,7 @@ regionToDf <- function(ldRegionId, colnames = c("chrom", "start", "end")) { #' @param regions A region string, character vector, or data.frame with #' chrom/start/end columns. #' @return A \code{GRanges} object. -#' @noRd +#' @export asGranges <- function(regions) { if (is.character(regions)) { df <- regionToDf(regions) diff --git a/man/CtwasResult.Rd b/man/CtwasResult.Rd new file mode 100644 index 00000000..38786c64 --- /dev/null +++ b/man/CtwasResult.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CtwasResult.R +\name{CtwasResult} +\alias{CtwasResult} +\title{Create a CtwasResult Collection} +\usage{ +CtwasResult( + gwasStudy, + study, + context, + method, + entry, + jointStudies = NULL, + jointContexts = NULL +) +} +\arguments{ +\item{gwasStudy}{Character vector of GWAS (disease) study identifiers.} + +\item{study}{Character vector of QTL study identifiers.} + +\item{context}{Character vector of context labels (scalar per row).} + +\item{method}{Character vector of weight-method names.} + +\item{entry}{List / \code{SimpleList} of \code{CtwasResultEntry} objects.} + +\item{jointStudies}{Optional character vector (length \code{length(gwasStudy)}) +listing the joined studies for a multi-study run, \code{NA_character_} +otherwise. \code{NULL} (default) omits the column.} + +\item{jointContexts}{Optional character vector for multi-context runs. Same +shape as \code{jointStudies}.} +} +\value{ +A \code{CtwasResult} object. +} +\description{ +Construct a \code{CtwasResult} DFrame-subclass collection of + cTWAS runs from per-tuple vectors and a list of + \code{\link{CtwasResultEntry}} payloads (one per run). A single-context run + is one row; a multi-context run emits per-context rows sharing the same + \code{jointContexts} set and jointly-estimated group priors. +} diff --git a/man/CtwasResultEntry.Rd b/man/CtwasResultEntry.Rd new file mode 100644 index 00000000..95e6728b --- /dev/null +++ b/man/CtwasResultEntry.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CtwasResultEntry.R +\name{CtwasResultEntry} +\alias{CtwasResultEntry} +\title{Create a CtwasResultEntry} +\usage{ +CtwasResultEntry( + finemap = NULL, + susieAlpha = NULL, + param = NULL, + regionInfo = NULL +) +} +\arguments{ +\item{finemap}{The per-gene (and, when SNPs are retained, per-SNP) posterior +summary table (\code{ctwas::finemap_regions} \code{finemap_res} shape), or +\code{NULL}.} + +\item{susieAlpha}{The per-effect susie alpha table +(\code{ctwas::finemap_regions} \code{susie_alpha_res} shape) — the fuller +cTWAS output retained so the raw run is reconstructable, or \code{NULL}.} + +\item{param}{The estimated \code{group_prior} / \code{group_prior_var} for +this run, or \code{NULL}.} + +\item{regionInfo}{Per-region metadata, or \code{NULL}.} +} +\value{ +A \code{CtwasResultEntry} object. +} +\description{ +Per-run cTWAS payload wrapping the fine-mapping posterior table, + the full per-effect susie alpha table, the estimated group prior(s), and + region metadata. Held in every row of a \code{\link{CtwasResult}} + collection. +} diff --git a/man/GwasFineMappingResult.Rd b/man/GwasFineMappingResult.Rd index 989589d3..f5c7c0b3 100644 --- a/man/GwasFineMappingResult.Rd +++ b/man/GwasFineMappingResult.Rd @@ -4,7 +4,14 @@ \alias{GwasFineMappingResult} \title{TWAS Weights Collection} \usage{ -GwasFineMappingResult(study, method, entry, region_id = NULL, ldSketch = NULL) +GwasFineMappingResult( + study, + method, + entry, + region_id = NULL, + region = NULL, + ldSketch = NULL +) } \arguments{ \item{study}{Character vector of study identifiers (per tuple).} @@ -20,6 +27,13 @@ id (\code{"region_1"}, \code{"region_2"}, ...) so the Supplying meaningful labels (e.g. \code{"chr22_10516173_17414263"}) is preferred for downstream consumers that join on region.} +\item{region}{Optional \code{GRanges} (length \code{length(study)}) genomic +anchor per row. When \code{NULL} (default) it is derived from +\code{region_id} (parsing \code{chr_start_end} / \code{chr:start-end}; a +0-width \code{chrUn} sentinel for ids that do not encode coordinates). +Carried forward as provenance (e.g. for cTWAS LD-block placement); not +part of the identity key.} + \item{ldSketch}{An optional \code{GenotypeHandle}.} } \value{ diff --git a/man/QtlFineMappingResult.Rd b/man/QtlFineMappingResult.Rd index fa2ca6a8..76e19286 100644 --- a/man/QtlFineMappingResult.Rd +++ b/man/QtlFineMappingResult.Rd @@ -13,6 +13,7 @@ QtlFineMappingResult( jointStudies = NULL, jointContexts = NULL, jointTraits = NULL, + region = NULL, ldSketch = NULL ) } @@ -42,6 +43,11 @@ Same shape as \code{jointStudies}.} \item{jointTraits}{Optional character vector for cross-trait joints. Same shape as \code{jointStudies}.} +\item{region}{Optional \code{GRanges} (length \code{length(study)}) giving +the genomic anchor of each row's trait (its own coordinates). Carried +forward as provenance (e.g. for cTWAS LD-block placement); not part of the +identity key. \code{NULL} (default) omits the column.} + \item{ldSketch}{An optional \code{GenotypeHandle} (the LD reference for RSS-derived fits), or \code{NULL} for individual-level fits.} } diff --git a/man/TwasWeights.Rd b/man/TwasWeights.Rd index 343bc918..a699d2ad 100644 --- a/man/TwasWeights.Rd +++ b/man/TwasWeights.Rd @@ -13,6 +13,7 @@ TwasWeights( jointStudies = NULL, jointContexts = NULL, jointTraits = NULL, + region = NULL, ldSketch = NULL ) } @@ -41,6 +42,13 @@ Same shape as \code{jointStudies}.} \item{jointTraits}{Optional character vector for cross-trait joints. Same shape as \code{jointStudies}.} +\item{region}{Optional \code{GRanges} (length \code{length(study)}) giving +the genomic anchor of each row's trait -- the trait's own coordinates when +built from a \code{QtlDataset}, or the summary-stat variant-span window +when built from a \code{QtlSumStats}. Carried forward as provenance (e.g. +for cTWAS LD-block placement); not part of the identity key. \code{NULL} +(default) omits the column.} + \item{ldSketch}{An optional \code{GenotypeHandle}, or \code{NULL} for individual-level fits.} } diff --git a/man/asCtwasResult.Rd b/man/asCtwasResult.Rd new file mode 100644 index 00000000..a3a9b155 --- /dev/null +++ b/man/asCtwasResult.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ctwasPipeline.R +\name{asCtwasResult} +\alias{asCtwasResult} +\title{Structure a granular cTWAS finemap result as a CtwasResult} +\usage{ +asCtwasResult(finemapResult, keepSnps = FALSE) +} +\arguments{ +\item{finemapResult}{A list from \code{\link{finemapCtwasRegions}} or +\code{\link{mergeCtwasBoundaryRegions}}.} + +\item{keepSnps}{Logical (length 1). Retain the context-agnostic SNP +background as one extra \code{study = context = "SNP"} row. Default +\code{FALSE}. See \code{\link{ctwasPipeline}}.} +} +\value{ +A \code{\link{CtwasResult}}. +} +\description{ +Decompose the raw list returned by + \code{\link{finemapCtwasRegions}} (optionally after + \code{\link{mergeCtwasBoundaryRegions}}) into the structured, per-(study, + context) \code{\link{CtwasResult}} — the same decomposition + \code{\link{ctwasPipeline}} applies to its one-shot output, exposed for the + granular \code{assembleCtwasInputs} \eqn{\to} \code{estCtwasParam} \eqn{\to} + \code{screenCtwasRegions} \eqn{\to} \code{finemapCtwasRegions} path. The + GWAS study is read from \code{z_snp$study} and the (single) weight method + from the gene ids. +} +\seealso{ +\code{\link{ctwasPipeline}}, \code{\link{finemapCtwasRegions}} +} diff --git a/man/asGranges.Rd b/man/asGranges.Rd new file mode 100644 index 00000000..2b98d404 --- /dev/null +++ b/man/asGranges.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/variantId.R +\name{asGranges} +\alias{asGranges} +\title{Convert region specifications to a GRanges object} +\usage{ +asGranges(regions) +} +\arguments{ +\item{regions}{A region string, character vector, or data.frame with +chrom/start/end columns.} +} +\value{ +A \code{GRanges} object. +} +\description{ +Accepts region strings ("chr1:100-200", "1_100_200"), character vectors of +such strings, or data.frames with chrom/start/end columns. Returns a +\code{\link[GenomicRanges]{GRanges}} object. +} diff --git a/man/assembleCtwasInputs.Rd b/man/assembleCtwasInputs.Rd index c5bdfc4f..68df5a6d 100644 --- a/man/assembleCtwasInputs.Rd +++ b/man/assembleCtwasInputs.Rd @@ -21,12 +21,19 @@ assembleCtwasInputs( by \code{region_id} (at least two entries). Each must have \code{getQcInfo()} non-empty.} -\item{twasWeights}{NAMED LIST of \code{\link{TwasWeights}} keyed by -\code{region_id}. Keys must be a SUBSET of \code{gwasSumStats}'s -keys: blocks without any TWAS weights still contribute their -SNP-level signal to ctwas's joint group prior estimate (matches -the legacy whole-chromosome pattern where only a few of many LD -blocks carry gene weights).} +\item{twasWeights}{The per-gene weight source. Either (a) a FLAT +\code{\link{TwasWeights}} / \code{QtlFineMappingResult} (or a +homogeneous list of them) carrying \code{region} provenance — each +gene is placed into its home LD block internally by +\code{start(region)} (matching cTWAS's \code{p0} assignment rule); +or (b) a pre-bucketed NAMED LIST keyed by \code{region_id} (keys a +SUBSET of \code{gwasSumStats}'s), used as-is. Blocks without any +TWAS weights still contribute their SNP-level signal to ctwas's +joint group prior estimate (the legacy whole-chromosome pattern +where only a few of many LD blocks carry gene weights). A gene whose +cis span straddles a block boundary is homed by its single anchor; +the cross-block signal is cTWAS's boundary-gene concern +(\code{\link{mergeCtwasBoundaryRegions}}), not placement.} \item{twasZ}{Optional \code{GRanges} of TWAS Z-scores (output of \code{\link{causalInferencePipeline}}). When supplied, the @@ -44,10 +51,11 @@ credible-set membership data used by the CS / PIP rescue filters \item{method}{Optional character (length 1). Picks which TWAS method's weights to feed into ctwas for each (study, context, trait) gene. When \code{NULL} (default): use \code{"ensemble"} if -that method is present across the TwasWeights; otherwise use the -sole method when only one is present; otherwise error. Passing -the name explicitly (e.g. \code{"mrash"}) overrides the default -resolution.} +that method is present across the weight sources; otherwise use the +sole method when only one is present; otherwise run \strong{every} +method as an independent cTWAS run (one \code{CtwasResult} row-set +per method). Passing the name explicitly (e.g. \code{"mrash"}) +restricts the run to that single method.} \item{twasWeightCutoff}{Numeric (length 1). Drop variants with \code{|weight| < twasWeightCutoff} from each gene's weight matrix diff --git a/man/combineFineMappingResults.Rd b/man/combineFineMappingResults.Rd new file mode 100644 index 00000000..2d6e570f --- /dev/null +++ b/man/combineFineMappingResults.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fineMappingPipeline.R +\name{combineFineMappingResults} +\alias{combineFineMappingResults} +\title{Combine FineMappingResult collections} +\usage{ +combineFineMappingResults(..., ldSketch = NULL) +} +\arguments{ +\item{...}{Two or more \code{FineMappingResultBase} objects, or a single +\code{list} of them.} + +\item{ldSketch}{Optional \code{\link{GenotypeHandle}} to attach to the +combined collection. Default \code{NULL}. Applied when combining two or +more inputs; a single input is returned unchanged.} +} +\value{ +A single combined fine-mapping result of the shared concrete class. +} +\description{ +Row-bind two or more fine-mapping result collections of the SAME concrete +class (all \code{\link{QtlFineMappingResult}} or all +\code{\link{GwasFineMappingResult}}) into one -- e.g. per-block GWAS results +into a genome-wide collection for cTWAS. Mixing the two concrete classes is +an error. +} +\seealso{ +\code{\link{combineTwasWeights}} +} diff --git a/man/combineTwasWeights.Rd b/man/combineTwasWeights.Rd new file mode 100644 index 00000000..957488e0 --- /dev/null +++ b/man/combineTwasWeights.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/twasWeightsPipeline.R +\name{combineTwasWeights} +\alias{combineTwasWeights} +\title{Combine TwasWeights collections} +\usage{ +combineTwasWeights(..., ldSketch = NULL) +} +\arguments{ +\item{...}{Two or more \code{TwasWeights} objects, or a single \code{list} +of them.} + +\item{ldSketch}{Optional \code{\link{GenotypeHandle}} to attach to the +combined collection. Default \code{NULL}. Applied when combining two or +more inputs; a single input is returned unchanged.} +} +\value{ +A single combined \code{TwasWeights}. +} +\description{ +Row-bind two or more \code{\link{TwasWeights}} collections into one -- e.g. +assembling per-gene weight sets into a single per-region collection for +cTWAS. Joint-specification metadata columns are carried through. +} +\seealso{ +\code{\link{combineFineMappingResults}} +} diff --git a/man/ctwasPipeline.Rd b/man/ctwasPipeline.Rd index f6c9f68d..9ffb134a 100644 --- a/man/ctwasPipeline.Rd +++ b/man/ctwasPipeline.Rd @@ -22,6 +22,11 @@ ctwasPipeline( minPipCutoff = 0, maxNumVariants = Inf, fallbackToPrefit = FALSE, + keepSnps = FALSE, + mergeBoundary = FALSE, + mergePipThresh = 0.5, + mergeFilterCs = FALSE, + mergeMaxSNP = Inf, ... ) } @@ -30,12 +35,19 @@ ctwasPipeline( by \code{region_id} (at least two entries). Each must have \code{getQcInfo()} non-empty.} -\item{twasWeights}{NAMED LIST of \code{\link{TwasWeights}} keyed by -\code{region_id}. Keys must be a SUBSET of \code{gwasSumStats}'s -keys: blocks without any TWAS weights still contribute their -SNP-level signal to ctwas's joint group prior estimate (matches -the legacy whole-chromosome pattern where only a few of many LD -blocks carry gene weights).} +\item{twasWeights}{The per-gene weight source. Either (a) a FLAT +\code{\link{TwasWeights}} / \code{QtlFineMappingResult} (or a +homogeneous list of them) carrying \code{region} provenance — each +gene is placed into its home LD block internally by +\code{start(region)} (matching cTWAS's \code{p0} assignment rule); +or (b) a pre-bucketed NAMED LIST keyed by \code{region_id} (keys a +SUBSET of \code{gwasSumStats}'s), used as-is. Blocks without any +TWAS weights still contribute their SNP-level signal to ctwas's +joint group prior estimate (the legacy whole-chromosome pattern +where only a few of many LD blocks carry gene weights). A gene whose +cis span straddles a block boundary is homed by its single anchor; +the cross-block signal is cTWAS's boundary-gene concern +(\code{\link{mergeCtwasBoundaryRegions}}), not placement.} \item{twasZ}{Optional \code{GRanges} of TWAS Z-scores (output of \code{\link{causalInferencePipeline}}). When supplied, the @@ -53,10 +65,11 @@ credible-set membership data used by the CS / PIP rescue filters \item{method}{Optional character (length 1). Picks which TWAS method's weights to feed into ctwas for each (study, context, trait) gene. When \code{NULL} (default): use \code{"ensemble"} if -that method is present across the TwasWeights; otherwise use the -sole method when only one is present; otherwise error. Passing -the name explicitly (e.g. \code{"mrash"}) overrides the default -resolution.} +that method is present across the weight sources; otherwise use the +sole method when only one is present; otherwise run \strong{every} +method as an independent cTWAS run (one \code{CtwasResult} row-set +per method). Passing the name explicitly (e.g. \code{"mrash"}) +restricts the run to that single method.} \item{thin, niterPrefit, niter, L}{Pass-throughs to \code{ctwas::ctwas_sumstats}.} @@ -94,12 +107,49 @@ NaN failure is recovered by falling back to the prefit estimates (mirrors the legacy ctwas_2 workaround on underpowered data). Default \code{FALSE}.} +\item{keepSnps}{Logical (length 1). When \code{TRUE}, retain the +context-agnostic SNP background of each run as one extra +\code{CtwasResult} row (\code{study = context = "SNP"}, mirroring +cTWAS's own \code{"SNP"} group) so the full ctwas output is +reconstructable from \code{\link{getFinemap}} / \code{getSusieAlpha}. +Default \code{FALSE} — the SNP rows are the null background and are +dropped from the structured gene-level result.} + +\item{mergeBoundary}{Logical (length 1). When \code{TRUE}, run +\code{\link{mergeCtwasBoundaryRegions}} after fine-mapping each +run: a high-PIP gene whose cis window straddles an LD-block boundary +has its adjacent regions merged and re-fine-mapped (the legacy +default-off \code{ctwas_3} post-processing). Default \code{FALSE}.} + +\item{mergePipThresh}{Numeric (length 1). PIP threshold for selecting +which boundary genes to merge (\code{\link{mergeCtwasBoundaryRegions}} +\code{pipThresh}). Default \code{0.5}. Ignored unless +\code{mergeBoundary = TRUE}.} + +\item{mergeFilterCs}{Logical (length 1). Require the boundary gene to +be in a credible set to be selected. Default \code{FALSE}. Ignored +unless \code{mergeBoundary = TRUE}.} + +\item{mergeMaxSNP}{Numeric (length 1). Per-merged-region SNP cap. +Default \code{Inf}. Ignored unless \code{mergeBoundary = TRUE}.} + \item{...}{Additional arguments forwarded to \code{ctwas::ctwas_sumstats}.} } \value{ -Whatever \code{ctwas::ctwas_sumstats} returns (a list with - \code{susie_alpha_res}, \code{param}, and other diagnostics). +A \code{\link{CtwasResult}} collection: one row per + \code{(gwasStudy, study, context, method)}. A single-context run is + one row per method; a multi-context (joint) run emits per-context + rows sharing the same \code{jointContexts} set and the + jointly-estimated group priors. Each row's + \code{\link{CtwasResultEntry}} payload carries that context's + per-gene fine-mapping posteriors (\code{finemap}), the run's + \code{param}, and its \code{regionInfo}. For the raw + \code{ctwas::finemap_regions} list (e.g. to feed + \code{\link{mergeCtwasBoundaryRegions}}), call the granular + \code{\link{assembleCtwasInputs}} \eqn{\to} \code{\link{estCtwasParam}} + \eqn{\to} \code{\link{screenCtwasRegions}} \eqn{\to} + \code{\link{finemapCtwasRegions}} path instead. } \description{ Pipeline that hands a per-block set of diff --git a/man/fineMappingPipeline.Rd b/man/fineMappingPipeline.Rd index a08bfea8..771e95ba 100644 --- a/man/fineMappingPipeline.Rd +++ b/man/fineMappingPipeline.Rd @@ -18,6 +18,13 @@ fineMappingPipeline(data, ...) traitId = NULL, region = NULL, cisWindow = NULL, + mafCutoff = NULL, + macCutoff = NULL, + xvarCutoff = NULL, + imissCutoff = NULL, + keepIndel = NULL, + keepSamples = NULL, + keepVariants = NULL, jointRegions = FALSE, jointSpecification = NULL, addSusieInf = TRUE, @@ -30,6 +37,7 @@ fineMappingPipeline(data, ...) medianAbsCorr = NULL, fineMappingResult = NULL, cvFolds = 0, + cvThreads = 1, samplePartition = NULL, pipCutoffToSkip = 0, usePCA = FALSE, @@ -66,6 +74,7 @@ fineMappingPipeline(data, ...) twasWeights = NULL, dataDrivenPriorWeightsCutoff = 1e-10, cvFolds = 0, + cvThreads = 1, samplePartition = NULL, pipCutoffToSkip = 0, seed = NULL, @@ -148,8 +157,10 @@ fineMappingPipeline(data, ...) \item{traitId}{Optional character vector of trait names to restrict processing to.} -\item{region}{Optional \code{GRanges} for QtlDataset trait -selection. Mutually exclusive with \code{traitId}.} +\item{region}{Optional variant window for QtlDataset trait selection: a +\code{GRanges}, a \code{"chr:start-end"} string, or a one-row data.frame +with \code{chrom}/\code{start}/\code{end}. Mutually exclusive with +\code{traitId}.} \item{cisWindow}{For QtlDataset: cis-window (bp) around each trait's genomic position when extracting variants. Required when @@ -209,6 +220,11 @@ this partition and feeds these predictions into the SR-TWAS ensemble. Individual-level (\code{QtlDataset} / \code{MultiStudyQtlDataset}) input only; ignored for sumstat inputs.} +\item{cvThreads}{Integer. Number of parallel workers for the +cross-validation fold refits (passed to the shared CV engine's +\code{numThreads}). Default \code{1} (serial); \code{-1} uses all cores. +Only consulted when \code{cvFolds > 1}.} + \item{samplePartition}{Optional pre-defined CV partition \code{data.frame} with columns \code{Sample} and \code{Fold}. When supplied (and \code{cvFolds > 1}), every method reuses this exact diff --git a/man/getContexts.Rd b/man/getContexts.Rd index a497ca47..f9a3247a 100644 --- a/man/getContexts.Rd +++ b/man/getContexts.Rd @@ -1,8 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AllGenerics.R, R/GwasFineMappingResult.R, -% R/QtlDataset.R, R/QtlFineMappingResult.R, R/qtlSumStats.R, R/twasWeights.R +% Please edit documentation in R/AllGenerics.R, R/CtwasResult.R, +% R/GwasFineMappingResult.R, R/QtlDataset.R, R/QtlFineMappingResult.R, +% R/qtlSumStats.R, R/twasWeights.R \name{getContexts} \alias{getContexts} +\alias{getContexts,CtwasResult-method} \alias{getContexts,GwasFineMappingResult-method} \alias{getContexts,QtlDataset-method} \alias{getContexts,QtlFineMappingResult-method} @@ -12,6 +14,8 @@ \usage{ getContexts(x) +\S4method{getContexts}{CtwasResult}(x) + \S4method{getContexts}{GwasFineMappingResult}(x) \S4method{getContexts}{QtlDataset}(x) diff --git a/man/getCtwasParam.Rd b/man/getCtwasParam.Rd new file mode 100644 index 00000000..671b8f10 --- /dev/null +++ b/man/getCtwasParam.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/CtwasResultEntry.R +\name{getCtwasParam} +\alias{getCtwasParam} +\alias{getCtwasParam,CtwasResultEntry-method} +\title{Get cTWAS Group Prior Parameters} +\usage{ +getCtwasParam(x, ...) + +\S4method{getCtwasParam}{CtwasResultEntry}(x, ...) +} +\arguments{ +\item{x}{A \code{CtwasResultEntry}.} + +\item{...}{Reserved for future use.} +} +\value{ +The stored parameter object (typically a list), or \code{NULL}. +} +\description{ +Return the estimated \code{group_prior} / \code{group_prior_var} + of a \code{\link{CtwasResultEntry}}. +} diff --git a/man/getFinemap.Rd b/man/getFinemap.Rd new file mode 100644 index 00000000..329cf79d --- /dev/null +++ b/man/getFinemap.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/CtwasResult.R, +% R/CtwasResultEntry.R +\name{getFinemap} +\alias{getFinemap} +\alias{getFinemap,CtwasResult-method} +\alias{getFinemap,CtwasResultEntry-method} +\title{Get cTWAS Fine-mapping Posteriors} +\usage{ +getFinemap(x, ...) + +\S4method{getFinemap}{CtwasResult}(x, ...) + +\S4method{getFinemap}{CtwasResultEntry}(x, ...) +} +\arguments{ +\item{x}{A \code{CtwasResultEntry} or \code{CtwasResult}.} + +\item{...}{Class-specific selection arguments.} +} +\value{ +A \code{data.frame} of posteriors (or \code{NULL} when absent). +} +\description{ +Return the per-gene (and per-SNP) fine-mapping posterior table + from a \code{\link{CtwasResultEntry}} or a \code{\link{CtwasResult}} + collection (aggregated across rows, tagged with run identity). +} diff --git a/man/getFits.Rd b/man/getFits.Rd index b584b0cd..6a8ae3d7 100644 --- a/man/getFits.Rd +++ b/man/getFits.Rd @@ -1,14 +1,17 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AllGenerics.R, R/TwasWeightsEntry.R, -% R/twasWeights.R +% Please edit documentation in R/AllGenerics.R, R/FineMappingEntry.R, +% R/TwasWeightsEntry.R, R/twasWeights.R \name{getFits} \alias{getFits} +\alias{getFits,FineMappingEntry-method} \alias{getFits,TwasWeightsEntry-method} \alias{getFits,TwasWeights-method} \title{Get Model Fits} \usage{ getFits(x, ...) +\S4method{getFits}{FineMappingEntry}(x, ...) + \S4method{getFits}{TwasWeightsEntry}(x, ...) \S4method{getFits}{TwasWeights}(x, study = NULL, context = NULL, trait = NULL, method = NULL, ...) diff --git a/man/getMethodNames.Rd b/man/getMethodNames.Rd index fac009c4..33762955 100644 --- a/man/getMethodNames.Rd +++ b/man/getMethodNames.Rd @@ -1,8 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AllGenerics.R, R/AllClasses.R, R/twasWeights.R +% Please edit documentation in R/AllGenerics.R, R/AllClasses.R, +% R/CtwasResult.R, R/twasWeights.R \name{getMethodNames} \alias{getMethodNames} \alias{getMethodNames,FineMappingResultBase-method} +\alias{getMethodNames,CtwasResult-method} \alias{getMethodNames,TwasWeights-method} \title{Get Method Names} \usage{ @@ -10,6 +12,8 @@ getMethodNames(x) \S4method{getMethodNames}{FineMappingResultBase}(x) +\S4method{getMethodNames}{CtwasResult}(x) + \S4method{getMethodNames}{TwasWeights}(x) } \arguments{ diff --git a/man/getRegion.Rd b/man/getRegion.Rd new file mode 100644 index 00000000..311da83a --- /dev/null +++ b/man/getRegion.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/AllClasses.R, R/twasWeights.R +\name{getRegion} +\alias{getRegion} +\alias{getRegion,FineMappingResultBase-method} +\alias{getRegion,TwasWeights-method} +\title{Get Per-Row Genomic Regions} +\usage{ +getRegion(x, ...) + +\S4method{getRegion}{FineMappingResultBase}(x, ...) + +\S4method{getRegion}{TwasWeights}(x, ...) +} +\arguments{ +\item{x}{The object.} + +\item{...}{Reserved for future use.} +} +\value{ +A \code{GRanges} with one range per row of \code{x}, or an empty + \code{GRanges} when the collection carries no region provenance. +} +\description{ +Return the genomic anchor of each row of a per-tuple collection + as a \code{GRanges} with one range per row -- the trait's own region for a + \code{TwasWeights}, or the fine-mapping window for a \code{FineMappingResult}. + This is location provenance (e.g. for cTWAS LD-block placement). +} diff --git a/man/getStandardized.Rd b/man/getStandardized.Rd index 63fe1d6d..6e425088 100644 --- a/man/getStandardized.Rd +++ b/man/getStandardized.Rd @@ -1,14 +1,17 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AllGenerics.R, R/TwasWeightsEntry.R, -% R/twasWeights.R +% Please edit documentation in R/AllGenerics.R, R/FineMappingEntry.R, +% R/TwasWeightsEntry.R, R/twasWeights.R \name{getStandardized} \alias{getStandardized} +\alias{getStandardized,FineMappingEntry-method} \alias{getStandardized,TwasWeightsEntry-method} \alias{getStandardized,TwasWeights-method} \title{Get Standardized Flag} \usage{ getStandardized(x, ...) +\S4method{getStandardized}{FineMappingEntry}(x, ...) + \S4method{getStandardized}{TwasWeightsEntry}(x, ...) \S4method{getStandardized}{TwasWeights}(x, study = NULL, context = NULL, trait = NULL, method = NULL) diff --git a/man/getStudy.Rd b/man/getStudy.Rd index 6cac6089..9b35d24b 100644 --- a/man/getStudy.Rd +++ b/man/getStudy.Rd @@ -1,10 +1,11 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/AllGenerics.R, R/AllClasses.R, R/QtlDataset.R, -% R/MultiStudyQtlDataset.R, R/twasWeights.R +% Please edit documentation in R/AllGenerics.R, R/AllClasses.R, +% R/CtwasResult.R, R/QtlDataset.R, R/MultiStudyQtlDataset.R, R/twasWeights.R \name{getStudy} \alias{getStudy} \alias{getStudy,SumStatsBase-method} \alias{getStudy,FineMappingResultBase-method} +\alias{getStudy,CtwasResult-method} \alias{getStudy,QtlDataset-method} \alias{getStudy,MultiStudyQtlDataset-method} \alias{getStudy,TwasWeights-method} @@ -16,6 +17,8 @@ getStudy(x) \S4method{getStudy}{FineMappingResultBase}(x) +\S4method{getStudy}{CtwasResult}(x) + \S4method{getStudy}{QtlDataset}(x) \S4method{getStudy}{MultiStudyQtlDataset}(x) diff --git a/man/getSusieAlpha.Rd b/man/getSusieAlpha.Rd new file mode 100644 index 00000000..9774600a --- /dev/null +++ b/man/getSusieAlpha.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/CtwasResult.R, +% R/CtwasResultEntry.R +\name{getSusieAlpha} +\alias{getSusieAlpha} +\alias{getSusieAlpha,CtwasResult-method} +\alias{getSusieAlpha,CtwasResultEntry-method} +\title{Get cTWAS Per-effect Susie Alpha Table} +\usage{ +getSusieAlpha(x, ...) + +\S4method{getSusieAlpha}{CtwasResult}(x, ...) + +\S4method{getSusieAlpha}{CtwasResultEntry}(x, ...) +} +\arguments{ +\item{x}{A \code{CtwasResultEntry} or \code{CtwasResult}.} + +\item{...}{Class-specific selection arguments.} +} +\value{ +A \code{data.frame} of per-effect alphas (or \code{NULL} when absent). +} +\description{ +Return the full per-effect susie alpha table + (\code{ctwas::finemap_regions} \code{susie_alpha_res} shape) from a + \code{\link{CtwasResultEntry}} or a \code{\link{CtwasResult}} collection + (aggregated across rows, tagged with run identity). This is the fuller + cTWAS output retained so the raw run is reconstructable. +} diff --git a/man/qtlSumStatsFromZMatrix.Rd b/man/qtlSumStatsFromZMatrix.Rd new file mode 100644 index 00000000..4bc1d80b --- /dev/null +++ b/man/qtlSumStatsFromZMatrix.Rd @@ -0,0 +1,62 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mashWrapper.R +\name{qtlSumStatsFromZMatrix} +\alias{qtlSumStatsFromZMatrix} +\title{Build a QtlSumStats from a Z-score matrix} +\usage{ +qtlSumStatsFromZMatrix( + z, + study, + ldSketch, + context = colnames(z), + trait = "mash", + genome = "GRCh38", + n = 1000L, + a1 = "A", + a2 = "G", + role = "mash" +) +} +\arguments{ +\item{z}{Numeric matrix (variants x conditions). \code{rownames(z)} are +variant ids (ideally \code{chr:pos:A2:A1}); \code{colnames(z)} label the +conditions.} + +\item{study}{Study identifier (recycled across conditions).} + +\item{ldSketch}{A \code{\link{GenotypeHandle}} embedded in the collection.} + +\item{context}{Condition context label(s): a single value recycled across +every column, or a length-\code{ncol(z)} vector (one per condition). +Defaults to \code{colnames(z)} -- one context per column.} + +\item{trait}{Condition trait label(s): a single value recycled across every +column, or a length-\code{ncol(z)} vector. Default \code{"mash"}. Pass +\code{colnames(z)} here (with a constant \code{context}) when the columns +are traits rather than contexts.} + +\item{genome}{Genome build. Default \code{"GRCh38"}.} + +\item{n}{Placeholder per-variant sample size. Default \code{1000}.} + +\item{a1, a2}{Placeholder alleles. Defaults \code{"A"} / \code{"G"}.} + +\item{role}{Tag stored in the \code{qcInfo} record. Default \code{"mash"}.} +} +\value{ +A \code{\link{QtlSumStats}} with one entry per condition (column). +} +\description{ +Assemble a per-condition \code{\link{QtlSumStats}} from a + \code{variants x conditions} Z-score matrix -- the input shape the mash + pipeline uses when only Z is available. Each column is one condition, and + conditions are distinguished by \code{context}, \code{trait}, or both: the + columns may be different cell types / tissues (contexts), different + molecular phenotypes (traits), or arbitrary context x trait pairs. + Chromosome / position are decoded from the row (variant) identifiers via + \code{\link{parseVariantId}} (with a synthetic-position fallback for ids + that do not encode coordinates); \code{A1} / \code{A2} / \code{N} are + placeholders because a Z-only input carries no alleles or sample sizes + (mash reads only Z). A pass-through \code{qcInfo} record is set so the + result clears the mash QC gate. +} diff --git a/man/resolveWeights.Rd b/man/resolveWeights.Rd new file mode 100644 index 00000000..bf145cfe --- /dev/null +++ b/man/resolveWeights.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AllGenerics.R, R/FineMappingEntry.R, +% R/TwasWeightsEntry.R +\name{resolveWeights} +\alias{resolveWeights} +\alias{resolveWeights,FineMappingEntry-method} +\alias{resolveWeights,TwasWeightsEntry-method} +\title{Resolve Per-Variant Weights From a Weight Source} +\usage{ +resolveWeights(x, ...) + +\S4method{resolveWeights}{FineMappingEntry}(x, ...) + +\S4method{resolveWeights}{TwasWeightsEntry}(x, ...) +} +\arguments{ +\item{x}{A \code{TwasWeightsEntry} or \code{FineMappingEntry}.} + +\item{...}{Reserved for future use.} +} +\value{ +A \code{list} with \code{variantIds} (character) and \code{weights} + (numeric) of equal length; both empty when no usable weights are present. +} +\description{ +Return an aligned \code{(variantIds, weights)} pair from a single + weight-source entry, so cTWAS and \code{\link{causalInferencePipeline}} + extract weights identically whether the source is a + \code{TwasWeightsEntry} (its learned weight vector) or a + \code{FineMappingEntry} (its topLoci posterior effect). +} diff --git a/man/sldscSubsetMeta.Rd b/man/sldscSubsetMeta.Rd new file mode 100644 index 00000000..6db9b363 --- /dev/null +++ b/man/sldscSubsetMeta.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sldscWrapper.R +\name{sldscSubsetMeta} +\alias{sldscSubsetMeta} +\title{Random-effects meta-analysis over a subset of sLDSC traits} +\usage{ +sldscSubsetMeta(postprocessResult, subsetTraits, targetCategories = NULL) +} +\arguments{ +\item{postprocessResult}{The list returned by +\code{\link{sldscPostprocessingPipeline}}. Must carry a \code{$per_trait} +element; \code{$params$target_categories} is used when +\code{targetCategories} is \code{NULL}.} + +\item{subsetTraits}{Character vector of trait ids to meta over; every id must +be present in \code{postprocessResult$per_trait}.} + +\item{targetCategories}{Optional character vector of target annotation names. +Defaults to \code{postprocessResult$params$target_categories}.} +} +\value{ +A list with \code{tau_star_single}, \code{tau_star_joint}, + \code{enrichment}, and \code{enrichstat}; each is a per-category named list + of \code{\link{metaSldscRandom}} results. +} +\description{ +Re-run the DerSimonian-Laird random-effects meta-analysis on a chosen subset +of the per-trait standardised tables produced by +\code{\link{sldscPostprocessingPipeline}} -- no regression is re-run, only the +already-standardised per-trait estimates are re-meta'd. Powers a "meta on a +subset of traits" workflow: pick traits, pick target annotation categories, +and get the per-category tau* / enrichment / enrichstat meta results back. +} +\seealso{ +\code{\link{sldscPostprocessingPipeline}}, \code{\link{metaSldscRandom}} +} diff --git a/man/twasWeightsPipeline.Rd b/man/twasWeightsPipeline.Rd index b7895800..0870dbd0 100644 --- a/man/twasWeightsPipeline.Rd +++ b/man/twasWeightsPipeline.Rd @@ -17,8 +17,13 @@ twasWeightsPipeline(data, ...) traitId = NULL, region = NULL, cisWindow = NULL, - minTwasMaf = NULL, - minTwasXvar = NULL, + mafCutoff = NULL, + macCutoff = NULL, + xvarCutoff = NULL, + imissCutoff = NULL, + keepIndel = NULL, + keepSamples = NULL, + keepVariants = NULL, jointRegions = FALSE, jointSpecification = NULL, fineMappingResult = NULL, @@ -108,23 +113,21 @@ processing to (QtlDataset / QtlSumStats inputs). Default \code{NULL} restrict processing to (QtlDataset / QtlSumStats inputs). Default \code{NULL}.} -\item{region}{Optional \code{GRanges} for QtlDataset trait selection. -Mutually exclusive with \code{traitId}.} +\item{region}{Optional variant window for QtlDataset trait selection: a +\code{GRanges}, a \code{"chr:start-end"} string, or a one-row data.frame +with \code{chrom}/\code{start}/\code{end}. Mutually exclusive with +\code{traitId}.} \item{cisWindow}{For QtlDataset: cis-window (bp) around each trait's genomic position when extracting variants. Required when \code{traitId} is supplied. Mutually exclusive with \code{region}.} -\item{minTwasMaf}{For QtlDataset: optional minimum minor-allele frequency -applied to the variant set used for TWAS weight learning, on top of the -dataset's construct-time \code{mafCutoff} (the effective cutoff is the -larger of the two). Lets the TWAS pass use a stricter MAF threshold than -fine mapping. \code{NULL} (default) leaves the construct-time cutoff in -place.} - -\item{minTwasXvar}{As \code{minTwasMaf} but for the per-variant genotype -variance cutoff (\code{xvarCutoff}). \code{NULL} (default) leaves the -construct-time cutoff in place.} +\item{mafCutoff, macCutoff, xvarCutoff, imissCutoff, keepIndel, keepSamples, keepVariants}{For QtlDataset: optional per-call genotype-filter overrides. Each replaces +the corresponding construct-time \code{\link{QtlDataset}} slot for this +call only (applied to a validated copy); \code{NULL} (default) leaves the +stored value in place. Variant QC is a property of the data, so these are +applied identically here and in \code{\link{fineMappingPipeline}} -- there +is deliberately no TWAS-specific variant filter.} \item{jointRegions}{For QtlDataset with a multi-range \code{region}: \code{FALSE} (default) learns weights for each range independently and diff --git a/tests/testthat/test_AllClasses.R b/tests/testthat/test_AllClasses.R index 865b0a78..ee7f5bd5 100644 --- a/tests/testthat/test_AllClasses.R +++ b/tests/testthat/test_AllClasses.R @@ -150,3 +150,42 @@ test_that("FineMappingResultBase: adjustPips on a zero-row collection returns th out <- adjustPips(empty, character(0)) expect_identical(out, empty) }) + +# =========================================================================== +# getTopLoci(type = "GRanges") + aggregate-view branches on FineMappingResultBase +# =========================================================================== +.alc_makeFmr2 <- function() { + QtlFineMappingResult( + study = c("Q1", "Q1"), context = c("c1", "c1"), trait = c("t1", "t2"), + method = c("susie", "susie"), + entry = list(.alc_makeFmEntry(), .alc_makeFmEntry()), + ldSketch = .alc_makeHandle()) +} + +test_that("getTopLoci(type='GRanges'): a pinned single entry returns a GRanges", { + fmr <- .alc_makeFmr2() + gr <- getTopLoci(fmr, type = "GRanges", + study = "Q1", context = "c1", trait = "t1", method = "susie") + expect_s4_class(gr, "GRanges") +}) + +test_that("getTopLoci(type='GRanges'): aggregating >1 entry is rejected", { + expect_error(getTopLoci(.alc_makeFmr2(), type = "GRanges"), + "requires type = 'data.frame'") +}) + +test_that("FineMappingResultBase aggregate view: empty per-entry views -> 0-row frame", { + # .alc_makeFmEntry has no credible sets, so getCs aggregates to nothing. + expect_equal(nrow(getCs(.alc_makeFmr2())), 0L) +}) + +test_that("FineMappingResultBase aggregate view: a no-match selector re-raises the selection error", { + expect_error(getCs(.alc_makeFmr2(), study = "nope"), "entries") +}) + +test_that("FineMappingResultBase aggregate view: an empty collection yields a 0-row frame", { + empty <- QtlFineMappingResult(study = character(0), context = character(0), + trait = character(0), method = character(0), entry = list(), + ldSketch = .alc_makeHandle()) + expect_equal(nrow(getCs(empty)), 0L) +}) diff --git a/tests/testthat/test_CtwasResult.R b/tests/testthat/test_CtwasResult.R new file mode 100644 index 00000000..737923dd --- /dev/null +++ b/tests/testthat/test_CtwasResult.R @@ -0,0 +1,115 @@ +context("CtwasResult / CtwasResultEntry") + +.cr_entry <- function(ids = c("g1", "g2"), pip = c(0.9, 0.1), prior = NULL) { + CtwasResultEntry( + finemap = data.frame(id = ids, susie_pip = pip, stringsAsFactors = FALSE), + susieAlpha = data.frame(id = ids, susie_alpha = pip, stringsAsFactors = FALSE), + param = prior, + regionInfo = data.frame(region_id = "r1", stringsAsFactors = FALSE)) +} + +test_that("CtwasResultEntry: constructor + accessors round-trip", { + e <- .cr_entry(prior = list(group_prior = c(brain = 0.01))) + expect_s4_class(e, "CtwasResultEntry") + expect_equal(nrow(getFinemap(e)), 2L) + expect_equal(nrow(getSusieAlpha(e)), 2L) + expect_equal(getCtwasParam(e)$group_prior, c(brain = 0.01)) + # empty entry is valid; accessors return NULL + expect_null(getFinemap(CtwasResultEntry())) + expect_null(getSusieAlpha(CtwasResultEntry())) +}) + +test_that("CtwasResult: constructor keyed by (gwasStudy, study, context, method)", { + cr <- CtwasResult( + gwasStudy = c("D1", "D1"), study = c("Q1", "Q1"), + context = c("brain", "liver"), method = c("susie", "susie"), + entry = list(.cr_entry(), .cr_entry(ids = "g3", pip = 0.5))) + expect_s4_class(cr, "CtwasResult") + expect_equal(nrow(cr), 2L) + expect_equal(getMethodNames(cr), "susie") + expect_setequal(getContexts(cr), c("brain", "liver")) + expect_equal(getStudy(cr), "Q1") +}) + +test_that("CtwasResult: getFinemap aggregates rows tagged with run identity", { + cr <- CtwasResult( + gwasStudy = c("D1", "D1"), study = c("Q1", "Q1"), + context = c("brain", "liver"), method = c("susie", "susie"), + entry = list(.cr_entry(ids = c("g1", "g2")), .cr_entry(ids = "g3", pip = 0.5))) + fm <- getFinemap(cr) + expect_true(all(c("gwasStudy", "study", "context", "method", "id", "susie_pip") + %in% names(fm))) + expect_equal(nrow(fm), 3L) # 2 + 1 gene rows + expect_equal(fm$context, c("brain", "brain", "liver")) +}) + +test_that("CtwasResult: getSusieAlpha aggregates the per-effect table with identity", { + cr <- CtwasResult( + gwasStudy = c("D1", "D1"), study = c("Q1", "Q1"), + context = c("brain", "liver"), method = c("susie", "susie"), + entry = list(.cr_entry(ids = c("g1", "g2")), .cr_entry(ids = "g3", pip = 0.5))) + sa <- getSusieAlpha(cr) + expect_equal(nrow(sa), 3L) + expect_true(all(c("gwasStudy", "study", "context", "method", "susie_alpha") + %in% names(sa))) + expect_equal(sa$context, c("brain", "brain", "liver")) +}) + +test_that("CtwasResult: uniqueness includes joint columns", { + # same (gwasStudy, study, context, method) collides without joint columns + expect_error( + CtwasResult(gwasStudy = c("D1", "D1"), study = c("Q1", "Q1"), + context = c("brain", "brain"), method = c("susie", "susie"), + entry = list(.cr_entry(), .cr_entry())), + "uniqueness violated") + # a single-context row and a multi-context (jointContexts) row over the same + # context are distinct once jointContexts joins the key + cr <- CtwasResult(gwasStudy = c("D1", "D1"), study = c("Q1", "Q1"), + context = c("brain", "brain"), method = c("susie", "susie"), + entry = list(.cr_entry(), .cr_entry()), + jointContexts = c(NA, "brain,liver")) + expect_equal(nrow(cr), 2L) +}) + +test_that("CtwasResult: rejects a non-CtwasResultEntry payload", { + expect_error( + CtwasResult(gwasStudy = "D1", study = "Q1", context = "brain", + method = "susie", entry = list("not an entry")), + "must be a CtwasResultEntry") +}) + +test_that("CtwasResult: constructor rejects mismatched core-vector lengths", { + expect_error( + CtwasResult(gwasStudy = c("D1", "D1"), study = "Q1", context = "brain", + method = "susie", entry = list(.cr_entry())), + "same length") +}) + +test_that("CtwasResult: a joint column of the wrong length is rejected", { + expect_error( + CtwasResult(gwasStudy = "D1", study = "Q1", context = "brain", + method = "susie", entry = list(.cr_entry()), + jointContexts = c("a", "b")), + "same length as") +}) + +test_that("CtwasResult: getFinemap/getSusieAlpha are NULL when empty or all-NULL", { + empty <- CtwasResult(gwasStudy = character(0), study = character(0), + context = character(0), method = character(0), + entry = list()) + expect_null(getFinemap(empty)) + expect_null(getSusieAlpha(empty)) + nullE <- CtwasResult(gwasStudy = "D1", study = "Q1", context = "brain", + method = "susie", + entry = list(CtwasResultEntry(finemap = NULL, + susieAlpha = NULL))) + expect_null(getFinemap(nullE)) + expect_null(getSusieAlpha(nullE)) +}) + +test_that("CtwasResult: show() prints a one-line-per-run summary", { + cr <- CtwasResult(gwasStudy = "D1", study = "Q1", context = "brain", + method = "susie", entry = list(.cr_entry())) + expect_output(show(cr), "CtwasResult: 1 run") + expect_output(show(cr), "susie") +}) diff --git a/tests/testthat/test_FineMappingEntry.R b/tests/testthat/test_FineMappingEntry.R index 1daa2e9a..fd02318c 100644 --- a/tests/testthat/test_FineMappingEntry.R +++ b/tests/testthat/test_FineMappingEntry.R @@ -164,6 +164,19 @@ test_that("FineMappingEntry: getPip returns named pip vector keyed by variant_id paste0("chr1:", 100 * 1:3, ":A:G")) }) +test_that("FineMappingEntry: resolveWeights returns topLoci posterior effect aligned to variant_id", { + entry <- .sc_makeFineMappingEntry(3) # topLoci posterior_mean = 0.05 for all + wr <- resolveWeights(entry) + expect_equal(wr$variantIds, paste0("chr1:", 100 * 1:3, ":A:G")) + expect_equal(wr$weights, rep(0.05, 3)) + expect_equal(length(wr$variantIds), length(wr$weights)) + # empty topLoci -> empty pair + empty <- FineMappingEntry(variantIds = character(0), susieFit = list(), + topLoci = data.frame(variant_id = character(0), pip = numeric(0), + stringsAsFactors = FALSE)) + expect_length(resolveWeights(empty)$variantIds, 0L) +}) + test_that("FineMappingEntry: getPip returns numeric(0) when topLoci is empty", { entry <- FineMappingEntry( diff --git a/tests/testthat/test_GwasFineMappingResult.R b/tests/testthat/test_GwasFineMappingResult.R index 813401c1..dcdeaaaa 100644 --- a/tests/testthat/test_GwasFineMappingResult.R +++ b/tests/testthat/test_GwasFineMappingResult.R @@ -13,6 +13,20 @@ test_that("GwasFineMappingResult: builds a collection keyed by 2-tuple", { expect_equal(nrow(res), 2L) }) +test_that("GwasFineMappingResult: region is auto-derived from region_id", { + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- GwasFineMappingResult( + study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("chr1_1_500", "chr2_600_900"), entry = list(e1, e2)) + reg <- getRegion(res) + expect_equal(as.character(GenomicRanges::seqnames(reg)), c("chr1", "chr2")) + expect_equal(GenomicRanges::start(reg), c(1L, 600L)) + expect_equal(GenomicRanges::end(reg), c(500L, 900L)) + # synthetic region_id (no coordinates) -> a chrUn 0-width sentinel + res2 <- GwasFineMappingResult(study = "g1", method = "susie", entry = list(e1)) + expect_equal(as.character(GenomicRanges::seqnames(getRegion(res2))), "chrUn") +}) + test_that("GwasFineMappingResult: validity does not recurse on key subset (#546)", { # The validity method builds a key-column data.frame to check tuple @@ -164,6 +178,50 @@ test_that("GwasFineMappingResult: getCs/getTopLoci/getSusieFit/getVariantIds dis expect_equal(length(getVariantIds(res)), 3L) }) +test_that("GwasFineMappingResult: getTopLoci aggregates per-block rows genome-wide", { + # A genome-wide collection: same (study, method) across two region blocks. + # With no selectors getTopLoci now stacks both blocks, tagging each variant + # with its region_id; context/trait are NA-filled (GWAS keys on region). + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- GwasFineMappingResult( + study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("chr1:1-100", "chr1:200-300"), entry = list(e1, e2)) + agg <- getTopLoci(res, signalCutoff = 0) + expect_equal(nrow(agg), 5L) + expect_equal(agg$region_id, c("chr1:1-100", "chr1:1-100", "chr1:1-100", + "chr1:200-300", "chr1:200-300")) + expect_true(all(is.na(agg$context))) + expect_true(all(is.na(agg$trait))) + expect_true("variant_id" %in% names(agg)) +}) + +test_that("GwasFineMappingResult: getTopLoci region= selects a single block", { + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- GwasFineMappingResult( + study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("chr1:1-100", "chr1:200-300"), entry = list(e1, e2)) + # region= pins one row, so this hits the single-entry fast path (bare table). + tl <- getTopLoci(res, study = "g1", method = "susie", + region = "chr1:200-300", signalCutoff = 0) + expect_equal(nrow(tl), 2L) + expect_false("region_id" %in% names(tl)) +}) + +test_that("GwasFineMappingResult: getCs aggregates per-block credible sets tagged by region_id", { + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- GwasFineMappingResult( + study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("chr1:1-100", "chr1:200-300"), entry = list(e1, e2)) + cs <- getCs(res) + expect_equal(nrow(cs), 4L) # 2 CS members per block + expect_equal(cs$region_id, + c("chr1:1-100", "chr1:1-100", "chr1:200-300", "chr1:200-300")) + expect_true(all(is.na(cs$context))) + # region= pins one block -> bare table + bare <- getCs(res, study = "g1", method = "susie", region = "chr1:1-100") + expect_false("region_id" %in% names(bare)) +}) + test_that("GwasFineMappingResult: .tupleSelectRowGwasFmr requires both selectors for multi-row", { e <- .ca_makeFmEntry(3) diff --git a/tests/testthat/test_QtlFineMappingResult.R b/tests/testthat/test_QtlFineMappingResult.R index bb8ce76b..009d8192 100644 --- a/tests/testthat/test_QtlFineMappingResult.R +++ b/tests/testthat/test_QtlFineMappingResult.R @@ -290,6 +290,94 @@ test_that("QtlFineMappingResult: getTopLoci returns the entry's topLoci (project expect_equal(tl$variant_id, .ca_makeTopLoci(3)$variant_id) }) +test_that("QtlFineMappingResult: getTopLoci aggregates entries when selectors do not pin one", { + # Two entries, no selectors: the old behaviour errored ("2 entries; pass + # study/context/trait/method"). Now it stacks the per-variant tables, + # prefixed with the row identity so variants stay attributable. + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- QtlFineMappingResult( + study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + entry = list(e1, e2)) + agg <- getTopLoci(res, signalCutoff = 0) + expect_equal(nrow(agg), 5L) + expect_true(all(c("study", "context", "trait", "region_id", "method") %in% + names(agg))) + expect_equal(agg$context, c("c1", "c1", "c1", "c2", "c2")) + expect_true("variant_id" %in% names(agg)) + # QTL results key on context/trait, so region_id is NA-filled. + expect_true(all(is.na(agg$region_id))) + # the stamped per-variant `method` must not duplicate the identity column + expect_equal(sum(names(agg) == "method"), 1L) +}) + +test_that("QtlFineMappingResult: getTopLoci with a full tuple keeps the bare (id-free) table", { + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- QtlFineMappingResult( + study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + entry = list(e1, e2)) + tl <- getTopLoci(res, study = "s1", context = "c1", trait = "t1", + method = "susie", signalCutoff = 0) + expect_equal(nrow(tl), 3L) + expect_false("context" %in% names(tl)) # single-entry fast path, unchanged +}) + +test_that("QtlFineMappingResult: getTopLoci aggregates only the matching subset", { + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- QtlFineMappingResult( + study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + entry = list(e1, e2)) + sub <- getTopLoci(res, context = "c2", signalCutoff = 0) + expect_equal(nrow(sub), 2L) + expect_equal(unique(sub$context), "c2") +}) + +test_that("QtlFineMappingResult: getTopLoci on a single-row collection still carries identity columns", { + # A no-selector call takes the aggregate path even for one row, so the PIP + # table has a uniform shape (identity columns present) whatever the row count + # -- callers don't special-case single-entry results. + res <- QtlFineMappingResult(study = "s1", context = "c1", trait = "t1", + method = "susie", + entry = list(.sc_makeFineMappingEntry(4))) + tl <- getTopLoci(res, signalCutoff = 0) + expect_equal(nrow(tl), 4L) + expect_true(all(c("study", "context", "trait", "method", "variant_id") %in% + names(tl))) + expect_equal(unique(tl$context), "c1") +}) + +test_that("QtlFineMappingResult: getCs aggregates every entry's credible sets with identity columns", { + # e1 (n=3) has cs_95 = susie_1/susie_1/susie_0 -> 2 CS members; e2 (n=2) -> 2. + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- QtlFineMappingResult( + study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + entry = list(e1, e2)) + cs <- getCs(res) + expect_equal(nrow(cs), 4L) + expect_true(all(c("study", "context", "trait", "region_id", "method", + "variant_id") %in% names(cs))) + expect_equal(cs$context, c("c1", "c1", "c2", "c2")) + # full tuple still returns the bare per-entry credible-set table + bare <- getCs(res, study = "s1", context = "c1", trait = "t1", method = "susie") + expect_false("context" %in% names(bare)) + expect_equal(nrow(bare), 2L) +}) + +test_that("QtlFineMappingResult: getMarginalEffects aggregates every entry with identity columns", { + e1 <- .sc_makeFineMappingEntry(3); e2 <- .sc_makeFineMappingEntry(2) + res <- QtlFineMappingResult( + study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + entry = list(e1, e2)) + me <- getMarginalEffects(res) + expect_equal(nrow(me), 5L) # all topLoci rows (no CS filter) + expect_true(all(c("study", "context", "trait", "method") %in% names(me))) + expect_equal(unique(me$context), c("c1", "c2")) +}) + test_that("QtlFineMappingResult: getSusieFit reads the entry's trimmedFit", { e <- .ca_makeFmEntry(3) diff --git a/tests/testthat/test_TwasWeightsEntry.R b/tests/testthat/test_TwasWeightsEntry.R index 7c4b2546..1627e387 100644 --- a/tests/testthat/test_TwasWeightsEntry.R +++ b/tests/testthat/test_TwasWeightsEntry.R @@ -19,6 +19,18 @@ test_that("TwasWeightsEntry: constructor and accessors round-trip", { expect_equal(getDataType(e), "expression") }) +test_that("TwasWeightsEntry: resolveWeights returns the aligned (variantIds, weights) pair", { + e <- TwasWeightsEntry(variantIds = c("v1", "v2", "v3"), + weights = c(0.1, -0.2, 0.05)) + wr <- resolveWeights(e) + expect_equal(wr$variantIds, c("v1", "v2", "v3")) + expect_equal(wr$weights, c(0.1, -0.2, 0.05)) + # length mismatch (defensive) -> empty + e2 <- TwasWeightsEntry(variantIds = c("v1", "v2"), weights = c(0.1, -0.2)) + e2@weights <- matrix(0, 2, 2) # flattens to length 4 != 2 ids + expect_length(resolveWeights(e2)$variantIds, 0L) +}) + test_that("TwasWeightsEntry: standardized is coerced via isTRUE() semantics", { # isTRUE() only returns TRUE for a length-1 logical TRUE. Non-TRUE diff --git a/tests/testthat/test_crossValidation.R b/tests/testthat/test_crossValidation.R new file mode 100644 index 00000000..dcf37d11 --- /dev/null +++ b/tests/testthat/test_crossValidation.R @@ -0,0 +1,176 @@ +context("crossValidation (shared CV engine)") + +# The engine is exercised directly through .crossValidateWeights() with a mock +# per-fold fit, so these tests cover the harness mechanics (partitioning, fold +# aggregation, prediction, metrics, key format, subsampling, threading) without +# any real weight-learning. twasWeightsCv() and .fmWeightsCv() are thin callers +# that supply a domain-specific fitFold; their integration is tested elsewhere. + +cv <- function(...) pecotmr:::.crossValidateWeights(...) + +# One "mock" method whose weights are all 1s over the training columns, so a +# held-out prediction is the row sum of that sample's (training-column) dosages. +mock_fit_fold <- function(Xtr, Ytr, j) { + list(weights = list(mock = matrix(1, ncol(Xtr), ncol(Ytr), + dimnames = list(colnames(Xtr), NULL))), + fits = list()) +} + +mk_xy <- function(n = 30, p = 6, k = 1, seed = 1) { + set.seed(seed) + X <- matrix(rnorm(n * p), n, p, + dimnames = list(paste0("s", seq_len(n)), paste0("v", seq_len(p)))) + Y <- matrix(rnorm(n * k), n, k, + dimnames = list(rownames(X), paste0("c", seq_len(k)))) + list(X = X, Y = Y) +} + +test_that("input is validated", { + d <- mk_xy() + expect_error(cv(d$X, d$Y, fold = 0, fitFold = mock_fit_fold), "positive integer") + expect_error(cv(d$X, d$Y, fold = "a", fitFold = mock_fit_fold), "positive integer") + expect_error(cv(as.data.frame(d$X), d$Y, fold = 2, fitFold = mock_fit_fold), + "X must be a matrix") + expect_error(cv(d$X, d$Y[1:5, , drop = FALSE], fold = 2, fitFold = mock_fit_fold), + "same") + expect_error(cv(d$X, d$Y, fitFold = mock_fit_fold), + "Either 'fold' or 'samplePartitions'") +}) + +test_that("Y as a vector is converted to a matrix with a message", { + d <- mk_xy(k = 1) + expect_message( + cv(d$X, as.numeric(d$Y), fold = 2, fitFold = mock_fit_fold, verbose = 1), + "Y converted to matrix") +}) + +test_that("output keys use the canonical _predicted / _performance", { + d <- mk_xy() + r <- suppressMessages(cv(d$X, d$Y, fold = 3, fitFold = mock_fit_fold)) + expect_equal(names(r$prediction), "mock_predicted") + expect_equal(names(r$performance), "mock_performance") +}) + +test_that("performance carries the six metric colnames and outcome rownames", { + d <- mk_xy(k = 2) + r <- suppressMessages(cv(d$X, d$Y, fold = 3, fitFold = mock_fit_fold)) + perf <- r$performance[["mock_performance"]] + expect_equal(colnames(perf), c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE")) + expect_equal(rownames(perf), colnames(d$Y)) +}) + +test_that("every sample is predicted exactly once across folds", { + d <- mk_xy(n = 30) + r <- suppressMessages(cv(d$X, d$Y, fold = 5, fitFold = mock_fit_fold)) + pred <- r$prediction[["mock_predicted"]] + expect_equal(nrow(pred), nrow(d$X)) + expect_false(any(is.na(pred))) +}) + +test_that("a provided samplePartition is reused; a fold mismatch warns", { + d <- mk_xy(n = 12) + sp <- data.frame(Sample = rownames(d$X), Fold = rep(1:3, each = 4), + stringsAsFactors = FALSE) + r <- suppressMessages(cv(d$X, d$Y, samplePartitions = sp, fitFold = mock_fit_fold)) + expect_equal(r$samplePartition, sp) + expect_message( + cv(d$X, d$Y, fold = 2, samplePartitions = sp, fitFold = mock_fit_fold, + verbose = 1), + "does not match") +}) + +test_that("a samplePartition with unknown samples errors", { + d <- mk_xy() + sp <- data.frame(Sample = paste0("zzz", 1:5), Fold = rep(1:2, length.out = 5)) + expect_error(cv(d$X, d$Y, samplePartitions = sp, fitFold = mock_fit_fold), + "do not match") +}) + +test_that("maxNumVariants subsamples variants with a message", { + d <- mk_xy(p = 20) + expect_message( + cv(d$X, d$Y, fold = 2, fitFold = mock_fit_fold, maxNumVariants = 8, + verbose = 1), + "Randomly selecting 8 out of 20") +}) + +test_that("maxNumVariants with variantsToKeep retains the specified variants", { + d <- mk_xy(p = 20) + expect_message( + cv(d$X, d$Y, fold = 2, fitFold = mock_fit_fold, maxNumVariants = 8, + variantsToKeep = c("v1", "v2", "v3"), verbose = 1), + "Including 3 specified variants") +}) + +test_that("a degenerate fold (empty train/test) is skipped, not errored", { + d <- mk_xy(n = 10) + sp <- data.frame(Sample = rownames(d$X), Fold = rep(1L, nrow(d$X))) + r <- suppressMessages(cv(d$X, d$Y, samplePartitions = sp, fitFold = mock_fit_fold)) + expect_true(all(is.na(r$prediction[["mock_predicted"]]))) +}) + +test_that("zero-variance predictions yield NA metrics with a message", { + d <- mk_xy() + zero_fit <- function(Xtr, Ytr, j) { + list(weights = list(mock = matrix(0, ncol(Xtr), ncol(Ytr), + dimnames = list(colnames(Xtr), NULL))), + fits = list()) + } + expect_message(r <- cv(d$X, d$Y, fold = 3, fitFold = zero_fit, verbose = 1), + "zero variance") + expect_true(all(is.na(r$performance[["mock_performance"]]))) +}) + +test_that("the parallel fold path matches the serial one", { + d <- mk_xy(n = 40, seed = 7) + set.seed(1) + r1 <- suppressMessages(cv(d$X, d$Y, fold = 4, fitFold = mock_fit_fold, + numThreads = 1)) + set.seed(1) + r2 <- suppressMessages(cv(d$X, d$Y, fold = 4, fitFold = mock_fit_fold, + numThreads = 2)) + expect_equal(r1$prediction, r2$prediction) +}) + +test_that("retainFits collects per-fold fits only when requested", { + d <- mk_xy() + fit_with_model <- function(Xtr, Ytr, j) { + list(weights = list(mock = matrix(1, ncol(Xtr), ncol(Ytr), + dimnames = list(colnames(Xtr), NULL))), + fits = list(mock = list(fold = j))) + } + r_off <- suppressMessages(cv(d$X, d$Y, fold = 3, fitFold = fit_with_model, + retainFits = FALSE)) + expect_true(all(vapply(r_off$foldFits, length, integer(1)) == 0L)) + r_on <- suppressMessages(cv(d$X, d$Y, fold = 3, fitFold = fit_with_model, + retainFits = TRUE)) + expect_equal(r_on$foldFits[["fold_1"]][["mock"]]$fold, 1L) +}) + +test_that("X without rownames inherits the sample names (from Y)", { + d <- mk_xy() + X2 <- d$X; rownames(X2) <- NULL + r <- suppressMessages(cv(X2, d$Y, fold = 3, fitFold = mock_fit_fold)) + expect_equal(rownames(r$prediction$mock_predicted), rownames(d$Y)) +}) + +test_that("maxNumVariants subsamples from variantsToKeep when it already exceeds the cap", { + d <- mk_xy(p = 6) + expect_message( + cv(d$X, d$Y, fold = 2, fitFold = mock_fit_fold, maxNumVariants = 2, + variantsToKeep = c("v1", "v2", "v3"), verbose = 1), + "Randomly selecting 2 out of 3") +}) + +test_that("a NULL per-method weight matrix yields an all-NA prediction, not an error", { + d <- mk_xy() + fit_with_null <- function(Xtr, Ytr, j) + list(weights = list( + mock = matrix(1, ncol(Xtr), ncol(Ytr), + dimnames = list(colnames(Xtr), NULL)), + empty = NULL), # NULL W -> skipped per fold + fits = list()) + r <- suppressMessages(cv(d$X, d$Y, fold = 3, fitFold = fit_with_null)) + expect_true(all(is.na(r$prediction$empty_predicted))) + expect_false(all(is.na(r$prediction$mock_predicted))) +}) diff --git a/tests/testthat/test_ctwasPipeline.R b/tests/testthat/test_ctwasPipeline.R index 46fb9ce0..cf724b13 100644 --- a/tests/testthat/test_ctwasPipeline.R +++ b/tests/testthat/test_ctwasPipeline.R @@ -89,6 +89,23 @@ context("ctwasPipeline") ldSketch = .ctp_makeHandle()) } +# The same per-gene weights expressed as a QtlFineMappingResult weight source: +# the topLoci posterior_mean carries the weight vector (what resolveWeights +# reads). Mirrors .ctp_makeTwasWeights so the two sources are comparable. +.ctp_makeFmrWeightSource <- function() { + vids <- vapply(1:5, .ctp_snpId, character(1)) + w <- c(0.1, 0.05, -0.2, 0.3, 0.0) + tl <- data.frame( + variant_id = vids, chrom = rep("1", 5), pos = 100L * (1:5), + A1 = rep("A", 5), A2 = rep("G", 5), N = rep(100, 5), af = rep(0.3, 5), + pip = rep(0.5, 5), posterior_mean = w, posterior_sd = rep(0.1, 5), + cs_95 = rep("susie_1", 5), stringsAsFactors = FALSE) + fe <- FineMappingEntry(variantIds = vids, susieFit = list(), topLoci = tl) + QtlFineMappingResult(study = "Q1", context = "c1", trait = "t1", + method = "susie", entry = list(fe), + ldSketch = .ctp_makeHandle()) +} + # =========================================================================== # Input validation # ---------------------------------------------------------------------------- @@ -164,6 +181,30 @@ test_that("assembleCtwasInputs: allows twasWeights keys to be a subset of gwasSu expect_true(all(grepl("^block1\\|", names(inputs$weights)))) }) +test_that("assembleCtwasInputs: accepts a QtlFineMappingResult weight source (topLoci posterior effect)", { + ss <- .ctp_makeGwasSumstats() + # The FMR topLoci posterior effect is on the STANDARDIZED scale, so it matches + # a standardized TwasWeights carrying the same effect vector (both skip cTWAS's + # variance scaling). Compare against that, not the default (unstandardized) one. + tw_std <- TwasWeights( + study = "Q1", context = "c1", trait = "t1", method = "susie", + entry = list(TwasWeightsEntry( + variantIds = vapply(1:5, .ctp_snpId, character(1)), + weights = c(0.1, 0.05, -0.2, 0.3, 0.0), standardized = TRUE)), + ldSketch = .ctp_makeHandle()) + fmr <- .ctp_makeFmrWeightSource() + local_mocked_bindings(extractBlockGenotypes = .ctp_mockExtractor(), + .package = "pecotmr") + itw <- assembleCtwasInputs(gwasSumStats = list(block1 = ss, block2 = ss), + twasWeights = list(block1 = tw_std, block2 = tw_std)) + ifmr <- assembleCtwasInputs(gwasSumStats = list(block1 = ss, block2 = ss), + twasWeights = list(block1 = fmr, block2 = fmr)) + # Same gene keys, and the resolved weights match the standardized TwasWeights. + expect_equal(names(ifmr$weights), names(itw$weights)) + expect_true(length(ifmr$weights) > 0L) + expect_equal(ifmr$weights[[1L]], itw$weights[[1L]]) +}) + test_that("assembleCtwasInputs: filters TwasWeights against UNION of all blocks' GWAS variants", { # Build two blocks with NON-OVERLAPPING GWAS variants. Block 1 covers # v1..v3, block 2 covers v4..v6. The gene's weight spans v2..v5 — i.e. @@ -204,13 +245,16 @@ test_that("assembleCtwasInputs: filters TwasWeights against UNION of all blocks' expect_setequal(rownames(wgt), vapply(2:5, .ctp_snpId, character(1))) }) -test_that("ctwasPipeline: rejects bare (non-list) TwasWeights", { +test_that("ctwasPipeline: a bare TwasWeights is a FLAT source, placed by region", { skip_if_not_installed("ctwas") + # A bare TwasWeights is now accepted as a flat weight source and placed into + # LD blocks by region; the fixture carries no region, so placement errors + # (rather than the old 'must be a NAMED LIST' shape rejection). expect_error( ctwasPipeline(gwasSumStats = list(block1 = .ctp_makeGwasSumstats(), block2 = .ctp_makeGwasSumstats()), twasWeights = .ctp_makeTwasWeights()), - "NAMED LIST of TwasWeights" + "no `region` provenance" ) }) @@ -400,6 +444,216 @@ test_that(".ctwasResolveMethod: multi-method + no ensemble + no caller method er "Supply a `method` argument") }) +# --------------------------------------------------------------------------- +# .ctwasResolveMethods (plural) — the pipeline's fan-out resolver +# --------------------------------------------------------------------------- +.ctp_multiMethodTw <- function(methods = c("mrash", "susie")) { + TwasWeights( + study = rep("Q1", length(methods)), context = rep("c1", length(methods)), + trait = rep("t1", length(methods)), method = methods, + entry = lapply(methods, function(m) + TwasWeightsEntry(variantIds = paste0("v", 1:3), weights = c(0.1, 0.2, 0.3))), + ldSketch = .ctp_makeHandle()) +} + +test_that(".ctwasResolveMethods: unspecified + multiple + no ensemble -> ALL methods", { + tw <- .ctp_multiMethodTw(c("mrash", "susie")) + expect_setequal(pecotmr:::.ctwasResolveMethods(list(r1 = tw)), c("mrash", "susie")) +}) + +test_that(".ctwasResolveMethods: ensemble wins, single auto-picks, explicit restricts", { + expect_equal(pecotmr:::.ctwasResolveMethods( + list(r1 = .ctp_multiMethodTw(c("mrash", "ensemble")))), "ensemble") + expect_equal(pecotmr:::.ctwasResolveMethods( + list(r1 = .ctp_multiMethodTw("mrash"))), "mrash") + expect_equal(pecotmr:::.ctwasResolveMethods( + list(r1 = .ctp_multiMethodTw(c("mrash", "susie")), r2 = NULL), method = "susie"), + "susie") + expect_error(pecotmr:::.ctwasResolveMethods( + list(r1 = .ctp_multiMethodTw("mrash")), method = "lasso"), "not present") +}) + +# --------------------------------------------------------------------------- +# .ctwasGwasStudy — single-disease invariant +# --------------------------------------------------------------------------- +test_that(".ctwasGwasStudy: unique study, errors on mixed studies", { + # .ctwasGwasStudy only reads `$study`, so plain data.frame stand-ins suffice + # (and dodge the S4 `$<-` coercion barrier on GwasSumStats). + mk <- function(s) data.frame(study = s, variant_id = "chr1:1:G:A") + expect_equal(pecotmr:::.ctwasGwasStudy(list(a = mk("D1"), b = mk("D1"))), "D1") + expect_error(pecotmr:::.ctwasGwasStudy(list(a = mk("D1"), b = mk("D2"))), + "multiple GWAS studies") +}) + +# --------------------------------------------------------------------------- +# Phase 5: internal LD-block placement (by start(region) == cTWAS p0 rule) +# --------------------------------------------------------------------------- +test_that(".ctwasPlaceByAnchor: homes each gene by start(region), half-open, chr-agnostic", { + region <- GenomicRanges::GRanges( + c("chr1", "chr1", "22", "chr2"), + IRanges::IRanges(start = c(500, 1500, 800, 100), + end = c(600, 1600, 900, 200))) + grid <- setNames(vector("list", 3), + c("chr1_1_1000", "chr1_1000_2000", "chr22_1_1000")) + home <- pecotmr:::.ctwasPlaceByAnchor(region, grid) + # chr1:500 -> [1,1000); chr1:1500 -> [1000,2000); "22":800 -> chr22 block + # (chr prefix normalized); chr2 -> no block -> NA. + expect_equal(home, c("chr1_1_1000", "chr1_1000_2000", "chr22_1_1000", NA)) +}) + +test_that(".ctwasPlaceByAnchor: block interval is half-open [start, stop)", { + region <- GenomicRanges::GRanges("chr1", + IRanges::IRanges(start = c(1000, 999), end = c(1000, 999))) + grid <- setNames(vector("list", 2), c("chr1_1_1000", "chr1_1000_2000")) + # start==1000 falls in the SECOND block (>= 1000), 999 in the first. + expect_equal(pecotmr:::.ctwasPlaceByAnchor(region, grid), + c("chr1_1000_2000", "chr1_1_1000")) +}) + +test_that(".ctwasIsPreBucketed / .ctwasCombineWeightSources: dispatch flat vs pre-bucketed", { + tw <- .ctp_makeTwasWeights() + grid <- list(block1 = 1, block2 = 2) + expect_true(pecotmr:::.ctwasIsPreBucketed(list(block1 = tw), grid)) # named list + expect_false(pecotmr:::.ctwasIsPreBucketed(tw, grid)) # flat S4 + expect_false(pecotmr:::.ctwasIsPreBucketed(list(tw), grid)) # unnamed + expect_s4_class(pecotmr:::.ctwasCombineWeightSources(tw), "TwasWeights") + expect_s4_class(pecotmr:::.ctwasCombineWeightSources(list(tw)), "TwasWeights") +}) + +test_that(".ctwasBucketWeights: places a flat 2-gene source into its home blocks", { + mkE <- function() TwasWeightsEntry( + variantIds = vapply(1:5, .ctp_snpId, character(1)), + weights = c(0.1, 0.05, -0.2, 0.3, 0.0)) + # gA anchor @chr1:100 -> block chr1_1_350; gB anchor @chr1:400 -> chr1_350_700. + tw <- TwasWeights( + study = c("Q1", "Q1"), context = c("c1", "c1"), + trait = c("gA", "gB"), method = c("susie", "susie"), + entry = list(mkE(), mkE()), + region = GenomicRanges::GRanges(c("chr1", "chr1"), + IRanges::IRanges(c(100, 400), c(150, 450))), + ldSketch = .ctp_makeHandle()) + grid <- list(chr1_1_350 = .ctp_makeGwasSumstats(), + chr1_350_700 = .ctp_makeGwasSumstats()) + bucketed <- pecotmr:::.ctwasBucketWeights(tw, grid) + expect_named(bucketed, c("chr1_1_350", "chr1_350_700"), ignore.order = TRUE) + expect_equal(as.character(bucketed[["chr1_1_350"]]$trait), "gA") + expect_equal(as.character(bucketed[["chr1_350_700"]]$trait), "gB") + expect_false(is.null(getLdSketch(bucketed[["chr1_1_350"]]))) +}) + +test_that(".ctwasBucketWeights: errors when the weight source has no region", { + expect_error( + pecotmr:::.ctwasBucketWeights(.ctp_makeTwasWeights(), + list(chr1_1_350 = .ctp_makeGwasSumstats(), + chr1_350_700 = .ctp_makeGwasSumstats())), + "no `region` provenance") +}) + +# --------------------------------------------------------------------------- +# .ctwasParseGeneIds — id -> identity components +# --------------------------------------------------------------------------- +test_that(".ctwasParseGeneIds: splits region|study|context|trait|method", { + p <- pecotmr:::.ctwasParseGeneIds(c("blk|Q1|brain|gA|susie", "blk2|Q1|liver|gB|lasso")) + expect_equal(p$rid, c("blk", "blk2")) + expect_equal(p$study, c("Q1", "Q1")) + expect_equal(p$context, c("brain", "liver")) + expect_equal(p$trait, c("gA", "gB")) + expect_equal(p$method, c("susie", "lasso")) + expect_error(pecotmr:::.ctwasParseGeneIds("too|few|fields"), "malformed") +}) + +# --------------------------------------------------------------------------- +# .ctwasRunToRows — decompose a cTWAS run into per-context CtwasResult rows +# --------------------------------------------------------------------------- +.ctp_runResult <- function(ids, pips, prior, snpIds = character(0), + region = data.frame(region_id = "blk")) { + geneFm <- data.frame(id = ids, type = "gene", susie_pip = pips, + stringsAsFactors = FALSE) + snpFm <- if (length(snpIds)) + data.frame(id = snpIds, type = "SNP", + susie_pip = seq(0.01, by = 0.01, length.out = length(snpIds)), + stringsAsFactors = FALSE) else NULL + fm <- rbind(geneFm, snpFm) + list( + weights = setNames(as.list(seq_along(ids)), ids), + finemap_res = fm, + susie_alpha_res = cbind(fm, susie_alpha = 0.5), + param = list(group_prior = prior), + region_info = region) +} + +test_that(".ctwasRunToRows: single-context run -> one row, no jointContexts", { + run <- .ctp_runResult(c("blk|Q1|c1|gA|susie", "blk|Q1|c1|gB|susie"), + c(0.9, 0.2), c(c1 = 0.01, SNP = 1e-4)) + rows <- pecotmr:::.ctwasRunToRows(run, gwasStudy = "D1", method = "susie") + expect_length(rows, 1L) + expect_equal(rows[[1L]]$context, "c1") + expect_equal(rows[[1L]]$study, "Q1") + expect_equal(rows[[1L]]$gwasStudy, "D1") + expect_true(is.na(rows[[1L]]$jointContexts)) + expect_equal(nrow(getFinemap(rows[[1L]]$entry)), 2L) + expect_equal(getCtwasParam(rows[[1L]]$entry)$group_prior[["c1"]], 0.01) +}) + +test_that(".ctwasRunToRows: multi-context run -> per-context rows sharing jointContexts + param", { + ids <- c("blk|Q1|brain|gA|susie", "blk|Q1|brain|gB|susie", + "blk|Q1|liver|gA|susie", "blk|Q1|liver|gB|susie") + run <- .ctp_runResult(ids, c(0.9, 0.1, 0.8, 0.2), + c(brain = 0.01, liver = 0.02, SNP = 1e-4)) + rows <- pecotmr:::.ctwasRunToRows(run, gwasStudy = "D1", method = "susie") + expect_length(rows, 2L) + expect_setequal(vapply(rows, function(r) r$context, ""), c("brain", "liver")) + expect_true(all(vapply(rows, function(r) r$jointContexts, "") == "brain,liver")) + # Each per-context row keeps only its own genes but shares the joint param. + expect_equal(nrow(getFinemap(rows[[1L]]$entry)), 2L) + expect_named(getCtwasParam(rows[[1L]]$entry)$group_prior, + c("brain", "liver", "SNP")) +}) + +test_that(".ctwasRunToRows: rejects multi-context runs with unshared gene sets", { + run <- .ctp_runResult(c("blk|Q1|brain|gA|susie", "blk|Q1|liver|gC|susie"), + c(0.9, 0.8), c(brain = 0.01, liver = 0.02, SNP = 1e-4)) + expect_error(pecotmr:::.ctwasRunToRows(run, "D1", "susie"), + "SAME gene set in every context") +}) + +test_that(".ctwasRunToRows: empty finemap_res still yields the modeled row", { + run <- list(weights = setNames(list(1), "blk|Q1|c1|gA|susie"), + finemap_res = NULL, param = list(), region_info = NULL) + rows <- pecotmr:::.ctwasRunToRows(run, "D1", "susie") + expect_length(rows, 1L) + expect_null(getFinemap(rows[[1L]]$entry)) + expect_null(getSusieAlpha(rows[[1L]]$entry)) +}) + +test_that(".ctwasRunToRows: susieAlpha subset mirrors finemap per context", { + run <- .ctp_runResult(c("blk|Q1|c1|gA|susie", "blk|Q1|c1|gB|susie"), + c(0.9, 0.2), c(c1 = 0.01, SNP = 1e-4)) + rows <- pecotmr:::.ctwasRunToRows(run, "D1", "susie") + sa <- getSusieAlpha(rows[[1L]]$entry) + expect_equal(nrow(sa), 2L) + expect_true("susie_alpha" %in% names(sa)) +}) + +test_that(".ctwasRunToRows: keepSnps adds a dedicated SNP row; default drops SNPs", { + ids <- c("blk|Q1|c1|gA|susie", "blk|Q1|c1|gB|susie") + run <- .ctp_runResult(ids, c(0.9, 0.2), c(c1 = 0.01, SNP = 1e-4), + snpIds = c("chr1:1:G:A", "chr1:2:C:T")) + # default: SNP background dropped, only the gene (context) row. + rowsOff <- pecotmr:::.ctwasRunToRows(run, "D1", "susie") + expect_length(rowsOff, 1L) + expect_true(all(getFinemap(rowsOff[[1L]]$entry)$type == "gene")) + # keepSnps: one extra study = context = "SNP" row carrying the SNP background. + rowsOn <- pecotmr:::.ctwasRunToRows(run, "D1", "susie", keepSnps = TRUE) + expect_length(rowsOn, 2L) + snpRow <- rowsOn[[2L]] + expect_equal(snpRow$study, "SNP") + expect_equal(snpRow$context, "SNP") + expect_equal(nrow(getFinemap(snpRow$entry)), 2L) + expect_true(all(getFinemap(snpRow$entry)$type == "SNP")) + expect_equal(nrow(getSusieAlpha(snpRow$entry)), 2L) +}) + test_that(".ctwasFilterMethod: subsets rows to the requested method", { tw <- TwasWeights( study = c("Q1", "Q1"), context = c("c1", "c1"), @@ -782,8 +1036,17 @@ test_that("ctwasPipeline: dispatches assemble → est → screen → finemap and }, finemap_regions = function(...) { capturedFinemap <<- list(...) - list(finemap_res = data.frame(id = "t1"), - susie_alpha_res = data.frame(susie_pip = 0.9)) + # Return per-gene finemap rows keyed by our cTWAS gene ids + # (region|study|context|trait|method) so the decomposition into + # CtwasResult rows finds a matching finemap payload. + list(finemap_res = data.frame( + id = c("block1|Q1|c1|t1|susie", "block2|Q1|c1|t1|susie"), + type = c("gene", "gene"), context = c("c1", "c1"), + susie_pip = c(0.9, 0.8), stringsAsFactors = FALSE), + susie_alpha_res = data.frame( + id = c("block1|Q1|c1|t1|susie", "block2|Q1|c1|t1|susie"), + type = c("gene", "gene"), + susie_alpha = c(0.9, 0.8), stringsAsFactors = FALSE)) }, .package = "ctwas") local_mocked_bindings(extractBlockGenotypes = .ctp_mockExtractor(), @@ -802,13 +1065,199 @@ test_that("ctwasPipeline: dispatches assemble → est → screen → finemap and expect_equal(unname(capturedScreen$group_prior), c(0.1, 0.0001)) # finemap consumes screen's screened_region_data. expect_named(capturedFinemap$region_data, "block1") - # Output mirrors the ctwas_sumstats top-level shape. - expect_setequal( - names(out), - c("z_gene", "param", "finemap_res", "susie_alpha_res", - "region_data", "boundary_genes", "screen_res", - "region_info", "z_snp", "weights", "snp_map", - "LD_map", "LD_loader_fun", "snpinfo_loader_fun")) + # Output is a CtwasResult: one (gwasStudy, study, context, method) row. + expect_s4_class(out, "CtwasResult") + expect_equal(nrow(out), 1L) + expect_equal(getMethodNames(out), "susie") + expect_equal(getStudy(out), "Q1") + expect_equal(getContexts(out), "c1") + expect_equal(as.character(out$gwasStudy), "G1") + # The run's jointly-estimated param is carried on the row. + expect_equal(unname(getCtwasParam(out$entry[[1L]])$group_prior), + c(0.1, 0.0001)) + # getFinemap aggregates the per-gene rows, tagged with run identity. + fm <- getFinemap(out) + expect_equal(nrow(fm), 2L) + expect_true(all(c("gwasStudy", "study", "context", "method", "susie_pip") + %in% names(fm))) + expect_setequal(fm$id, c("block1|Q1|c1|t1|susie", "block2|Q1|c1|t1|susie")) + # getSusieAlpha aggregates the fuller per-effect table, likewise tagged. + sa <- getSusieAlpha(out) + expect_equal(nrow(sa), 2L) + expect_true(all(c("gwasStudy", "context", "susie_alpha") %in% names(sa))) +}) + +test_that("ctwasPipeline: keepSnps retains the SNP background as a dedicated row", { + skip_if_not_installed("ctwas") + inp <- .ctp_makeMultiBlockInputs() + local_mocked_bindings( + assemble_region_data = function(...) list(block1 = list(stub = TRUE)), + get_boundary_genes = function(...) data.frame(id = "t1", n_regions = 2L), + est_param = function(...) list(group_prior = c(c1 = 0.1, SNP = 1e-4), + group_prior_var = c(c1 = 5, SNP = 5)), + screen_regions = function(...) list( + screened_region_data = list(block1 = list(stub = TRUE)), screen_res_meta = "m"), + finemap_regions = function(...) list( + finemap_res = data.frame( + id = c("block1|Q1|c1|t1|susie", "chr1:100:G:A"), + type = c("gene", "SNP"), susie_pip = c(0.9, 0.02), + stringsAsFactors = FALSE), + susie_alpha_res = data.frame( + id = c("block1|Q1|c1|t1|susie", "chr1:100:G:A"), + type = c("gene", "SNP"), susie_alpha = c(0.9, 0.02), + stringsAsFactors = FALSE)), + .package = "ctwas") + local_mocked_bindings(extractBlockGenotypes = .ctp_mockExtractor(), + .package = "pecotmr") + out <- ctwasPipeline(inp$gwasSumStats, inp$twasWeights, keepSnps = TRUE) + # A gene (c1) row plus the dedicated SNP row. + expect_equal(nrow(out), 2L) + expect_setequal(getContexts(out), c("c1", "SNP")) + fm <- getFinemap(out) + expect_setequal(fm$type, c("gene", "SNP")) + expect_true("chr1:100:G:A" %in% fm$id) +}) + +# =========================================================================== +# mergeCtwasBoundaryRegions + ctwasPipeline(mergeBoundary = TRUE) +# =========================================================================== +.ctp_finemapResult <- function(hasLd = TRUE) { + fm <- data.frame(id = "block1|Q1|c1|t1|susie", type = "gene", + susie_pip = 0.90, stringsAsFactors = FALSE) + base <- list( + finemap_res = fm, + susie_alpha_res = cbind(fm, susie_alpha = 0.90), + region_data = list(block1 = list(stub = TRUE)), + region_info = data.frame(region_id = "block1", chrom = 1, + start = 1, stop = 1000), + z_snp = data.frame(id = "chr1:100:G:A", z = 1.0), + z_gene = data.frame(id = "block1|Q1|c1|t1|susie", z = 2.0), + weights = setNames(list(list(wgt = 1)), "block1|Q1|c1|t1|susie"), + snp_map = list(block1 = data.frame(id = "chr1:100:G:A")), + param = list(group_prior = c(c1 = 0.1, SNP = 1e-4), + group_prior_var = c(c1 = 5, SNP = 5)), + LD_map = data.frame(region_id = "block1", LD_file = "f", SNP_file = "f", + stringsAsFactors = FALSE)) + if (hasLd) { + base$LD_loader_fun <- function(...) NULL + base$snpinfo_loader_fun <- function(...) NULL + } + base +} + +.ctp_mergeReturn <- function() list( + updated_finemap_res = data.frame(id = "block1|Q1|c1|t1|susie", + type = "gene", susie_pip = 0.99, + stringsAsFactors = FALSE), + updated_susie_alpha_res = data.frame(id = "block1|Q1|c1|t1|susie", + susie_alpha = 0.99), + updated_region_data = list(merged = TRUE), + updated_region_info = data.frame(region_id = "block1_merged"), + updated_LD_map = data.frame(region_id = "block1_merged"), + updated_snp_map = list(block1_merged = 1)) + +test_that("mergeCtwasBoundaryRegions: LD path splices updated_* back + merge_res", { + skip_if_not_installed("ctwas") + captured <- NULL + local_mocked_bindings( + postprocess_region_merging = function(...) { captured <<- list(...); .ctp_mergeReturn() }, + .package = "ctwas") + out <- mergeCtwasBoundaryRegions(.ctp_finemapResult(hasLd = TRUE), + pipThresh = 0.5) + # LD path passed the loader closures + pip_thresh through. + expect_equal(captured$pip_thresh, 0.5) + expect_true(is.function(captured$LD_loader_fun)) + # updated_* spliced onto the result; merge_res carries the full postprocess out. + expect_equal(out$finemap_res$susie_pip, 0.99) + expect_equal(out$susie_alpha_res$susie_alpha, 0.99) + expect_equal(out$region_info$region_id, "block1_merged") + expect_true(!is.null(out$merge_res)) +}) + +test_that("mergeCtwasBoundaryRegions: no-LD path uses postprocess_region_merging_noLD", { + skip_if_not_installed("ctwas") + usedNoLd <- FALSE + local_mocked_bindings( + postprocess_region_merging = function(...) { .ctp_mergeReturn() }, + postprocess_region_merging_noLD = function(...) { usedNoLd <<- TRUE; .ctp_mergeReturn() }, + .package = "ctwas") + out <- mergeCtwasBoundaryRegions(.ctp_finemapResult(hasLd = FALSE)) + expect_true(usedNoLd) + expect_equal(out$finemap_res$susie_pip, 0.99) +}) + +test_that("mergeCtwasBoundaryRegions: empty first-pass finemap returns unchanged", { + skip_if_not_installed("ctwas") + fmr <- .ctp_finemapResult(); fmr$finemap_res <- fmr$finemap_res[0, ] + expect_message(out <- mergeCtwasBoundaryRegions(fmr), "no first-pass finemap") + expect_null(out$merge_res) +}) + +test_that("ctwasPipeline: mergeBoundary = TRUE re-fine-maps and flows into CtwasResult", { + skip_if_not_installed("ctwas") + inp <- .ctp_makeMultiBlockInputs() + merged <- FALSE + local_mocked_bindings( + assemble_region_data = function(...) list(block1 = list(stub = TRUE)), + get_boundary_genes = function(...) data.frame(id = "t1", n_regions = 2L), + est_param = function(...) list(group_prior = c(c1 = 0.1, SNP = 1e-4), + group_prior_var = c(c1 = 5, SNP = 5)), + screen_regions = function(...) list( + screened_region_data = list(block1 = list(stub = TRUE)), screen_res_meta = "m"), + finemap_regions = function(...) list( + finemap_res = data.frame(id = "block1|Q1|c1|t1|susie", type = "gene", + susie_pip = 0.90, stringsAsFactors = FALSE), + susie_alpha_res = data.frame(id = "block1|Q1|c1|t1|susie", + susie_alpha = 0.90)), + postprocess_region_merging = function(...) { + merged <<- TRUE + list(updated_finemap_res = data.frame(id = "block1|Q1|c1|t1|susie", + type = "gene", susie_pip = 0.99, stringsAsFactors = FALSE), + updated_susie_alpha_res = data.frame(id = "block1|Q1|c1|t1|susie", + susie_alpha = 0.99)) + }, + .package = "ctwas") + local_mocked_bindings(extractBlockGenotypes = .ctp_mockExtractor(), + .package = "pecotmr") + out <- ctwasPipeline(inp$gwasSumStats, inp$twasWeights, mergeBoundary = TRUE) + expect_true(merged) # merging ran + expect_s4_class(out, "CtwasResult") + # The post-merge PIP (0.99) is what the decomposition carries. + expect_equal(getFinemap(out)$susie_pip, 0.99) +}) + +# =========================================================================== +# asCtwasResult — structure a granular finemap result (the wrapper path) +# =========================================================================== +test_that("asCtwasResult: structures a granular finemap result into a CtwasResult", { + fmr <- .ctp_finemapResult() + fmr$z_snp$study <- "D1" # GWAS study read from z_snp + cr <- asCtwasResult(fmr) + expect_s4_class(cr, "CtwasResult") + expect_equal(nrow(cr), 1L) + expect_equal(as.character(cr$gwasStudy), "D1") + expect_equal(getStudy(cr), "Q1") # derived from the gene ids + expect_equal(getContexts(cr), "c1") + expect_equal(getMethodNames(cr), "susie") + expect_equal(getFinemap(cr)$susie_pip, 0.90) + expect_false(is.null(getSusieAlpha(cr))) +}) + +test_that("asCtwasResult: keepSnps adds the dedicated SNP row", { + fmr <- .ctp_finemapResult() + fmr$z_snp$study <- "D1" + fmr$finemap_res <- rbind(fmr$finemap_res, + data.frame(id = "chr1:100:G:A", type = "SNP", susie_pip = 0.02, + stringsAsFactors = FALSE)) + cr <- asCtwasResult(fmr, keepSnps = TRUE) + expect_setequal(getContexts(cr), c("c1", "SNP")) +}) + +test_that("asCtwasResult: errors when the weights mix methods", { + fmr <- .ctp_finemapResult() + fmr$weights <- setNames(list(1, 2), + c("block1|Q1|c1|t1|susie", "block1|Q1|c1|t2|lasso")) + expect_error(asCtwasResult(fmr), "mixes weight methods") }) test_that("estCtwasParam: fallbackToPrefit recovers from accurate-EM NaN divergence", { @@ -929,16 +1378,26 @@ test_that("ctwasPipeline: real-engine end-to-end on the bundled example panel", min_p_single_effect = 0, filter_L = FALSE))) - # ctwas_sumstats returns these 7 elements on success. - expect_named(res, c("z_gene", "param", "finemap_res", "susie_alpha_res", - "region_data", "boundary_genes", "screen_res", - "region_info", "z_snp", "weights", "snp_map", - "LD_map", "LD_loader_fun", "snpinfo_loader_fun"), - ignore.order = TRUE) - # The gene we passed in came through. - expect_true(any(grepl("study1\\|brain\\|ENSG_example\\|susie", res$z_gene$id))) - expect_true(all(c("susie_pip", "susie_alpha", "region_id") - %in% colnames(res$susie_alpha_res))) + # The one-shot pipeline returns a CtwasResult: a single (brain, study1, + # susie) row (single-context run, so no jointContexts). + expect_s4_class(res, "CtwasResult") + expect_equal(nrow(res), 1L) + expect_equal(getContexts(res), "brain") + expect_equal(getStudy(res), "study1") + expect_equal(getMethodNames(res), "susie") + expect_false("jointContexts" %in% names(res)) + # The gene we passed in was fine-mapped and shows up in the aggregated + # finemap table, tagged with the run identity. + fm <- getFinemap(res) + expect_true(!is.null(fm) && nrow(fm) > 0L) + expect_true(all(c("gwasStudy", "study", "context", "method", "susie_pip") + %in% names(fm))) + expect_true(any(grepl("study1\\|brain\\|ENSG_example\\|susie", fm$id))) + expect_true(all(fm$context == "brain")) + # The fuller per-effect table is retained and reconstructable. + sa <- getSusieAlpha(res) + expect_true(!is.null(sa) && nrow(sa) > 0L) + expect_true(all(c("susie_pip", "susie_alpha", "region_id") %in% names(sa))) }) # =========================================================================== @@ -1138,14 +1597,16 @@ test_that("assembleCtwasInputs: rejects an unnamed gwasSumStats list", { "named list keyed by region_id") }) -test_that("assembleCtwasInputs: rejects an unnamed twasWeights list", { +test_that("assembleCtwasInputs: an unnamed twasWeights list is a FLAT source (needs region)", { skip_if_not_installed("ctwas") ss <- .ctp_makeGwasSumstats() - tw <- .ctp_makeTwasWeights() + tw <- .ctp_makeTwasWeights() # no region provenance + # An unnamed list is now treated as a flat weight source to place by region; + # without region provenance that placement can't happen. expect_error( assembleCtwasInputs(gwasSumStats = list(block1 = ss, block2 = ss), - twasWeights = list(tw)), # unnamed - "named list keyed by region_id") + twasWeights = list(tw)), # unnamed -> flat + "no `region` provenance") }) test_that("assembleCtwasInputs: rejects a non-GwasSumStats list entry", { @@ -1159,14 +1620,14 @@ test_that("assembleCtwasInputs: rejects a non-GwasSumStats list entry", { "is not a GwasSumStats") }) -test_that("assembleCtwasInputs: rejects a non-TwasWeights list entry", { +test_that("assembleCtwasInputs: rejects a weight source that is neither TwasWeights nor QtlFineMappingResult", { skip_if_not_installed("ctwas") ss <- .ctp_makeGwasSumstats() expect_error( assembleCtwasInputs( gwasSumStats = list(block1 = ss, block2 = ss), - twasWeights = list(block1 = "not a TwasWeights")), - "is not a TwasWeights") + twasWeights = list(block1 = "not a weight source")), + "must be a TwasWeights or QtlFineMappingResult") }) test_that("assembleCtwasInputs: rejects a non-FineMappingResultBase fineMappingResult", { @@ -1412,3 +1873,111 @@ test_that(".ctwasSnpInfoForGwasBlock: returns an empty frame when the block has expect_equal(nrow(out), 0L) expect_setequal(colnames(out), colnames(panelInfo)) }) + +# =========================================================================== +# Coverage: internal error / edge branches +# =========================================================================== +test_that(".ctwasCombineWeightSources: type + emptiness guards; FMR list", { + expect_error(pecotmr:::.ctwasCombineWeightSources(42), "must be a TwasWeights") + expect_error(pecotmr:::.ctwasCombineWeightSources(list(NULL)), "empty weight source") + fmr <- .ctp_makeFmrWeightSource() # a QtlFineMappingResult + expect_s4_class(pecotmr:::.ctwasCombineWeightSources(list(fmr)), + "QtlFineMappingResult") +}) + +test_that(".ctwasGwasStudy / .ctwasGwasStudyFromZSnp: NA when no study present", { + expect_true(is.na(pecotmr:::.ctwasGwasStudy(list(a = data.frame(study = NA_character_))))) + expect_true(is.na(pecotmr:::.ctwasGwasStudyFromZSnp(NULL))) + expect_true(is.na(pecotmr:::.ctwasGwasStudyFromZSnp(data.frame(study = character(0))))) +}) + +test_that(".ctwasSubsetById / .ctwasSubsetSnp: NULL on empty / no type column", { + df <- data.frame(id = c("a", "b"), type = c("gene", "SNP"), stringsAsFactors = FALSE) + expect_null(pecotmr:::.ctwasSubsetById(df, "zzz")) + expect_null(pecotmr:::.ctwasSubsetById(NULL, "a")) + expect_null(pecotmr:::.ctwasSubsetSnp(data.frame(id = "a"))) # no type column + expect_equal(nrow(pecotmr:::.ctwasSubsetSnp(df)), 1L) +}) + +test_that(".ctwasRowsToResult / .ctwasMethodFromWeights: empty-input guards", { + expect_error(pecotmr:::.ctwasRowsToResult(list()), "no genes were modeled") + expect_error(pecotmr:::.ctwasMethodFromWeights(NULL), "carries no weights") + expect_error(pecotmr:::.ctwasMethodFromWeights(setNames(list(1, 2), + c("blk|Q1|c1|gA|susie", "blk|Q1|c1|gB|lasso"))), + "mixes weight methods") +}) + +test_that(".ctwasRunToRows: empty weights -> no rows; a context mixing studies errors", { + expect_equal( + pecotmr:::.ctwasRunToRows(list(weights = setNames(list(), character(0))), + "D1", "susie"), list()) + run <- list(weights = setNames(list(1, 2), + c("blk|Q1|c1|gA|susie", "blk|Q2|c1|gB|susie")), + finemap_res = NULL, susie_alpha_res = NULL, param = list(), + region_info = NULL) + expect_error(pecotmr:::.ctwasRunToRows(run, "D1", "susie"), + "mixes multiple QTL studies") +}) + +test_that(".ctwasBucketWeights: unplaced genes warn + drop; empty blocks skipped", { + mkE <- function() TwasWeightsEntry( + variantIds = vapply(1:5, .ctp_snpId, character(1)), + weights = c(0.1, 0.05, -0.2, 0.3, 0.0)) + tw <- TwasWeights( + study = c("Q1", "Q1"), context = c("c1", "c1"), trait = c("gA", "gX"), + method = c("susie", "susie"), entry = list(mkE(), mkE()), + region = GenomicRanges::GRanges(c("chr1", "chr9"), + IRanges::IRanges(c(100, 100), c(150, 150))), + ldSketch = .ctp_makeHandle()) + grid <- list(chr1_1_350 = .ctp_makeGwasSumstats(), + chr1_350_700 = .ctp_makeGwasSumstats()) + b <- expect_warning(pecotmr:::.ctwasBucketWeights(tw, grid), "fell in no LD block") + expect_named(b, "chr1_1_350") # 2nd block empty -> skipped + expect_equal(as.character(b[["chr1_1_350"]]$trait), "gA") +}) + +test_that(".ctwasBucketWeights: errors when no gene lands in any block", { + mkE <- function() TwasWeightsEntry( + variantIds = vapply(1:5, .ctp_snpId, character(1)), + weights = c(0.1, 0.05, -0.2, 0.3, 0.0)) + twOff <- TwasWeights( + study = "Q1", context = "c1", trait = "gX", method = "susie", + entry = list(mkE()), + region = GenomicRanges::GRanges("chr9", IRanges::IRanges(100, 150)), + ldSketch = .ctp_makeHandle()) + grid <- list(chr1_1_350 = .ctp_makeGwasSumstats(), + chr1_350_700 = .ctp_makeGwasSumstats()) + expect_error(suppressWarnings(pecotmr:::.ctwasBucketWeights(twOff, grid)), + "no gene placed") +}) + +test_that("assembleCtwasInputs: bare gwasSumStats + NULL twasWeights guards", { + skip_if_not_installed("ctwas") + ss <- .ctp_makeGwasSumstats() + expect_error(assembleCtwasInputs(ss, .ctp_makeTwasWeights()), + "NAMED LIST of GwasSumStats") + expect_error(assembleCtwasInputs(list(b1 = ss, b2 = ss), NULL), + "twasWeights` is required") +}) + +test_that("weight-source / method / gwas-study guards reject degenerate input", { + expect_error(pecotmr:::.ctwasRequireNamedLists(list(a = 1, b = 2), 42), + "must be a TwasWeights") + expect_error(pecotmr:::.ctwasResolveMethods(list()), "no method entries") + expect_error(pecotmr:::.ctwasGwasStudyFromZSnp(data.frame(study = c("D1", "D2"))), + "multiple GWAS studies") +}) + +test_that("finemapCtwasRegions: an empty screened-region set returns a NULL finemap", { + skip_if_not_installed("ctwas") + screenStub <- list( + screened_region_data = list(), + z_gene = NULL, region_data = list(), boundary_genes = NULL, screen_res = NULL, + param = list(group_prior = c(g = 0.1), group_prior_var = c(g = 5)), + region_info = data.frame(), z_snp = data.frame(), weights = list(), + snp_map = list(), LD_map = data.frame(), + LD_loader_fun = NULL, snpinfo_loader_fun = NULL) + out <- finemapCtwasRegions(screenStub) + expect_null(out$finemap_res) + expect_null(out$susie_alpha_res) +}) diff --git a/tests/testthat/test_fineMappingPipeline.R b/tests/testthat/test_fineMappingPipeline.R index bb341f0b..a67275c4 100644 --- a/tests/testthat/test_fineMappingPipeline.R +++ b/tests/testthat/test_fineMappingPipeline.R @@ -436,6 +436,24 @@ test_that(".rbindFineMappingResult: concatenates two GwasFineMappingResult colle expect_equal(nrow(out), 2L) }) +test_that("combineFineMappingResults: row-binds same-class collections; rejects mixed class/empty", { + e <- FineMappingEntry( + variantIds = "v1", susieFit = list(token = "susie"), + topLoci = data.frame(variant_id = "v1", pip = 0.5, stringsAsFactors = FALSE)) + g1 <- GwasFineMappingResult(study = "g1", method = "susie", + region_id = "r1", entry = list(e)) + g2 <- GwasFineMappingResult(study = "g1", method = "susie", + region_id = "r2", entry = list(e)) + out <- combineFineMappingResults(g1, g2) # variadic + expect_s4_class(out, "GwasFineMappingResult") + expect_equal(nrow(out), 2L) + expect_equal(nrow(combineFineMappingResults(list(g1, g2))), 2L) # list form + q <- QtlFineMappingResult(study = "s1", context = "c1", trait = "t1", + method = "susie", entry = list(e)) + expect_error(combineFineMappingResults(g1, q), "same concrete class") + expect_error(combineFineMappingResults(), "nothing to combine") +}) + # =========================================================================== # .fmExtractZN # =========================================================================== @@ -829,19 +847,19 @@ test_that("fineMappingPipeline(QtlDataset, cvFolds=0): leaves cvResult NULL", { }) # =========================================================================== -# Cross-validation internals: .fmMakeSamplePartition / .fmCrossValidate / +# Cross-validation internals: .cvSamplePartition / .fmWeightsCv / # .fmSliceCv / .fmAttachCv (unit-level counterparts to the cvFolds end-to-end # tests above). # =========================================================================== -test_that(".fmMakeSamplePartition partitions every sample into the requested folds", { - part <- pecotmr:::.fmMakeSamplePartition(paste0("s", 1:20), fold = 4L) +test_that(".cvSamplePartition partitions every sample into the requested folds", { + part <- pecotmr:::.cvSamplePartition(paste0("s", 1:20), fold = 4L) expect_setequal(part$Sample, paste0("s", 1:20)) expect_setequal(sort(unique(part$Fold)), 1:4) expect_equal(nrow(part), 20L) }) -test_that(".fmCrossValidate returns twasWeightsCv-shaped output keyed by snake method", { +test_that(".fmWeightsCv returns twasWeightsCv-shaped output keyed by snake method", { skip_if_not_installed("susieR") set.seed(42) n <- 60L; p <- 12L @@ -849,7 +867,7 @@ test_that(".fmCrossValidate returns twasWeightsCv-shaped output keyed by snake m dimnames = list(paste0("s", seq_len(n)), paste0("v", seq_len(p)))) y <- X[, 2] * 1.5 + rnorm(n, sd = 0.5) names(y) <- rownames(X) - cv <- pecotmr:::.fmCrossValidate( + cv <- pecotmr:::.fmWeightsCv( X, y, tokens = "susie", methodArgs = list(susie = list()), fold = 3L, coverage = 0.95, verbose = 0) @@ -868,15 +886,15 @@ test_that(".fmCrossValidate returns twasWeightsCv-shaped output keyed by snake m expect_gt(perf[1, "corr"], 0) }) -test_that(".fmCrossValidate reuses a supplied samplePartition verbatim", { +test_that(".fmWeightsCv reuses a supplied samplePartition verbatim", { skip_if_not_installed("susieR") set.seed(7) n <- 40L; p <- 8L X <- matrix(rnorm(n * p), n, p, dimnames = list(paste0("s", seq_len(n)), paste0("v", seq_len(p)))) y <- X[, 1] + rnorm(n, sd = 0.5); names(y) <- rownames(X) - part <- pecotmr:::.fmMakeSamplePartition(rownames(X), fold = 4L) - cv <- pecotmr:::.fmCrossValidate(X, y, tokens = "susie", + part <- pecotmr:::.cvSamplePartition(rownames(X), fold = 4L) + cv <- pecotmr:::.fmWeightsCv(X, y, tokens = "susie", methodArgs = list(susie = list()), fold = 4L, samplePartition = part, coverage = 0.95, verbose = 0) @@ -2297,9 +2315,9 @@ test_that(".fmTwasMethodKey: bare token without adapter returned unchanged", { expect_equal(pecotmr:::.fmTwasMethodKey("susie"), "susie") # adapter -> stripped }) -test_that(".fmCvMetricRow: < 3 usable predictions -> all-NA row", { - expect_true(all(is.na(pecotmr:::.fmCvMetricRow(c(1, 2), c(1, 2))))) # < 3 (1182) - ok <- pecotmr:::.fmCvMetricRow(c(1, 2, 3, 4, 5), c(1.1, 2, 2.9, 4, 5)) +test_that(".cvMetricRow: < 3 usable predictions -> all-NA row", { + expect_true(all(is.na(pecotmr:::.cvMetricRow(c(1, 2), c(1, 2))))) # < 3 (1182) + ok <- pecotmr:::.cvMetricRow(c(1, 2, 3, 4, 5), c(1.1, 2, 2.9, 4, 5)) expect_false(is.na(ok[["rsq"]])) }) @@ -2391,10 +2409,10 @@ test_that(".fmRelabelCs returns non-matching labels unchanged", { expect_equal(out, c("susie_3", "nomatch", "susie_0")) }) -test_that(".fmCrossValidate + .fmFoldWeights cover the mvSuSiE CV path (mocked fitter)", { +test_that(".fmWeightsCv + .fmFoldWeights cover the mvSuSiE CV path (mocked fitter)", { # Mock one level below the orchestration: fitMvsusie/mvsusieWeights return # canned outputs (sized to the per-fold training columns), so the real - # .fmFoldWeights mvsusie branch + .fmCrossValidate fold loop run at ~no cost. + # .fmFoldWeights mvsusie branch + .fmWeightsCv fold loop run at ~no cost. local_mocked_bindings( fitMvsusie = function(X, Y, ...) list(vn = colnames(X), R = ncol(as.matrix(Y))), mvsusieWeights = function(mvsusieFit = NULL, ...) @@ -2407,7 +2425,7 @@ test_that(".fmCrossValidate + .fmFoldWeights cover the mvSuSiE CV path (mocked f dimnames = list(paste0("s", 1:n), paste0("v", 1:p))) Y <- matrix(rnorm(n * R), n, R, dimnames = list(rownames(X), c("c1", "c2"))) - cv <- pecotmr:::.fmCrossValidate( + cv <- pecotmr:::.fmWeightsCv( X, Y, tokens = "mvsusie", methodArgs = list(mvsusie = list()), fold = 3L, coverage = 0.95, verbose = 0) expect_named(cv, c("samplePartition", "prediction", "performance")) @@ -2480,7 +2498,7 @@ test_that(".fmPostprocessOne errors when output carries no FineMappingEntry", { "FineMappingEntry payload") }) -test_that(".fmCrossValidate covers per-fold prior, NULL-weights, and no-overlap branches", { +test_that(".fmWeightsCv covers per-fold prior, NULL-weights, and no-overlap branches", { # mvPriorCv supplies a prior for fold "1" only: fold 1 takes the per-fold # prior (else-branch) and returns weights whose rownames don't overlap the # test columns (no-common `next`); fold 2 has no prior, so .fmFoldWeights @@ -2495,7 +2513,7 @@ test_that(".fmCrossValidate covers per-fold prior, NULL-weights, and no-overlap X <- matrix(rbinom(40, 2, 0.4), 20, 2, dimnames = list(paste0("s", 1:20), c("v1", "v2"))) Y <- matrix(rnorm(40), 20, 2, dimnames = list(rownames(X), c("c1", "c2"))) - cv <- pecotmr:::.fmCrossValidate( + cv <- pecotmr:::.fmWeightsCv( X, Y, tokens = "mvsusie", methodArgs = list(mvsusie = list()), fold = 2L, coverage = 0.95, verbose = 0, mvPriorCv = list("1" = list(priorVariance = diag(2)))) @@ -2549,31 +2567,31 @@ test_that(".fmFoldWeights covers mvPrior residual var, missing rownames, unknown expect_null(pecotmr:::.fmFoldWeights("bogus", X, Y, 0.95, list(), NULL)) }) -test_that(".fmCrossValidate returns NULL for empty tokens", { +test_that(".fmWeightsCv returns NULL for empty tokens", { X <- matrix(0, 10, 2, dimnames = list(paste0("s", 1:10), c("v1", "v2"))) - expect_null(pecotmr:::.fmCrossValidate( + expect_null(pecotmr:::.fmWeightsCv( X, matrix(0, 10, 1), tokens = character(0), methodArgs = list(), fold = 2L)) }) -test_that(".fmCrossValidate fills Y rownames and reports per-fold fit failures", { +test_that(".fmWeightsCv fills Y rownames and reports per-fold fit failures", { local_mocked_bindings(.fmFoldWeights = function(...) stop("boom"), .package = "pecotmr") X <- matrix(rbinom(40, 2, 0.4), 20, 2, dimnames = list(paste0("s", 1:20), c("v1", "v2"))) Y <- matrix(rnorm(20), 20, 1) # no rownames -> filled from X (1258) expect_message( - cv <- suppressWarnings(pecotmr:::.fmCrossValidate( + cv <- suppressWarnings(pecotmr:::.fmWeightsCv( X, Y, tokens = "susie", methodArgs = list(susie = list()), fold = 2L, verbose = 1)), "CV fold .* failed") # 1291-1294 }) -test_that(".fmCrossValidate skips a fold that holds out every sample", { +test_that(".fmWeightsCv skips a fold that holds out every sample", { X <- matrix(rbinom(40, 2, 0.4), 20, 2, dimnames = list(paste0("s", 1:20), c("v1", "v2"))) Y <- matrix(rnorm(20), 20, 1, dimnames = list(rownames(X), NULL)) sp <- data.frame(Sample = rownames(X), Fold = 1L) # single fold = all test -> 1273 - cv <- pecotmr:::.fmCrossValidate( + cv <- pecotmr:::.fmWeightsCv( X, Y, tokens = "susie", methodArgs = list(susie = list()), fold = 1L, samplePartition = sp, verbose = 0) expect_true(all(is.na(cv$prediction[["susie_predicted"]]))) diff --git a/tests/testthat/test_genotypeHandle.R b/tests/testthat/test_genotypeHandle.R index 9634cc8d..1482f048 100644 --- a/tests/testthat/test_genotypeHandle.R +++ b/tests/testthat/test_genotypeHandle.R @@ -356,6 +356,40 @@ test_that("genoMeta meta-file resolves payloads relative to its own directory", expect_equal(nrow(h@snpInfo), 2L * 349L) }) +test_that(".parseChromMeta matches chrom/path columns by name (order- and extra-column-tolerant)", { + mk <- function(lines) { f <- tempfile(fileext = ".tsv"); writeLines(lines, f); f } + # extra legacy start/end columns + path not in position 2 + m1 <- pecotmr:::.parseChromMeta( + mk(c("#chrom\tstart\tend\tpath", "21\t1\t9\t/abs/a.bed", "22\t1\t9\t/abs/b.bed"))) + expect_equal(m1$chrom, c("21", "22")) + expect_equal(m1$path, c("/abs/a.bed", "/abs/b.bed")) + # reordered header (path first) + alias column names + m2 <- pecotmr:::.parseChromMeta( + mk(c("genotype\tchr", "/abs/a.bed\t21", "/abs/b.bed\t22"))) + expect_equal(m2$chrom, c("21", "22")) + expect_equal(m2$path, c("/abs/a.bed", "/abs/b.bed")) +}) + +test_that(".parseChromMeta falls back to the first two columns when no alias header matches", { + mk <- function(lines) { f <- tempfile(fileext = ".tsv"); writeLines(lines, f); f } + m <- pecotmr:::.parseChromMeta( + mk(c("foo\tbar", "21\t/abs/a.bed", "22\t/abs/b.bed"))) + expect_equal(m$chrom, c("21", "22")) + expect_equal(m$path, c("/abs/a.bed", "/abs/b.bed")) +}) + +test_that(".parseChromMeta errors when a chromosome maps to multiple payloads", { + mk <- function(lines) { f <- tempfile(fileext = ".tsv"); writeLines(lines, f); f } + expect_error( + pecotmr:::.parseChromMeta( + mk(c("#chr\tpath", "21\t/abs/a.bed", "21\t/abs/b.bed"))), + "map to multiple genotype payloads") + # same guard on the named-vector path + expect_error( + pecotmr:::.parseChromMeta(c("21" = "/abs/a.bed", "21" = "/abs/b.bed")), + "map to multiple genotype payloads") +}) + test_that("genoMeta errors on mismatched samples across shards", { skip_if_not_installed("snpStats") # protocol_example.genotype is a different sample panel. diff --git a/tests/testthat/test_jointEngine.R b/tests/testthat/test_jointEngine.R index a1e84a4b..e258218a 100644 --- a/tests/testthat/test_jointEngine.R +++ b/tests/testthat/test_jointEngine.R @@ -83,10 +83,10 @@ test_that(".runJointCell: cross-context FM uses the per-fold mr.mash CV prior", rescaleCovW0 = function(w0) c(K = 1), fitMvsusie = function(...) list(), .fmPostprocessOne = .je_mockPostprocess, - .fmCrossValidate = function(X, Y, tokens, methodArgs, fold, + .fmWeightsCv = function(X, Y, tokens, methodArgs, fold, samplePartition = NULL, coverage = 0.95, pos = NULL, verbose = 1, mvPrior = NULL, - mvPriorCv = NULL) { + mvPriorCv = NULL, numThreads = 1) { captured <<- list(mvPriorCv = mvPriorCv, samplePartition = samplePartition) list(samplePartition = samplePartition, prediction = list(), performance = list()) @@ -845,7 +845,7 @@ test_that("fitJointGroup(Individual, Fm): fsusie honest per-fold CV is attached" fitFsusie = function(...) list(), fsusieWeights = function(fsusieFit, variantIds) NULL, .fmPostprocessOne = .je_mockPostprocess, - .fmCrossValidate = function(X, Y, token, methodArgs, fold, ...) { + .fmWeightsCv = function(X, Y, token, methodArgs, fold, ...) { cvCalled <<- TRUE; list(samplePartition = NULL) }, .fmSliceCv = function(cv, token) list(prediction = NULL), .fmAttachCv = function(e, cv) e, diff --git a/tests/testthat/test_mashWrapper.R b/tests/testthat/test_mashWrapper.R index f2bbbffa..f5f33d0c 100644 --- a/tests/testthat/test_mashWrapper.R +++ b/tests/testthat/test_mashWrapper.R @@ -950,3 +950,94 @@ test_that("mergeMashData returns resData when oneData is an empty list", { d1 <- list(random = data.frame(a = 1:3)) expect_equal(mergeMashData(d1, list()), d1) }) + +# =========================================================================== +# qtlSumStatsFromZMatrix +# =========================================================================== + +.qszm_gh <- function() { + new("GenotypeHandle", + path = "/tmp/sketch.gds", format = "gds", + snpInfo = data.frame(SNP = "v1", CHR = "1", BP = 100L, + A1 = "A", A2 = "G", stringsAsFactors = FALSE), + nSamples = 10L, sampleIds = paste0("s", seq_len(10L)), pgenPtr = NULL) +} + +test_that("qtlSumStatsFromZMatrix: one row per context, Z preserved verbatim", { + z <- matrix(c(1.1, -2.2, 0.3, 0.4, -0.5, 0.6), nrow = 3, + dimnames = list(c("chr1:100:A:G", "chr1:200:C:T", "chr2:50:A:T"), + c("brain", "liver"))) + qss <- qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh()) + expect_s4_class(qss, "QtlSumStats") + expect_equal(nrow(qss), 2L) + expect_equal(as.character(qss$context), c("brain", "liver")) + expect_equal(unique(as.character(qss$study)), "s1") + expect_equal(unique(as.character(qss$trait)), "mash") + # Z column for each context matches the input matrix column (values only; + # the mcols column is unnamed whereas z[, j] carries the row ids as names). + expect_equal(S4Vectors::mcols(qss$entry[[1]])$Z, unname(z[, 1])) + expect_equal(S4Vectors::mcols(qss$entry[[2]])$Z, unname(z[, 2])) +}) + +test_that("qtlSumStatsFromZMatrix: decodes chrom/pos from ids, synthesises when they don't parse", { + z <- matrix(rnorm(2), ncol = 1, + dimnames = list(c("chr1:250:A:G", "not_a_variant"), "ctx")) + qss <- qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh()) + e <- qss$entry[[1]] + expect_equal(GenomicRanges::start(e)[1], 250L) # decoded + expect_true(GenomicRanges::start(e)[2] >= 1L) # synthetic fallback + # un-parseable chrom falls back to chr1 (never NA) + expect_false(any(is.na(as.character(GenomicRanges::seqnames(e))))) + # SNP ids are carried through unchanged + expect_equal(S4Vectors::mcols(e)$SNP, rownames(z)) +}) + +test_that("qtlSumStatsFromZMatrix: NULL rownames get synthetic variant ids", { + z <- matrix(rnorm(4), nrow = 2, dimnames = list(NULL, c("a", "b"))) + qss <- qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh()) + expect_equal(S4Vectors::mcols(qss$entry[[1]])$SNP, c("var1", "var2")) +}) + +test_that("qtlSumStatsFromZMatrix: placeholders and pass-through qcInfo are set", { + z <- matrix(rnorm(6), nrow = 3, dimnames = list(NULL, c("x", "y"))) + qss <- qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh(), + n = 500L, a1 = "T", a2 = "C", role = "strong") + mc <- S4Vectors::mcols(qss$entry[[1]]) + expect_equal(unique(mc$A1), "T") + expect_equal(unique(mc$A2), "C") + expect_equal(unique(mc$N), 500L) + expect_equal(qss@qcInfo$role, "strong") + expect_equal(length(qss@qcInfo$entryAudit), 2L) # one slot per context +}) + +test_that("qtlSumStatsFromZMatrix: columns can map to traits or context x trait pairs", { + z <- matrix(rnorm(6), nrow = 3, dimnames = list(NULL, c("geneA", "geneB"))) + # columns as traits: constant context, one trait per column + qss <- qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh(), + context = "brain", trait = colnames(z)) + expect_equal(as.character(qss$context), c("brain", "brain")) + expect_equal(as.character(qss$trait), c("geneA", "geneB")) + # columns as (context, trait) pairs + qss2 <- qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh(), + context = c("brain", "liver"), + trait = c("geneA", "geneA")) + expect_equal(as.character(qss2$context), c("brain", "liver")) + expect_equal(as.character(qss2$trait), c("geneA", "geneA")) +}) + +test_that("qtlSumStatsFromZMatrix: a condition label of the wrong length errors", { + z <- matrix(rnorm(6), nrow = 3, dimnames = list(NULL, c("a", "b"))) + expect_error( + qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh(), + trait = c("t1", "t2", "t3")), + "must be length 1 or ncol") +}) + +test_that("qtlSumStatsFromZMatrix: rejects non-matrix input and unlabelled conditions", { + expect_error(qtlSumStatsFromZMatrix(1:5, study = "s1", ldSketch = .qszm_gh()), + "variants x conditions matrix") + # no colnames -> the default context = colnames(z) is NULL + z <- matrix(rnorm(4), nrow = 2) + expect_error(qtlSumStatsFromZMatrix(z, study = "s1", ldSketch = .qszm_gh()), + "column names") +}) diff --git a/tests/testthat/test_sldscWrapper.R b/tests/testthat/test_sldscWrapper.R index 3598749f..a16cb6e8 100644 --- a/tests/testthat/test_sldscWrapper.R +++ b/tests/testthat/test_sldscWrapper.R @@ -448,6 +448,44 @@ test_that("metaSldscRandom works for enrichstat", { expect_true(is.finite(result$mean)) }) +# A per-trait $summary carrying the wide (Single/Joint) columns that +# sldscPostprocessingPipeline emits and .sldscViewForMeta consumes. +.sms_makeWideTrait <- function() list(summary = data.frame( + target = c("A_0", "B_0"), + tauStarSingle = c(2, 3), tauStarSeSingle = c(.5, .5), + enrichmentSingle = c(1.5, 2), enrichmentSeSingle = c(.3, .3), + enrichmentPSingle = c(.01, .02), + enrichstatSingle = c(1, 1.2), enrichstatSeSingle = c(.2, .2), + tauStarJoint = c(1.8, 2.5), tauStarSeJoint = c(.4, .4), + enrichmentJoint = c(1.4, 1.9), enrichmentSeJoint = c(.3, .3), + enrichmentPJoint = c(.02, .03), + enrichstatJoint = c(.9, 1.1), enrichstatSeJoint = c(.2, .2), + stringsAsFactors = FALSE)) + +test_that("sldscSubsetMeta metas a chosen trait/category subset with the four blocks", { + res <- list( + per_trait = list(t1 = .sms_makeWideTrait(), t2 = .sms_makeWideTrait(), + t3 = .sms_makeWideTrait()), + params = list(target_categories = c("A_0", "B_0"))) + out <- sldscSubsetMeta(res, c("t1", "t2"), targetCategories = "A_0") + expect_setequal(names(out), + c("tau_star_single", "tau_star_joint", "enrichment", "enrichstat")) + expect_equal(names(out$tau_star_single), "A_0") # only the chosen category + expect_equal(out$tau_star_single[["A_0"]]$nTraits, 2L) # only the chosen traits + expect_true(is.finite(out$tau_star_joint[["A_0"]]$mean)) +}) + +test_that("sldscSubsetMeta defaults categories from params and validates traits", { + res <- list( + per_trait = list(t1 = .sms_makeWideTrait(), t2 = .sms_makeWideTrait()), + params = list(target_categories = c("A_0", "B_0"))) + out <- sldscSubsetMeta(res, c("t1", "t2")) # targetCategories = NULL + expect_setequal(names(out$enrichment), c("A_0", "B_0")) + expect_error(sldscSubsetMeta(res, c("t1", "ghost")), "absent from") + expect_error(sldscSubsetMeta(list(per_trait = res$per_trait), c("t1")), + "no `targetCategories`") +}) + test_that("metaSldscRandom returns NA with < 2 traits", { pt <- .make_per_trait_meta(nTraits = 1) result <- metaSldscRandom(pt, "A_0", "tauStar") diff --git a/tests/testthat/test_tupleSelectors.R b/tests/testthat/test_tupleSelectors.R index 7764b00b..9d9a7d09 100644 --- a/tests/testthat/test_tupleSelectors.R +++ b/tests/testthat/test_tupleSelectors.R @@ -167,3 +167,137 @@ test_that(".tupleSelectRowGwasFmr: ambiguous multi-match (no region) lists candi pecotmr:::.tupleSelectRowGwasFmr(multi, study = "g1", method = "susie"), "pass `region` to disambiguate") }) + +# =========================================================================== +# .fmrRowsMatching (aggregate row selector -- never errors on ambiguity) +# =========================================================================== + +test_that(".fmrRowsMatching: no selectors returns every row", { + df <- data.frame(study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + stringsAsFactors = FALSE) + expect_equal(pecotmr:::.fmrRowsMatching(df), c(1L, 2L)) +}) + +test_that(".fmrRowsMatching: matches a subset without erroring on ambiguity", { + df <- data.frame(study = c("s1", "s1", "s2"), context = c("c1", "c2", "c1"), + trait = c("t1", "t1", "t1"), method = c("susie", "susie", "susie"), + stringsAsFactors = FALSE) + expect_equal(pecotmr:::.fmrRowsMatching(df, study = "s1"), c(1L, 2L)) + expect_equal(pecotmr:::.fmrRowsMatching(df, context = "c1"), c(1L, 3L)) +}) + +test_that(".fmrRowsMatching: selectors on absent columns are ignored", { + # GWAS-shaped frame has no context/trait column; passing context must not + # error and must not constrain the result. + gwas <- data.frame(study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("r1", "r2"), stringsAsFactors = FALSE) + expect_equal(pecotmr:::.fmrRowsMatching(gwas, context = "c1"), c(1L, 2L)) +}) + +test_that(".fmrRowsMatching: `region` matches the region_id column", { + gwas <- data.frame(study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("r1", "r2"), stringsAsFactors = FALSE) + expect_equal(pecotmr:::.fmrRowsMatching(gwas, region = "r2"), 2L) +}) + +test_that(".fmrRowsMatching: a vector selector matches any listed value", { + df <- data.frame(study = c("s1", "s2", "s3"), context = c("c1", "c2", "c3"), + trait = c("t1", "t1", "t1"), method = c("susie", "susie", "susie"), + stringsAsFactors = FALSE) + expect_equal(pecotmr:::.fmrRowsMatching(df, study = c("s1", "s3")), c(1L, 3L)) +}) + +# =========================================================================== +# .fmrRowMetadata (stable 5-column identity frame) +# =========================================================================== + +test_that(".fmrRowMetadata: emits all five identity columns, NA-filling absent ones", { + qtl <- data.frame(study = c("s1", "s1"), context = c("c1", "c2"), + trait = c("t1", "t1"), method = c("susie", "susie"), + stringsAsFactors = FALSE) + m <- pecotmr:::.fmrRowMetadata(qtl) + expect_equal(names(m), + c("study", "context", "trait", "region_id", "method")) + expect_equal(m$context, c("c1", "c2")) + expect_true(all(is.na(m$region_id))) # QTL frame has no region_id +}) + +test_that(".fmrRowMetadata: GWAS frame NA-fills context/trait, keeps region_id", { + gwas <- data.frame(study = c("g1", "g1"), method = c("susie", "susie"), + region_id = c("r1", "r2"), stringsAsFactors = FALSE) + m <- pecotmr:::.fmrRowMetadata(gwas) + expect_equal(m$region_id, c("r1", "r2")) + expect_true(all(is.na(m$context))) + expect_true(all(is.na(m$trait))) +}) + +test_that(".fmrRowMetadata: zero-row input yields a zero-row 5-column frame", { + empty <- data.frame(study = character(0), method = character(0), + region_id = character(0), stringsAsFactors = FALSE) + m <- pecotmr:::.fmrRowMetadata(empty) + expect_equal(nrow(m), 0L) + expect_equal(names(m), + c("study", "context", "trait", "region_id", "method")) +}) + +# =========================================================================== +# .rbindAligned (union-of-columns rbind) +# =========================================================================== + +test_that(".rbindAligned: a single part is returned unchanged", { + df <- data.frame(a = 1:2, b = c("x", "y"), stringsAsFactors = FALSE) + expect_identical(pecotmr:::.rbindAligned(list(df)), df) +}) + +test_that(".rbindAligned: identical columns rbind and preserve type", { + p1 <- data.frame(a = 1L, b = "x", stringsAsFactors = FALSE) + p2 <- data.frame(a = 2L, b = "y", stringsAsFactors = FALSE) + out <- pecotmr:::.rbindAligned(list(p1, p2)) + expect_equal(nrow(out), 2L) + expect_type(out$a, "integer") + expect_equal(out$b, c("x", "y")) +}) + +test_that(".rbindAligned: differing columns align on the union, NA-filling gaps", { + p1 <- data.frame(a = 1L, b = "x", stringsAsFactors = FALSE) + p2 <- data.frame(a = 2L, c = "z", stringsAsFactors = FALSE) + out <- pecotmr:::.rbindAligned(list(p1, p2)) + expect_equal(sort(names(out)), c("a", "b", "c")) + expect_equal(out$a, c(1L, 2L)) + expect_equal(out$b, c("x", NA)) + expect_equal(out$c, c(NA, "z")) +}) + +# =========================================================================== +# Coverage: region-column + rbind helpers +# =========================================================================== +test_that(".naLikeColumn: a non-character/non-GRanges exemplar pads with logical NA", { + expect_identical(pecotmr:::.naLikeColumn(TRUE, 3L), rep(NA, 3L)) +}) + +test_that(".rbindCollections: all-NULL input returns NULL", { + expect_null(pecotmr:::.rbindCollections(list(NULL, NULL))) +}) + +test_that(".getRegionColumn: an absent region column yields an empty GRanges", { + gr <- pecotmr:::.getRegionColumn(S4Vectors::DataFrame(a = 1:2)) + expect_s4_class(gr, "GRanges") + expect_length(gr, 0L) +}) + +test_that(".appendRegionCol: region must be a GRanges of matching length", { + expect_error(pecotmr:::.appendRegionCol(list(), "notgranges", 1L), + "must be a GRanges") + expect_error( + pecotmr:::.appendRegionCol(list(), + GenomicRanges::GRanges(c("chr1", "chr1"), IRanges::IRanges(1:2, 1:2)), 1L), + "same length as") +}) + +test_that(".validateRegionColumn: reports a non-GRanges region column", { + expect_equal( + pecotmr:::.validateRegionColumn(S4Vectors::DataFrame(region = c("x", "y"))), + "'region' column must be a GRanges") + expect_length(pecotmr:::.validateRegionColumn(S4Vectors::DataFrame(a = 1L)), 0L) +}) diff --git a/tests/testthat/test_twasWeights.R b/tests/testthat/test_twasWeights.R index 50e0feff..9e02d7e5 100644 --- a/tests/testthat/test_twasWeights.R +++ b/tests/testthat/test_twasWeights.R @@ -362,75 +362,6 @@ test_that("twasWeights: enetWeights produces correct structure with real glmnet" # # =========================================================================== -test_that("twasWeightsCv: fold must be positive integer", { - d <- make_data() - expect_error( - twasWeightsCv(d$X, d$Y, fold = -1), - "Invalid value for 'fold'" - ) - expect_error( - twasWeightsCv(d$X, d$Y, fold = 0), - "Invalid value for 'fold'" - ) -}) - -test_that("twasWeightsCv: fold as string errors", { - d <- make_data() - expect_error( - twasWeightsCv(d$X, d$Y, fold = "abc"), - "Invalid value for 'fold'" - ) -}) - -test_that("twasWeightsCv: X must be a matrix", { - d <- make_data() - expect_error( - twasWeightsCv(as.data.frame(d$X), d$Y, fold = 5), - "X must be a matrix" - ) -}) - -test_that("twasWeightsCv: Y must be a matrix or vector", { - d <- make_data() - # In R, is.vector(list(...)) returns TRUE, so a list passes the initial - # type check and gets converted via matrix(). The resulting 3-row matrix - # mismatches X's 50 rows, triggering the row count error. - expect_error( - twasWeightsCv(d$X, list(1, 2, 3), fold = 5), - "The number of rows in X and Y must be the same" - ) -}) - -test_that("twasWeightsCv: fold or sample_partitions must be provided", { - d <- make_data() - expect_error( - twasWeightsCv(d$X, d$Y), - "Either 'fold' or 'samplePartitions' must be provided" - ) -}) - -test_that("twasWeightsCv: row count mismatch between X and Y errors", { - d <- make_data() - Y_wrong <- d$Y[1:20, , drop = FALSE] - expect_error( - twasWeightsCv(d$X, Y_wrong, fold = 5), - "The number of rows in X and Y must be the same" - ) -}) - -test_that("twasWeightsCv: Y as vector is accepted and converted", { - d <- make_data() - y_vec <- as.numeric(d$Y) - - # With NULL weight_methods, should return just samplePartition - expect_message( - result <- twasWeightsCv(d$X, y_vec, fold = 3, weightMethods = NULL), - "Y converted to matrix" - ) - expect_true(is.list(result)) - expect_true("samplePartition" %in% names(result)) -}) - test_that("twasWeightsCv: NULL weight_methods returns only samplePartition", { d <- make_data() result <- twasWeightsCv(d$X, d$Y, fold = 3, weightMethods = NULL) @@ -465,173 +396,6 @@ test_that("twasWeightsCv: character weight_methods are accepted", { expect_true("prediction" %in% names(result)) }) -test_that("twasWeightsCv: max_num_variants subsets X columns", { - d <- make_data(n = 50, p = 20) - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - expect_message( - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list()), - maxNumVariants = 5 - ), - "Randomly selecting 5 out of 20" - ) - expect_true(is.list(result)) - # The result should have prediction and performance entries for the one method - expect_true("prediction" %in% names(result)) - expect_true("performance" %in% names(result)) - expect_true("lassoPredicted" %in% names(result$prediction)) - # Weight method returned zero weights, so predictions should exist but all be zero - pred <- result$prediction[["lassoPredicted"]] - expect_equal(nrow(pred), nrow(d$X)) - expect_true(all(pred == 0)) -}) - -test_that("twasWeightsCv: max_num_variants with variants_to_keep", { - d <- make_data(n = 50, p = 20) - keep_vars <- colnames(d$X)[1:3] - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - expect_message( - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list()), - maxNumVariants = 8, - variantsToKeep = keep_vars - ), - "Including 3 specified variants" - ) - expect_true(is.list(result)) -}) - -test_that("twasWeightsCv: sample_partitions with mismatched samples errors", { - d <- make_data() - bad_partitions <- data.frame( - Sample = c("nonexistent_1", "nonexistent_2"), - Fold = c(1, 2), - stringsAsFactors = FALSE - ) - expect_error( - twasWeightsCv(d$X, d$Y, samplePartitions = bad_partitions), - "Some samples in 'samplePartitions' do not match" - ) -}) - -test_that("twasWeightsCv: provided sample_partitions are used", { - d <- make_data(n = 20, p = 5) - sp <- data.frame( - Sample = rownames(d$X), - Fold = rep(1:2, each = 10), - stringsAsFactors = FALSE - ) - result <- twasWeightsCv(d$X, d$Y, samplePartitions = sp, weightMethods = NULL) - expect_equal(result$samplePartition, sp) -}) - -test_that("twasWeightsCv: rownames are auto-generated when missing", { - set.seed(42) - n <- 30 - p <- 5 - X <- matrix(rnorm(n * p), nrow = n, ncol = p) - Y <- matrix(rnorm(n), ncol = 1) - # No row names on X or Y - - result <- twasWeightsCv(X, Y, fold = 2, weightMethods = NULL) - sp <- result$samplePartition - - # Should have auto-generated sample names - expect_true(all(grepl("^sample_", sp$Sample))) -}) - -test_that("twasWeightsCv: colnames are auto-generated when missing", { - set.seed(42) - n <- 30 - p <- 5 - X <- matrix(rnorm(n * p), nrow = n, ncol = p) - rownames(X) <- paste0("s", 1:n) - Y <- matrix(rnorm(n), ncol = 1) - rownames(Y) <- paste0("s", 1:n) - # No col names on X or Y - - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - result <- twasWeightsCv(X, Y, fold = 2, weightMethods = list(lassoWeights = list())) - # Should not error; column names are auto-generated - expect_true(!is.null(result$prediction)) -}) - -test_that("twasWeightsCv: fold and sample_partitions mismatch prints message", { - d <- make_data(n = 20, p = 5) - sp <- data.frame( - Sample = rownames(d$X), - Fold = rep(1:4, each = 5), - stringsAsFactors = FALSE - ) - # fold=2 but sample_partitions has 4 folds - expect_message( - twasWeightsCv(d$X, d$Y, fold = 2, samplePartitions = sp, weightMethods = NULL), - "fold number provided does not match" - ) -}) - -test_that("twasWeightsCv: zero-variance predictions yield NA metrics with message", { - d <- make_data() - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - expect_message( - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list()) - ), - "zero variance" - ) - perf <- result$performance[["lassoPerformance"]] - expect_true(all(is.na(perf))) -}) - -test_that("twasWeightsCv: performance names use _performance suffix", { - d <- make_data() - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)), - enetWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list(), enetWeights = list()) - ) - expect_equal( - sort(names(result$performance)), - sort(c("lassoPerformance", "enetPerformance")) - ) -}) - -test_that("twasWeightsCv: prediction names use _predicted suffix", { - d <- make_data() - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)), - enetWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list(), enetWeights = list()) - ) - expect_equal( - sort(names(result$prediction)), - sort(c("lassoPredicted", "enetPredicted")) - ) -}) - # --------------------------------------------------------------------------- # CV with real lassoWeights (integration test) # --------------------------------------------------------------------------- @@ -655,14 +419,14 @@ test_that("twasWeightsCv: basic CV with lassoWeights produces correct metrics st expect_true("timeElapsed" %in% names(result)) # Prediction name transformation - expect_equal(names(result$prediction), "lassoPredicted") + expect_equal(names(result$prediction), "lasso_predicted") # Prediction dimensions should match Y - pred <- result$prediction[["lassoPredicted"]] + pred <- result$prediction[["lasso_predicted"]] expect_equal(dim(pred), dim(d$Y)) # Performance table structure - perf <- result$performance[["lassoPerformance"]] + perf <- result$performance[["lasso_performance"]] expect_true(is.matrix(perf)) expect_equal(colnames(perf), c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE")) expect_equal(nrow(perf), ncol(d$Y)) @@ -671,23 +435,6 @@ test_that("twasWeightsCv: basic CV with lassoWeights produces correct metrics st expect_true(perf[1, "corr"] > 0) }) -test_that("twasWeightsCv: metrics table has correct column names (mocked)", { - d <- make_data() - local_mocked_bindings( - lassoWeights = function(X, y, ...) { - # Return weights that produce non-zero-variance predictions - rep(0.1, ncol(X)) - } - ) - set.seed(42) - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list()) - ) - perf <- result$performance[["lassoPerformance"]] - expect_equal(colnames(perf), c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE")) -}) - test_that("twasWeightsCv: multiple real methods produce per-method metrics", { skip_if_not_installed("glmnet") d <- make_data(n = 50, p = 10) @@ -704,27 +451,10 @@ test_that("twasWeightsCv: multiple real methods produce per-method metrics", { expect_equal(length(result$prediction), 2) expect_equal(length(result$performance), 2) - expect_true("lassoPredicted" %in% names(result$prediction)) - expect_true("enetPredicted" %in% names(result$prediction)) - expect_true("lassoPerformance" %in% names(result$performance)) - expect_true("enetPerformance" %in% names(result$performance)) -}) - -test_that("twasWeightsCv: all samples appear exactly once in predictions", { - skip_if_not_installed("glmnet") - d <- make_data(n = 50, p = 10) - - set.seed(77) - result <- twasWeightsCv( - d$X, d$Y, - fold = 5, - weightMethods = list(lassoWeights = list()) - ) - - pred <- result$prediction[["lassoPredicted"]] - # No NAs -- every sample was predicted in exactly one fold - expect_false(any(is.na(pred))) - expect_equal(nrow(pred), nrow(d$X)) + expect_true("lasso_predicted" %in% names(result$prediction)) + expect_true("enet_predicted" %in% names(result$prediction)) + expect_true("lasso_performance" %in% names(result$performance)) + expect_true("enet_performance" %in% names(result$performance)) }) # =========================================================================== @@ -750,11 +480,11 @@ test_that("twasWeightsCv: multivariate Y with multiple columns", { weightMethods = list(lassoWeights = list()) ) - pred <- result$prediction[["lassoPredicted"]] + pred <- result$prediction[["lasso_predicted"]] expect_equal(ncol(pred), 2) expect_equal(nrow(pred), 50) - perf <- result$performance[["lassoPerformance"]] + perf <- result$performance[["lasso_performance"]] expect_equal(nrow(perf), 2) expect_equal(rownames(perf), c("outcome_1", "outcome_2")) }) @@ -824,69 +554,6 @@ test_that("twasWeights: SuSiE-inf is fitted before and initializes ordinary SuSi # twasWeightsCv: extra split_data / sample-name / variant-selection branches # =========================================================================== -test_that("twasWeightsCv: split_data errors when a fold leaves train or test empty", { - d <- make_data(n = 10, p = 5) - # All samples in fold 1 -> with fold = 1, every sample is a "test" row and - # the train set has zero rows, hitting the split_data zero-row stop. - sp <- data.frame( - Sample = rownames(d$X), - Fold = rep(1L, nrow(d$X)), - stringsAsFactors = FALSE - ) - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)) - ) - expect_error( - suppressMessages(twasWeightsCv( - d$X, d$Y, - samplePartitions = sp, - weightMethods = list(lassoWeights = list()) - )), - "One of the datasets \\(train or test\\) has zero rows" - ) -}) - -test_that("twasWeightsCv: rownames(X) get reassigned to rownames(Y) when they differ", { - d <- make_data(n = 20, p = 5) - rownames(d$X) <- paste0("xname_", seq_len(nrow(d$X))) # differ from rownames(Y) - set.seed(42) - result <- twasWeightsCv(d$X, d$Y, fold = 2, weightMethods = NULL) - # samplePartition$Sample should now use rownames(Y), not rownames(X) - expect_true(all(result$samplePartition$Sample %in% rownames(d$Y))) - expect_false(any(grepl("^xname_", result$samplePartition$Sample))) -}) - -test_that("twasWeightsCv: sampleNames taken from Y when only Y has rownames", { - set.seed(42) - n <- 20; p <- 5 - X <- matrix(rnorm(n * p), nrow = n, ncol = p) - rownames(X) <- NULL - Y <- matrix(rnorm(n), ncol = 1) - rownames(Y) <- paste0("yonly_", seq_len(n)) - result <- twasWeightsCv(X, Y, fold = 2, weightMethods = NULL) - expect_true(all(grepl("^yonly_", result$samplePartition$Sample))) -}) - -test_that("twasWeightsCv: variants_to_keep >= max_num_variants samples from variants_to_keep only", { - d <- make_data(n = 50, p = 20) - # 10 keep variants, maxNumVariants = 5 => length(variants_to_keep) >= max_num_variants - keep_vars <- colnames(d$X)[1:10] - local_mocked_bindings( - lassoWeights = function(X, y, ...) rep(0, ncol(X)) - ) - set.seed(42) - expect_message( - result <- twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list()), - maxNumVariants = 5, - variantsToKeep = keep_vars - ), - "Randomly selecting 5 out of 10 input variants" - ) - expect_true("prediction" %in% names(result)) -}) - test_that("twasWeightsCv: NA values in Y trigger NA-removal branch in metrics", { set.seed(42) n <- 30; p <- 5 @@ -908,7 +575,7 @@ test_that("twasWeightsCv: NA values in Y trigger NA-removal branch in metrics", X, Y, fold = 2, weightMethods = list(lassoWeights = list()) ) - perf <- result$performance[["lassoPerformance"]] + perf <- result$performance[["lasso_performance"]] # NA-removal branch ran; metrics should be finite (not all-NA) expect_true(is.finite(perf[1, "rsq"])) }) @@ -1507,21 +1174,7 @@ test_that("twasWeightsCv: univariate fitter runs under verbose=2 (no quiet wrapp weightMethods = list(lassoWeights = list()), verbose = 2)) expect_true("prediction" %in% names(result)) - expect_equal(nrow(result$prediction[["lassoPredicted"]]), nrow(d$X)) -}) - -test_that("twasWeightsCv: parallel fold path (numThreads = 2)", { - d <- make_data(n = 30, p = 6) - local_mocked_bindings( - lassoWeights = function(X, y, ...) { w <- rep(0, ncol(X)); w[1] <- 0.4; w } - ) - set.seed(1) - result <- suppressMessages(twasWeightsCv( - d$X, d$Y, fold = 2, - weightMethods = list(lassoWeights = list()), - numThreads = 2)) - expect_true("prediction" %in% names(result)) - expect_equal(nrow(result$prediction[["lassoPredicted"]]), nrow(d$X)) + expect_equal(nrow(result$prediction[["lasso_predicted"]]), nrow(d$X)) }) # =========================================================================== diff --git a/tests/testthat/test_twasWeightsPipeline.R b/tests/testthat/test_twasWeightsPipeline.R index 04765e5f..60483db0 100644 --- a/tests/testthat/test_twasWeightsPipeline.R +++ b/tests/testthat/test_twasWeightsPipeline.R @@ -177,9 +177,16 @@ test_that("twasWeightsPipeline(QtlDataset): runs end-to-end with mocked solvers" expect_equal(nrow(res), 4L) expect_setequal(getMethodNames(res), c("lasso", "enet")) expect_setequal(getTraits(res), c("ENSG_A", "ENSG_B")) + # region provenance is populated from each trait's rowRanges and aligned per + # row (.tp_makeSe: ENSG_A -> chr1:1000-1499, ENSG_B -> chr1:2000-2499). + reg <- getRegion(res) + expect_s4_class(reg, "GRanges") + expect_equal(length(reg), 4L) + expect_true(all(GenomicRanges::start(reg)[res$trait == "ENSG_A"] == 1000L)) + expect_true(all(GenomicRanges::start(reg)[res$trait == "ENSG_B"] == 2000L)) }) -test_that("twasWeightsPipeline(QtlDataset): minTwasMaf/minTwasXvar tighten the variant set", { +test_that("twasWeightsPipeline(QtlDataset): mafCutoff/xvarCutoff overrides tighten the variant set", { qd <- .tp_makeQtlDataset(contexts = "brain", traits = "ENSG_A") do.call(local_mocked_bindings, c(list(extractBlockGenotypes = .tp_mockExtractor()), @@ -193,12 +200,13 @@ test_that("twasWeightsPipeline(QtlDataset): minTwasMaf/minTwasXvar tighten the v # Baseline: all 20 mock variants are retained. expect_equal(vcount(call()), 20L) # Mock MAF spans 0.225..0.388; a 0.25 cutoff drops the 3 lowest-MAF variants - # while leaving the rest, so the TWAS variant set strictly shrinks. - tight <- vcount(call(minTwasMaf = 0.25)) + # while leaving the rest, so the variant set strictly shrinks. The filter is + # the same one fine-mapping uses (unified QtlDataset QC), not TWAS-specific. + tight <- vcount(call(mafCutoff = 0.25)) expect_lt(tight, 20L) expect_gt(tight, 0L) - # minTwasXvar is wired through the same elevation; a no-op value is accepted. - expect_s4_class(call(minTwasXvar = 0), "TwasWeights") + # xvarCutoff overrides the per-variant variance slot; a no-op value is accepted. + expect_s4_class(call(xvarCutoff = 0), "TwasWeights") }) test_that("twasWeightsPipeline(QtlDataset): contexts filter restricts the per-context loop", { @@ -1982,6 +1990,16 @@ test_that(".rbindTwasWeights: concatenates two collections and rejects non-TwasW "expects two TwasWeights") }) +test_that("combineTwasWeights: variadic and list forms row-bind; guards empty/mixed", { + tw1 <- .tp_tw(trait = "g1"); tw2 <- .tp_tw(trait = "g2"); tw3 <- .tp_tw(trait = "g3") + expect_equal(nrow(combineTwasWeights(tw1, tw2)), 2L) # variadic + expect_equal(nrow(combineTwasWeights(list(tw1, tw2, tw3))), 3L) # single list + expect_equal(nrow(combineTwasWeights(tw1)), 1L) # single -> as-is + expect_s4_class(combineTwasWeights(tw1, tw2), "TwasWeights") + expect_error(combineTwasWeights(), "nothing to combine") + expect_error(combineTwasWeights(tw1, "nope"), "must be a TwasWeights") +}) + # ----------------------------------------------------------------------------- # .twasMergeRegions / .twasMergeRegionEntries / .twasFitsForRegion (dead-ish # multi-region helpers and per-region fit selection) diff --git a/tests/testthat/test_variantId.R b/tests/testthat/test_variantId.R index 5551aaf7..b0a84231 100644 --- a/tests/testthat/test_variantId.R +++ b/tests/testthat/test_variantId.R @@ -839,6 +839,17 @@ test_that("asGranges errors on malformed input", { expect_error(pecotmr:::asGranges(list(1)), "character vector or data.frame") }) +test_that("asGranges (exported) converts region strings and data.frames", { + gr <- asGranges("chr1:100-200") + expect_s4_class(gr, "GRanges") + expect_equal(as.character(GenomicRanges::seqnames(gr)), "chr1") + expect_equal(GenomicRanges::start(gr), 100L) + expect_equal(GenomicRanges::end(gr), 200L) + # numeric chrom gets a chr prefix; data.frame input is accepted too + gr2 <- asGranges(data.frame(chrom = "2", start = 5, end = 9)) + expect_equal(as.character(GenomicRanges::seqnames(gr2)), "chr2") +}) + test_that("matchVariants allowFlip = FALSE accepts data.frame input", { a <- .vid_df("1", 100, "A", "G") b <- .vid_df("1", 100, "A", "G")