diff --git a/NAMESPACE b/NAMESPACE index 736321cc..88b6b056 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ S3method(postprocessFinemappingFit,susie) S3method(postprocessFinemappingFit,susieInf) S3method(postprocessFinemappingFit,susieRss) export(.fullFitColumns) +export(.overlapPrefixNonKey) export(AnnotationMatrix) export(CtwasResult) export(CtwasResultEntry) @@ -157,10 +158,10 @@ export(getRefVariantInfo) export(getRegion) export(getResidualizedGenotypes) export(getResidualizedPhenotypes) -export(getSE) export(getSampleIds) export(getScaleResiduals) export(getScoreStats) +export(getSe) export(getSignificantQtls) export(getSnpIdx) export(getSnpInfo) @@ -206,6 +207,7 @@ export(ldPruneByCorrelation) export(learnTwasWeights) export(loadGenotypeRegion) export(loadGwasSumStatsFromManifest) +export(loadLdBlock) export(loadLdMatrix) export(loadLdSketch) export(loadMultiStudyQtlDatasetFromManifest) @@ -253,7 +255,6 @@ export(mvsusieWeights) export(nSignificantScore) export(nSnps) export(normalizeVariantId) -export(overlapTopLoci) export(parseCsCorr) export(parseRegion) export(parseVariantId) @@ -389,10 +390,10 @@ exportMethods(getRefPanel) exportMethods(getRegion) exportMethods(getResidualizedGenotypes) exportMethods(getResidualizedPhenotypes) -exportMethods(getSE) exportMethods(getSampleIds) exportMethods(getScaleResiduals) exportMethods(getScoreStats) +exportMethods(getSe) exportMethods(getSignificantQtls) exportMethods(getSnpIdx) exportMethods(getSnpInfo) diff --git a/R/AllClasses.R b/R/AllClasses.R index bec58385..6b8dc737 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -101,7 +101,7 @@ setMethod("getZ", "SumStatsBase", function(x, ...) mcols(getSumStats(x, ...))$Z) #' @export setMethod("getN", "SumStatsBase", function(x, ...) mcols(getSumStats(x, ...))$N) -# getP / getBeta / getSE are first-class alongside getZ / getN: they read the +# getP / getBeta / getSe are first-class alongside getZ / getN: they read the # optional P / BETA / SE mcols and return NULL when the entry does not carry # them (DataFrame `$` semantics), so a p-value-primary sumstats (e.g. TensorQTL # cis output) is an equal citizen to a Z-primary GWAS sumstats. @@ -113,9 +113,9 @@ setMethod("getP", "SumStatsBase", function(x, ...) mcols(getSumStats(x, ...))$P) #' @export setMethod("getBeta", "SumStatsBase", function(x, ...) mcols(getSumStats(x, ...))$BETA) -#' @rdname getSE +#' @rdname getSe #' @export -setMethod("getSE", "SumStatsBase", function(x, ...) mcols(getSumStats(x, ...))$SE) +setMethod("getSe", "SumStatsBase", function(x, ...) mcols(getSumStats(x, ...))$SE) #' @rdname getMaf #' @export diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 722243da..bdd1c3fa 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -137,7 +137,7 @@ setGeneric("getN", function(x, ...) standardGeneric("getN")) #' @description Extract the association p-value vector from a #' \code{GwasSumStats} or \code{QtlSumStats} entry, selected by its identity #' tuple. Part of the first-class summary-statistic column set alongside -#' \code{\link{getZ}} / \code{\link{getBeta}} / \code{\link{getSE}}. +#' \code{\link{getZ}} / \code{\link{getBeta}} / \code{\link{getSe}}. #' @param x A \code{GwasSumStats} or \code{QtlSumStats} object. #' @param ... Class-specific selection arguments. #' @return Numeric vector of p-values, or \code{NULL} if not available. @@ -162,7 +162,7 @@ setGeneric("getBeta", function(x, ...) standardGeneric("getBeta")) #' @param ... Class-specific selection arguments. #' @return Numeric vector of standard errors, or \code{NULL} if not available. #' @export -setGeneric("getSE", function(x, ...) standardGeneric("getSE")) +setGeneric("getSe", function(x, ...) standardGeneric("getSe")) #' @title Get Minor Allele Frequencies #' @description Extract MAF vector from a GwasSumStats object. @@ -1000,12 +1000,13 @@ setGeneric("getH2", function(x) standardGeneric("getH2")) setGeneric("fitJointGroup", function(group, pipeline, token, args) standardGeneric("fitJointGroup")) -# construct(pipeline, rows) -- assemble the per-pipeline result collection -# (QtlFineMappingResult vs TwasWeights) from accumulated joint rows. The joint -# row identity (which axes collapse to "joint" + jointStudies/Contexts/Traits) -# is derived from each group's `conditions` by the rows accumulator. +# construct(pipeline, records) -- assemble the per-pipeline result collection +# (QtlFineMappingResult vs TwasWeights) from the driver's list of per-row +# records. The joint row identity (which axes collapse to "joint" + +# jointStudies/Contexts/Traits) is carried on each record, derived from each +# group's `conditions`. setGeneric("construct", - function(pipeline, rows, ...) standardGeneric("construct")) + function(pipeline, records, ...) standardGeneric("construct")) # ---- SldscData accessors ---- #' @title Get the annotation table from an SldscData diff --git a/R/QtlDataset.R b/R/QtlDataset.R index 36f50fbf..04befb0a 100644 --- a/R/QtlDataset.R +++ b/R/QtlDataset.R @@ -722,7 +722,7 @@ setMethod("getPhenotypeCovariates", "QtlDataset", # the redundant columns automatically. Optionally rescales each residual # column to unit standard deviation; constant-valued columns are left # unchanged. -.qtlResidualizeQR <- function(Y, C, scaleResiduals = TRUE) { +.qtlResidualizeQr <- function(Y, C, scaleResiduals = TRUE) { X <- if (is.null(C) || ncol(C) == 0L) { matrix(1, nrow = nrow(Y), ncol = 1L, dimnames = list(rownames(Y), "intercept")) @@ -751,6 +751,23 @@ setMethod("getPhenotypeCovariates", "QtlDataset", res } +# Resolve the phenotype-covariate selection for one context: NULL requested -> +# all available covariates; otherwise validate the requested names are present. +# @noRd +.qtlResolveOne <- function(ctx, requested, x) { + se <- x@phenotypes[[ctx]] + avail <- colnames(SummarizedExperiment::colData(se)) + if (is.null(requested)) return(avail) + keep <- intersect(requested, avail) + if (length(keep) != length(requested)) { + missingNames <- setdiff(requested, avail) + stop(sprintf( + "phenotypeCovariatesToResidualize: context '%s' has no covariate(s) named: %s", + ctx, paste(missingNames, collapse = ", "))) + } + keep +} + # Internal: validate and resolve the `*ToResidualize` argument against a # set of contexts and the covariates actually present in those contexts' # colData. Accepts either NULL (use all), a character vector (apply to all @@ -762,21 +779,8 @@ setMethod("getPhenotypeCovariates", "QtlDataset", # (per the rule: named-list keys must equal `contexts`) # - an explicitly requested name matches no actual covariate .qtlResolvePhenoSelection <- function(x, contexts, toResidualize) { - resolveOne <- function(ctx, requested) { - se <- x@phenotypes[[ctx]] - avail <- colnames(SummarizedExperiment::colData(se)) - if (is.null(requested)) return(avail) - keep <- intersect(requested, avail) - if (length(keep) != length(requested)) { - missingNames <- setdiff(requested, avail) - stop(sprintf( - "phenotypeCovariatesToResidualize: context '%s' has no covariate(s) named: %s", - ctx, paste(missingNames, collapse = ", "))) - } - keep - } if (is.null(toResidualize)) { - out <- lapply(contexts, resolveOne, requested = NULL) + out <- lapply(contexts, .qtlResolveOne, requested = NULL, x = x) names(out) <- contexts return(out) } @@ -798,12 +802,12 @@ setMethod("getPhenotypeCovariates", "QtlDataset", "context set as `contexts`. Missing keys: ", paste(missingKeys, collapse = ", ")) } - out <- lapply(contexts, function(ctx) resolveOne(ctx, toResidualize[[ctx]])) + out <- lapply(contexts, function(ctx) .qtlResolveOne(ctx, toResidualize[[ctx]], x)) names(out) <- contexts return(out) } if (is.character(toResidualize)) { - out <- lapply(contexts, resolveOne, requested = toResidualize) + out <- lapply(contexts, .qtlResolveOne, requested = toResidualize, x = x) names(out) <- contexts return(out) } @@ -1012,7 +1016,7 @@ setMethod("getResidualizedGenotypes", "QtlDataset", G <- G[common, , drop = FALSE] C <- C[common, , drop = FALSE] } - .qtlResidualizeQR(G, C, scaleResiduals = x@scaleResiduals) + .qtlResidualizeQr(G, C, scaleResiduals = x@scaleResiduals) }) #' @rdname getResidualizedPhenotypes @@ -1091,7 +1095,7 @@ setMethod("getResidualizedPhenotypes", "QtlDataset", } else { Cctx <- NULL } - Yres <- .qtlResidualizeQR(Y, Cctx, scaleResiduals = x@scaleResiduals) + Yres <- .qtlResidualizeQr(Y, Cctx, scaleResiduals = x@scaleResiduals) # Outlier detection on the residualized scale (samples whose # residualized phenotype is unusual *given* their covariates). if (outlierAction != "keep") { diff --git a/R/crossValidation.R b/R/crossValidation.R index 8d9c3320..8d4c3bb0 100644 --- a/R/crossValidation.R +++ b/R/crossValidation.R @@ -48,9 +48,42 @@ out } +# One CV fold: split train/test by fold `j`, drop zero-variance training columns, +# fit via `fitFold(Xtr, Ytr, j, fitFoldCtx)`, and predict the held-out samples. +# @noRd +.cvRunFold <- function(j, cv) { + X <- cv$X; Y <- cv$Y; samplePartition <- cv$samplePartition + foldIds <- cv$foldIds; fitFold <- cv$fitFold; fitFoldCtx <- cv$fitFoldCtx + retainFits <- cv$retainFits; verbose <- cv$verbose + 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, fitFoldCtx) + 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()) +} + +# No-op fold fitter: used when the caller only wants the fold partition. +# @noRd +.cvNoopFitFold <- function(Xtr, Ytr, j, fitFoldCtx) list(weights = list(), fits = list()) + # Shared K-fold cross-validation engine. # -# `fitFold(Xtrain, Ytrain, foldIndex)` must return +# `fitFold(Xtrain, Ytrain, foldIndex, fitFoldCtx)` must return # list(weights = (variants x outcomes) weight matrix, # rownames indexing colnames(Xtrain)>, # fits = fitted model or NULL>) @@ -64,9 +97,9 @@ #' @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) { + fitFold, fitFoldCtx = NULL, 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.") } @@ -136,36 +169,16 @@ 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())) + cvState <- list(X = X, Y = Y, samplePartition = samplePartition, + foldIds = foldIds, fitFold = fitFold, fitFoldCtx = fitFoldCtx, + retainFits = retainFits, verbose = verbose) foldResults <- if (numCores >= 2) { - bplapply(foldIds, runFold, + bplapply(foldIds, .cvRunFold, cv = cvState, BPPARAM = MulticoreParam(workers = numCores, RNGseed = 1L)) } else { - lapply(foldIds, runFold) + lapply(foldIds, .cvRunFold, cv = cvState) } metricNames <- c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE") diff --git a/R/ctwasPipeline.R b/R/ctwasPipeline.R index e411bf52..75ebc7be 100644 --- a/R/ctwasPipeline.R +++ b/R/ctwasPipeline.R @@ -936,6 +936,11 @@ mergeCtwasBoundaryRegions <- function(finemapResult, .ctwasBucketWeights(twasWeights, gwasSumStats) } +# Extract the character vector of method names carried by a weight source +# (NULL-safe: a NULL source contributes no methods). +# @noRd +.ctwasMethodsOf <- function(tw) if (is.null(tw)) NULL else as.character(tw$method) + # 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). @@ -946,12 +951,11 @@ mergeCtwasBoundaryRegions <- function(finemapResult, # `.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 + .ctwasMethodsOf(twasWeightsList) # a flat weight source + else unlist(lapply(twasWeightsList, .ctwasMethodsOf))) # a list of them if (length(available) == 0L) stop("ctwasPipeline: weight sources carry no method entries.") if (!is.null(method) && nzchar(method)) { @@ -979,6 +983,11 @@ mergeCtwasBoundaryRegions <- function(finemapResult, studies } +# Extract field `i` (as character) from each `region|study|context|trait|method` +# split in `parts`. +# @noRd +.ctwasPickField <- function(i, parts) vapply(parts, function(p) p[[i]], character(1)) + # 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 @@ -991,12 +1000,11 @@ mergeCtwasBoundaryRegions <- function(finemapResult, 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), + rid = .ctwasPickField(1L, parts), + study = .ctwasPickField(2L, parts), + context = .ctwasPickField(3L, parts), 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), @@ -1057,6 +1065,13 @@ mergeCtwasBoundaryRegions <- function(finemapResult, if (nrow(sub) == 0L) NULL else `rownames<-`(sub, NULL) } +# Build a CtwasResultEntry from a finemap + susieAlpha slice, stamping the run's +# param + region_info. +# @noRd +.ctwasMkEntry <- function(fm, sa, runResult) CtwasResultEntry( + finemap = fm, susieAlpha = sa, param = runResult$param, + regionInfo = runResult$region_info) + # 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 @@ -1082,10 +1097,6 @@ mergeCtwasBoundaryRegions <- function(finemapResult, 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]) @@ -1095,8 +1106,8 @@ mergeCtwasBoundaryRegions <- function(finemapResult, idsCx <- parsed$id[inCx] list(gwasStudy = gwasStudy, study = studyCx, context = cx, method = method, jointContexts = jointStr, - entry = mkEntry(.ctwasSubsetById(fmDf, idsCx), - .ctwasSubsetById(saDf, idsCx))) + entry = .ctwasMkEntry(.ctwasSubsetById(fmDf, idsCx), + .ctwasSubsetById(saDf, idsCx), runResult)) }) if (keepSnps) { @@ -1105,7 +1116,7 @@ mergeCtwasBoundaryRegions <- function(finemapResult, 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)))) + jointContexts = jointStr, entry = .ctwasMkEntry(snpFm, snpSa, runResult)))) } rows } diff --git a/R/fineMappingPipeline.R b/R/fineMappingPipeline.R index 2c6c7f25..471c3aef 100644 --- a/R/fineMappingPipeline.R +++ b/R/fineMappingPipeline.R @@ -503,6 +503,14 @@ setGeneric("fineMappingPipeline", } } +# TRUE if method token `tk` is unknown (kept; validated elsewhere) or its +# capability advertises a non-NULL `capField`. +# @noRd +.fmMethodOk <- function(tk, capField, caps) { + info <- caps[[tk]] + is.null(info) || !is.null(info[[capField]]) +} + # Keep only the tokens in `methods` whose capability has a non-NULL `capField` # (individualImpl / sumstatImpl), so a sumstat-only method (e.g. ser) is dropped # from the individual-level recursion and an individual-only method from the @@ -510,9 +518,8 @@ setGeneric("fineMappingPipeline", # of per-token args; unknown tokens pass through (handled elsewhere). .fmFilterMethodsForKind <- function(methods, capField) { caps <- .fineMappingMethodCapabilities - ok <- function(tk) { info <- caps[[tk]]; is.null(info) || !is.null(info[[capField]]) } - if (is.character(methods)) methods[vapply(methods, ok, logical(1))] - else if (is.list(methods)) methods[vapply(names(methods), ok, logical(1))] + if (is.character(methods)) methods[vapply(methods, .fmMethodOk, logical(1), capField, caps)] + else if (is.list(methods)) methods[vapply(names(methods), .fmMethodOk, logical(1), capField, caps)] else methods } @@ -579,6 +586,26 @@ setGeneric("fineMappingPipeline", } +# Append one QTL-side result row to the accumulator env `acc` (holds the parallel +# rowStudy/rowContext/rowTrait/rowMethod vectors + rowEntries list). +# @noRd +.fmPushQtlRow <- function(acc, st, ctx, tr, mt, ent) { + acc$rowStudy <- c(acc$rowStudy, st) + acc$rowContext <- c(acc$rowContext, ctx) + acc$rowTrait <- c(acc$rowTrait, tr) + acc$rowMethod <- c(acc$rowMethod, mt) + acc$rowEntries[[length(acc$rowEntries) + 1L]] <- ent +} + +# Append one GWAS-side result row (region-keyed) to the accumulator env `acc`. +# @noRd +.fmPushGwasRow <- function(acc, st, mt, rg, ent) { + acc$rowStudy <- c(acc$rowStudy, st) + acc$rowMethod <- c(acc$rowMethod, mt) + acc$rowRegion <- c(acc$rowRegion, rg) + acc$rowEntries[[length(acc$rowEntries) + 1L]] <- ent +} + # Build a QtlFineMappingResult collection from per-tuple parallel vectors. # `jointStudies`, `jointContexts`, `jointTraits` are optional character # vectors (length matches `studies`) describing semicolon-joined joint @@ -824,6 +851,14 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { } } +# The canonical (non-reweighted) mvSuSiE mixture prior for residual variance `V`: +# create_mixture_prior(R) restricted to the group's conditions. +# @noRd +.fmCanonicalPrior <- function(V, conditionNames, R) list( + priorVariance = mvsusieR::create_mixture_prior( + R = R, include_indices = conditionNames), + residualVariance = V) + # Rebuild the mvSuSiE data-driven *reweighted* mixture prior + residual variance # from a stored mr.mash fit -- the lean payload # (list(dataDrivenPriorMatrices, w0, V)) that mrmashWeights(retainFit = TRUE) @@ -843,20 +878,16 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { .buildMvsusieReweightedPrior <- function(fitParts, conditionNames, weightsTol = 1e-10, overrideU = NULL) { R <- length(conditionNames) - canonical <- function(V) list( - priorVariance = mvsusieR::create_mixture_prior( - R = R, include_indices = conditionNames), - residualVariance = V) - if (is.null(fitParts)) return(canonical(NULL)) + if (is.null(fitParts)) return(.fmCanonicalPrior(NULL, conditionNames, R)) # `overrideU` (mode C / hybrid): reuse this fit's reweighted mixture weights # (w0) and residual variance (V) but swap in a different set of data-driven # covariance matrices -- the per-fold mash prior U. Components are matched to # w0 by name, so the override U must share component names with the fit. ddpm <- if (!is.null(overrideU)) overrideU else fitParts$dataDrivenPriorMatrices - if (is.null(ddpm) || is.null(ddpm$U)) return(canonical(fitParts$V)) + if (is.null(ddpm) || is.null(ddpm$U)) return(.fmCanonicalPrior(fitParts$V, conditionNames, R)) w0Updated <- rescaleCovW0(fitParts$w0) w0Updated <- w0Updated[names(w0Updated) %in% names(ddpm$U)] - if (length(w0Updated) == 0L) return(canonical(fitParts$V)) + if (length(w0Updated) == 0L) return(.fmCanonicalPrior(fitParts$V, conditionNames, R)) mixture <- list(matrices = ddpm$U[names(w0Updated)], weights = w0Updated) list( priorVariance = mvsusieR::create_mixture_prior( @@ -1118,7 +1149,7 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { # GwasSumStats entry GRanges. Errors when Z or N is missing. Wraps the # shared `.entryToSumstatDf` helper (R/sumstatsQc.R). # @noRd -.fmExtractZN <- function(gr, label) { +.fmExtractZn <- function(gr, label) { df <- .entryToSumstatDf(gr, require = c("SNP", "Z", "N"), label = label) @@ -1172,6 +1203,14 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { sub("_weights$", "", adapter$methodKey) } +# Coerce a weight vector to a single-column matrix (rows named by the vector's +# names); pass matrices through unchanged. +# @noRd +.fmAsMat <- function(w) { + if (is.matrix(w)) return(w) + matrix(w, ncol = 1L, dimnames = list(names(w), NULL)) +} + # 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 @@ -1179,10 +1218,6 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { # @noRd .fmFoldWeights <- function(token, Xtr, Ytr, coverage, userArgs, pos, mvPrior = NULL) { - asMat <- function(w) { - if (is.matrix(w)) return(w) - matrix(w, ncol = 1L, dimnames = list(names(w), NULL)) - } if (token %in% c("susie", "susieInf", "susieAsh")) { y <- if (is.matrix(Ytr)) Ytr[, 1L] else Ytr fit <- .fmFitSusieIndiv(Xtr, y, token, coverage = coverage, @@ -1193,7 +1228,7 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { susieAsh = susieAshWeights(susieAshFit = fit)) w <- as.numeric(w) names(w) <- colnames(Xtr) - return(asMat(w)) + return(.fmAsMat(w)) } if (token == "mvsusie") { # Reuse the data-driven reweighted prior + residual covariance from the @@ -1222,6 +1257,34 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { NULL } +# Per-fold fine-mapping fit for the CV engine. `ctx` carries mvPrior, mvPriorCv, +# tokens, coverage, methodArgs, pos, verbose. Weights keyed by canonical method key. +# @noRd +.fmFitFold <- function(Xtr, Ytr, j, ctx) { + mvPrior <- ctx$mvPrior; mvPriorCv <- ctx$mvPriorCv; tokens <- ctx$tokens + coverage <- ctx$coverage; methodArgs <- ctx$methodArgs; pos <- ctx$pos + verbose <- ctx$verbose + # 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 + weights <- list() + for (tk in tokens) { + weights[[.fmTwasMethodKey(tk)]] <- tryCatch( + .fmFoldWeights(tk, Xtr, Ytr, coverage, methodArgs[[tk]], pos, + mvPriorThisFold), + error = function(e) { + if (verbose >= 1) + message(sprintf(" CV fold %s, method %s failed: %s", + j, tk, conditionMessage(e))) + NULL + }) + } + list(weights = weights, fits = list()) +} + # 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 @@ -1235,33 +1298,16 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { pos = NULL, verbose = 1, mvPrior = NULL, mvPriorCv = NULL, numThreads = 1) { if (length(tokens) == 0L) return(NULL) - # 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 - weights <- list() - for (tk in tokens) { - weights[[.fmTwasMethodKey(tk)]] <- tryCatch( - .fmFoldWeights(tk, Xtr, Ytr, coverage, methodArgs[[tk]], pos, - mvPriorThisFold), - error = function(e) { - if (verbose >= 1) - message(sprintf(" CV fold %s, method %s failed: %s", - j, tk, conditionMessage(e))) - NULL - }) - } - list(weights = weights, fits = list()) - } + # Per-fold fit context passed to the shared engine's top-level fitter + # (.fmFitFold). Weights are keyed by the canonical method key so + # _predicted / _performance line up with the TwasWeights method column. + cvFitCtx <- list(mvPrior = mvPrior, mvPriorCv = mvPriorCv, tokens = tokens, + coverage = coverage, methodArgs = methodArgs, pos = pos, + verbose = verbose) res <- .crossValidateWeights( X, Y, fold = fold, samplePartitions = samplePartition, - fitFold = fitFold, numThreads = numThreads, verbose = verbose) + fitFold = .fmFitFold, fitFoldCtx = cvFitCtx, + numThreads = numThreads, verbose = verbose) list(samplePartition = res$samplePartition, prediction = res$prediction, performance = res$performance) } @@ -1454,19 +1500,12 @@ setMethod("fineMappingPipeline", "QtlDataset", chain <- .fmResolveSusieChain(univTokens, addSusieInf) - rowStudy <- character(0) - rowContext <- character(0) - rowTrait <- character(0) - rowMethod <- character(0) - rowEntries <- list() - - pushRow <- function(st, ctx, tr, mt, ent) { - rowStudy <<- c(rowStudy, st) - rowContext <<- c(rowContext, ctx) - rowTrait <<- c(rowTrait, tr) - rowMethod <<- c(rowMethod, mt) - rowEntries[[length(rowEntries) + 1L]] <<- ent - } + acc <- new.env(parent = emptyenv()) + acc$rowStudy <- character(0) + acc$rowContext <- character(0) + acc$rowTrait <- character(0) + acc$rowMethod <- character(0) + acc$rowEntries <- list() # ---- Univariate dispatch: per (context, trait), per method. # X is drawn from each window in `xRegions` (cis = one trait-derived block; @@ -1481,7 +1520,7 @@ setMethod("fineMappingPipeline", "QtlDataset", for (tk in univTokens) { cached <- .fmCacheLookup(fineMappingResult, study, ctx, tid, tk) if (!is.null(cached)) { - pushRow(study, ctx, tid, tk, cached) + .fmPushQtlRow(acc, study, ctx, tid, tk, cached) } else { toRun <- c(toRun, tk) } @@ -1531,7 +1570,7 @@ setMethod("fineMappingPipeline", "QtlDataset", ents <- lapply(blockEntries, function(be) be[[tk]]) if (any(vapply(ents, is.null, logical(1)))) next entry <- if (length(ents) == 1L) ents[[1L]] else .fmMergeEntries(ents) - pushRow(study, ctx, tid, tk, entry) + .fmPushQtlRow(acc, study, ctx, tid, tk, entry) } } } @@ -1555,7 +1594,7 @@ setMethod("fineMappingPipeline", "QtlDataset", ncol(scores), ctx, length(traits))) for (pcName in colnames(scores)) { cached <- .fmCacheLookup(fineMappingResult, study, ctx, pcName, "susie") - if (!is.null(cached)) { pushRow(study, ctx, pcName, "susie", cached); next } + if (!is.null(cached)) { .fmPushQtlRow(acc, study, ctx, pcName, "susie", cached); next } pcY <- scores[, pcName] blockEntries <- lapply(xRegions, function(rg) { X <- if (is.null(rg)) { @@ -1581,7 +1620,7 @@ setMethod("fineMappingPipeline", "QtlDataset", ents <- lapply(blockEntries, function(be) be[["susie"]]) if (any(vapply(ents, is.null, logical(1)))) next entry <- if (length(ents) == 1L) ents[[1L]] else .fmMergeEntries(ents) - pushRow(study, ctx, pcName, "susie", entry) + .fmPushQtlRow(acc, study, ctx, pcName, "susie", entry) } } } @@ -1614,6 +1653,9 @@ setMethod("fineMappingPipeline", "QtlDataset", ldSketch = NULL) } + rowStudy <- acc$rowStudy; rowContext <- acc$rowContext + rowTrait <- acc$rowTrait; rowMethod <- acc$rowMethod + rowEntries <- acc$rowEntries perTupleResult <- if (length(rowEntries) > 0L) .fmBuildQtlResult(rowStudy, rowContext, rowTrait, rowMethod, rowEntries, region = tryCatch( @@ -1641,6 +1683,38 @@ setMethod("fineMappingPipeline", "QtlDataset", # MultiStudyQtlDataset method # ============================================================================= +# Per-embedded-study fine-mapping worker for .multiStudyPipelineDriver: recurse +# fineMappingPipeline on one QtlDataset with the individual-capable methods. +# `cfg` bundles the parent call's forwarded arguments. +# @noRd +.fmPerStudy <- function(qd, cfg) { + m <- .fmFilterMethodsForKind(cfg$methods, "individualImpl") + if (length(m) == 0L) return(NULL) + do.call(fineMappingPipeline, c(list( + data = qd, methods = m, contexts = cfg$contexts, traitId = cfg$traitId, + region = cfg$region, cisWindow = cfg$cisWindow, jointRegions = cfg$jointRegions, + jointSpecification = NULL, addSusieInf = cfg$addSusieInf, coverage = cfg$coverage, + secondaryCoverage = cfg$secondaryCoverage, signalCutoff = cfg$signalCutoff, + minAbsCorr = cfg$minAbsCorr, fineMappingResult = cfg$fineMappingResult, + cvFolds = cfg$cvFolds, cvThreads = cfg$cvThreads, + samplePartition = cfg$samplePartition, pipCutoffToSkip = cfg$pipCutoffToSkip, + seed = cfg$seed, naAction = cfg$naAction, verbose = cfg$verbose), cfg$dotArgs)) +} + +# Embedded-sumstats fine-mapping worker for .multiStudyPipelineDriver: recurse +# fineMappingPipeline on the QtlSumStats with the sumstat-capable methods. +# @noRd +.fmSumStats <- function(ss, cfg) { + m <- .fmFilterMethodsForKind(cfg$methods, "sumstatImpl") + if (length(m) == 0L) return(NULL) + do.call(fineMappingPipeline, c(list( + data = ss, methods = m, contexts = cfg$contexts, traitId = cfg$traitId, + jointSpecification = NULL, addSusieInf = cfg$addSusieInf, coverage = cfg$coverage, + secondaryCoverage = cfg$secondaryCoverage, signalCutoff = cfg$signalCutoff, + minAbsCorr = cfg$minAbsCorr, fineMappingResult = cfg$fineMappingResult, + verbose = cfg$verbose), cfg$dotArgs)) +} + #' @rdname fineMappingPipeline #' @export setMethod("fineMappingPipeline", "MultiStudyQtlDataset", @@ -1719,31 +1793,16 @@ setMethod("fineMappingPipeline", "MultiStudyQtlDataset", # Route each method to the components it supports: individual-capable # methods to the per-study QtlDatasets, sumstat-capable methods (incl. the # sumstat-only `ser`) to the embedded QtlSumStats. - perStudyFn <- function(qd) { - m <- .fmFilterMethodsForKind(methods, "individualImpl") - if (length(m) == 0L) return(NULL) - do.call(fineMappingPipeline, c(list( - data = qd, methods = m, contexts = contexts, traitId = traitId, - region = region, cisWindow = cisWindow, jointRegions = jointRegions, - jointSpecification = NULL, addSusieInf = addSusieInf, coverage = coverage, - secondaryCoverage = secondaryCoverage, signalCutoff = signalCutoff, - minAbsCorr = minAbsCorr, fineMappingResult = fineMappingResult, - cvFolds = cvFolds, cvThreads = cvThreads, samplePartition = samplePartition, - pipCutoffToSkip = pipCutoffToSkip, seed = seed, naAction = naAction, - verbose = verbose), dotArgs)) - } - sumStatsFn <- function(ss) { - m <- .fmFilterMethodsForKind(methods, "sumstatImpl") - if (length(m) == 0L) return(NULL) - do.call(fineMappingPipeline, c(list( - data = ss, methods = m, contexts = contexts, traitId = traitId, - jointSpecification = NULL, addSusieInf = addSusieInf, coverage = coverage, - secondaryCoverage = secondaryCoverage, signalCutoff = signalCutoff, - minAbsCorr = minAbsCorr, fineMappingResult = fineMappingResult, - verbose = verbose), dotArgs)) - } + cfg <- list(methods = methods, contexts = contexts, traitId = traitId, + region = region, cisWindow = cisWindow, jointRegions = jointRegions, + addSusieInf = addSusieInf, coverage = coverage, + secondaryCoverage = secondaryCoverage, signalCutoff = signalCutoff, + minAbsCorr = minAbsCorr, fineMappingResult = fineMappingResult, + cvFolds = cvFolds, cvThreads = cvThreads, + samplePartition = samplePartition, pipCutoffToSkip = pipCutoffToSkip, + seed = seed, naAction = naAction, verbose = verbose, dotArgs = dotArgs) .multiStudyPipelineDriver( - data, jointResult, perStudyFn, sumStatsFn, + data, jointResult, .fmPerStudy, .fmSumStats, cfg, .rbindFineMappingResult, QtlFineMappingResult, "fineMappingPipeline") }) @@ -1846,19 +1905,13 @@ setMethod("fineMappingPipeline", "QtlSumStats", (isTRUE(serFallback) || !identical(rMismatch, "none"))) getNSamples(ldSketch) else rFinite - rowStudy <- character(0) - rowContext <- character(0) - rowTrait <- character(0) - rowMethod <- character(0) - rowEntries <- list() + acc <- new.env(parent = emptyenv()) + acc$rowStudy <- character(0) + acc$rowContext <- character(0) + acc$rowTrait <- character(0) + acc$rowMethod <- character(0) + acc$rowEntries <- list() nSkipped <- 0L - pushRow <- function(st, ctx, tr, mt, ent) { - rowStudy <<- c(rowStudy, st) - rowContext <<- c(rowContext, ctx) - rowTrait <<- c(rowTrait, tr) - rowMethod <<- c(rowMethod, mt) - rowEntries[[length(rowEntries) + 1L]] <<- ent - } # ---- Univariate dispatch: per (study, context, trait), per method. if (length(univTokens) > 0L) { @@ -1870,7 +1923,7 @@ setMethod("fineMappingPipeline", "QtlSumStats", for (tk in univTokens) { cached <- .fmCacheLookup(fineMappingResult, st, ctx, tr, tk) if (!is.null(cached)) { - pushRow(st, ctx, tr, tk, cached) + .fmPushQtlRow(acc, st, ctx, tr, tk, cached) } else { toRun <- c(toRun, tk) } @@ -1878,7 +1931,7 @@ setMethod("fineMappingPipeline", "QtlSumStats", if (length(toRun) == 0L) next # A trait screened out by summaryStatsQc(pipCutoffToSkip) is empty here; - # skip it gracefully (no row, no error) rather than tripping .fmExtractZN. + # skip it gracefully (no row, no error) rather than tripping .fmExtractZn. skip <- .fmEntrySkipInfo(data, i) if (isTRUE(skip$skipped)) { nSkipped <- nSkipped + 1L @@ -1889,7 +1942,7 @@ setMethod("fineMappingPipeline", "QtlSumStats", next } entry <- data$entry[[i]] - zn <- .fmExtractZN(entry, + zn <- .fmExtractZn(entry, sprintf("fineMappingPipeline(QtlSumStats): entry %d (study='%s', context='%s', trait='%s')", i, st, ctx, tr)) variantIds <- zn$variantIds z <- zn$z @@ -1913,7 +1966,7 @@ setMethod("fineMappingPipeline", "QtlSumStats", keepFullFit = keepFullFit) # The method column carries the bare token, independent of the # postprocess class. - for (tk in names(ents)) pushRow(st, ctx, tr, tk, ents[[tk]]) + for (tk in names(ents)) .fmPushQtlRow(acc, st, ctx, tr, tk, ents[[tk]]) } } @@ -1941,6 +1994,9 @@ setMethod("fineMappingPipeline", "QtlSumStats", ldSketch = ldSketch) } + rowStudy <- acc$rowStudy; rowContext <- acc$rowContext + rowTrait <- acc$rowTrait; rowMethod <- acc$rowMethod + rowEntries <- acc$rowEntries perTupleResult <- if (length(rowEntries) > 0L) .fmBuildQtlResult(rowStudy, rowContext, rowTrait, rowMethod, rowEntries, # QtlSumStats: region = the entry's variant span (no @@ -2021,23 +2077,18 @@ setMethod("fineMappingPipeline", "GwasSumStats", getNSamples(ldSketch) else rFinite studyCol <- as.character(data$study) - rowStudy <- character(0) - rowMethod <- character(0) - rowRegion <- character(0) - rowEntries <- list() + acc <- new.env(parent = emptyenv()) + acc$rowStudy <- character(0) + acc$rowMethod <- character(0) + acc$rowRegion <- character(0) + acc$rowEntries <- list() nSkipped <- 0L - pushRow <- function(st, mt, rg, ent) { - rowStudy <<- c(rowStudy, st) - rowMethod <<- c(rowMethod, mt) - rowRegion <<- c(rowRegion, rg) - rowEntries[[length(rowEntries) + 1L]] <<- ent - } for (i in seq_len(nrow(data))) { st <- studyCol[[i]] gr <- data$entry[[i]] # A region screened out by summaryStatsQc(pipCutoffToSkip) is empty here; - # skip it gracefully (no row, no error) rather than tripping .fmExtractZN. + # skip it gracefully (no row, no error) rather than tripping .fmExtractZn. skip <- .fmEntrySkipInfo(data, i) if (isTRUE(skip$skipped)) { nSkipped <- nSkipped + 1L @@ -2047,7 +2098,7 @@ setMethod("fineMappingPipeline", "GwasSumStats", st, skip$reason)) next } - zn <- .fmExtractZN(gr, + zn <- .fmExtractZn(gr, sprintf("fineMappingPipeline(GwasSumStats): study='%s'", st)) variantIds <- zn$variantIds z <- zn$z @@ -2077,7 +2128,7 @@ setMethod("fineMappingPipeline", "GwasSumStats", .fmCacheLookupGwas(fineMappingResult, st, tk, region_id) } else NULL if (!is.null(cached)) { - pushRow(st, tk, region_id, cached) + .fmPushGwasRow(acc, st, tk, region_id, cached) } else { toRun <- c(toRun, tk) } @@ -2095,11 +2146,13 @@ setMethod("fineMappingPipeline", "GwasSumStats", serFallback = serFallback, rFinite = rFiniteResolved, rMismatch = rMismatch, rssControl = rssControl, keepFullFit = keepFullFit) - for (tk in names(ents)) pushRow(st, tk, region_id, ents[[tk]]) + for (tk in names(ents)) .fmPushGwasRow(acc, st, tk, region_id, ents[[tk]]) } # An all-screened (or empty-input) collection legitimately yields a 0-row # result -- allow it instead of erroring "no ... tuples produced a result". + rowStudy <- acc$rowStudy; rowMethod <- acc$rowMethod + rowRegion <- acc$rowRegion; rowEntries <- acc$rowEntries .fmBuildGwasResult(rowStudy, rowMethod, rowEntries, region_ids = rowRegion, ldSketch = ldSketch, diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index c663fe79..6827379e 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -637,6 +637,47 @@ computeCsTable <- function(fit, dataX, coverage, csInput = c("X", "Xcorr", "fsus cols } +# Slice a susie posterior array to the active condition (3-D fit) or coerce a +# 2-D array to matrix; NULL for a 3-D fit with no conditionIdx. +# @noRd +.fmSliceCond <- function(arr, conditionIdx) { + if (is.null(arr)) return(NULL) + if (length(dim(arr)) == 3L) { + if (is.null(conditionIdx)) return(NULL) + return(as.matrix(arr[, , conditionIdx])) + } + as.matrix(arr) +} + +# Per-variant CS index at coverage `targetCov` (0 = not in any CS; on overlap the +# smallest cs_idx wins). +# @noRd +.fmCsIdxAtCoverage <- function(targetCov, coverageValues, csTables, nV) { + out <- integer(nV) + hit <- which(abs(coverageValues - targetCov) < 1e-12) + if (length(hit) == 0L) return(out) + sets <- csTables[[hit[1L]]]$sets$cs + if (is.null(sets) || length(sets) == 0L) return(out) + for (csIdx in seq_along(sets)) { + vi <- as.integer(sets[[csIdx]]) + vi <- vi[vi >= 1L & vi <= nV & out[vi] == 0L] + out[vi] <- csIdx + } + out +} + +# Per-variant CS purity (min.abs.corr) at coverage `targetCov`; 0 for non-CS +# variants. +# @noRd +.fmPurityAtCoverage <- function(targetCov, idxVec, coverageValues, csTables) { + h <- which(abs(coverageValues - targetCov) < 1e-12) + pv <- if (length(h) > 0L) .csPurityVec(csTables[[h[1L]]]) else numeric() + vapply(idxVec, function(i) { + if (i <= 0L || i > length(pv)) return(0) + v <- pv[i]; if (is.na(v)) 0 else as.numeric(v) + }, numeric(1)) +} + buildTopLoci <- function(fit, csTables, variantNames, sumstats = NULL, af = NULL, method, signalCutoff = 0, dataX = NULL, dataY = NULL, @@ -670,16 +711,8 @@ buildTopLoci <- function(fit, csTables, variantNames, sumstats = NULL, # (the per-context-row representation). conditionIdx = NULL keeps the 2-D path # (univariate); a 3-D fit without a conditionIdx leaves posterior NA. alpha <- as.matrix(fit$alpha) - sliceCond <- function(arr) { - if (is.null(arr)) return(NULL) - if (length(dim(arr)) == 3L) { - if (is.null(conditionIdx)) return(NULL) - return(as.matrix(arr[, , conditionIdx])) - } - as.matrix(arr) - } - mu <- sliceCond(fit$mu) - mu2 <- sliceCond(fit$mu2) + mu <- .fmSliceCond(fit$mu, conditionIdx) + mu2 <- .fmSliceCond(fit$mu2, conditionIdx) postMean <- if (!is.null(mu) && all(dim(alpha) == dim(mu))) { colSums(alpha * mu) } else rep(NA_real_, length(variantNames)) @@ -732,31 +765,10 @@ buildTopLoci <- function(fit, csTables, variantNames, sumstats = NULL, # Per-coverage CS membership: for each variant, which CS at each # coverage level (cs_idx, or 0 if not in any). If a variant belongs # to multiple CSs at a given coverage, the smallest cs_idx wins. - csIdxAtCoverage <- function(targetCov) { - out <- integer(nV) - hit <- which(abs(coverageValues - targetCov) < 1e-12) - if (length(hit) == 0L) return(out) - sets <- csTables[[hit[1L]]]$sets$cs - if (is.null(sets) || length(sets) == 0L) return(out) - for (csIdx in seq_along(sets)) { - vi <- as.integer(sets[[csIdx]]) - vi <- vi[vi >= 1L & vi <= nV & out[vi] == 0L] - out[vi] <- csIdx - } - out - } # Per-variant CS purity at a coverage (0 for non-CS variants). Purity # (min.abs.corr) is a CS-quality measure, independent of the coverage # (confidence) level; the accessors expose an independent `minPurity` filter # over these columns. - purityAtCoverage <- function(targetCov, idxVec) { - h <- which(abs(coverageValues - targetCov) < 1e-12) - pv <- if (length(h) > 0L) .csPurityVec(csTables[[h[1L]]]) else numeric() - vapply(idxVec, function(i) { - if (i <= 0L || i > length(pv)) return(0) - v <- pv[i]; if (is.na(v)) 0 else as.numeric(v) - }, numeric(1)) - } # Derive CS membership + purity for EVERY coverage the pipeline actually # produced (attr(csTables, "coverage")), rather than assuming fixed # 0.95/0.70/0.50 levels — otherwise a non-default secondaryCoverage would be @@ -764,8 +776,11 @@ buildTopLoci <- function(fit, csTables, variantNames, sumstats = NULL, # names match getCs's `cs_` lookup. covSorted <- sort(unique(coverageValues[is.finite(coverageValues)]), decreasing = TRUE) - csIdxByCov <- lapply(covSorted, csIdxAtCoverage) - csPurityByCov <- Map(purityAtCoverage, covSorted, csIdxByCov) + csIdxByCov <- lapply(covSorted, .fmCsIdxAtCoverage, coverageValues, csTables, nV) + csPurityByCov <- mapply(.fmPurityAtCoverage, covSorted, csIdxByCov, + MoreArgs = list(coverageValues = coverageValues, + csTables = csTables), + SIMPLIFY = FALSE) csColNames <- paste0("cs_", covSorted * 100) # Per-condition posterior conditional effect + local false sign rate for @@ -1555,6 +1570,18 @@ mvsusieWeights <- function(mvsusieFit = NULL, X = NULL, Y = NULL, return(mvsusieR::coef.mvsusie(mvsusieFit)[-1, ]) } +# One wavelet basis row: inverse-DWT (wr) of the unit coefficient vector e_k, +# using the fit's template DWT object. +# @noRd +.fmReconstructUnit <- function(k, nWac, scaleCols, template) { + coeffRow <- numeric(nWac) + coeffRow[k] <- 1 + temp <- template + temp$D <- coeffRow[-scaleCols] + temp$C[length(temp$C)] <- sum(coeffRow[scaleCols]) + as.numeric(wavethresh::wr(temp)) +} + # Build the wavelet synthesis (inverse-DWT) matrix S (n_wac x nFeat) for the # basis fSuSiE uses, by reconstructing each unit wavelet coefficient through the # SAME $D / $C assignment as out_prep.susiF (detail columns -> $D, the coarsest @@ -1566,15 +1593,7 @@ mvsusieWeights <- function(mvsusieFit = NULL, X = NULL, Y = NULL, # @noRd .fsusieSynthesisMatrix <- function(nWac, scaleCols) { template <- wavethresh::wd(rep(0, nWac)) - reconstructUnit <- function(k) { - coeffRow <- numeric(nWac) - coeffRow[k] <- 1 - temp <- template - temp$D <- coeffRow[-scaleCols] - temp$C[length(temp$C)] <- sum(coeffRow[scaleCols]) - as.numeric(wavethresh::wr(temp)) - } - do.call(rbind, lapply(seq_len(nWac), reconstructUnit)) + do.call(rbind, lapply(seq_len(nWac), .fmReconstructUnit, nWac, scaleCols, template)) } #' Compute fSuSiE feature-level TWAS weights @@ -1769,6 +1788,141 @@ mvsusieRssWeights <- function(stat, LD, mvsusieRssFit = NULL, # Cross-condition credible-set merging # ============================================================================= +# Identify variant IDs that are associated with more than one credible set. +# @noRd +.identifyOverlapSets <- function(variantsSetsAndPipsList) { + overlapSets <- list() + for (variantId in names(variantsSetsAndPipsList)) { + sets <- variantsSetsAndPipsList[[variantId]][["sets"]] + if (length(sets) > 1) { + overlapSets[[variantId]] <- sets + } + } + return(overlapSets) +} + +# Union-find root of `x` following the `parent` map. +# @noRd +.ufFindRoot <- function(x, parent) { + while (!identical(parent[[x]], x)) x <- parent[[x]] + x +} + +# Union-find merge of `a` and `b` in `parent`; returns the updated parent map. +# @noRd +.ufUnion <- function(a, b, parent) { + rootA <- .ufFindRoot(a, parent) + rootB <- .ufFindRoot(b, parent) + if (!identical(rootA, rootB)) parent[[rootB]] <- rootA + parent +} + +# Merge overlapping credible sets using connected components (union-find). +# @noRd +.mergeAndUpdateOverlapSets <- function(variantsSetsAndPipsList, overlapSets) { + allSets <- unique(unlist(overlapSets)) + if (length(allSets) == 0) return(list()) + + parent <- setNames(allSets, allSets) + for (sets in overlapSets) { + if (length(sets) > 1) { + for (s in sets[-1]) parent <- .ufUnion(sets[[1]], s, parent) + } + } + + components <- split(names(parent), + vapply(names(parent), .ufFindRoot, character(1), parent)) + setNameMap <- list() + for (members in components) { + label <- paste(sort(members), collapse = ",") + for (s in members) { + setNameMap[[s]] <- label + } + } + + # Update each variant's credible set names + updatedCredibleSets <- lapply( + setNames(names(variantsSetsAndPipsList), names(variantsSetsAndPipsList)), + function(variantId) { + currentSets <- variantsSetsAndPipsList[[variantId]][["sets"]] + mapped <- intersect(currentSets, names(setNameMap)) + if (length(mapped) > 0) { + setNameMap[[mapped[1]]] + } else { + paste(sort(unique(currentSets)), collapse = ",") + } + } + ) + return(updatedCredibleSets) +} + +# Collapse the per-variant extracted-CS map into a top-loci data frame: merge +# overlapping credible sets, then one row per variant with its merged CS label, +# max PIP and median PIP. +# @noRd +.combineTopLoci <- function(extractedResult) { + if (length(extractedResult) == 0) return(NULL) + + overlapSets <- .identifyOverlapSets(extractedResult) + hasOverlaps <- length(overlapSets) != 0 + mergedSets <- if (hasOverlaps) { + .mergeAndUpdateOverlapSets(extractedResult, overlapSets = overlapSets) + } else { + NULL + } + + topLociDf <- do.call(rbind, lapply(names(extractedResult), function(variantId) { + maxPip <- max(unlist(extractedResult[[variantId]]$pips)) + medianPip <- median(unlist(extractedResult[[variantId]]$pips)) + credibleSetNames <- if (hasOverlaps) { + mergedSets[[variantId]] + } else { + paste(sort(unique(unlist(extractedResult[[variantId]]$sets))), collapse = ",") + } + data.frame( + variant_id = variantId, credibleSetNames = credibleSetNames, + maxPip = maxPip, medianPip = medianPip, stringsAsFactors = FALSE + ) + })) + return(topLociDf) +} + +# Build the per-variant extracted-CS map from a fine-mapping result: for each +# entry, one record per (variant, credible set) labelled cs__, +# aggregated by variant preserving first-seen order. +# @noRd +.fmExtractTopLoci <- function(fineMappingResult, csCol) { + entries <- fineMappingResult$entry + rows <- map_dfr(seq_along(entries), function(i) { + topLoci <- .translateLegacyTopLociCsColumns(getTopLoci(entries[[i]])) + if (is.null(topLoci) || nrow(topLoci) == 0 || !(csCol %in% names(topLoci))) + return(NULL) + pipCol <- resolvePipColumn(topLoci) + if (is.null(pipCol)) return(NULL) + csIdx <- .fmCsIdx(topLoci[[csCol]]) + setNum <- unique(csIdx) + setNum <- setNum[!is.na(setNum) & setNum != 0] + if (length(setNum) == 0) return(NULL) + + map_dfr(setNum, function(sn) { + keep <- !is.na(csIdx) & csIdx == sn + df <- topLoci[keep, c("variant_id", pipCol), drop = FALSE] + names(df)[names(df) == pipCol] <- "pip" + df$set_name <- paste0("cs_", i, "_", sn) + df + }) + }) + + if (is.null(rows) || nrow(rows) == 0) return(list()) + + # Aggregate by variant_id preserving first-seen order. + seenOrder <- unique(rows$variant_id) + splitRows <- split(rows, factor(rows$variant_id, levels = seenOrder)) + lapply(splitRows, function(df) { + list(sets = df$set_name, pips = df$pip) + }) +} + #' Merge SuSiE credible sets across conditions #' #' Reconciles per-condition (univariate) SuSiE fine-mapping into a single set of @@ -1802,129 +1956,12 @@ mergeSusieCs <- function(fineMappingResult, coverage = 0.95) { } csCol <- paste0("cs_", as.integer(round(coverage * 100))) - # Identify variant IDs that are associated with more than one credible set. - identifyOverlapSets <- function(variantsSetsAndPipsList) { - overlapSets <- list() - for (variantId in names(variantsSetsAndPipsList)) { - sets <- variantsSetsAndPipsList[[variantId]][["sets"]] - if (length(sets) > 1) { - overlapSets[[variantId]] <- sets - } - } - return(overlapSets) - } - # Merge overlapping credible sets using connected components. - mergeAndUpdateOverlapSets <- function(variantsSetsAndPipsList, overlapSets) { - allSets <- unique(unlist(overlapSets)) - if (length(allSets) == 0) return(list()) - - parent <- setNames(allSets, allSets) - findRoot <- function(x) { - while (!identical(parent[[x]], x)) x <- parent[[x]] - x - } - unionSets <- function(a, b) { - rootA <- findRoot(a) - rootB <- findRoot(b) - if (!identical(rootA, rootB)) parent[[rootB]] <<- rootA - } - - for (sets in overlapSets) { - if (length(sets) > 1) { - for (s in sets[-1]) unionSets(sets[[1]], s) - } - } - - components <- split(names(parent), vapply(names(parent), findRoot, character(1))) - setNameMap <- list() - for (members in components) { - label <- paste(sort(members), collapse = ",") - for (s in members) { - setNameMap[[s]] <- label - } - } - - # Update each variant's credible set names - updatedCredibleSets <- lapply( - setNames(names(variantsSetsAndPipsList), names(variantsSetsAndPipsList)), - function(variantId) { - currentSets <- variantsSetsAndPipsList[[variantId]][["sets"]] - mapped <- intersect(currentSets, names(setNameMap)) - if (length(mapped) > 0) { - setNameMap[[mapped[1]]] - } else { - paste(sort(unique(currentSets)), collapse = ",") - } - } - ) - return(updatedCredibleSets) - } - # Each row (entry) of the fine-mapping result is one condition. Build a flat # data frame of (variant_id, pip, set_name) across conditions, giving each # condition's credible sets a unique "cs__" label. - extractTopLoci <- function() { - entries <- fineMappingResult$entry - rows <- map_dfr(seq_along(entries), function(i) { - topLoci <- .translateLegacyTopLociCsColumns(getTopLoci(entries[[i]])) - if (is.null(topLoci) || nrow(topLoci) == 0 || !(csCol %in% names(topLoci))) - return(NULL) - pipCol <- resolvePipColumn(topLoci) - if (is.null(pipCol)) return(NULL) - csIdx <- .fmCsIdx(topLoci[[csCol]]) - setNum <- unique(csIdx) - setNum <- setNum[!is.na(setNum) & setNum != 0] - if (length(setNum) == 0) return(NULL) - - map_dfr(setNum, function(sn) { - keep <- !is.na(csIdx) & csIdx == sn - df <- topLoci[keep, c("variant_id", pipCol), drop = FALSE] - names(df)[names(df) == pipCol] <- "pip" - df$set_name <- paste0("cs_", i, "_", sn) - df - }) - }) - - if (is.null(rows) || nrow(rows) == 0) return(list()) - - # Aggregate by variant_id preserving first-seen order. - seenOrder <- unique(rows$variant_id) - splitRows <- split(rows, factor(rows$variant_id, levels = seenOrder)) - lapply(splitRows, function(df) { - list(sets = df$set_name, pips = df$pip) - }) - } - - combineTopLoci <- function(extractedResult) { - if (length(extractedResult) == 0) return(NULL) - - overlapSets <- identifyOverlapSets(extractedResult) - hasOverlaps <- length(overlapSets) != 0 - mergedSets <- if (hasOverlaps) { - mergeAndUpdateOverlapSets(extractedResult, overlapSets = overlapSets) - } else { - NULL - } - - topLociDf <- do.call(rbind, lapply(names(extractedResult), function(variantId) { - maxPip <- max(unlist(extractedResult[[variantId]]$pips)) - medianPip <- median(unlist(extractedResult[[variantId]]$pips)) - credibleSetNames <- if (hasOverlaps) { - mergedSets[[variantId]] - } else { - paste(sort(unique(unlist(extractedResult[[variantId]]$sets))), collapse = ",") - } - data.frame( - variant_id = variantId, credibleSetNames = credibleSetNames, - maxPip = maxPip, medianPip = medianPip, stringsAsFactors = FALSE - ) - })) - return(topLociDf) - } - - extractedTopLoci <- extractTopLoci() + extractedTopLoci <- .fmExtractTopLoci(fineMappingResult, csCol) if (length(extractedTopLoci) == 0) return(NULL) - combinedTopLociDf <- combineTopLoci(extractedTopLoci) + combinedTopLociDf <- .combineTopLoci(extractedTopLoci) if (is.null(combinedTopLociDf) || nrow(combinedTopLociDf) == 0) return(NULL) combinedTopLociDf <- combinedTopLociDf[!duplicated(combinedTopLociDf$variant_id), ] rownames(combinedTopLociDf) <- NULL diff --git a/R/genotypeIo.R b/R/genotypeIo.R index fdb8995c..a9605bc8 100644 --- a/R/genotypeIo.R +++ b/R/genotypeIo.R @@ -38,6 +38,37 @@ setMethod("readGenotypes", # Handle constructors — read metadata, defer genotype loading # ============================================================================= +# Record each variant's original 1-based position in the genotype file. This +# lets @snpInfo be row-subset (e.g. to the range of a study's summary stats) +# while genotype reads that index BY FILE POSITION still resolve correctly: +# PLINK2's ReadList(variant_subset=) and the per-chromosome PLINK2 view inside +# sharded routing read fileIdx[snpIdx]; the by-id backends (plink1/gds/vcf) look +# up snpInfo$SNP[snpIdx] / $BP[snpIdx] and ignore fileIdx entirely. For a full +# (unsubset) handle fileIdx == seq_len(nrow), so reads are unchanged. +# @noRd +.withFileIdx <- function(snpInfo) { + snpInfo$fileIdx <- seq_len(nrow(snpInfo)) + snpInfo +} + +# Restrict a GenotypeHandle's @snpInfo to `keep` (a logical mask or integer row +# indices into @snpInfo). Genotype reads stay correct because fileIdx carries +# each kept variant's original file position; everything else (path/format/ +# pgenPtr/chromPaths/samples) is preserved. NULL-safe; a no-op when nothing is +# dropped. Handles built before the fileIdx column existed are NOT subset (the +# read path would be positional) -- return them unchanged. +# @noRd +.subsetGenotypeHandle <- function(handle, keep) { + if (is.null(handle)) return(NULL) + si <- handle@snpInfo + keepIdx <- if (is.logical(keep)) which(keep) else as.integer(keep) + if (length(keepIdx) >= nrow(si)) return(handle) # nothing dropped + if (!"fileIdx" %in% names(si)) return(handle) # legacy handle: unsafe + handle@snpInfo <- si[keepIdx, , drop = FALSE] + rownames(handle@snpInfo) <- NULL + handle +} + #' @keywords internal .makeGdsHandle <- function(path) { # nocov start @@ -58,7 +89,7 @@ setMethod("readGenotypes", new("GenotypeHandle", path = path, format = "gds", - snpInfo = snpInfo, + snpInfo = .withFileIdx(snpInfo), nSamples = as.integer(nSamples), sampleIds = sampleIds, pgenPtr = NULL @@ -96,7 +127,7 @@ setMethod("readGenotypes", new("GenotypeHandle", path = normalizePath(path), format = "vcf", - snpInfo = snpInfo, + snpInfo = .withFileIdx(snpInfo), nSamples = as.integer(nSamples), sampleIds = sampleIds, pgenPtr = NULL @@ -145,7 +176,7 @@ setMethod("readGenotypes", new("GenotypeHandle", path = stem, format = "plink1", - snpInfo = snpInfo, + snpInfo = .withFileIdx(snpInfo), nSamples = as.integer(nSamples), sampleIds = sampleIds, pgenPtr = NULL @@ -188,7 +219,7 @@ setMethod("readGenotypes", new("GenotypeHandle", path = stem, format = "plink2", - snpInfo = snpInfo, + snpInfo = .withFileIdx(snpInfo), nSamples = as.integer(nSamples), sampleIds = sampleIds, pgenPtr = pgen @@ -267,6 +298,24 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { ) } +# Extract one chromosome's block from a sharded handle: reslice the handle to +# that chrom's file + SNP subset and delegate to extractBlockGenotypes. +# @noRd +.extractBlockForChrom <- function(chrom, posInReq, handle, snpIdx, unifiedChr, meanImpute) { + if (!chrom %in% names(handle@chromPaths)) + stop("extractBlockGenotypes: no per-chromosome file for chromosome '", + chrom, "' (have: ", + paste(names(handle@chromPaths), collapse = ", "), ").") + blockGlobal <- which(unifiedChr == chrom) # file-order global indices + localIdx <- match(snpIdx[posInReq], blockGlobal) + th <- handle + th@path <- handle@chromPaths[[chrom]] + th@snpInfo <- handle@snpInfo[blockGlobal, , drop = FALSE] + th@pgenPtr <- NULL + th@chromPaths <- character(0) # treat as single-file + extractBlockGenotypes(th, localIdx, meanImpute = meanImpute) +} + # Extract a block from a one-file-per-chromosome (sharded) handle. The global # snpIdx index the unified @snpInfo; we group them by chromosome, route each # group to its per-chromosome payload via a transient single-file view (with @@ -290,22 +339,10 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { reqChr <- unifiedChr[snpIdx] groups <- split(seq_along(snpIdx), reqChr) - buildForChrom <- function(chrom, posInReq) { - if (!chrom %in% names(handle@chromPaths)) - stop("extractBlockGenotypes: no per-chromosome file for chromosome '", - chrom, "' (have: ", - paste(names(handle@chromPaths), collapse = ", "), ").") - blockGlobal <- which(unifiedChr == chrom) # file-order global indices - localIdx <- match(snpIdx[posInReq], blockGlobal) - th <- handle - th@path <- handle@chromPaths[[chrom]] - th@snpInfo <- handle@snpInfo[blockGlobal, , drop = FALSE] - th@pgenPtr <- NULL - th@chromPaths <- character(0) # treat as single-file - extractBlockGenotypes(th, localIdx, meanImpute = meanImpute) - } - - ses <- Map(buildForChrom, names(groups), groups) + ses <- mapply(.extractBlockForChrom, names(groups), groups, + MoreArgs = list(handle = handle, snpIdx = snpIdx, + unifiedChr = unifiedChr, meanImpute = meanImpute), + SIMPLIFY = FALSE) if (length(ses) == 1L) return(ses[[1L]]) # Combine at the assay/rowRanges level (rather than rbind-ing the SEs, which @@ -387,15 +424,21 @@ extractBlockGenotypes <- function(handle, snpIdx, meanImpute = TRUE) { # pointer errors out. Opening is cheap relative to dosage extraction. ptr <- getPgenPtr(handle) paths <- resolvePlink2Paths(.genotypeReadPath(handle)) + # `variant_subset` indexes the .pgen by FILE position. `snpIdx` is a position + # into @snpInfo, which may have been row-subset; translate through fileIdx to + # recover the true .pgen index. For a full handle fileIdx == seq_len(nrow), so + # this is a no-op. Older RDS handles predate the column -> fall back to snpIdx. + fileIdx <- getSnpInfo(handle)$fileIdx + variantSubset <- if (is.null(fileIdx)) snpIdx else fileIdx[snpIdx] # A sharded handle routes through a transient view with pgenPtr = NULL (one # pgen per chromosome), and a deserialized pointer is stale; open a fresh # pgen up front in those cases rather than provoking a caught read error. if (is.null(ptr)) ptr <- pgenlibr::NewPgen(paths$pgen) geno <- tryCatch( - pgenlibr::ReadList(ptr, variant_subset = snpIdx, meanimpute = FALSE), + pgenlibr::ReadList(ptr, variant_subset = variantSubset, meanimpute = FALSE), error = function(e) { reopened <- pgenlibr::NewPgen(paths$pgen) - pgenlibr::ReadList(reopened, variant_subset = snpIdx, meanimpute = FALSE) + pgenlibr::ReadList(reopened, variant_subset = variantSubset, meanimpute = FALSE) }) storage.mode(geno) <- "double" geno @@ -983,8 +1026,8 @@ matchVariantsToKeep <- function(variantInfo, keepVariantsPath) { } } -NoSNPsError <- function(message) { - structure(list(message = message), class = c("NoSNPsError", "error", "condition")) +NoSnpsError <- function(message) { + structure(list(message = message), class = c("NoSnpsError", "error", "condition")) } @@ -1038,7 +1081,7 @@ loadGenotypeRegion <- function(genotype, region = NULL, keepIndel = TRUE, if (!is.null(region)) { snpIdx <- .regionToSnpIdx(handleSnpInfo, region) if (length(snpIdx) == 0) { - stop(NoSNPsError(paste("No SNPs found in the specified region", region))) + stop(NoSnpsError(paste("No SNPs found in the specified region", region))) } } else { snpIdx <- seq_len(nrow(handleSnpInfo)) diff --git a/R/gwasSumStats.R b/R/gwasSumStats.R index acbb84c7..6e3faceb 100644 --- a/R/gwasSumStats.R +++ b/R/gwasSumStats.R @@ -71,6 +71,16 @@ NULL # Constructor # ============================================================================= +# Recycle a length-1 per-study scalar to one value per study (or validate an +# already-per-study vector), coerced to numeric. +# @noRd +.recyclePerStudy <- function(v, nm, study) { + if (length(v) == 1L && length(study) > 1L) v <- rep(v, length(study)) + if (length(v) != length(study)) + stop("`", nm, "` must have length 1 or length(study).") + as.numeric(v) +} + #' @title Create a GwasSumStats Collection Object #' @description Construct a \code{GwasSumStats} S4 DFrame-subclass #' collection from per-study tuple vectors and a list of \code{GRanges} @@ -126,14 +136,7 @@ GwasSumStats <- function(study, entry, genome, ldSketch = NULL, stop("length(entry) (", length(entry), ") must equal length(study) (", length(study), ").") } - # Per-study scalars: recycle a single value to one-per-study. - recycle <- function(v, nm) { - if (length(v) == 1L && length(study) > 1L) v <- rep(v, length(study)) - if (length(v) != length(study)) - stop("`", nm, "` must have length 1 or length(study).") - as.numeric(v) - } - varY <- recycle(varY, "varY") + varY <- .recyclePerStudy(varY, "varY", study) cols <- list( study = as.character(study), @@ -144,9 +147,9 @@ GwasSumStats <- function(study, entry, genome, ldSketch = NULL, # supplied (default NULL), so quantitative-trait GwasSumStats keep the # original schema. Provide a vector (length = n studies) with NA for the # non-case/control studies in a mixed collection. - if (!is.null(nCase)) cols$nCase <- recycle(nCase, "nCase") - if (!is.null(nControl)) cols$nControl <- recycle(nControl, "nControl") - if (!is.null(nSample)) cols$nSample <- recycle(nSample, "nSample") + if (!is.null(nCase)) cols$nCase <- .recyclePerStudy(nCase, "nCase", study) + if (!is.null(nControl)) cols$nControl <- .recyclePerStudy(nControl, "nControl", study) + if (!is.null(nSample)) cols$nSample <- .recyclePerStudy(nSample, "nSample", study) extras <- list(...) for (nm in names(extras)) cols[[nm]] <- extras[[nm]] df <- do.call(S4Vectors::DataFrame, c(cols, list(check.names = FALSE))) diff --git a/R/h2EstimationWrappers.R b/R/h2EstimationWrappers.R index aa491a85..dc11c1db 100644 --- a/R/h2EstimationWrappers.R +++ b/R/h2EstimationWrappers.R @@ -1129,6 +1129,72 @@ NULL # Univariate HDL # ============================================================================= +# Per-eigenvalue variance sigma2_i = n/M * sum_a(tau_a * d_i * ldAnnot_{a,i}) + 1 +# (scalar-tau fallback when the block carries no annotation matrix). +# @noRd +.hdlComputeSigma2 <- function(tau, bd, n, M) { + if (!is.null(bd$ldAnnot)) { + n / M * bd$d * as.vector(bd$ldAnnot %*% tau) + 1 + } else { + n * tau[1] * bd$d / M + 1 + } +} + +# HDL stratified negative log-likelihood as a function of the tau vector (with +# optional L2 penalty `lambda`); `...` absorbs the gradient's extra optim args. +# @noRd +.hdlNll <- function(tau, blockData, n, M, lambda, ...) { + val <- 0 + for (bd in blockData) { + sigma2 <- .hdlComputeSigma2(tau, bd, n, M) + if (any(sigma2 <= 0)) return(1e10) + val <- val + 0.5 * sum(log(sigma2) + bd$zRot^2 / sigma2) + } + if (lambda > 0) val <- val + lambda * sum(tau^2) + val +} + +# Gradient of the HDL stratified negative log-likelihood w.r.t. tau. +# @noRd +.hdlNllGrad <- function(tau, blockData, n, M, lambda, nTau) { + grad <- numeric(nTau) + for (bd in blockData) { + sigma2 <- .hdlComputeSigma2(tau, bd, n, M) + if (any(sigma2 <= 0)) return(rep(0, nTau)) + # dsigma2/dtau_a = n/M * d_i * ld_annot_{a,i} + dsig <- n / M * bd$d * bd$ldAnnot # (nEigen x nTau) + # dNLL/dsigma2_i = 0.5 * (1/sigma2_i - z_rot_i^2/sigma2_i^2) + dNLL_dsig <- 0.5 * (1 / sigma2 - bd$zRot^2 / sigma2^2) + grad <- grad + as.vector(crossprod(dsig, dNLL_dsig)) + } + if (lambda > 0) grad <- grad + 2 * lambda * tau + grad +} + +# Single-parameter HDL negative log-likelihood (unstratified scalar h2). +# @noRd +.hdlNllScalar <- function(h2, blockData, n, M) { + val <- 0 + for (bd in blockData) { + sigma2 <- n * h2 * bd$d / M + 1 + val <- val + 0.5 * sum(log(sigma2) + bd$zRot^2 / sigma2) + } + val +} + +# Local HDL negative log-likelihood: a local deviation deltaH2 on top of the +# block's baseline variance. +# @noRd +.hdlNllLocal <- function(deltaH2, sigma2Baseline, bd, n, M) { + if (!is.null(sigma2Baseline)) { + sigma2 <- sigma2Baseline + n * deltaH2 * bd$d / M + } else { + sigma2 <- n * deltaH2 * bd$d / M + 1 + } + if (any(sigma2 <= 0)) return(1e10) + 0.5 * sum(log(sigma2) + bd$zRot^2 / sigma2) +} + #' @title Univariate HDL #' @description Estimate h2 via HDL likelihood. #' @param z Numeric vector of z-scores. @@ -1173,67 +1239,23 @@ hdlUnivariate <- function(z, n, eigenRef, annotations = NULL, snpIdx = idx, p = length(idx)) }) - # Compute per-eigenvalue variance: sigma2_i = n/M * sum_a(tau_a * d_i * ld_annot_{a,i}) + 1 - .computeSigma2 <- function(tau, bd) { - if (!is.null(bd$ldAnnot)) { - n / M * bd$d * as.vector(bd$ldAnnot %*% tau) + 1 - } else { - n * tau[1] * bd$d / M + 1 - } - } - if (!is.null(baselineMat)) { nTau <- ncol(baselineMat) - # Negative log-likelihood as function of tau vector (with optional L2 penalty) - nll <- function(tau) { - val <- 0 - for (bd in blockData) { - sigma2 <- .computeSigma2(tau, bd) - if (any(sigma2 <= 0)) return(1e10) - val <- val + 0.5 * sum(log(sigma2) + bd$zRot^2 / sigma2) - } - if (lambda > 0) val <- val + lambda * sum(tau^2) - val - } - - # Gradient of negative log-likelihood - nllGrad <- function(tau) { - grad <- numeric(nTau) - for (bd in blockData) { - sigma2 <- .computeSigma2(tau, bd) - if (any(sigma2 <= 0)) return(rep(0, nTau)) - # dsigma2/dtau_a = n/M * d_i * ld_annot_{a,i} - dsig <- n / M * bd$d * bd$ldAnnot # (nEigen x nTau) - # dNLL/dsigma2_i = 0.5 * (1/sigma2_i - z_rot_i^2/sigma2_i^2) - dNLL_dsig <- 0.5 * (1 / sigma2 - bd$zRot^2 / sigma2^2) - grad <- grad + as.vector(crossprod(dsig, dNLL_dsig)) - } - if (lambda > 0) grad <- grad + 2 * lambda * tau - grad - } - # Initialize with uniform h2 across annotations tauInit <- rep(0.5 / nTau, nTau) - opt <- optim(tauInit, nll, gr = nllGrad, method = "BFGS", - control = list(maxit = 200, reltol = 1e-8)) + opt <- optim(tauInit, .hdlNll, gr = .hdlNllGrad, + blockData = blockData, n = n, M = M, lambda = lambda, nTau = nTau, + method = "BFGS", control = list(maxit = 200, reltol = 1e-8)) tau <- opt$par # h2 = sum_a tau_a * M_a h2 <- sum(tau * colSums(baselineMat)) } else { # Single-parameter case: use optimize (current behavior) - nll <- function(h2) { - val <- 0 - for (bd in blockData) { - sigma2 <- n * h2 * bd$d / M + 1 - val <- val + 0.5 * sum(log(sigma2) + bd$zRot^2 / sigma2) - } - val - } - - opt <- optimize(nll, interval = c(-0.5, 1.5), tol = 1e-8) + opt <- optimize(.hdlNllScalar, interval = c(-0.5, 1.5), + blockData = blockData, n = n, M = M, tol = 1e-8) h2 <- opt$minimum tau <- h2 # scalar, used by downstream functions } @@ -1407,16 +1429,8 @@ hdlUnivariate <- function(z, n, eigenRef, annotations = NULL, # Local likelihood: optimize a local deviation deltaH2 # sigma2_i = sigma2_baseline_i + n * delta_h2 * d_i / M - nllLocal <- function(deltaH2) { - if (!is.null(sigma2Baseline)) { - sigma2 <- sigma2Baseline + n * deltaH2 * bd$d / M - } else { - sigma2 <- n * deltaH2 * bd$d / M + 1 - } - if (any(sigma2 <= 0)) return(1e10) - 0.5 * sum(log(sigma2) + bd$zRot^2 / sigma2) - } - opt <- optimize(nllLocal, interval = c(-0.5, 0.5), tol = 1e-8) + opt <- optimize(.hdlNllLocal, interval = c(-0.5, 0.5), + sigma2Baseline = sigma2Baseline, bd = bd, n = n, M = M, tol = 1e-8) deltaH2 <- opt$minimum # Total local h2 = global baseline contribution + local deviation diff --git a/R/jointEngine.R b/R/jointEngine.R index a517fe50..e150db51 100644 --- a/R/jointEngine.R +++ b/R/jointEngine.R @@ -20,18 +20,23 @@ NULL # ---- identity derivation ---------------------------------------------------- +# A group's constant value on axis `ax` (study/context/trait), or NULL when the +# axis varies (jointed) so the prior lookup matches any value there. +# @noRd +.jointAxisVal <- function(ax, conditions) { + u <- unique(as.character(conditions[[ax]])) + if (length(u) > 1L) NULL else u[[1L]] +} + # The data-driven-prior LOOKUP key for a group's conditions: a varying (jointed) # axis -> NULL (match-any, because the shared joint mr.mash fit lives on every # per-context row), a constant axis -> its single value. Used only to find the # mr.mash fit; the OUTPUT rows carry each condition's REAL (study, context, # trait). .jointPriorKey <- function(conditions) { - axisVal <- function(ax) { - u <- unique(as.character(conditions[[ax]])) - if (length(u) > 1L) NULL else u[[1L]] - } - list(study = axisVal("study"), context = axisVal("context"), - trait = axisVal("trait")) + list(study = .jointAxisVal("study", conditions), + context = .jointAxisVal("context", conditions), + trait = .jointAxisVal("trait", conditions)) } # The ";"-joined distinct members of a varying axis (the per-row provenance tag @@ -69,17 +74,11 @@ NULL foldFits = cvRes$foldFits) } -# 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$records <- list() - e$add <- function(...) e$records[[length(e$records) + 1L]] <- list(...) - e -} +# The engine assembles one per-row RECORD per fitted entry: a named list whose +# names match the target collection constructor's parameters (see +# .jointEntryRecords()). The driver collects these into a plain list and +# construct() folds them into a collection via .buildJointResult(). Nothing here +# enumerates columns, so adding a column needs only that it be put in the record. # --- Trait-position and fine-mapping-region provenance ---------------------- # Two DISTINCT per-row anchors, kept separate because they mean different things @@ -323,25 +322,29 @@ setMethod("fitJointGroup", signature("SumStatsJointGroup", "FmJointPipeline"), includeAllCs = cfg$includeAllCs)) }) +# Select the list element whose name (stripped of a _predicted/_performance +# suffix) matches `token`; NULL if none. +# @noRd +.jointPickByBase <- function(lst, token) { + if (is.null(lst) || length(lst) == 0L) return(NULL) + bare <- sub("(_predicted|Predicted|_performance|Performance)$", "", names(lst)) + hit <- which(bare == token) + if (length(hit) == 0L) NULL else lst[[hit[[1L]]]] +} + # Reshape a twasWeightsCv() result into the single joint entry's cvResult: the # out-of-fold prediction matrix, the per-condition metric rows, and the per-fold # mr.mash fits (named fold_) that fineMappingPipeline's mvSuSiE path consumes. .jointTwasCvResult <- function(cv, token) { if (is.null(cv)) return(NULL) - pickByBase <- function(lst) { - if (is.null(lst) || length(lst) == 0L) return(NULL) - bare <- sub("(_predicted|Predicted|_performance|Performance)$", "", names(lst)) - hit <- which(bare == token) - if (length(hit) == 0L) NULL else lst[[hit[[1L]]]] - } ffKey <- paste0(token, "_weights") foldFits <- if (!is.null(cv$foldFits)) { ff <- lapply(cv$foldFits, function(f) f[[ffKey]]) if (all(vapply(ff, is.null, logical(1)))) NULL else ff } else NULL list(samplePartition = cv$samplePartition, - predictions = pickByBase(cv$prediction), - metrics = pickByBase(cv$performance), + predictions = .jointPickByBase(cv$prediction, token), + metrics = .jointPickByBase(cv$performance, token), foldFits = foldFits) } @@ -520,13 +523,13 @@ setMethod("fitJointGroup", signature("SumStatsJointGroup", "TwasJointPipeline"), } setMethod("construct", "FmJointPipeline", - function(pipeline, rows, ...) - .buildJointResult(QtlFineMappingResult, rows$records, + function(pipeline, records, ...) + .buildJointResult(QtlFineMappingResult, records, pipeline@config$ldSketch)) setMethod("construct", "TwasJointPipeline", - function(pipeline, rows, ...) - .buildJointResult(TwasWeights, rows$records, pipeline@config$ldSketch)) + function(pipeline, records, ...) + .buildJointResult(TwasWeights, records, pipeline@config$ldSketch)) # ---- enumerators (pattern x dataForm -> list) -------------------- @@ -541,7 +544,7 @@ setMethod("construct", "TwasJointPipeline", verbose <- if (is.null(args$verbose)) 1 else args$verbose groups <- list() for (tid in scopedTraits) { - xy <- .buildIndividualCrossContextXY( + xy <- .buildIndividualCrossContextXy( data, tid, scopedContexts, args$cisWindow, verbose, label = "jointCrossContext", region = args$region) if (is.null(xy)) next @@ -591,7 +594,7 @@ setMethod("construct", "TwasJointPipeline", verbose <- if (is.null(args$verbose)) 1 else args$verbose groups <- list() for (cx in scopedContexts) { - xy <- .buildIndividualCrossTraitXY( + xy <- .buildIndividualCrossTraitXy( data, cx, scopedTraits, args$cisWindow, verbose, label = "jointCrossTrait", study = study, region = args$region) if (is.null(xy)) next @@ -677,7 +680,7 @@ setMethod("construct", "TwasJointPipeline", study <- getStudy(data) if (!(study %in% scope$studies)) return(list()) verbose <- if (is.null(args$verbose)) 1 else args$verbose - xy <- .buildComposedIndividualXY( + xy <- .buildComposedIndividualXy( data, scope, study, args$cisWindow, verbose, label = "composed", region = args$region) if (is.null(xy)) return(list()) @@ -783,6 +786,32 @@ setMethod("construct", "TwasJointPipeline", out } +# Append one output record per (condition, method) to the joint-rows accumulator +# `rows` (mutated by reference), resolving each condition's fine-mapping region + +# trait position. `grp` bundles the per-group invariants list(cond, js, jc, jt). +# @noRd +.jointEntryRecords <- function(entries, method, grp, data, cisWindow) { + cond <- grp$cond + recs <- lapply(seq_len(min(length(entries), nrow(cond))), function(i) { + e <- entries[[i]] + if (is.null(e)) return(NULL) + ctx <- as.character(cond$context[[i]]) + tid <- as.character(cond$trait[[i]]) + # region = the fine-mapping window (cis-window-expanded for a QtlDataset, + # variant span for a QtlSumStats); traitPos = the bare trait position + # (NULL when a QtlSumStats caller supplied none -> the column is omitted + # and getTraitPosition() reports NA). + reg <- .fitRegionFor(data, ctx, tid, cisWindow) + tpos <- .traitPosFor(data, ctx, tid) + list(study = as.character(cond$study[[i]]), + context = ctx, trait = tid, + method = method, entry = e, + jointStudies = grp$js, jointContexts = grp$jc, jointTraits = grp$jt, + region = reg, traitPos = tpos) + }) + Filter(Negate(is.null), recs) +} + # Run one dispatch cell: enumerate joint groups, fit each method (S4 dispatch on # the group x pipeline pair) per group, accumulate per-context rows, build the # per-pipeline result. The loop is GROUP-outer / token-inner so the twas ensemble @@ -795,40 +824,22 @@ setMethod("construct", "TwasJointPipeline", if (length(groups) == 0L) return(NULL) doEnsemble <- is(pipeline, "TwasJointPipeline") && isTRUE(pipeline@config$ensemble) - rows <- .jointRows() + records <- list() for (g in groups) { cond <- g@conditions # Provenance: the ";"-joined members of each varying axis, identical on every # per-context row of this joint group. - js <- .jointAxisMembers(cond, "study") - jc <- .jointAxisMembers(cond, "context") - jt <- .jointAxisMembers(cond, "trait") - addEntries <- function(entries, method) { - 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]]) - # region = the fine-mapping window (cis-window-expanded for a QtlDataset, - # variant span for a QtlSumStats); traitPos = the bare trait position - # (NULL when a QtlSumStats caller supplied none -> the column is omitted - # and getTraitPosition() reports NA). - reg <- .fitRegionFor(data, ctx, tid, args$cisWindow) - tpos <- .traitPosFor(data, ctx, tid) - rows$add(study = as.character(cond$study[[i]]), - context = ctx, trait = tid, - method = method, entry = e, - jointStudies = js, jointContexts = jc, jointTraits = jt, - region = reg, traitPos = tpos) - } - } + grp <- list(cond = cond, + js = .jointAxisMembers(cond, "study"), + jc = .jointAxisMembers(cond, "context"), + jt = .jointAxisMembers(cond, "trait")) # Twas: resolve this group's fine-mapping fits + CV (keyed on its first # condition; the joint fit is shared across conditions) and fix ONE fold # partition up front, so every method's out-of-fold CV predictions are # aligned for the ensemble layer. FM leaves args untouched. fitArgs <- .twasGroupArgs(g, pipeline, args) - # Per-method fit -> per-condition entries -> rows (shared FM + twas). Retain - # each method's entries so the twas ensemble layer can combine them. + # Per-method fit -> per-condition entries -> records (shared FM + twas). + # Retain each method's entries so the twas ensemble layer can combine them. perTokenEntries <- list() for (token in tokens) { # Resume cache: if every condition of this group is already present in the @@ -848,17 +859,19 @@ setMethod("construct", "TwasJointPipeline", if (is.null(entries)) entries <- fitJointGroup(g, pipeline, token, fitArgs) if (is.null(entries) || length(entries) == 0L) next perTokenEntries[[token]] <- entries - addEntries(entries, token) + records <- c(records, + .jointEntryRecords(entries, token, grp, data, args$cisWindow)) } # SR-TWAS ensemble layer: combine the group's per-method per-condition fits # (CV predictions + weights) into ensemble per-context rows -- built ON TOP # of the shared per-method fitting above, never inside it. if (doEnsemble && length(perTokenEntries) >= 2L) { - addEntries(.twasEnsembleLayer(g, perTokenEntries, pipeline@config), - "ensemble") + records <- c(records, .jointEntryRecords( + .twasEnsembleLayer(g, perTokenEntries, pipeline@config), + "ensemble", grp, data, args$cisWindow)) } } - construct(pipeline, rows) + construct(pipeline, records) } # SR-TWAS ensemble LAYER (twas only): combine a group's per-method per-condition diff --git a/R/jointSpecification.R b/R/jointSpecification.R index f24b4773..c743a5d6 100644 --- a/R/jointSpecification.R +++ b/R/jointSpecification.R @@ -404,6 +404,24 @@ parseTraitIds <- function(traitId, data) { out } +# Validate one leaf method vector: non-empty character, all tokens known (in +# `caps`), and none in `rejectedAtUser`. +# @noRd +.jointValidateLeafVec <- function(vec, label, caps, rejectedAtUser) { + if (!is.character(vec) || length(vec) == 0L) + stop(label, ": method vector must be a non-empty character vector") + bad <- setdiff(vec, names(caps)) + if (length(bad) > 0L) + stop(label, ": unknown method token(s): ", + paste(bad, collapse = ", "), + ". Known tokens: ", paste(names(caps), collapse = ", ")) + rejected <- intersect(vec, rejectedAtUser) + if (length(rejected) > 0L) + stop(label, ": method(s) cannot be user-requested on this pipeline: ", + paste(rejected, collapse = ", ")) + invisible(NULL) +} + # @noRd parseMethods <- function(methods, sumStatsMethods = NULL, @@ -430,24 +448,9 @@ parseMethods <- function(methods, stop("`qtlDatasetMethods` must be a non-empty character vector.") } - validateLeafVec <- function(vec, label) { - if (!is.character(vec) || length(vec) == 0L) - stop(label, ": method vector must be a non-empty character vector") - bad <- setdiff(vec, names(caps)) - if (length(bad) > 0L) - stop(label, ": unknown method token(s): ", - paste(bad, collapse = ", "), - ". Known tokens: ", paste(names(caps), collapse = ", ")) - rejected <- intersect(vec, rejectedAtUser) - if (length(rejected) > 0L) - stop(label, ": method(s) cannot be user-requested on this pipeline: ", - paste(rejected, collapse = ", ")) - invisible(NULL) - } - if (splitGiven) { - validateLeafVec(sumStatsMethods, "sumStatsMethods") - validateLeafVec(qtlDatasetMethods, "qtlDatasetMethods") + .jointValidateLeafVec(sumStatsMethods, "sumStatsMethods", caps, rejectedAtUser) + .jointValidateLeafVec(qtlDatasetMethods, "qtlDatasetMethods", caps, rejectedAtUser) } else { walked <- .spWalkMethods(methods, label = "methods", maxDepth = 3L) studyNames <- .spListStudies(data) @@ -455,7 +458,7 @@ parseMethods <- function(methods, lab <- sprintf("methods[[%s]]", paste0("'", leaf$path, "'", collapse = "$")) if (length(leaf$path) == 0L) lab <- "methods" - validateLeafVec(leaf$methods, lab) + .jointValidateLeafVec(leaf$methods, lab, caps, rejectedAtUser) # Multi-axis methods may not appear at per-context or per-trait levels. if (leaf$depth >= 2L) { bad <- intersect(leaf$methods, multivariateMethods) @@ -638,7 +641,7 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { # NULL when fewer than 2 contexts carry `tid` or the sample / complete-Y # subset is too small to fit. # @noRd -.buildIndividualCrossContextXY <- function(data, tid, scopedContexts, +.buildIndividualCrossContextXy <- function(data, tid, scopedContexts, cisWindow, verbose, label, region = NULL) { perTraitContexts <- character(0) @@ -707,7 +710,7 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { # when fewer than 2 traits live in the context or the sample / complete-Y # subset is too small. # @noRd -.buildIndividualCrossTraitXY <- function(data, cx, scopedTraits, +.buildIndividualCrossTraitXy <- function(data, cx, scopedTraits, cisWindow, verbose, label, study, region = NULL) { se <- getPhenotypes(data, contexts = cx) @@ -742,7 +745,7 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { # Build a composed-axes (context, trait) X/Y for individual-level # QtlDataset. Returns list(X, Y, tuples) or NULL. # @noRd -.buildComposedIndividualXY <- function(data, scope, study, cisWindow, +.buildComposedIndividualXy <- function(data, scope, study, cisWindow, verbose, label, region = NULL) { scopedContexts <- scope$contexts[[study]] scopedTraits <- scope$traits[[study]] @@ -832,6 +835,13 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { # Fine-mapping dispatchers # ============================================================================= +# Identity-tuple key (study/context/trait/method joined by "\r"), used to align +# per-region result entries when merging. Shared by the fm/twas mergers. +# @noRd +.mergeResultKeyOf <- function(r) paste(as.character(r$study), as.character(r$context), + as.character(r$trait), as.character(r$method), + sep = "\r") + # Top-level joint dispatcher for fineMappingPipeline(QtlDataset). # @noRd # Merge per-region QtlFineMappingResult collections (same keys across regions) @@ -842,13 +852,10 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { base <- results[[1L]] n <- nrow(base) if (n == 0L) return(base) - keyOf <- function(r) paste(as.character(r$study), as.character(r$context), - as.character(r$trait), as.character(r$method), - sep = "\r") - baseKeys <- keyOf(base) + baseKeys <- .mergeResultKeyOf(base) mergedEntries <- lapply(seq_len(n), function(i) { perRegion <- lapply(results, function(r) { - hit <- which(keyOf(r) == baseKeys[[i]]) + hit <- which(.mergeResultKeyOf(r) == baseKeys[[i]]) if (length(hit)) r$entry[[hit[[1L]]]] else NULL }) .fmMergeEntries(Filter(Negate(is.null), perRegion)) @@ -861,22 +868,30 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { list(ldSketch = NULL))) } +# One passthrough column of a joint result row as a character vector (the joint- +# key columns), or NULL when absent. +# @noRd +.jointStrCol <- function(nm, df) if (nm %in% names(df)) as.character(df[[nm]]) else NULL + +# One passthrough column carried through UNCOERCED (GRanges provenance columns +# region / traitPos), or NULL when absent. +# @noRd +.jointRawCol <- function(nm, df) if (nm %in% names(df)) df[[nm]] else NULL + # The optional passthrough columns of a per-tuple result row, as a named list # (NULL for any absent column): the three joint-key columns (jointStudies / # jointContexts / jointTraits) plus the GRanges provenance columns (region / # traitPos). Spliced into the QtlFineMappingResult / TwasWeights constructors so # by-key / cross-study rebuilds preserve them. .jointCols <- function(df) { - pick <- function(nm) if (nm %in% names(df)) as.character(df[[nm]]) else NULL # region / traitPos are GRanges provenance columns: carry them through # uncoerced so by-key / cross-study rebuilds keep the fine-mapping window and # trait position instead of silently dropping them. - pickCol <- function(nm) if (nm %in% names(df)) df[[nm]] else NULL - list(jointStudies = pick("jointStudies"), - jointContexts = pick("jointContexts"), - jointTraits = pick("jointTraits"), - region = pickCol("region"), - traitPos = pickCol("traitPos")) + list(jointStudies = .jointStrCol("jointStudies", df), + jointContexts = .jointStrCol("jointContexts", df), + jointTraits = .jointStrCol("jointTraits", df), + region = .jointRawCol("region", df), + traitPos = .jointRawCol("traitPos", df)) } # Shared tail of the MultiStudyQtlDataset fineMapping / twasWeights pipeline @@ -889,19 +904,19 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { # `jointResult` in. # @noRd .multiStudyPipelineDriver <- function(data, jointResult, perStudyFn, sumStatsFn, - rbindFn, resultCtor, pipelineName, + cfg, rbindFn, resultCtor, pipelineName, noun = "a result") { qtlDatasets <- getQtlDatasets(data) sumStats <- getSumStats(data) out <- NULL embeddedLd <- NULL for (qdName in names(qtlDatasets)) { - res <- perStudyFn(qtlDatasets[[qdName]]) + res <- perStudyFn(qtlDatasets[[qdName]], cfg) if (!is.null(res)) out <- if (is.null(out)) res else rbindFn(out, res, ldSketch = NULL) } if (!is.null(sumStats)) { - ssRes <- sumStatsFn(sumStats) + ssRes <- sumStatsFn(sumStats, cfg) if (!is.null(ssRes)) { embeddedLd <- getLdSketch(ssRes) out <- if (is.null(out)) ssRes else rbindFn(out, ssRes, ldSketch = embeddedLd) @@ -1139,13 +1154,10 @@ validateMethodsVsJointSpec <- function(methodsParsed, jointSpecParsed) { base <- results[[1L]] n <- length(base$method) if (n == 0L) return(base) - keyOf <- function(r) paste(as.character(r$study), as.character(r$context), - as.character(r$trait), as.character(r$method), - sep = "\r") - baseKeys <- keyOf(base) + baseKeys <- .mergeResultKeyOf(base) mergedEntries <- lapply(seq_len(n), function(i) { perRegion <- lapply(results, function(r) { - hit <- which(keyOf(r) == baseKeys[[i]]) + hit <- which(.mergeResultKeyOf(r) == baseKeys[[i]]) if (length(hit)) r$entry[[hit[[1L]]]] else NULL }) keep <- !vapply(perRegion, is.null, logical(1)) diff --git a/R/ld.R b/R/ld.R index 725d42ad..46dd783e 100644 --- a/R/ld.R +++ b/R/ld.R @@ -217,23 +217,25 @@ extractLdForRegion <- function(ldMatrix, variants, region, extractCoordinates) { list(extractedLdMatrix = mat, extractedLdVariants = extracted) } +# Concatenate per-block variant-id lists into one deduplicated vector, dropping a +# repeated boundary variant shared between adjacent blocks. +# @noRd +.ldMergeVariants <- function(variantList) { + merged <- character(0) + for (v in variantList) { + ids <- if (is.list(v) && !is.null(v$variants)) v$variants else v + if (length(ids) == 0) next + if (length(merged) > 0 && tail(merged, 1) == ids[1]) ids <- ids[-1] + merged <- c(merged, ids) + } + merged +} + #' Combine multiple block-level LD matrices into one, handling boundary overlaps. #' @importFrom utils tail #' @noRd createLdMatrix <- function(ldMatrices, variants) { - # Merge variant lists, deduplicating boundary overlaps - mergeVariants <- function(variantList) { - merged <- character(0) - for (v in variantList) { - ids <- if (is.list(v) && !is.null(v$variants)) v$variants else v - if (length(ids) == 0) next - if (length(merged) > 0 && tail(merged, 1) == ids[1]) ids <- ids[-1] - merged <- c(merged, ids) - } - merged - } - - allVariants <- mergeVariants(variants) + allVariants <- .ldMergeVariants(variants) combined <- matrix(0, nrow = length(allVariants), ncol = length(allVariants), dimnames = list(allVariants, allVariants)) @@ -1387,6 +1389,14 @@ dropCollinearColumns <- function(X, problematicCols, X[, !(colnames(X) %in% colToRemove), drop = FALSE] } +# Design matrix [1 | X | C] with the intercept + X columns named. +# @noRd +.ldBuildDesign <- function(X, C) { + XD <- cbind(1, X, C) + colnames(XD)[seq_len(ncol(X) + 1L)] <- c("Intercept", colnames(X)) + XD +} + #' Iteratively enforce full column rank on a design matrix #' #' Given a candidate predictor matrix \code{X} and an optional unnamed @@ -1437,13 +1447,7 @@ enforceDesignFullRank <- function(X, C, initialNcol <- ncol(X) iteration <- 0L - buildDesign <- function(X) { - XD <- cbind(1, X, C) - colnames(XD)[seq_len(ncol(X) + 1L)] <- c("Intercept", colnames(X)) - XD - } - - Xdesign <- buildDesign(X) + Xdesign <- .ldBuildDesign(X, C) matrixRank <- qr(Xdesign)$rank if (verbose) { message("enforceDesignFullRank: initial rank ", matrixRank, @@ -1461,7 +1465,7 @@ enforceDesignFullRank <- function(X, C, if (length(problematicColnames) > 0) { Xtemp <- X[, !(colnames(X) %in% problematicColnames), drop = FALSE] - if (qr(buildDesign(Xtemp))$rank == ncol(buildDesign(Xtemp))) { + if (qr(.ldBuildDesign(Xtemp, C))$rank == ncol(.ldBuildDesign(Xtemp, C))) { if (verbose) { message("enforceDesignFullRank: full rank after batch-removing ", length(problematicColnames), " column(s)") @@ -1489,7 +1493,7 @@ enforceDesignFullRank <- function(X, C, X <- dropCollinearColumns(X, problematicColnames, strategy = strategy, response = response, verbose = verbose) - Xdesign <- buildDesign(X) + Xdesign <- .ldBuildDesign(X, C) matrixRank <- qr(Xdesign)$rank iteration <- iteration + 1L if (verbose) { @@ -1504,7 +1508,7 @@ enforceDesignFullRank <- function(X, C, } # Correlation-threshold fallback. - Xdesign <- buildDesign(X) + Xdesign <- .ldBuildDesign(X, C) matrixRank <- qr(Xdesign)$rank if (matrixRank < ncol(Xdesign)) { if (verbose) { @@ -1514,7 +1518,7 @@ enforceDesignFullRank <- function(X, C, filterResult <- ldPruneByCorrelation(X, corThres = threshold, verbose = verbose) X <- filterResult$X.new - Xdesign <- buildDesign(X) + Xdesign <- .ldBuildDesign(X, C) matrixRank <- qr(Xdesign)$rank if (verbose) { message("enforceDesignFullRank: threshold ", threshold, @@ -1686,16 +1690,17 @@ extractLdMatrix <- function(ld, wantGenotype = FALSE) { #' @param maxVariants Integer or \code{NULL}. If set, randomly subsample #' blocks larger than this to control memory usage. #' -#' @return A function \code{loader(g)} that, given a block index \code{g}, -#' returns the corresponding LD matrix or genotype matrix. +#' @return An \code{ldLoaderSpec} object (an opaque list describing the source). +#' Pass it with a block index to \code{\link{loadLdBlock}} to load one block. #' +#' @seealso \code{\link{loadLdBlock}} #' @examples #' # List mode with pre-computed LD #' R1 <- diag(10) #' R2 <- diag(15) -#' loader <- ldLoader(rList = list(R1, R2)) -#' loader(1) # returns R1 -#' loader(2) # returns R2 +#' spec <- ldLoader(rList = list(R1, R2)) +#' loadLdBlock(spec, 1) # returns R1 +#' loadLdBlock(spec, 2) # returns R2 #' #' @export ldLoader <- function(rList = NULL, xList = NULL, @@ -1703,90 +1708,131 @@ ldLoader <- function(rList = NULL, xList = NULL, ldInfo = NULL, returnGenotype = FALSE, maxVariants = NULL) { - # Validate: exactly one source + # Validate eagerly, at spec construction: exactly one source, plus the + # per-mode requirements the branch loaders used to check lazily. nSources <- sum(!is.null(rList), !is.null(xList), !is.null(ldMetaPath), !is.null(ldInfo)) if (nSources != 1) stop("Provide exactly one of rList, xList, ldMetaPath, or ldInfo.") + if (!is.null(ldMetaPath) && is.null(regions)) + stop("'regions' is required when using ldMetaPath.") + if (!is.null(ldInfo) && (!is.data.frame(ldInfo) || + !"LD_file" %in% colnames(ldInfo))) + stop("ldInfo must be a data.frame with column 'LD_file'.") + + mode <- if (!is.null(rList)) "rList" + else if (!is.null(xList)) "xList" + else if (!is.null(ldMetaPath)) "meta" + else "info" + structure( + list(mode = mode, rList = rList, xList = xList, + ldMetaPath = ldMetaPath, regions = regions, ldInfo = ldInfo, + returnGenotype = returnGenotype, maxVariants = maxVariants), + class = "ldLoaderSpec") +} - if (!is.null(rList)) { - # List mode (R matrices) - loader <- function(g) { - R <- rList[[g]] - if (!is.null(maxVariants) && ncol(R) > maxVariants) { - keep <- sort(sample(ncol(R), maxVariants)) - R <- R[keep, keep] - } - R - } - } else if (!is.null(xList)) { - # List mode (genotype matrices) - loader <- function(g) { - X <- xList[[g]] - if (!is.null(maxVariants) && ncol(X) > maxVariants) { - keep <- sort(sample(ncol(X), maxVariants)) - X <- X[, keep] - } - X - } - } else if (!is.null(ldMetaPath)) { - # Region mode: load on the fly via loadLdMatrix() - if (is.null(regions)) - stop("'regions' is required when using ldMetaPath.") - - loader <- function(g) { - ld <- loadLdMatrix(ldMetaPath, region = regions[g], - returnGenotype = returnGenotype) - mat <- extractLdMatrix(ld, wantGenotype = returnGenotype) - if (!is.null(maxVariants) && ncol(mat) > maxVariants) { - keep <- sort(sample(ncol(mat), maxVariants)) - if (returnGenotype || nrow(mat) > ncol(mat)) { - mat <- mat[, keep] - } else { - mat <- mat[keep, keep] - } - } - # Center and scale genotype matrices - if (returnGenotype || nrow(mat) > ncol(mat)) { - mat <- scale(mat) - mat[is.na(mat)] <- 0 - } - mat + +# ---- Per-block loaders (one per ldLoader source mode) ----------------------- +# Top-level workers dispatched by .ldLoadBlock() on the spec's $mode; formerly +# branch closures inside ldLoader(). + +# @noRd +.ldLoadRList <- function(g, rList, maxVariants) { + R <- rList[[g]] + if (!is.null(maxVariants) && ncol(R) > maxVariants) { + keep <- sort(sample(ncol(R), maxVariants)) + R <- R[keep, keep] + } + R +} + +# @noRd +.ldLoadXList <- function(g, xList, maxVariants) { + X <- xList[[g]] + if (!is.null(maxVariants) && ncol(X) > maxVariants) { + keep <- sort(sample(ncol(X), maxVariants)) + X <- X[, keep] + } + X +} + +# @noRd +.ldLoadRegionMeta <- function(g, ldMetaPath, regions, returnGenotype, + maxVariants) { + ld <- loadLdMatrix(ldMetaPath, region = regions[g], + returnGenotype = returnGenotype) + mat <- extractLdMatrix(ld, wantGenotype = returnGenotype) + if (!is.null(maxVariants) && ncol(mat) > maxVariants) { + keep <- sort(sample(ncol(mat), maxVariants)) + if (returnGenotype || nrow(mat) > ncol(mat)) { + mat <- mat[, keep] + } else { + mat <- mat[keep, keep] } - } else { - # ldInfo mode: load LD blocks by index from file paths - # Supports all genotype formats (PLINK2, PLINK1, VCF, GDS) and - # pre-computed .cor.xz + .bim/.pvar blocks - if (!is.data.frame(ldInfo) || !"LD_file" %in% colnames(ldInfo)) - stop("ldInfo must be a data.frame with column 'LD_file'.") - - loader <- function(g) { - ldPath <- ldInfo$LD_file[g] - - # Auto-detect format: genotype source or pre-computed block - if (isGenotypeSource(ldPath)) { - geno <- loadGenotypeRegion(ldPath) - mat <- computeLd(geno) - } else { - # Pre-computed .cor.xz block - snpFile <- if ("SNP_file" %in% colnames(ldInfo)) { - ldInfo$SNP_file[g] - } else { - NULL # let processLdMatrix auto-detect .bim/.pvar/.pvar.zst - } - ld <- processLdMatrix(ldPath, snpFile) - mat <- extractLdMatrix(ld) - } + } + # Center and scale genotype matrices + if (returnGenotype || nrow(mat) > ncol(mat)) { + mat <- scale(mat) + mat[is.na(mat)] <- 0 + } + mat +} - if (!is.null(maxVariants) && ncol(mat) > maxVariants) { - keep <- sort(sample(ncol(mat), maxVariants)) - mat <- mat[keep, keep] - } - mat +# @noRd +.ldLoadIdInfo <- function(g, ldInfo, maxVariants) { + ldPath <- ldInfo$LD_file[g] + + # Auto-detect format: genotype source or pre-computed block + if (isGenotypeSource(ldPath)) { + geno <- loadGenotypeRegion(ldPath) + mat <- computeLd(geno) + } else { + # Pre-computed .cor.xz block + snpFile <- if ("SNP_file" %in% colnames(ldInfo)) { + ldInfo$SNP_file[g] + } else { + NULL # let processLdMatrix auto-detect .bim/.pvar/.pvar.zst } + ld <- processLdMatrix(ldPath, snpFile) + mat <- extractLdMatrix(ld) } - loader + if (!is.null(maxVariants) && ncol(mat) > maxVariants) { + keep <- sort(sample(ncol(mat), maxVariants)) + mat <- mat[keep, keep] + } + mat +} + +# Dispatch a single block load by the spec's source mode. +# @noRd +.ldLoadBlock <- function(spec, g) { + switch(spec$mode, + rList = .ldLoadRList(g, spec$rList, spec$maxVariants), + xList = .ldLoadXList(g, spec$xList, spec$maxVariants), + meta = .ldLoadRegionMeta(g, spec$ldMetaPath, spec$regions, + spec$returnGenotype, spec$maxVariants), + info = .ldLoadIdInfo(g, spec$ldInfo, spec$maxVariants)) +} + +#' Load one LD block from an ldLoader spec +#' +#' Given an \code{ldLoaderSpec} (from \code{\link{ldLoader}}) and a block index +#' \code{g}, load the corresponding LD correlation matrix (or the genotype +#' matrix, in region mode with \code{returnGenotype = TRUE}). +#' +#' @param spec An \code{ldLoaderSpec} object returned by \code{\link{ldLoader}}. +#' @param g Integer block index (1-based). +#' @return The LD correlation matrix or genotype matrix for block \code{g}. +#' @seealso \code{\link{ldLoader}} +#' @examples +#' spec <- ldLoader(rList = list(diag(10), diag(15))) +#' loadLdBlock(spec, 1) +#' @export +loadLdBlock <- function(spec, g) { + if (!inherits(spec, "ldLoaderSpec")) + stop("`spec` must be an ldLoaderSpec (from ldLoader()).") + .ldLoadBlock(spec, g) } diff --git a/R/manifestLoaders.R b/R/manifestLoaders.R index e8192158..8d1f55d4 100644 --- a/R/manifestLoaders.R +++ b/R/manifestLoaders.R @@ -384,6 +384,28 @@ NULL } } +# Resolve one canonical sumstats field to its source column: honor an explicit +# columnMapping (under any accepted key spelling), else fall back to known column +# aliases; NA when unresolved. +# @noRd +.resolveSumstatKey <- function(key, df, mapping, label) { + # Look the field up in the mapping under any accepted standard-key spelling + # (`z`/`Z`, `n_sample`/`N`, ...). `intersect` avoids the "subscript out of + # bounds" a named character vector throws for an absent `[[` key. + if (!is.null(mapping)) { + cand <- intersect(.sumstatMappingKeys[[key]], names(mapping)) + if (length(cand) >= 1L) { + src <- mapping[[cand[[1L]]]] + if (!(src %in% names(df))) { + stop(label, ": columnMapping['", cand[[1L]], "'] = '", src, + "' is not a column in the sumstats file.") + } + return(src) + } + } + intersect(.sumstatColumnAliases[[key]], names(df))[1L] +} + # Standardise the columns of a raw sumstats data.frame into the canonical # schema expected by .dfToEntryGranges: chrom, pos, SNP, A1, A2, Z, N (+ # optional N_CASE, N_CONTROL, BETA, SE, P, MAF, INFO). N is required UNLESS @@ -399,25 +421,8 @@ NULL # read paths, and so explicit columnMapping values match the bare name. if (ncol(df) > 0L) names(df)[1L] <- sub("^#", "", names(df)[1L]) mapping <- .readColumnMapping(columnMapping) - resolveKey <- function(key) { - # Look the field up in the mapping under any accepted standard-key spelling - # (`z`/`Z`, `n_sample`/`N`, ...). `intersect` avoids the "subscript out of - # bounds" a named character vector throws for an absent `[[` key. - if (!is.null(mapping)) { - cand <- intersect(.sumstatMappingKeys[[key]], names(mapping)) - if (length(cand) >= 1L) { - src <- mapping[[cand[[1L]]]] - if (!(src %in% names(df))) { - stop(label, ": columnMapping['", cand[[1L]], "'] = '", src, - "' is not a column in the sumstats file.") - } - return(src) - } - } - intersect(.sumstatColumnAliases[[key]], names(df))[1L] - } required <- c("chrom", "pos", "variant_id", "A1", "A2") - resolved <- setNames(lapply(required, resolveKey), required) + resolved <- setNames(lapply(required, .resolveSumstatKey, df, mapping, label), required) missingKeys <- required[vapply(resolved, function(x) is.na(x) || is.null(x), logical(1))] if (length(missingKeys) > 0L) { stop(label, ": sumstats file is missing required field(s): ", @@ -427,9 +432,9 @@ NULL # z-score: a Z column, OR BETA + SE (from which the Wald z = beta/se is # derived, matching susie_rss()'s z_method="wald"). At least one is required; # when both a Z column and BETA/SE are present, the Z column takes precedence. - zSrc <- resolveKey("Z") - betaSrc <- resolveKey("BETA") - seSrc <- resolveKey("SE") + zSrc <- .resolveSumstatKey("Z", df, mapping, label) + betaSrc <- .resolveSumstatKey("BETA", df, mapping, label) + seSrc <- .resolveSumstatKey("SE", df, mapping, label) hasZ <- !is.na(zSrc) && !is.null(zSrc) hasBetaSe <- !is.na(betaSrc) && !is.null(betaSrc) && !is.na(seSrc) && !is.null(seSrc) @@ -442,9 +447,9 @@ NULL # which summaryStatsQc(effectiveN=) derives the effective sample size). At # least one is required; both may be present (N_CASE/N_CONTROL take priority # downstream when effectiveN is on). - nSrc <- resolveKey("N") - ncaseSrc <- resolveKey("N_CASE") - ncontrolSrc <- resolveKey("N_CONTROL") + nSrc <- .resolveSumstatKey("N", df, mapping, label) + ncaseSrc <- .resolveSumstatKey("N_CASE", df, mapping, label) + ncontrolSrc <- .resolveSumstatKey("N_CONTROL", df, mapping, label) hasN <- !is.na(nSrc) && !is.null(nSrc) hasCounts <- !is.na(ncaseSrc) && !is.null(ncaseSrc) && !is.na(ncontrolSrc) && !is.null(ncontrolSrc) @@ -472,7 +477,7 @@ NULL out$N_CONTROL <- as.numeric(df[[ncontrolSrc]]) } for (key in c("BETA", "SE", "P", "MAF", "INFO")) { - src <- resolveKey(key) + src <- .resolveSumstatKey(key, df, mapping, label) if (!is.na(src) && !is.null(src)) out[[key]] <- as.numeric(df[[src]]) } out @@ -499,6 +504,17 @@ NULL "); pass `sampleSelect` to choose the study column.") } +# Read column `tag` from the GWAS-VCF geno matrix for one sample as numeric, or +# NULL when the tag is absent/unmapped. +# @noRd +.vcfGetField <- function(tag, geno, sample) { + if (!is.null(tag) && tag %in% names(geno)) { + as.numeric(geno[[tag]][, sample]) + } else { + NULL + } +} + # Convert a VCF (GWAS-VCF format) into the canonical sumstats data.frame. The # effect allele (A1) is ALT; the other allele (A2) is REF. Stats come from the # per-study FORMAT fields ES/SE/LP/SS/EAF, with Z = ES / SE. @@ -511,18 +527,11 @@ NULL ref <- as.character(VariantAnnotation::ref(vcf)) geno <- VariantAnnotation::geno(vcf) sample <- .pickVcfSample(colnames(vcf), sampleSelect, label) - getField <- function(tag) { - if (!is.null(tag) && tag %in% names(geno)) { - as.numeric(geno[[tag]][, sample]) - } else { - NULL - } - } - es <- getField(fmap$BETA) - se <- getField(fmap$SE) - lp <- getField(fmap$P) - ss <- getField(fmap$N) - eaf <- getField(fmap$MAF) + es <- .vcfGetField(fmap$BETA, geno, sample) + se <- .vcfGetField(fmap$SE, geno, sample) + lp <- .vcfGetField(fmap$P, geno, sample) + ss <- .vcfGetField(fmap$N, geno, sample) + eaf <- .vcfGetField(fmap$MAF, geno, sample) if (is.null(es) || is.null(se)) { stop(label, ": GWAS-VCF must carry ES and SE FORMAT fields to derive Z ", "(configure via `formatMapping`).") @@ -843,6 +852,16 @@ loadQtlDatasetFromManifest <- function(manifest, study = NULL, genome = c("genome"), ldSketchPath = c("ldSketchPath", "ld_sketch_path", "ld_sketch")) +# TRUE if study `i` has usable case/control counts OR a total-N scalar, so a +# missing per-variant N is acceptable for that study. +# @noRd +.manifestHasStudyScalar <- function(i, nCaseCol, nControlCol, nSampleCol) { + ccOk <- !is.null(nCaseCol) && !is.null(nControlCol) && + is.finite(nCaseCol[[i]]) && is.finite(nControlCol[[i]]) + nOk <- !is.null(nSampleCol) && is.finite(nSampleCol[[i]]) + isTRUE(ccOk) || isTRUE(nOk) +} + #' @title Load a GwasSumStats collection from a manifest #' @description Build a \code{\link{GwasSumStats}} from a manifest with one row #' per study. No QC is run (the result carries \code{qcInfo = list()}). Each @@ -897,13 +916,6 @@ loadGwasSumStatsFromManifest <- function(manifest, genome = NULL, nCaseCol <- if ("nCase" %in% names(df)) as.numeric(df$nCase) else NULL nControlCol <- if ("nControl" %in% names(df)) as.numeric(df$nControl) else NULL nSampleCol <- if ("nSample" %in% names(df)) as.numeric(df$nSample) else NULL - hasStudyScalar <- function(i) { - ccOk <- !is.null(nCaseCol) && !is.null(nControlCol) && - is.finite(nCaseCol[[i]]) && is.finite(nControlCol[[i]]) - nOk <- !is.null(nSampleCol) && is.finite(nSampleCol[[i]]) - isTRUE(ccOk) || isTRUE(nOk) - } - entries <- lapply(seq_len(nrow(df)), function(i) { label <- paste0("GwasSumStats[study=", df$study[[i]], "]") mapping <- if ("columnMapping" %in% names(df) && @@ -915,11 +927,15 @@ loadGwasSumStatsFromManifest <- function(manifest, genome = NULL, } gr <- .loadSumStatsEntry(.resolveRel(as.character(df$sumStatsPath[[i]]), base), region, mapping, sampleSelect, formatMapping, label, - allowNoN = hasStudyScalar(i)) + allowNoN = .manifestHasStudyScalar(i, nCaseCol, nControlCol, nSampleCol)) .checkLdContainment(ldSketch, gr, minLdOverlapWarn, label) gr }) + # Trim a genome-wide LD sketch to the summary stats' per-chromosome position + # span so the object doesn't carry a full-genome snpInfo. + ldSketch <- .subsetSketchToRange(ldSketch, entries) + args <- list(study = as.character(df$study), entry = entries, genome = genome, ldSketch = ldSketch) if (!is.null(nCaseCol)) args$nCase <- nCaseCol @@ -1004,6 +1020,9 @@ loadQtlSumStatsFromManifest <- function(manifest, genome = NULL, minLdOverlapWarn, columnMapping, sampleSelect, formatMapping, allowNoN = allowNoN) + # Trim a genome-wide LD sketch to the summary stats' per-chromosome position + # span so the object doesn't carry a full-genome snpInfo. + ldSketch <- .subsetSketchToRange(ldSketch, entries) args <- list(study = as.character(df$study), context = as.character(df$context), trait = as.character(df$trait), diff --git a/R/mashWrapper.R b/R/mashWrapper.R index 30bc29c5..accbbff2 100644 --- a/R/mashWrapper.R +++ b/R/mashWrapper.R @@ -6,20 +6,37 @@ filterBySignificance <- function(zMatrix, sigPCutoff) { which(apply(zMatrix, 1, function(row) any(abs(row) >= zThreshold))) } +# Coerce every column to numeric, then replace NaN/Inf/NA with `replaceWith`. +# @noRd +.mashReplaceValues <- function(df, replaceWith) { + df <- df %>% + mutate(across(everything(), as.numeric)) %>% + mutate(across(everything(), ~ replace(., is.nan(.) | is.infinite(.) | is.na(.), replaceWith))) +} + +# Coerce z-scores to a matrix (NaN/Inf/NA -> 0) and, when a missing-rate +# threshold is given, drop rows falling below it. +# @noRd +.mashProcessZ <- function(zData, filterByMissingRate) { + zData <- as.matrix(.mashReplaceValues(zData, 0)) + + if (!is.null(filterByMissingRate)) { + proportionNonzero <- apply(zData, 1, function(row) mean(row != 0)) + zData <- zData[proportionNonzero >= filterByMissingRate, , drop = FALSE] + } + + return(zData) +} + #' @importFrom vroom vroom #' @export filterInvalidSummaryStat <- function(datList, bhat = NULL, sbhat = NULL, z = NULL, btoz = FALSE, sigPCutoff = 1E-6, filterByMissingRate = 0.2) { - replaceValues <- function(df, replaceWith) { - df <- df %>% - mutate(across(everything(), as.numeric)) %>% - mutate(across(everything(), ~ replace(., is.nan(.) | is.infinite(.) | is.na(.), replaceWith))) - } # Function to process bhat, sbhat if (!is.null(bhat) && !is.null(sbhat) && all(c(bhat, sbhat) %in% names(datList))) { # If the element is a list with 'bhat' and 'sbhat' if (!is.null(datList[[bhat]]) && !is.null(datList[[sbhat]])) { - datList[[bhat]] <- as.matrix(replaceValues(datList[[bhat]], 0)) - datList[[sbhat]] <- as.matrix(replaceValues(datList[[sbhat]], 1000)) + datList[[bhat]] <- as.matrix(.mashReplaceValues(datList[[bhat]], 0)) + datList[[sbhat]] <- as.matrix(.mashReplaceValues(datList[[sbhat]], 1000)) if (("null.b" %in% names(datList)) || ("random.b" %in% names(datList))) { if (!is.null(filterByMissingRate)) { proportionNonzero <- apply(datList[[bhat]], 1, function(row) { @@ -58,21 +75,10 @@ filterInvalidSummaryStat <- function(datList, bhat = NULL, sbhat = NULL, z = NUL } # Function to process z-scores and filter directly if (!is.null(z)) { - processZ <- function(zData) { - zData <- as.matrix(replaceValues(zData, 0)) - - if (!is.null(filterByMissingRate)) { - proportionNonzero <- apply(zData, 1, function(row) mean(row != 0)) - zData <- zData[proportionNonzero >= filterByMissingRate, , drop = FALSE] - } - - return(zData) - } - # Process each component if it exists for (comp in c("strong", "random", "null")) { if (!is.null(datList[[comp]]) && !is.null(datList[[comp]]$z)) { - datList[[comp]]$z <- processZ(datList[[comp]]$z) + datList[[comp]]$z <- .mashProcessZ(datList[[comp]]$z, filterByMissingRate) } } @@ -137,58 +143,59 @@ filterMixtureComponents <- function(conditionsToKeep, U, w = NULL, wCutoff = 1e- } -#' @export -mashRandNullSample <- function(dat, nRandom, nNull, excludeCondition, seed = NULL) { - # Function to extract one data set - extractOneData <- function(dat, nRandom, nNull) { - if (is.null(dat)) { - return(NULL) - } +# Draw the random + null sub-samples used to estimate the null correlation. +# @noRd +.mashExtractOneData <- function(dat, nRandom, nNull) { + if (is.null(dat)) { + return(NULL) + } - if ("z" %in% names(dat)) { - absZ <- abs(dat$z) - zData <- dat$z - } else { - absZ <- abs(dat$bhat / dat$sbhat) - zData <- NULL - } + if ("z" %in% names(dat)) { + absZ <- abs(dat$z) + zData <- dat$z + } else { + absZ <- abs(dat$bhat / dat$sbhat) + zData <- NULL + } - sampleIdx <- 1:nrow(absZ) - randomIdx <- sample(sampleIdx, min(nRandom, length(sampleIdx)), replace = FALSE) + sampleIdx <- 1:nrow(absZ) + randomIdx <- sample(sampleIdx, min(nRandom, length(sampleIdx)), replace = FALSE) - if (!is.null(zData)) { - random <- list(z = zData[randomIdx, , drop = FALSE]) - } else { - random <- list( - bhat = dat$bhat[randomIdx, , drop = FALSE], - sbhat = dat$sbhat[randomIdx, , drop = FALSE] - ) - } + if (!is.null(zData)) { + random <- list(z = zData[randomIdx, , drop = FALSE]) + } else { + random <- list( + bhat = dat$bhat[randomIdx, , drop = FALSE], + sbhat = dat$sbhat[randomIdx, , drop = FALSE] + ) + } - null.id <- which(apply(absZ, 1, max) < 2) - if (length(null.id) == 0) { - warning(paste("no variants are included in the null dataset because absZ > 2 for all variants in", dat$region)) + null.id <- which(apply(absZ, 1, max) < 2) + if (length(null.id) == 0) { + warning(paste("no variants are included in the null dataset because absZ > 2 for all variants in", dat$region)) + null <- list() + } else { + if (length(null.id) < ncol(absZ)) { + warning(paste("not enough null data to estimate null correlation in", dat$region)) null <- list() } else { - if (length(null.id) < ncol(absZ)) { - warning(paste("not enough null data to estimate null correlation in", dat$region)) - null <- list() + nullIdx <- sample(null.id, min(nNull, length(null.id)), replace = FALSE) + if (!is.null(zData)) { + null <- list(z = zData[nullIdx, , drop = FALSE]) } else { - nullIdx <- sample(null.id, min(nNull, length(null.id)), replace = FALSE) - if (!is.null(zData)) { - null <- list(z = zData[nullIdx, , drop = FALSE]) - } else { - null <- list( - bhat = dat$bhat[nullIdx, , drop = FALSE], - sbhat = dat$sbhat[nullIdx, , drop = FALSE] - ) - } + null <- list( + bhat = dat$bhat[nullIdx, , drop = FALSE], + sbhat = dat$sbhat[nullIdx, , drop = FALSE] + ) } } - dat <- list(random = random, null = null) - return(dat) } + dat <- list(random = random, null = null) + return(dat) +} +#' @export +mashRandNullSample <- function(dat, nRandom, nNull, excludeCondition, seed = NULL) { if (!is.null(seed)) { set.seed(seed) } @@ -204,7 +211,7 @@ mashRandNullSample <- function(dat, nRandom, nNull, excludeCondition, seed = NUL } } - result <- extractOneData(dat, nRandom, nNull) + result <- .mashExtractOneData(dat, nRandom, nNull) return(result) } diff --git a/R/overlapTopLoci.R b/R/overlapTopLoci.R index 4c18326b..70279c29 100644 --- a/R/overlapTopLoci.R +++ b/R/overlapTopLoci.R @@ -24,6 +24,32 @@ #' @seealso \code{\link{getTopLoci}}, \code{\link{matchVariants}} #' @include AllGenerics.R AllClasses.R QtlFineMappingResult.R GwasFineMappingResult.R #' @export +# Prefix a top-loci frame's NON-key columns with `pfx` (key columns unchanged). +# @noRd +.overlapPrefixNonKey <- function(df, pfx, keyCols) + stats::setNames(df, ifelse(names(df) %in% keyCols, names(df), + paste0(pfx, names(df)))) + +# The zero-row merged frame returned when either side has no signal: an +# identity-only join on variant_id (GWAS contributes only its non-coord cols). +# @noRd +.overlapEmptyMerge <- function(qtlTl, gwasTl, coordCols, keyCols) { + qp <- .overlapPrefixNonKey(qtlTl[0, , drop = FALSE], "qtl_", keyCols) + gp <- .overlapPrefixNonKey( + gwasTl[0, setdiff(names(gwasTl), coordCols), drop = FALSE], "gwas_", keyCols) + # A collection with no signal above the cutoff yields an identity-only + # getTopLoci frame (study/context/trait/method) that carries no + # variant_id; restore the join key so the empty merge still resolves + # instead of erroring in merge()'s `by` check. + if (!"variant_id" %in% names(qp)) qp$variant_id <- character(0) + if (!"variant_id" %in% names(gp)) gp$variant_id <- character(0) + merge(qp, gp, by = "variant_id") +} + +# Convert the merged frame to the requested output type. +# @noRd +.overlapFinish <- function(df, type) if (type == "GRanges") .overlapToGRanges(df) else df + setGeneric("overlapTopLoci", function(qtl, gwas, ...) standardGeneric("overlapTopLoci")) @@ -40,33 +66,16 @@ setMethod("overlapTopLoci", qtlTl <- as.data.frame(getTopLoci(qtl, signalCutoff = signalCutoff)) gwasTl <- as.data.frame(getTopLoci(gwas, signalCutoff = signalCutoff)) - prefixNonKey <- function(df, pfx) - stats::setNames(df, ifelse(names(df) %in% keyCols, names(df), - paste0(pfx, names(df)))) - # The GWAS side contributes only variant_id (the join key) + its non-key - # columns; the shared coordinate columns come from the QTL side. - emptyMerge <- function() { - qp <- prefixNonKey(qtlTl[0, , drop = FALSE], "qtl_") - gp <- prefixNonKey( - gwasTl[0, setdiff(names(gwasTl), coordCols), drop = FALSE], "gwas_") - # A collection with no signal above the cutoff yields an identity-only - # getTopLoci frame (study/context/trait/method) that carries no - # variant_id; restore the join key so the empty merge still resolves - # instead of erroring in merge()'s `by` check. - if (!"variant_id" %in% names(qp)) qp$variant_id <- character(0) - if (!"variant_id" %in% names(gp)) gp$variant_id <- character(0) - merge(qp, gp, by = "variant_id") - } - finish <- function(df) if (type == "GRanges") .overlapToGRanges(df) else df - - if (nrow(qtlTl) == 0L || nrow(gwasTl) == 0L) return(finish(emptyMerge())) + if (nrow(qtlTl) == 0L || nrow(gwasTl) == 0L) + return(.overlapFinish(.overlapEmptyMerge(qtlTl, gwasTl, coordCols, keyCols), type)) # Allele-aware correspondence between the unique QTL and GWAS variant sets. # target = GWAS, ref = QTL, so `sign` is the flip applied to the GWAS side. uq <- unique(qtlTl$variant_id) ug <- unique(gwasTl$variant_id) m <- matchVariants(ug, uq, allowFlip = TRUE) - if (length(m$idxA) == 0L) return(finish(emptyMerge())) + if (length(m$idxA) == 0L) + return(.overlapFinish(.overlapEmptyMerge(qtlTl, gwasTl, coordCols, keyCols), type)) map <- data.frame(gwas_vid = ug[m$idxA], canon_vid = uq[m$idxB], .sign = as.numeric(m$sign), stringsAsFactors = FALSE) @@ -75,7 +84,8 @@ setMethod("overlapTopLoci", g <- merge(gwasTl, map, by.x = "variant_id", by.y = "gwas_vid") # Unreachable defensive guard: map$gwas_vid is drawn from gwasTl$variant_id, # so this inner merge can never drop every row (length(m$idxA) > 0 here). - if (nrow(g) == 0L) return(finish(emptyMerge())) # nocov + if (nrow(g) == 0L) + return(.overlapFinish(.overlapEmptyMerge(qtlTl, gwasTl, coordCols, keyCols), type)) # nocov for (cc in intersect(c("beta", "z", "conditional_effect"), names(g))) g[[cc]] <- g[[cc]] * g$.sign if ("af" %in% names(g)) @@ -84,12 +94,12 @@ setMethod("overlapTopLoci", g <- g[, setdiff(names(g), c(coordCols, "canon_vid", ".sign")), drop = FALSE] # Wide cross-product per shared variant, variant key kept once (QTL side). - merged <- merge(prefixNonKey(qtlTl, "qtl_"), prefixNonKey(g, "gwas_"), - by = "variant_id") + merged <- merge(.overlapPrefixNonKey(qtlTl, "qtl_", keyCols), + .overlapPrefixNonKey(g, "gwas_", keyCols), by = "variant_id") # Restore key-column order (merge puts the `by` column first). ordered <- c(intersect(keyCols, names(merged)), setdiff(names(merged), keyCols)) - finish(merged[, ordered, drop = FALSE]) + .overlapFinish(merged[, ordered, drop = FALSE], type) }) # Build a GRanges from an overlap table: variants as width-1 ranges, all other diff --git a/R/qtlEnrichmentPipeline.R b/R/qtlEnrichmentPipeline.R index adf2e40e..357060f5 100644 --- a/R/qtlEnrichmentPipeline.R +++ b/R/qtlEnrichmentPipeline.R @@ -114,23 +114,8 @@ qtlEnrichmentPipeline <- function(gwasFineMappingResult, # empty union means no GWAS study has usable PIPs, so every study is skipped # before this is ever reached (and the length guard skips relabeling entirely # when the union GWAS panel is empty). - alignedByTuple <- vector("list", nrow(qtlTuples)) - alignTuple <- function(k) { - if (!is.null(alignedByTuple[[k]])) return(alignedByTuple[[k]]) - aligned <- lapply(qtlRegionsByTuple[[k]], function(x) { - if (!is.null(names(x$pip)) && length(unionGwasNames) > 0L) { - # Relabel matched pip names to the GWAS convention via the shared - # matcher (tuple match; unmatched names kept as-is). - mm <- matchVariants(names(x$pip), unionGwasNames) - nm <- names(x$pip) - nm[mm$idxA] <- unionGwasNames[mm$idxB] - names(x$pip) <- nm - } - x - }) - alignedByTuple[[k]] <<- aligned - aligned - } + alignCache <- new.env(parent = emptyenv()) + alignCache$byTuple <- vector("list", nrow(qtlTuples)) results <- list() for (gi in seq_along(gwasStudies)) { @@ -154,7 +139,7 @@ qtlEnrichmentPipeline <- function(gwasFineMappingResult, enr <- tryCatch( qtlEnrichment( gwasPip = gwasPip, - susieQtlRegions = alignTuple(k), + susieQtlRegions = .enrAlignTuple(k, alignCache, qtlRegionsByTuple, unionGwasNames), numGwas = numGwas, piQtl = piQtl, lambda = lambda, @@ -199,6 +184,26 @@ qtlEnrichmentPipeline <- function(gwasFineMappingResult, # Internal helpers # ============================================================================= +# Align (memoized) the QTL regions of tuple `k` to the GWAS naming convention, +# caching the result in the `cache` environment's `byTuple` list. +# @noRd +.enrAlignTuple <- function(k, cache, qtlRegionsByTuple, unionGwasNames) { + if (!is.null(cache$byTuple[[k]])) return(cache$byTuple[[k]]) + aligned <- lapply(qtlRegionsByTuple[[k]], function(x) { + if (!is.null(names(x$pip)) && length(unionGwasNames) > 0L) { + # Relabel matched pip names to the GWAS convention via the shared + # matcher (tuple match; unmatched names kept as-is). + mm <- matchVariants(names(x$pip), unionGwasNames) + nm <- names(x$pip) + nm[mm$idxA] <- unionGwasNames[mm$idxB] + names(x$pip) <- nm + } + x + }) + cache$byTuple[[k]] <- aligned + aligned +} + # Build a named GWAS PIP vector for one study. Walks every row of the # GwasFineMappingResult tagged with that study, extracts the per-row # pip from each FineMappingEntry, and concatenates with variant-id @@ -269,6 +274,15 @@ qtlEnrichmentPipeline <- function(gwasFineMappingResult, out } +# Pull one enrichment field from qtlEnrichment's list output as a scalar numeric +# (NA when absent). +# @noRd +.enrPickScalar <- function(field, enr) { + v <- enr[[field]] + if (is.null(v)) NA_real_ + else as.numeric(v[[1L]]) +} + # Coerce qtlEnrichment's variable-shape output into a single-row # named list with the canonical columns the caller documents. The # underlying estimator returns either a list with named numeric scalars @@ -277,15 +291,10 @@ qtlEnrichmentPipeline <- function(gwasFineMappingResult, # @noRd .enrFlattenEnrichment <- function(enr) { if (is.list(enr) && is.null(dim(enr))) { - pickScalar <- function(field) { - v <- enr[[field]] - if (is.null(v)) NA_real_ - else as.numeric(v[[1L]]) - } list( - enrichment = pickScalar("enrichment"), - enrichmentSe = pickScalar("enrichmentSe"), - enrichmentLogOdds = pickScalar("enrichmentLogOdds")) + enrichment = .enrPickScalar("enrichment", enr), + enrichmentSe = .enrPickScalar("enrichmentSe", enr), + enrichmentLogOdds = .enrPickScalar("enrichmentLogOdds", enr)) } else if (is.matrix(enr) || is.data.frame(enr)) { df <- as.data.frame(enr, stringsAsFactors = FALSE) if (nrow(df) == 0L) { diff --git a/R/regularizedRegressionWrappers.R b/R/regularizedRegressionWrappers.R index 2a142784..1862064e 100644 --- a/R/regularizedRegressionWrappers.R +++ b/R/regularizedRegressionWrappers.R @@ -1006,14 +1006,71 @@ lassosumRss <- function(bhat, R, n, result } +# Per-`s` fit for one RSS method (`method` selects the solver + which `config` +# fields apply). Returns list(beta, meta); l0learn's inner lambda0 sweep is here. +# @noRd +.rssFitOne <- function(method, solverInput, LDs, n, sVal, config) { + switch(method, + lassosum = { + model <- do.call(lassosumRss, + c(list(bhat = solverInput, R = LDs, n = n), config$dotArgs)) + list(beta = model$beta, + meta = data.frame(s = rep(sVal, length(model$lambda)), + lambda = model$lambda, fbeta = model$fbeta, + stringsAsFactors = FALSE)) + }, + penalized = { + model <- do.call(penalizedRss, + c(list(bhat = solverInput, R = LDs, n = n, + penalty = config$penalty, gamma = config$gamma, + alpha = config$alpha, lambda0 = config$lambda0, + lambda2 = config$lambda2), config$dotArgs)) + list(beta = model$beta, + meta = data.frame(s = rep(sVal, length(model$lambda)), + lambda = model$lambda, fbeta = model$fbeta, + stringsAsFactors = FALSE)) + }, + l0learn = { + beta <- NULL + meta <- list() + for (l0Val in config$lambda0) { + model <- do.call(penalizedRss, + c(list(bhat = solverInput, R = LDs, n = n, + penalty = config$penalty, lambda = config$lambda, + lambda0 = l0Val, lambda2 = config$lambda2, + maxSwaps = config$maxSwaps), config$dotArgs)) + beta <- cbind(beta, model$beta) + meta[[length(meta) + 1L]] <- data.frame( + s = rep(sVal, length(model$lambda)), + lambda0 = rep(l0Val, length(model$lambda)), + lambda = model$lambda, fbeta = model$fbeta, + stringsAsFactors = FALSE) + } + list(beta = beta, meta = do.call(rbind, meta)) + }) +} + +# Stamp the method-specific selection attribute onto the chosen coefficient vector. +# @noRd +.rssFinalize <- function(method, bestBeta, sel, meta, config) { + base <- c(mode = sel$mode, index = sel$index) + attr(bestBeta, if (method == "lassosum") "lassosum_selection" + else "penalized_rss_selection") <- switch(method, + lassosum = c(base, s = meta$s[sel$index], lambda = meta$lambda[sel$index]), + penalized = c(base, penalty = config$penalty, + s = meta$s[sel$index], lambda = meta$lambda[sel$index]), + l0learn = c(base, penalty = config$penalty, + s = meta$s[sel$index], lambda0 = meta$lambda0[sel$index], + lambda = meta$lambda[sel$index])) + bestBeta +} + # Shared scaffold for the RSS shrinkage-grid weight functions # (lassosumRssWeights / .penalizedRssWeights / l0learnRssWeights). Standardizes # the stat -> solverInput conversion, the outer LD-shrinkage grid over `s`, the -# candidate accumulation, and the ld_quadratic / min_fbeta selection. `fitOne` -# supplies the per-`s` fit as list(beta, meta); any inner grid (e.g. l0learn's -# lambda0 sweep) lives inside it. `finalize` stamps the function-specific -# selection attribute onto the returned coefficient vector. -.rssShrinkGridWeights <- function(stat, LD, s, fitOne, finalize, +# candidate accumulation, and the ld_quadratic / min_fbeta selection. `method` + +# `config` pick the per-`s` solver (.rssFitOne) and the finalizer (.rssFinalize). +.rssShrinkGridWeights <- function(stat, LD, s, method, config, selection = c("ld_quadratic", "min_fbeta")) { selection <- match.arg(selection) n <- median(stat$n) @@ -1024,7 +1081,7 @@ lassosumRss <- function(bhat, R, n, candidateMeta <- list() for (sVal in s) { LDs <- (1 - sVal) * LD + sVal * diag(p) - one <- fitOne(solverInput, LDs, n, sVal) + one <- .rssFitOne(method, solverInput, LDs, n, sVal, config) candidateBeta <- cbind(candidateBeta, one$beta) candidateMeta[[length(candidateMeta) + 1L]] <- one$meta } @@ -1035,7 +1092,7 @@ lassosumRss <- function(bhat, R, n, .lassosumSelectMinFbeta(candidateBeta, candidateMeta) } bestBeta <- as.numeric(selectorResult$beta) - finalize(bestBeta, selectorResult, candidateMeta) + .rssFinalize(method, bestBeta, selectorResult, candidateMeta, config) } #' Extract weights from lassosumRss with shrinkage grid search @@ -1073,23 +1130,8 @@ lassosumRssWeights <- function(stat, LD, s = c(0.2, 0.5, 0.9, 1.0), selection = c("ld_quadratic", "min_fbeta"), ...) { selection <- match.arg(selection) - dotArgs <- list(...) - fitOne <- function(solverInput, LDs, n, sVal) { - model <- do.call(lassosumRss, - c(list(bhat = solverInput, R = LDs, n = n), - dotArgs)) - list(beta = model$beta, - meta = data.frame(s = rep(sVal, length(model$lambda)), - lambda = model$lambda, fbeta = model$fbeta, - stringsAsFactors = FALSE)) - } - finalize <- function(bestBeta, sel, meta) { - attr(bestBeta, "lassosum_selection") <- c( - mode = sel$mode, index = sel$index, - s = meta$s[sel$index], lambda = meta$lambda[sel$index]) - bestBeta - } - .rssShrinkGridWeights(stat, LD, s, fitOne, finalize, selection) + .rssShrinkGridWeights(stat, LD, s, "lassosum", + list(dotArgs = list(...)), selection) } #' Penalized Regression on RSS (Summary Statistics) Objective @@ -1191,24 +1233,10 @@ penalizedRss <- function(bhat, R, n, selection = c("ld_quadratic", "min_fbeta"), ...) { selection <- match.arg(selection) - dotArgs <- list(...) - fitOne <- function(solverInput, LDs, n, sVal) { - model <- do.call(penalizedRss, - c(list(bhat = solverInput, R = LDs, n = n, - penalty = penalty, gamma = gamma, alpha = alpha, - lambda0 = lambda0, lambda2 = lambda2), dotArgs)) - list(beta = model$beta, - meta = data.frame(s = rep(sVal, length(model$lambda)), - lambda = model$lambda, fbeta = model$fbeta, - stringsAsFactors = FALSE)) - } - finalize <- function(bestBeta, sel, meta) { - attr(bestBeta, "penalized_rss_selection") <- c( - mode = sel$mode, index = sel$index, penalty = penalty, - s = meta$s[sel$index], lambda = meta$lambda[sel$index]) - bestBeta - } - .rssShrinkGridWeights(stat, LD, s, fitOne, finalize, selection) + .rssShrinkGridWeights(stat, LD, s, "penalized", + list(penalty = penalty, gamma = gamma, alpha = alpha, + lambda0 = lambda0, lambda2 = lambda2, + dotArgs = list(...)), selection) } #' Compute SCAD-Penalized Weights from Summary Statistics @@ -1307,35 +1335,10 @@ l0learnRssWeights <- function(stat, LD, } } - dotArgs <- list(...) - # Inner sweep over the lambda0 (L0) path for each shrinkage level; the outer - # `s` grid, selection, and bestBeta live in .rssShrinkGridWeights. - fitOne <- function(solverInput, LDs, n, sVal) { - beta <- NULL - meta <- list() - for (l0Val in lambda0) { - model <- do.call(penalizedRss, - c(list(bhat = solverInput, R = LDs, n = n, - penalty = penalty, lambda = lambda, - lambda0 = l0Val, lambda2 = lambda2, - maxSwaps = maxSwaps), dotArgs)) - beta <- cbind(beta, model$beta) - meta[[length(meta) + 1L]] <- data.frame( - s = rep(sVal, length(model$lambda)), - lambda0 = rep(l0Val, length(model$lambda)), - lambda = model$lambda, fbeta = model$fbeta, - stringsAsFactors = FALSE) - } - list(beta = beta, meta = do.call(rbind, meta)) - } - finalize <- function(bestBeta, sel, meta) { - attr(bestBeta, "penalized_rss_selection") <- c( - mode = sel$mode, index = sel$index, penalty = penalty, - s = meta$s[sel$index], lambda0 = meta$lambda0[sel$index], - lambda = meta$lambda[sel$index]) - bestBeta - } - .rssShrinkGridWeights(stat, LD, s, fitOne, finalize, selection) + .rssShrinkGridWeights(stat, LD, s, "l0learn", + list(penalty = penalty, lambda = lambda, lambda0 = lambda0, + lambda2 = lambda2, maxSwaps = maxSwaps, + dotArgs = list(...)), selection) } #' Compute Weights Using ncvreg with SCAD or MCP Penalty @@ -1823,34 +1826,37 @@ computeCoefficientsGlasso <- function(X, Y, standardize, nthreads, Xnew = NULL) } -### Function to compute coefficients for univariate glmnet -computeCoefficientsUnivGlmnet <- function(X, Y, alpha, standardize, nthreads, Xnew = NULL) { - r <- ncol(Y) +# Fit a cv.glmnet for outcome column `i` on its non-missing rows and return the +# lambda.min coefficients (plus predictions on `Xnew` when supplied). +# @noRd +.linreg <- function(i, X, Y, alpha, standardize, nthreads, Xnew) { + samplesKept <- which(!is.na(Y[, i])) + Ynomiss <- Y[samplesKept, i, drop = FALSE] + Xnomiss <- X[samplesKept, , drop = FALSE] - linreg <- function(i, X, Y, alpha, standardize, nthreads, Xnew) { - samplesKept <- which(!is.na(Y[, i])) - Ynomiss <- Y[samplesKept, i, drop = FALSE] - Xnomiss <- X[samplesKept, , drop = FALSE] + cvfit <- glmnet::cv.glmnet( + x = Xnomiss, y = Ynomiss, family = "gaussian", alpha = alpha, + standardize = standardize, parallel = FALSE + ) + coeffic <- as.vector(coef(cvfit, s = "lambda.min")) + lambdaSeq <- cvfit$lambda - cvfit <- glmnet::cv.glmnet( - x = Xnomiss, y = Ynomiss, family = "gaussian", alpha = alpha, - standardize = standardize, parallel = FALSE - ) - coeffic <- as.vector(coef(cvfit, s = "lambda.min")) - lambdaSeq <- cvfit$lambda + # Make predictions if requested + if (!is.null(Xnew)) { + yhatGlmnet <- drop(predict(cvfit, newx = Xnew, s = "lambda.min")) + res <- list(bhat = coeffic, lambda_seq = lambdaSeq, yhat_new = yhatGlmnet) + } else { + res <- list(bhat = coeffic, lambda_seq = lambdaSeq) + } - # Make predictions if requested - if (!is.null(Xnew)) { - yhatGlmnet <- drop(predict(cvfit, newx = Xnew, s = "lambda.min")) - res <- list(bhat = coeffic, lambda_seq = lambdaSeq, yhat_new = yhatGlmnet) - } else { - res <- list(bhat = coeffic, lambda_seq = lambdaSeq) - } + return(res) +} - return(res) - } +### Function to compute coefficients for univariate glmnet +computeCoefficientsUnivGlmnet <- function(X, Y, alpha, standardize, nthreads, Xnew = NULL) { + r <- ncol(Y) - out <- lapply(1:r, linreg, X, Y, alpha, standardize, nthreads, Xnew) + out <- lapply(1:r, .linreg, X, Y, alpha, standardize, nthreads, Xnew) Bhat <- sapply(out, "[[", "bhat") diff --git a/R/relatednessQc.R b/R/relatednessQc.R index b53afccd..d59d1f3f 100644 --- a/R/relatednessQc.R +++ b/R/relatednessQc.R @@ -108,26 +108,22 @@ filterRelatedness <- function( !(relatedness[[relatednessIid2]] %in% highRelatedIndiv), ] # --- Phase 2: plinkQC-based filtering ---- - runPlinkqc <- function(relDf) { - plinkQC::relatednessFilter( - relatedness = relDf, - otherCriterion = otherCriterion, - relatednessTh = relatednessThreshold, - relatednessIID1 = relatednessIid1, - relatednessIID2 = relatednessIid2, - otherCriterionTh = otherCriterionThreshold, - otherCriterionThDirection = otherCriterionDirection, - relatednessFID1 = relatednessFid1, - relatednessFID2 = relatednessFid2, - relatednessRelatedness = relatednessValue, - otherCriterionIID = otherCriterionIid, - otherCriterionMeasure = otherCriterionMeasure, - verbose = verbose - )$failIDs - } + plinkqcArgs <- list( + otherCriterion = otherCriterion, + relatednessTh = relatednessThreshold, + relatednessIID1 = relatednessIid1, + relatednessIID2 = relatednessIid2, + otherCriterionTh = otherCriterionThreshold, + otherCriterionThDirection = otherCriterionDirection, + relatednessFID1 = relatednessFid1, + relatednessFID2 = relatednessFid2, + relatednessRelatedness = relatednessValue, + otherCriterionIID = otherCriterionIid, + otherCriterionMeasure = otherCriterionMeasure, + verbose = verbose) if (analysisType == "maximize_unrelated") { - rel <- runPlinkqc(kin) + rel <- .relatednessRunPlinkqc(kin, plinkqcArgs) allExclude <- rel$IID } else { @@ -149,7 +145,7 @@ filterRelatedness <- function( caseKin <- kin[ kin[[relatednessIid1]] %in% relatedCases & kin[[relatednessIid2]] %in% relatedCases, ] - relCases <- runPlinkqc(caseKin) + relCases <- .relatednessRunPlinkqc(caseKin, plinkqcArgs) casesKeep <- setdiff(relatedCases, relCases$IID) # Step 2: Remove controls related to retained cases @@ -169,7 +165,7 @@ filterRelatedness <- function( controlKin <- kin[ kin[[relatednessIid1]] %in% controlsKeep & kin[[relatednessIid2]] %in% controlsKeep, ] - relControls <- runPlinkqc(controlKin) + relControls <- .relatednessRunPlinkqc(controlKin, plinkqcArgs) allExclude <- c(relCases$IID, controlsExclude, relControls$IID) } @@ -184,7 +180,7 @@ filterRelatedness <- function( while (nrow(remaining) > 0 && iter < maxIterations) { if (verbose) message("Iteration ", iter + 1L, ": ", nrow(remaining), " related pairs remaining.") - additional <- runPlinkqc(remaining) + additional <- .relatednessRunPlinkqc(remaining, plinkqcArgs) allExclude <- c(allExclude, additional$IID) remaining <- kin[ !(kin[[relatednessIid1]] %in% allExclude) & @@ -206,3 +202,11 @@ filterRelatedness <- function( allExclude } + +# Run plinkQC::relatednessFilter with the pre-bound column names + thresholds +# (`args`), returning its $failIDs. +# @noRd +.relatednessRunPlinkqc <- function(relDf, args) { + do.call(plinkQC::relatednessFilter, + c(list(relatedness = relDf), args))$failIDs +} diff --git a/R/sldscPostprocessingPipeline.R b/R/sldscPostprocessingPipeline.R index ceafdc36..afaea359 100644 --- a/R/sldscPostprocessingPipeline.R +++ b/R/sldscPostprocessingPipeline.R @@ -182,33 +182,10 @@ sldscPostprocessingPipeline <- function(sldscData, ptViewSingle <- .sldscViewForMeta(perTrait, "single") ptViewJoint <- .sldscViewForMeta(perTrait, "joint") - buildTable <- function(quantity, view, label) { - rows <- list() - for (cat in targetCategories) { - m <- metaSldscRandom(view, cat, quantity) - rows[[cat]] <- data.frame( - target = cat, - isBinary = unname(isBinary[cat]), - mean = m$mean, - se = m$se, - p = m$p, - nTraits = m$nTraits, - stringsAsFactors = FALSE - ) - } - df <- do.call(rbind, rows) - rownames(df) <- NULL - nmOld <- c("mean", "se", "p") - nmNew <- paste0(label, toupper(substring(nmOld, 1, 1)), - substring(nmOld, 2)) - names(df)[names(df) %in% nmOld] <- nmNew - df - } - - metaTauStarSingle <- buildTable("tauStar", ptViewSingle, "single") - metaTauStarJoint <- buildTable("tauStar", ptViewJoint, "joint") - metaESingle <- buildTable("enrichment", ptViewSingle, "single") - metaEsSingle <- buildTable("enrichstat", ptViewSingle, "single") + metaTauStarSingle <- .sldscBuildTable("tauStar", ptViewSingle, "single", isBinary, targetCategories) + metaTauStarJoint <- .sldscBuildTable("tauStar", ptViewJoint, "joint", isBinary, targetCategories) + metaESingle <- .sldscBuildTable("enrichment", ptViewSingle, "single", isBinary, targetCategories) + metaEsSingle <- .sldscBuildTable("enrichstat", ptViewSingle, "single", isBinary, targetCategories) # Combine tauStar single + joint into one wide frame. metaTauStar <- metaTauStarSingle @@ -252,21 +229,20 @@ sldscPostprocessingPipeline <- function(sldscData, length(targetLabels), length(targetCategories), paste(targetCategories, collapse = ", "))) relab <- setNames(targetLabels, targetCategories) - relabVec <- function(x) { y <- unname(relab[x]); y[is.na(y)] <- x[is.na(y)]; y } for (t in names(res$per_trait)) { pt <- res$per_trait[[t]] if (!is.null(pt$summary) && "target" %in% names(pt$summary)) - res$per_trait[[t]]$summary$target <- relabVec(pt$summary$target) + res$per_trait[[t]]$summary$target <- .sldscRelabVec(pt$summary$target, relab) for (bn in c("tau_star_blocks_single", "tau_star_blocks_joint")) { b <- pt[[bn]] if (!is.null(b) && !is.null(colnames(b))) - colnames(res$per_trait[[t]][[bn]]) <- relabVec(colnames(b)) + colnames(res$per_trait[[t]][[bn]]) <- .sldscRelabVec(colnames(b), relab) } } for (mn in names(res$meta)) { if (!is.null(res$meta[[mn]]) && "target" %in% names(res$meta[[mn]])) - res$meta[[mn]]$target <- relabVec(res$meta[[mn]]$target) + res$meta[[mn]]$target <- .sldscRelabVec(res$meta[[mn]]$target, relab) } res$params$target_categories_orig <- res$params$target_categories res$params$target_categories <- unname(relab[targetCategories]) @@ -277,3 +253,34 @@ sldscPostprocessingPipeline <- function(sldscData, res } + +# Build a per-category meta table for one quantity/view, with `label`-prefixed +# mean/se/p columns and an isBinary flag. +# @noRd +.sldscBuildTable <- function(quantity, view, label, isBinary, targetCategories) { + rows <- list() + for (cat in targetCategories) { + m <- metaSldscRandom(view, cat, quantity) + rows[[cat]] <- data.frame( + target = cat, + isBinary = unname(isBinary[cat]), + mean = m$mean, + se = m$se, + p = m$p, + nTraits = m$nTraits, + stringsAsFactors = FALSE + ) + } + df <- do.call(rbind, rows) + rownames(df) <- NULL + nmOld <- c("mean", "se", "p") + nmNew <- paste0(label, toupper(substring(nmOld, 1, 1)), + substring(nmOld, 2)) + names(df)[names(df) %in% nmOld] <- nmNew + df +} + +# Relabel a target-category vector via the `relab` map, leaving unmapped values +# unchanged. +# @noRd +.sldscRelabVec <- function(x, relab) { y <- unname(relab[x]); y[is.na(y)] <- x[is.na(y)]; y } diff --git a/R/sldscWrapper.R b/R/sldscWrapper.R index 09da46c0..d2c1d686 100644 --- a/R/sldscWrapper.R +++ b/R/sldscWrapper.R @@ -464,6 +464,26 @@ metaSldscRandom <- function(perTraitEstimates, category, } +# Append the single/joint enrichment columns (suffix-capitalized) to `out`, +# aligned to out$target; missing sources fill NA. +# @noRd +.sldscAddCols <- function(out, src, suffix) { + colsToAdd <- c("tau", "tauSe", "tauStar", "tauStarSe", + "enrichment", "enrichmentSe", "enrichmentP", + "enrichstat", "enrichstatSe") + suffixCap <- paste0(toupper(substring(suffix, 1, 1)), + substring(suffix, 2)) + for (c in colsToAdd) { + newcol <- paste0(c, suffixCap) + if (!is.null(src) && c %in% names(src)) { + out[[newcol]] <- src[[c]][match(out$target, src$target)] + } else { + out[[newcol]] <- NA_real_ + } + } + out +} + # Internal helper: assemble a wide per-trait summary frame with single + joint # columns side by side. .sldscAssembleTraitSummary <- function(singleDf, jointDf, targetCategories, @@ -474,24 +494,8 @@ metaSldscRandom <- function(perTraitEstimates, category, isBinary = unname(isBinaryVec[rows]), stringsAsFactors = FALSE) - addCols <- function(out, src, suffix) { - colsToAdd <- c("tau", "tauSe", "tauStar", "tauStarSe", - "enrichment", "enrichmentSe", "enrichmentP", - "enrichstat", "enrichstatSe") - suffixCap <- paste0(toupper(substring(suffix, 1, 1)), - substring(suffix, 2)) - for (c in colsToAdd) { - newcol <- paste0(c, suffixCap) - if (!is.null(src) && c %in% names(src)) { - out[[newcol]] <- src[[c]][match(out$target, src$target)] - } else { - out[[newcol]] <- NA_real_ - } - } - out - } - out <- addCols(out, singleDf, "single") - out <- addCols(out, jointDf, "joint") + out <- .sldscAddCols(out, singleDf, "single") + out <- .sldscAddCols(out, jointDf, "joint") out } @@ -518,6 +522,14 @@ metaSldscRandom <- function(perTraitEstimates, category, }) } +# Meta-analyze `quantity` across all target categories for one view, returning a +# per-category named list of metaSldscRandom results. +# @noRd +.sldscPerCategory <- function(view, quantity, targetCategories) + setNames(lapply(targetCategories, + function(cat) metaSldscRandom(view, cat, quantity)), + targetCategories) + #' Random-effects meta-analysis over a subset of sLDSC traits #' #' Re-run the random-effects meta-analysis (DerSimonian-Laird, via @@ -558,13 +570,9 @@ sldscSubsetMeta <- function(postprocessResult, subsetTraits, 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")) + tau_star_single = .sldscPerCategory(viewSingle, "tauStar", targetCategories), + tau_star_joint = .sldscPerCategory(viewJoint, "tauStar", targetCategories), + enrichment = .sldscPerCategory(viewSingle, "enrichment", targetCategories), + enrichstat = .sldscPerCategory(viewSingle, "enrichstat", targetCategories)) } diff --git a/R/sumstatsQc.R b/R/sumstatsQc.R index a3686190..c379042d 100644 --- a/R/sumstatsQc.R +++ b/R/sumstatsQc.R @@ -26,6 +26,27 @@ NULL # package-internal and called from here (.matchAgainstSketch), ctwasPipeline, # and the pipeline join sites. +# Standardize a variant set (GRanges or data.frame) to a chrom/pos/alt/ref frame. +# @noRd +.variantsToDf <- function(x) { + if (is(x, "GRanges")) { + mc <- as.data.frame(mcols(x)) + mc$chrom <- as.character(seqnames(x)) + mc$pos <- start(x) + mc[, c("chrom", "pos", "alt", "ref")] + } else { + as.data.frame(x)[, c("chrom", "pos", "alt", "ref")] + } +} + +# Canonical per-variant key from sorted alleles, so strand/allele flips collide. +# @noRd +.canonicalAlleleKey <- function(df) { + aMin <- pmin(df$alt, df$ref) + aMax <- pmax(df$alt, df$ref) + paste(df$chrom, df$pos, aMin, aMax) +} + #' Merge variant info from two sources with allele-flip-aware matching #' #' Merges variant metadata (chromosome, position, ref, alt) from two sources, @@ -43,30 +64,11 @@ NULL #' \code{ref}, deduplicated by position and alleles. #' @export mergeVariantInfo <- function(variants1, variants2, all = TRUE) { - # Convert GRanges to data.frame if needed - toDf <- function(x) { - if (is(x, "GRanges")) { - mc <- as.data.frame(mcols(x)) - mc$chrom <- as.character(seqnames(x)) - mc$pos <- start(x) - mc[, c("chrom", "pos", "alt", "ref")] - } else { - as.data.frame(x)[, c("chrom", "pos", "alt", "ref")] - } - } + df1 <- .variantsToDf(variants1) + df2 <- .variantsToDf(variants2) - df1 <- toDf(variants1) - df2 <- toDf(variants2) - - # Create a canonical key from sorted alleles so flipped pairs match - makeKey <- function(df) { - aMin <- pmin(df$alt, df$ref) - aMax <- pmax(df$alt, df$ref) - paste(df$chrom, df$pos, aMin, aMax) - } - - key1 <- makeKey(df1) - key2 <- makeKey(df2) + key1 <- .canonicalAlleleKey(df1) + key2 <- .canonicalAlleleKey(df2) # Detect flips: where df2's alt matches df1's ref at the same key matchIdx <- match(key2, key1) @@ -281,6 +283,22 @@ dentist <- function(sumStat, R = NULL, X = NULL, nSample = NULL, return(dentistResult) } +# Module-scoped accumulator for the rsq_eigen warnings raised (per call) by the +# C++ dentistIterativeImpute(); dentistSingleWindow resets $msgs before the call +# and reads it back after. dentist windows are processed serially, so the shared +# env carries no cross-call state. +.dentistRsqAcc <- new.env(parent = emptyenv()) + +# withCallingHandlers warning handler: collect the "rsq_eigen exceeding 1" +# warnings into .dentistRsqAcc and muffle them (they are summarized afterwards). +# @noRd +.dentistRsqHandler <- function(w) { + if (grepl("Adjusted rsq_eigen value exceeding 1", w$message)) { + .dentistRsqAcc$msgs <- c(.dentistRsqAcc$msgs, w$message) + invokeRestart("muffleWarning") + } +} + #' Perform DENTIST on a single window #' #' Detect outliers in GWAS summary statistics using LD-based iterative imputation. @@ -348,14 +366,9 @@ dentistSingleWindow <- function(zScore, R = NULL, X = NULL, nSample = NULL, ldMat <- dedupRes$filteredLD } - # Run C++ iterative imputation (collect rsq warnings) - rsqWarnings <- character(0) - warningHandler <- function(w) { - if (grepl("Adjusted rsq_eigen value exceeding 1", w$message)) { - rsqWarnings <<- c(rsqWarnings, w$message) - invokeRestart("muffleWarning") - } - } + # Run C++ iterative imputation; rsq_eigen warnings are collected into the + # module-level .dentistRsqAcc (reset here, read back after the call). + .dentistRsqAcc$msgs <- character(0) verboseIter <- getOption("pecotmr.dentist.verbose", FALSE) res <- withCallingHandlers( # cpp11 requires exact integer types for int parameters @@ -365,8 +378,9 @@ dentistSingleWindow <- function(zScore, R = NULL, X = NULL, nSample = NULL, gPvalueThreshold, as.integer(ncpus), correctChenEtAlBug, verboseIter ), - warning = warningHandler + warning = .dentistRsqHandler ) + rsqWarnings <- .dentistRsqAcc$msgs if (length(rsqWarnings) > 0) { warning(sprintf("%d rsq_eigen values exceeded 1 (capped at 1.0). Max reported: %s", length(rsqWarnings), rsqWarnings[length(rsqWarnings)])) @@ -655,6 +669,14 @@ slidingWindowLoop <- function(allGaps, n, buildSegmentResult(startList, endList, fillStartList, fillEndList, n, verbose) } +# Apply the quarter-distance index map `quaterIdx` `n` times to `x` (n = 1..4 +# gives the 1st..4th quarter boundary from x). Used by segmentByDist. +# @noRd +.nthQuaterIdx <- function(x, n, quaterIdx) { + for (i in seq_len(n)) x <- quaterIdx[x] + x +} + #' Segment Genomic Region by Distance (Original DENTIST Algorithm) #' #' Implements the same windowing/segmentation algorithm as the original DENTIST C++ binary's @@ -729,10 +751,6 @@ segmentByDist <- function(pos, maxDist = 2000000, minDim = 2000, verbose = FALSE quaterIdx <- pmax(quaterIdx, 1L) # Helper to chain quaterIdx lookups (equivalent to quaterIdx[quaterIdx[x]] in C++) - q1 <- function(x) quaterIdx[x] - q2 <- function(x) quaterIdx[quaterIdx[x]] - q3 <- function(x) quaterIdx[quaterIdx[quaterIdx[x]]] - q4 <- function(x) quaterIdx[quaterIdx[quaterIdx[quaterIdx[x]]]] # Find gaps > cutoff/4 allGaps <- detectGaps(pos, gapThreshold = cutoff / 4, verbose = verbose) @@ -743,21 +761,21 @@ segmentByDist <- function(pos, maxDist = 2000000, minDim = 2000, verbose = FALSE blockSize >= minBlockSize / 2 && (blockSize - minDim) >= 0 }, initEndFn = function(startIdx, blockEnd) { - min(q4(startIdx) + 1, blockEnd) + min(.nthQuaterIdx(startIdx, 4, quaterIdx) + 1, blockEnd) }, fillFn = function(startIdx, endIdx, notStartInterval, notLastInterval) { # Distance mode: fill is always q1 to q3 (inner 50% by distance); # first/last corrections are handled by fix_block_fills in the loop - list(start = q1(startIdx), end = q3(startIdx)) + list(start = .nthQuaterIdx(startIdx, 1, quaterIdx), end = .nthQuaterIdx(startIdx, 3, quaterIdx)) }, stepFn = function(startIdx, blockEnd) { - nextStart <- q2(startIdx) - list(startIdx = nextStart, endIdx = min(q4(nextStart) + 1, blockEnd)) + nextStart <- .nthQuaterIdx(startIdx, 2, quaterIdx) + list(startIdx = nextStart, endIdx = min(.nthQuaterIdx(nextStart, 4, quaterIdx) + 1, blockEnd)) }, adjustLastFn = function(startIdx, oldStartIdx, endIdx, blockEnd) { # If last interval is small, go back one step - if (as.numeric(pos[min(endIdx - 1, n)]) - as.numeric(pos[q1(oldStartIdx)]) < cutoff) { - q1(oldStartIdx) + if (as.numeric(pos[min(endIdx - 1, n)]) - as.numeric(pos[.nthQuaterIdx(oldStartIdx, 1, quaterIdx)]) < cutoff) { + .nthQuaterIdx(oldStartIdx, 1, quaterIdx) } else { startIdx } @@ -873,6 +891,37 @@ mergeWindows <- function(dentistResultByWindow, windowDividedRes) { # SLALoM: Approximate Bayes Factor single-causal-variant outlier detection # ============================================================================= +# Numerically stable log-sum-exp. +# @noRd +.logSumExp <- function(x) { + maxX <- max(x, na.rm = TRUE) + sumExp <- sum(exp(x - maxX), na.rm = TRUE) + return(maxX + log(sumExp)) +} + +# Approximate (Wakefield) Bayes factors from z-scores and standard errors: +# per-variant log-BF and normalized posterior probabilities. +# @noRd +.approxBayesFactor <- function(z, se, W = 0.04) { + V <- se^2 + r <- W / (W + V) + lbf <- 0.5 * (log(1 - r) + (r * z^2)) + denom <- .logSumExp(lbf) + prob <- exp(lbf - denom) + return(list(lbf = lbf, prob = prob)) +} + +# Greedy credible set: indices (ordered by decreasing prob) whose cumulative +# posterior first exceeds `coverage`. +# @noRd +.slalomCredibleSet <- function(prob, coverage = 0.95) { + ordering <- order(prob, decreasing = TRUE) + cumprob <- cumsum(prob[ordering]) + idx <- which(cumprob > coverage)[1] + cs <- ordering[1:idx] + return(cs) +} + #' Slalom Function for Summary Statistics QC for Fine-Mapping Analysis #' #' Performs Approximate Bayesian Factor (ABF) analysis, identifies credible sets, @@ -922,35 +971,12 @@ slalom <- function(zScore, R = NULL, X = NULL, standardError = rep(1, length(zSc # This selects the most negative z-score as lead when leadVariantChoice == "pvalue". pvalue <- pnorm(zScore) - logSumExp <- function(x) { - maxX <- max(x, na.rm = TRUE) - sumExp <- sum(exp(x - maxX), na.rm = TRUE) - return(maxX + log(sumExp)) - } - - abf <- function(z, se, W = 0.04) { - V <- se^2 - r <- W / (W + V) - lbf <- 0.5 * (log(1 - r) + (r * z^2)) - denom <- logSumExp(lbf) - prob <- exp(lbf - denom) - return(list(lbf = lbf, prob = prob)) - } - - abfResults <- abf(zScore, standardError, W = abfPriorVariance) + abfResults <- .approxBayesFactor(zScore, standardError, W = abfPriorVariance) lbf <- abfResults$lbf prob <- abfResults$prob - getCs <- function(prob, coverage = 0.95) { - ordering <- order(prob, decreasing = TRUE) - cumprob <- cumsum(prob[ordering]) - idx <- which(cumprob > coverage)[1] - cs <- ordering[1:idx] - return(cs) - } - - cs <- getCs(prob, coverage = 0.95) - cs99 <- getCs(prob, coverage = 0.99) + cs <- .slalomCredibleSet(prob, coverage = 0.95) + cs99 <- .slalomCredibleSet(prob, coverage = 0.99) leadIdx <- if (leadVariantChoice == "pvalue") { which.min(pvalue) @@ -1152,6 +1178,27 @@ extractTopPipInfo <- function(conData) { ) } +# Parse one comma-joined cs_corr string into its values + max/min |corr| (self- +# correlations equal to 1 excluded); invalid/empty input yields empties + NA. +# @noRd +.extractCorrelations <- function(x) { + # Early return if x is invalid + if(is.na(x) || x == "" || is.null(x) || !grepl(",", as.character(x))) { + return(list(values = numeric(0), max_corr = NA_real_, min_corr = NA_real_)) + } + + # Convert and filter values + values <- as.numeric(unlist(strsplit(x, ","))) + valuesFiltered <- abs(values[values != 1]) + + # Return list with NA if no valid correlations + list( + values = values, + max_corr = if(length(valuesFiltered) > 0) max(abs(valuesFiltered), na.rm = TRUE) else NA_real_, + min_corr = if(length(valuesFiltered) > 0) min(abs(valuesFiltered), na.rm = TRUE) else NA_real_ + ) +} + #' Parse Credible Set Correlations from extractCsInfo() Output #' #' This function takes the output from `extractCsInfo()` and expands the `cs_corr` column @@ -1183,25 +1230,8 @@ parseCsCorr <- function(df) { # Ensure we work with a data frame df <- as.data.frame(df) - extractCorrelations <- function(x) { - # Early return if x is invalid - if(is.na(x) || x == "" || is.null(x) || !grepl(",", as.character(x))) { - return(list(values = numeric(0), max_corr = NA_real_, min_corr = NA_real_)) - } - - # Convert and filter values - values <- as.numeric(unlist(strsplit(x, ","))) - valuesFiltered <- abs(values[values != 1]) - - # Return list with NA if no valid correlations - list( - values = values, - max_corr = if(length(valuesFiltered) > 0) max(abs(valuesFiltered), na.rm = TRUE) else NA_real_, - min_corr = if(length(valuesFiltered) > 0) min(abs(valuesFiltered), na.rm = TRUE) else NA_real_ - ) - } # Process correlations - processedResults <- lapply(df$cs_corr, extractCorrelations) + processedResults <- lapply(df$cs_corr, .extractCorrelations) # If no valid results, add NA columns and return if(all(sapply(processedResults, function(x) length(x$values) == 0))) { df$cs_corr_max <- NA_real_ @@ -1516,6 +1546,49 @@ raissSingleMatrixFromX <- function(refPanel, knownZscores, X, lamb = 0.01, )) } +# Sequentially append RAISS block results, resolving the shared boundary variant +# (duplicated across adjacent blocks) by keeping the higher-R2 imputation. +# @noRd +.combineWithBoundaryCheck <- function(combinedResult, newResult) { + # If either is empty, simply return the non-empty one or empty data frame + if (is.null(combinedResult)) { + return(newResult) + } + if (is.null(newResult)) { + return(combinedResult) + } + + # Check if the last variant of combined matches the first of new + lastVar <- combinedResult$variant_id[nrow(combinedResult)] + firstVar <- newResult$variant_id[1] + + if (lastVar == firstVar) { + newR2 <- newResult$raissR2[1] + oldR2 <- combinedResult$raissR2[nrow(combinedResult)] + if (is.na(newR2) && is.na(oldR2)) { + # Both are NA - keep the existing one + } else if (is.na(oldR2)) { + # Old is NA but new is not - use new + combinedResult[nrow(combinedResult), ] <- newResult[1, ] + } else if (is.na(newR2)) { + # New is NA but old is not - keep old + } else if (newR2 > oldR2) { + # Both are non-NA and new is better - use new + combinedResult[nrow(combinedResult), ] <- newResult[1, ] + } + + # Add remaining rows from new (excluding first) + if (nrow(newResult) > 1) { + combinedResult <- bind_rows(combinedResult, newResult[-1, ]) + } + } else { + # No overlap - combine all rows + combinedResult <- bind_rows(combinedResult, newResult) + } + + return(combinedResult) +} + #' Impute Summary Statistics Using LD (RAISS) #' #' This function is a part of the statistical library for SNP imputation from: @@ -1623,46 +1696,6 @@ raiss <- function(refPanel, knownZscores, ldMatrix = NULL, # For list of matrices, process each block if (verbose) message("Processing multiple LD blocks...") - combineWithBoundaryCheck <- function(combinedResult, newResult) { - # If either is empty, simply return the non-empty one or empty data frame - if (is.null(combinedResult)) { - return(newResult) - } - if (is.null(newResult)) { - return(combinedResult) - } - - # Check if the last variant of combined matches the first of new - lastVar <- combinedResult$variant_id[nrow(combinedResult)] - firstVar <- newResult$variant_id[1] - - if (lastVar == firstVar) { - newR2 <- newResult$raissR2[1] - oldR2 <- combinedResult$raissR2[nrow(combinedResult)] - if (is.na(newR2) && is.na(oldR2)) { - # Both are NA - keep the existing one - } else if (is.na(oldR2)) { - # Old is NA but new is not - use new - combinedResult[nrow(combinedResult), ] <- newResult[1, ] - } else if (is.na(newR2)) { - # New is NA but old is not - keep old - } else if (newR2 > oldR2) { - # Both are non-NA and new is better - use new - combinedResult[nrow(combinedResult), ] <- newResult[1, ] - } - - # Add remaining rows from new (excluding first) - if (nrow(newResult) > 1) { - combinedResult <- bind_rows(combinedResult, newResult[-1, ]) - } - } else { - # No overlap - combine all rows - combinedResult <- bind_rows(combinedResult, newResult) - } - - return(combinedResult) - } - resultsList <- list() variantIndices <- ldMatrix$variantIndices blockIds <- unique(variantIndices$blockId) @@ -1704,12 +1737,12 @@ raiss <- function(refPanel, knownZscores, ldMatrix = NULL, if (length(resultsList) > 1) { for (i in 2:length(resultsList)) { - combinedNofilter <- combineWithBoundaryCheck( + combinedNofilter <- .combineWithBoundaryCheck( combinedNofilter, resultsList[[i]]$resultNofilter ) - combinedFilter <- combineWithBoundaryCheck( + combinedFilter <- .combineWithBoundaryCheck( combinedFilter, resultsList[[i]]$resultFilter ) @@ -1820,6 +1853,13 @@ mergeRaissDf <- function(raissDf, knownZscores) { return(mergedDf) } +# Format one aligned "label: value" report line (label left-padded to +# maxLabelLength). +# @noRd +.formatRaissLine <- function(label, value, maxLabelLength) { + sprintf("%-*s %d", maxLabelLength, paste0(label, ":"), value) +} + filterRaissOutput <- function(zscores, r2Threshold = 0.6, minimumLd = 5, verbose = TRUE) { # Reset the index and subset the data frame zscores <- zscores[, c("chrom", "pos", "variant_id", "A1", "A2", "z", "Var", "raissLdScore")] @@ -1848,17 +1888,13 @@ filterRaissOutput <- function(zscores, r2Threshold = 0.6, minimumLd = 5, verbose "Remaining variants after filter:" ))) - formatLine <- function(label, value) { - sprintf("%-*s %d", maxLabelLength, paste0(label, ":"), value) - } - message("IMPUTATION REPORT\n") - message(formatLine("Variants before filter", nSnpsBfFilt)) - message(formatLine("Non-imputed variants", nSnpsInitial)) - message(formatLine("Imputed variants", nSnpsImputed)) - message(formatLine("Variants filtered because of low LD score", nSnpsLdFilt)) - message(formatLine("Variants filtered because of low R2", nSnpsR2Filt)) - message(formatLine("Remaining variants after filter", nSnpsAfFilt)) + message(.formatRaissLine("Variants before filter", nSnpsBfFilt, maxLabelLength)) + message(.formatRaissLine("Non-imputed variants", nSnpsInitial, maxLabelLength)) + message(.formatRaissLine("Imputed variants", nSnpsImputed, maxLabelLength)) + message(.formatRaissLine("Variants filtered because of low LD score", nSnpsLdFilt, maxLabelLength)) + message(.formatRaissLine("Variants filtered because of low R2", nSnpsR2Filt, maxLabelLength)) + message(.formatRaissLine("Remaining variants after filter", nSnpsAfFilt, maxLabelLength)) } return(zscore_list = list(zscoresNofilter = zscoresNofilter, zscores = zscores)) } @@ -2159,7 +2195,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, # # `require` character vector of mcol names that MUST be present; errors # when any is missing. Use this for the strict callers -# (e.g. `.fmExtractZN` needs SNP + Z + N to proceed). +# (e.g. `.fmExtractZn` needs SNP + Z + N to proceed). # `derive` when "zFromBetaSe" and `z` is absent but BETA + SE are # present, set z := BETA/SE. Default "none". # `label` error-message prefix for missing-`require` errors. @@ -2591,15 +2627,22 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, list(df = df, skipped = FALSE) } +# Prefix QC-track log lines with the entry label `lbl` (as `[lbl] ...`), or emit +# them bare when `lbl` is NA. +# @noRd +.qcEmit <- function(lbl, ...) { + if (is.na(lbl)) message(...) else message("[", lbl, "] ", ...) +} + # Internal: canonicalize the working per-variant `N`. The N source is resolved # by a four-level priority: (1) per-variant N_CASE / N_CONTROL columns, (2) study # -level nCase / nControl scalars, (3) a per-variant N column, (4) a study-level # nSample scalar (total N). Levels 1-2 give the effective sample size (default) # or the raw total (escape hatch). Returns list(df=, nSource=), where nSource is # "effective" | "column" | "total" | "study-n" | NA_character_ (no source). -# `emit` logs the counts-win override. A study-level scalar fills `df$N` even -# when the entry has no per-variant N column. -.resolveEffectiveN <- function(df, opts, emit) { +# The entry label `lbl` prefixes the counts-win override log. A study-level +# scalar fills `df$N` even when the entry has no per-variant N column. +.resolveEffectiveN <- function(df, opts, lbl) { hasCols <- all(c("N_CASE", "N_CONTROL") %in% colnames(df)) hasScalar <- !is.null(opts$nCase) && !is.null(opts$nControl) && length(opts$nCase) == 1L && length(opts$nControl) == 1L && @@ -2631,13 +2674,13 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, # Default: per-variant c/c -> study c/c -> per-variant N -> study nSample. if (hasCols) { - if (hasN) emit("QC track: N overridden by effective N from per-variant ", + if (hasN) .qcEmit(lbl, "QC track: N overridden by effective N from per-variant ", "n_case/n_control.") df$N <- effectiveN(df$N_CASE, df$N_CONTROL) return(list(df = df, nSource = "effective")) } if (hasScalar) { - if (hasN) emit("QC track: N overridden by effective N from study ", + if (hasN) .qcEmit(lbl, "QC track: N overridden by effective N from study ", "nCase/nControl.") df$N <- rep(effectiveN(opts$nCase, opts$nControl), nRow) return(list(df = df, nSource = "effective")) @@ -2664,10 +2707,6 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, imputeBefore = NA_integer_, imputeAfter = NA_integer_) lbl <- if (!is.null(entryLabel) && nzchar(entryLabel)) entryLabel else NA_character_ - emit <- function(...) { - if (is.na(lbl)) message(...) else message("[", lbl, "] ", ...) - } - df <- .entryGrangesToDf(gr) entryAudit$variantsIn <- nrow(df) nStudyIn <- nrow(df) @@ -2690,14 +2729,14 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, df <- sanity$df if (length(sanity$audit) > 0L) entryAudit$sanityChecks <- sanity$audit if (nSanIn > 0L && nrow(df) != nSanIn) { - emit("QC track: sanity checks kept ", nrow(df), " of ", nSanIn, + .qcEmit(lbl, "QC track: sanity checks kept ", nrow(df), " of ", nSanIn, " variant(s).") } # Canonicalize N to the effective sample size for case/control input, # BEFORE the N-cutoff filter so the filter, kriging (median df$N), and the # rebuilt entry all consume N_eff. No-op for quantitative traits (no counts). - nRes <- .resolveEffectiveN(df, opts, emit) + nRes <- .resolveEffectiveN(df, opts, lbl) df <- nRes$df entryAudit$nSource <- nRes$nSource # Keep the PIP-screen / kriging sample size consistent with the canonicalized @@ -2719,7 +2758,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, if (length(contentFiltered$audit) > 0L) entryAudit$contentFilters <- contentFiltered$audit if (nFiltIn > 0L && nrow(df) != nFiltIn) { - emit("QC track: MAF/INFO/N filters kept ", nrow(df), " of ", nFiltIn, + .qcEmit(lbl, "QC track: MAF/INFO/N filters kept ", nrow(df), " of ", nFiltIn, " variant(s).") } @@ -2790,13 +2829,13 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, qcCount$harmCorrSign <- harmCounts$signFlip qcCount$harmCorrStrand <- harmCounts$strandFlip qcCount$harmDropped <- nHarmIn - nrow(df) - emit("QC track: harmonization kept ", nrow(df), " of ", nHarmIn, + .qcEmit(lbl, "QC track: harmonization kept ", nrow(df), " of ", nHarmIn, " variant(s) (corrected: sign-flipped ", harmCounts$signFlip, ", strand-flipped ", harmCounts$strandFlip, "; dropped ", qcCount$harmDropped, ").") } else { qcCount$harmDropped <- nHarmIn - nrow(df) - emit("QC track: harmonization kept ", nrow(df), " of ", nHarmIn, + .qcEmit(lbl, "QC track: harmonization kept ", nrow(df), " of ", nHarmIn, " variant(s).") } @@ -2826,7 +2865,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, entryAudit$krigingFlipped <- nKr entryAudit$krigingDiagnostics <- kr$diagnostics qcCount$krigingFlipped <- nKr - emit("QC track: kriging sign-flipped ", nKr, " of ", nKrIn, + .qcEmit(lbl, "QC track: kriging sign-flipped ", nKr, " of ", nKrIn, " LD-inconsistent variant(s).") } @@ -2844,7 +2883,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, if (!is.null(ldQc$diagnostics)) entryAudit$ldMismatchDiagnostics <- ldQc$diagnostics qcCount$mismatchRemoved <- ldQc$outliers - emit("QC track: ", opts$zMismatchQc, " removed ", ldQc$outliers, " of ", + .qcEmit(lbl, "QC track: ", opts$zMismatchQc, " removed ", ldQc$outliers, " of ", nMmIn, " LD-mismatch outlier(s).") } @@ -2854,22 +2893,24 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, # Scope the reference panel + dosage to the analysis-region window. RAISS # imputation is local (each missing variant is filled from its LD # neighbours), so a per-chromosome / genome-wide sketch must NOT materialize - # its full dosage here -- restrict to [min(pos) - flank, max(pos) + flank] on - # the region chromosome, mirroring the region-scoping every other QC step + # its full dosage here -- restrict to [min(pos) - flank, max(pos) + flank] + # computed PER CHROMOSOME (a multi-chromosome entry must not impute across a + # cross-chromosome span), mirroring the region-scoping every other QC step # (harmonization, kriging, LD-mismatch) already does via df$SNP. flank <- if (is.null(opts$imputeOpts$flank)) 0L else as.integer(opts$imputeOpts$flank) - regChrom <- unique(sub("^chr", "", as.character(df$chrom), ignore.case = TRUE)) - lo <- min(as.integer(df$pos), na.rm = TRUE) - flank - hi <- max(as.integer(df$pos), na.rm = TRUE) + flank + dfChrom <- sub("^chr", "", as.character(df$chrom), ignore.case = TRUE) + dfPos <- as.integer(df$pos) + loByChr <- tapply(dfPos, dfChrom, min, na.rm = TRUE) - flank + hiByChr <- tapply(dfPos, dfChrom, max, na.rm = TRUE) + flank sketchSnpInfo <- getSnpInfo(ldSketch) skChrom <- sub("^chr", "", as.character(sketchSnpInfo$CHR), ignore.case = TRUE) - windowIdx <- which(skChrom %in% regChrom & - as.integer(sketchSnpInfo$BP) >= lo & - as.integer(sketchSnpInfo$BP) <= hi) + skBp <- as.integer(sketchSnpInfo$BP) + windowIdx <- which(skChrom %in% names(loByChr) & + skBp >= loByChr[skChrom] & skBp <= hiByChr[skChrom]) if (length(windowIdx) == 0L) { - emit("QC track: RAISS imputation skipped (no LD-panel variants in the ", + .qcEmit(lbl, "QC track: RAISS imputation skipped (no LD-panel variants in the ", "region window).") entryAudit$raissImputedVariants <- 0L qcCount$imputeAfter <- nrow(df) @@ -2940,7 +2981,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, entryAudit$raissImputedVariants <- 0L } qcCount$imputeAfter <- nrow(df) - emit("QC track: RAISS imputation ", qcCount$imputeBefore, " -> ", + .qcEmit(lbl, "QC track: RAISS imputation ", qcCount$imputeBefore, " -> ", qcCount$imputeAfter, " variant(s) (net ", sprintf("%+d", qcCount$imputeAfter - qcCount$imputeBefore), ").") } @@ -2993,7 +3034,7 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, paste0(" | imputed ", sprintf("%+d", qcCount$imputeAfter - qcCount$imputeBefore)) } else "" - emit("QC summary: ", nStudyIn, " in -> ", nrow(df), " out", + .qcEmit(lbl, "QC summary: ", nStudyIn, " in -> ", nrow(df), " out", " | corrected: ", correctedSeg, if (length(removedSegs) > 0L) paste0(" | removed: ", paste(removedSegs, collapse = ", ")) @@ -3004,6 +3045,53 @@ krigingOutlierQc <- function(zScore, R, n, variantIds = NULL, list(gr = .dfToEntryGranges(df), audit = entryAudit) } +# Shrink an LD-sketch GenotypeHandle to the panel variants inside the summary +# statistics' per-chromosome position span. `entries` is a list/SimpleList of +# per-study (or per-tuple) GRanges. A genome-wide sketch otherwise carries a +# full-genome snpInfo; only variants inside [min,max] BP of each represented +# chromosome are reachable by harmonization or within-range imputation, so the +# rest is dropped at load time. NULL-safe; a no-op when the span already covers +# the panel. See [[.subsetGenotypeHandle]] for why this is read-safe. +# @noRd +.subsetSketchToRange <- function(ldSketch, entries) { + if (is.null(ldSketch)) return(NULL) + chrom <- unlist(lapply(entries, function(gr) + canonChrom(as.character(GenomicRanges::seqnames(gr)))), use.names = FALSE) + pos <- unlist(lapply(entries, function(gr) + as.integer(GenomicRanges::start(gr))), use.names = FALSE) + ok <- !is.na(chrom) & !is.na(pos) + chrom <- chrom[ok]; pos <- pos[ok] + if (length(pos) == 0L) return(ldSketch) + lo <- tapply(pos, chrom, min) + hi <- tapply(pos, chrom, max) + si <- getSnpInfo(ldSketch) + siChrom <- canonChrom(as.character(si$CHR)) + siBp <- as.integer(si$BP) + keep <- siChrom %in% names(lo) & + siBp >= lo[siChrom] & siBp <= hi[siChrom] + keep[is.na(keep)] <- FALSE + .subsetGenotypeHandle(ldSketch, keep) +} + +# Shrink an LD-sketch GenotypeHandle to EXACTLY the variants present across the +# QC'd `entries` (imputation may have added variants; QC may have dropped some), +# matched by canonical variant id. Applied at the end of summaryStatsQc so the +# retained sketch mirrors the object's final variant set. NULL-safe. +# @noRd +.subsetSketchToIds <- function(ldSketch, entries) { + if (is.null(ldSketch)) return(NULL) + ids <- unlist(lapply(entries, function(gr) { + if (is.null(gr)) return(character(0)) + snp <- S4Vectors::mcols(gr)$SNP + if (is.null(snp)) character(0) else as.character(snp) + }), use.names = FALSE) + if (length(ids) == 0L) return(ldSketch) + si <- getSnpInfo(ldSketch) + keep <- normalizeVariantId(as.character(si$SNP)) %in% + normalizeVariantId(unique(ids)) + .subsetGenotypeHandle(ldSketch, keep) +} + #' Run QC on a SumStats Collection #' #' Applies a single QC pass to a \code{QtlSumStats} or \code{GwasSumStats} @@ -3239,13 +3327,17 @@ summaryStatsQc <- function(sumstats, dropNonpositiveSe = dropNonpositiveSe), entryAudit = entryAudits) + # Trim the LD sketch to the variants actually present after QC (imputation may + # have added some; filters removed others) so it mirrors the final object. + newLdSketch <- .subsetSketchToIds(getLdSketch(sumstats), newEntries) + # Rebuild the SumStats with new entries and qcInfo. if (methods::is(sumstats, "GwasSumStats")) { GwasSumStats( study = as.character(sumstats$study), entry = newEntries, genome = getGenome(sumstats), - ldSketch = getLdSketch(sumstats), + ldSketch = newLdSketch, varY = as.numeric(sumstats$varY), # Preserve the optional per-study case/control counts through QC. nCase = if ("nCase" %in% names(sumstats)) @@ -3262,7 +3354,7 @@ summaryStatsQc <- function(sumstats, trait = as.character(sumstats$trait), entry = newEntries, genome = getGenome(sumstats), - ldSketch = getLdSketch(sumstats), + ldSketch = newLdSketch, varY = as.numeric(sumstats$varY), # Preserve the optional per-tuple study-level total N through QC. nSample = if ("nSample" %in% names(sumstats)) diff --git a/R/tupleSelectors.R b/R/tupleSelectors.R index 1d2c6979..eb80c4c7 100644 --- a/R/tupleSelectors.R +++ b/R/tupleSelectors.R @@ -213,6 +213,14 @@ rep(NA, n) } +# First existing value for column `cn` across `parts` (a type exemplar used for +# NA-filling); NULL if no part has it. +# @noRd +.exemplarColumn <- function(cn, parts) { + for (p in parts) if (cn %in% names(p)) return(p[[cn]]) + NULL # nocov (unreachable: cn is always drawn from allCols, so some part has it) +} + # 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 @@ -225,13 +233,9 @@ 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 # nocov (unreachable: cn is always drawn from allCols, so some part has it) - } combined <- lapply(allCols, function(cn) { pieces <- lapply(parts, function(p) - if (cn %in% names(p)) p[[cn]] else .naLikeColumn(exemplar(cn), nrow(p))) + if (cn %in% names(p)) p[[cn]] else .naLikeColumn(.exemplarColumn(cn, parts), nrow(p))) suppressWarnings(do.call(c, pieces)) # GRanges concat may warn on seqinfo }) names(combined) <- allCols diff --git a/R/twasWeights.R b/R/twasWeights.R index e18f35b4..2e321bef 100644 --- a/R/twasWeights.R +++ b/R/twasWeights.R @@ -422,6 +422,13 @@ setMethod("show", "TwasWeights", function(object) { result } +# TRUE if `name` is a function visible in the search path or in namespace `ns`. +# @noRd +.functionExistsInNs <- function(name, ns) { + exists(name, mode = "function") || + exists(name, mode = "function", envir = ns, inherits = FALSE) +} + # Resolve the actual function name for a method key. Honors an "impl" attribute # on the per-method args list (set by .twasMethodLookup), and otherwise applies # a snake_case -> camelCase transformation as a fallback for user-supplied @@ -430,25 +437,46 @@ setMethod("show", "TwasWeights", function(object) { # Search pecotmr's namespace explicitly so this works equally well when the # function is called either from inside the package or from a user session. ns <- asNamespace("pecotmr") - fnExists <- function(name) { - exists(name, mode = "function") || - exists(name, mode = "function", envir = ns, inherits = FALSE) - } impl <- if (!is.null(methodArgs)) attr(methodArgs, "impl") else NULL - if (!is.null(impl) && nzchar(impl) && fnExists(impl)) { + if (!is.null(impl) && nzchar(impl) && .functionExistsInNs(impl, ns)) { return(impl) } # Direct match (e.g. caller already passed camelCase) - if (fnExists(methodKey)) return(methodKey) + if (.functionExistsInNs(methodKey, ns)) return(methodKey) # snake_case_weights -> camelCaseWeights parts <- strsplit(methodKey, "_", fixed = TRUE)[[1]] capRest <- paste0(toupper(substring(parts[-1], 1, 1)), substring(parts[-1], 2)) candidate <- paste0(parts[1], paste0(capRest, collapse = "")) - if (fnExists(candidate)) return(candidate) + if (.functionExistsInNs(candidate, ns)) return(candidate) methodKey } +# Validate a Sample/Fold partition data frame: required columns, no sample in +# two folds, and (when `sampleNames` given) exact coverage of all samples. +# @noRd +.validateFoldPartition <- function(df, sampleNames) { + if (!all(c("Sample", "Fold") %in% names(df))) + stop("samplePartition must have columns `Sample` and `Fold`.") + df$Sample <- as.character(df$Sample) + dup <- unique(df$Sample[duplicated(df$Sample)]) + if (length(dup) > 0L) + stop("Fold partition assigns sample(s) to more than one fold: ", + paste(dup, collapse = ", ")) + if (!is.null(sampleNames)) { + unknown <- setdiff(df$Sample, sampleNames) + if (length(unknown) > 0L) + stop("Fold partition references unknown sample(s): ", + paste(unknown, collapse = ", ")) + uncovered <- setdiff(sampleNames, df$Sample) + if (length(uncovered) > 0L) + stop("Fold partition does not cover ", length(uncovered), + " sample(s) (folds must partition all samples), e.g. ", + paste(utils::head(uncovered, 5L), collapse = ", ")) + } + df +} + # Normalize a cross-validation fold specification into the canonical # samplePartition data.frame(Sample, Fold) used throughout the CV machinery, so # callers can pass folds in any of three forms and downstream code has a single @@ -469,31 +497,10 @@ setMethod("show", "TwasWeights", function(object) { stop("Provide either a list-form `cvFolds` or an explicit ", "`samplePartition`, not both.") - validatePartition <- function(df) { - if (!all(c("Sample", "Fold") %in% names(df))) - stop("samplePartition must have columns `Sample` and `Fold`.") - df$Sample <- as.character(df$Sample) - dup <- unique(df$Sample[duplicated(df$Sample)]) - if (length(dup) > 0L) - stop("Fold partition assigns sample(s) to more than one fold: ", - paste(dup, collapse = ", ")) - if (!is.null(sampleNames)) { - unknown <- setdiff(df$Sample, sampleNames) - if (length(unknown) > 0L) - stop("Fold partition references unknown sample(s): ", - paste(unknown, collapse = ", ")) - uncovered <- setdiff(sampleNames, df$Sample) - if (length(uncovered) > 0L) - stop("Fold partition does not cover ", length(uncovered), - " sample(s) (folds must partition all samples), e.g. ", - paste(utils::head(uncovered, 5L), collapse = ", ")) - } - df - } if (!is.null(samplePartition)) { - df <- validatePartition(as.data.frame(samplePartition, - stringsAsFactors = FALSE)) + df <- .validateFoldPartition(as.data.frame(samplePartition, + stringsAsFactors = FALSE), sampleNames) return(list(samplePartition = df, nFolds = length(unique(df$Fold)))) } @@ -514,7 +521,7 @@ setMethod("show", "TwasWeights", function(object) { } data.frame(Sample = ids, Fold = k, stringsAsFactors = FALSE) }) - df <- validatePartition(do.call(rbind, rows)) + df <- .validateFoldPartition(do.call(rbind, rows), sampleNames) return(list(samplePartition = df, nFolds = length(cvFolds))) } @@ -596,6 +603,61 @@ setMethod("show", "TwasWeights", function(object) { weightMethods } +# Per-fold TWAS weight fit for the CV engine. `ctx` carries weightMethods, +# multivariateWeightMethods, cvArgs, retainFits, verbose. Weights are keyed by the +# canonical method key; captured fits keep the full method name. +# @noRd +.weightFitFold <- function(Xtr, Ytr, j, ctx) { + weightMethods <- ctx$weightMethods + multivariateWeightMethods <- ctx$multivariateWeightMethods + cvArgs <- ctx$cvArgs; retainFits <- ctx$retainFits; verbose <- ctx$verbose + 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]] + } + if (!is.null(cvArgs$reweightedMixturePriorCv) && + method %in% c("mvsusie_weights", "mvsusieWeights")) { + args$prior_variance <- cvArgs$reweightedMixturePriorCv[[j]] + } + 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 { + 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 + } + list(weights = weights, fits = fits) +} + #' Cross-Validation for weights selection in Transcriptome-Wide Association Studies (TWAS) #' #' Performs cross-validation for TWAS, supporting both univariate and multivariate methods. @@ -652,62 +714,18 @@ twasWeightsCv <- function(X, Y, fold = NULL, samplePartitions = NULL, weightMeth 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 + # Per-fold fit context passed to the shared engine's top-level fitter + # (.weightFitFold). Weights are keyed by the canonical method key; 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]] - } - if (!is.null(cvArgs$reweightedMixturePriorCv) && - method %in% c("mvsusie_weights", "mvsusieWeights")) { - args$prior_variance <- cvArgs$reweightedMixturePriorCv[[j]] - } - 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 { - 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 - } - list(weights = weights, fits = fits) - } + cvFitCtx <- list(weightMethods = weightMethods, + multivariateWeightMethods = multivariateWeightMethods, + cvArgs = cvArgs, retainFits = retainFits, verbose = verbose) # 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()), + fitFold = .cvNoopFitFold, numThreads = numThreads, maxNumVariants = maxNumVariants, variantsToKeep = variantsToKeep, retainFits = retainFits, verbose = verbose) return(list(samplePartition = res$samplePartition)) @@ -715,11 +733,145 @@ twasWeightsCv <- function(X, Y, fold = NULL, samplePartitions = NULL, weightMeth .crossValidateWeights( X, Y, fold = fold, samplePartitions = samplePartitions, - fitFold = weightFitFold, numThreads = numThreads, + fitFold = .weightFitFold, fitFoldCtx = cvFitCtx, numThreads = numThreads, maxNumVariants = maxNumVariants, variantsToKeep = variantsToKeep, retainFits = retainFits, verbose = verbose) } +# Fit one TWAS weight method by name against the filtered design matrix, embedding +# the fitted weights back into the full variant space. `ctx` carries the shared +# fit state (X, Y, Xfiltered, validColumns, retainFits, retainFitDetail, verbose). +# @noRd +.computeMethodWeights <- function(methodName, weightMethods, ctx) { + X <- ctx$X; Y <- ctx$Y; Xfiltered <- ctx$Xfiltered + validColumns <- ctx$validColumns; retainFits <- ctx$retainFits + retainFitDetail <- ctx$retainFitDetail; verbose <- ctx$verbose + shortName <- sub("_weights$", "", methodName) + if (verbose >= 1) { + message(sprintf(" Fitting %s ...", shortName)) + tic() + } + + # Hardcoded vector of multivariate methods (accept both snake and camel). + # fSuSiE is multivariate (variants x features weight matrix) but is never + # refit here — fsusieWeights extracts from the supplied fsusieFit. + multivariateWeightMethods <- c("mrmash_weights", "mvsusie_weights", + "fsusie_weights", + "mrmashWeights", "mvsusieWeights", + "fsusieWeights") + args <- weightMethods[[methodName]] + fnName <- .resolveMethodFunction(methodName, args) + + # Only pass retainFit (or its legacy snake_case alias) to functions that accept it + if (retainFits) { + fnFormals <- names(formals(fnName)) + if ("retainFit" %in% fnFormals) { + args$retainFit <- TRUE + } else if ("retain_fit" %in% fnFormals) { + args$retain_fit <- TRUE + } + # Propagate the slim/full payload choice to producers that support it + # (mr.mash individual + RSS), unless the caller already set it per-method. + if ("fitDetail" %in% fnFormals && is.null(args$fitDetail)) { + args$fitDetail <- retainFitDetail + } + } + + methodFit <- NULL + if (methodName %in% multivariateWeightMethods) { + # Apply multivariate method + weightsMatrix <- if (verbose < 2) { + .quietEval(do.call(fnName, c(list(X = Xfiltered, Y = Y), args))) + } else { + do.call(fnName, c(list(X = Xfiltered, Y = Y), args)) + } + if (retainFits) methodFit <- attr(weightsMatrix, "fit") + if (nrow(weightsMatrix) != length(validColumns)) weightsMatrix <- weightsMatrix[names(validColumns), , drop = FALSE] + } else { + # Apply univariate method to each column of Y + # Initialize it with zeros to avoid NA + weightsMatrix <- matrix(0, nrow = ncol(Xfiltered), ncol = ncol(Y)) + + for (k in 1:ncol(Y)) { + weightsVector <- if (verbose < 2) { + .quietEval(do.call(fnName, c(list(X = Xfiltered, y = Y[, k]), args))) + } else { + do.call(fnName, c(list(X = Xfiltered, y = Y[, k]), args)) + } + if (retainFits && is.null(methodFit)) { + methodFit <- attr(weightsVector, "fit") + } + if (is.matrix(weightsVector)) weightsVector <- weightsVector[, k] + weightsMatrix[, k] <- weightsVector + } + } + + result <- .embedWeights(weightsMatrix, validColumns, ncol(X), ncol(Y), colnames(X), colnames(Y)) + if (!is.null(methodFit)) attr(result, "fit") <- methodFit + if (verbose >= 1) { + elapsed <- toc(quiet = TRUE) + message(sprintf(" Fitting %s done in %.1fs", shortName, elapsed$toc - elapsed$tic)) + } + return(result) +} + +# Assemble the (study, context, trait, method, entry) row vectors for the +# TwasWeights collection from the fitted `weightsList`. `ctx` carries the shared +# identity + flags (study, context, trait, Y, retainFits, standardized, dataType). +# @noRd +.buildTwasWeightEntries <- function(weightsList, variantIds, ctx) { + Y <- ctx$Y; study <- ctx$study; context <- ctx$context; trait <- ctx$trait + retainFits <- ctx$retainFits; standardized <- ctx$standardized; dataType <- ctx$dataType + studies <- character(0) + contexts <- character(0) + traits <- character(0) + methodsV <- character(0) + entries <- list() + for (m in names(weightsList)) { + wMat <- weightsList[[m]] + fitVal <- attr(wMat, "fit") + attr(wMat, "fit") <- NULL + shortMethod <- sub("(_weights|Weights)$", "", m) + # When trait/context were supplied per-row (length == ncol(Y)), emit + # one row per (method, outcome). Otherwise emit one row per method + # and carry the (possibly multi-column) weights matrix as-is. + perOutcome <- length(trait) == ncol(Y) && + length(context) %in% c(1L, ncol(Y)) + if (perOutcome) { + contextV <- if (length(context) == 1L) rep(context, ncol(Y)) else context + studyV <- if (length(study) == 1L) rep(study, ncol(Y)) else study + for (k in seq_len(ncol(Y))) { + studies <- c(studies, studyV[k]) + contexts <- c(contexts, contextV[k]) + traits <- c(traits, trait[k]) + methodsV <- c(methodsV, shortMethod) + entries[[length(entries) + 1L]] <- TwasWeightsEntry( + variantIds = variantIds, + weights = wMat[, k], + fits = if (retainFits) fitVal else NULL, + cvResult = NULL, + standardized = isTRUE(standardized), + dataType = dataType) + } + } else { + studies <- c(studies, study[1L]) + contexts <- c(contexts, context[1L]) + traits <- c(traits, trait[1L]) + methodsV <- c(methodsV, shortMethod) + wPayload <- if (ncol(wMat) == 1L) drop(wMat) else wMat + entries[[length(entries) + 1L]] <- TwasWeightsEntry( + variantIds = variantIds, + weights = wPayload, + fits = if (retainFits) fitVal else NULL, + cvResult = NULL, + standardized = isTRUE(standardized), + dataType = dataType) + } + } + list(study = studies, context = contexts, trait = traits, + method = methodsV, entry = entries) +} + #' Run multiple TWAS weight methods #' #' Applies specified weight methods to the datasets X and Y, returning weight matrices for each method. @@ -786,83 +938,19 @@ learnTwasWeights <- function(X, Y, weightMethods, Xfiltered, Y, weightMethods, fittedModels ) - computeMethodWeights <- function(methodName, weightMethods) { - shortName <- sub("_weights$", "", methodName) - if (verbose >= 1) { - message(sprintf(" Fitting %s ...", shortName)) - tic() - } - - # Hardcoded vector of multivariate methods (accept both snake and camel). - # fSuSiE is multivariate (variants x features weight matrix) but is never - # refit here — fsusieWeights extracts from the supplied fsusieFit. - multivariateWeightMethods <- c("mrmash_weights", "mvsusie_weights", - "fsusie_weights", - "mrmashWeights", "mvsusieWeights", - "fsusieWeights") - args <- weightMethods[[methodName]] - fnName <- .resolveMethodFunction(methodName, args) - - # Only pass retainFit (or its legacy snake_case alias) to functions that accept it - if (retainFits) { - fnFormals <- names(formals(fnName)) - if ("retainFit" %in% fnFormals) { - args$retainFit <- TRUE - } else if ("retain_fit" %in% fnFormals) { - args$retain_fit <- TRUE - } - # Propagate the slim/full payload choice to producers that support it - # (mr.mash individual + RSS), unless the caller already set it per-method. - if ("fitDetail" %in% fnFormals && is.null(args$fitDetail)) { - args$fitDetail <- retainFitDetail - } - } - - methodFit <- NULL - if (methodName %in% multivariateWeightMethods) { - # Apply multivariate method - weightsMatrix <- if (verbose < 2) { - .quietEval(do.call(fnName, c(list(X = Xfiltered, Y = Y), args))) - } else { - do.call(fnName, c(list(X = Xfiltered, Y = Y), args)) - } - if (retainFits) methodFit <- attr(weightsMatrix, "fit") - if (nrow(weightsMatrix) != length(validColumns)) weightsMatrix <- weightsMatrix[names(validColumns), , drop = FALSE] - } else { - # Apply univariate method to each column of Y - # Initialize it with zeros to avoid NA - weightsMatrix <- matrix(0, nrow = ncol(Xfiltered), ncol = ncol(Y)) - - for (k in 1:ncol(Y)) { - weightsVector <- if (verbose < 2) { - .quietEval(do.call(fnName, c(list(X = Xfiltered, y = Y[, k]), args))) - } else { - do.call(fnName, c(list(X = Xfiltered, y = Y[, k]), args)) - } - if (retainFits && is.null(methodFit)) { - methodFit <- attr(weightsVector, "fit") - } - if (is.matrix(weightsVector)) weightsVector <- weightsVector[, k] - weightsMatrix[, k] <- weightsVector - } - } - - result <- .embedWeights(weightsMatrix, validColumns, ncol(X), ncol(Y), colnames(X), colnames(Y)) - if (!is.null(methodFit)) attr(result, "fit") <- methodFit - if (verbose >= 1) { - elapsed <- toc(quiet = TRUE) - message(sprintf(" Fitting %s done in %.1fs", shortName, elapsed$toc - elapsed$tic)) - } - return(result) - } + # Shared fit context threaded to the top-level per-method worker + entry builder. + ctx <- list(X = X, Y = Y, Xfiltered = Xfiltered, validColumns = validColumns, + study = study, context = context, trait = trait, + retainFits = retainFits, retainFitDetail = retainFitDetail, + standardized = standardized, dataType = dataType, verbose = verbose) if (numCores >= 2) { bpParam <- MulticoreParam(workers = numCores, RNGseed = 1L) weightsList <- bplapply(names(weightMethods), - computeMethodWeights, weightMethods, BPPARAM = bpParam) + .computeMethodWeights, weightMethods, ctx, BPPARAM = bpParam) } else { - weightsList <- names(weightMethods) %>% map(computeMethodWeights, weightMethods) + weightsList <- names(weightMethods) %>% map(.computeMethodWeights, weightMethods, ctx) } names(weightsList) <- names(weightMethods) @@ -884,58 +972,8 @@ learnTwasWeights <- function(X, Y, weightMethods, # carries the matrix for that method across outcomes via a single # `trait` value taken from the input `trait` arg (when length 1) or the # corresponding Y column name when `trait` matches `colnames(Y)`. - buildEntries <- function() { - studies <- character(0) - contexts <- character(0) - traits <- character(0) - methodsV <- character(0) - entries <- list() - for (m in names(weightsList)) { - wMat <- weightsList[[m]] - fitVal <- attr(wMat, "fit") - attr(wMat, "fit") <- NULL - shortMethod <- sub("(_weights|Weights)$", "", m) - # When trait/context were supplied per-row (length == ncol(Y)), emit - # one row per (method, outcome). Otherwise emit one row per method - # and carry the (possibly multi-column) weights matrix as-is. - perOutcome <- length(trait) == ncol(Y) && - length(context) %in% c(1L, ncol(Y)) - if (perOutcome) { - contextV <- if (length(context) == 1L) rep(context, ncol(Y)) else context - studyV <- if (length(study) == 1L) rep(study, ncol(Y)) else study - for (k in seq_len(ncol(Y))) { - studies <- c(studies, studyV[k]) - contexts <- c(contexts, contextV[k]) - traits <- c(traits, trait[k]) - methodsV <- c(methodsV, shortMethod) - entries[[length(entries) + 1L]] <- TwasWeightsEntry( - variantIds = variantIds, - weights = wMat[, k], - fits = if (retainFits) fitVal else NULL, - cvResult = NULL, - standardized = isTRUE(standardized), - dataType = dataType) - } - } else { - studies <- c(studies, study[1L]) - contexts <- c(contexts, context[1L]) - traits <- c(traits, trait[1L]) - methodsV <- c(methodsV, shortMethod) - wPayload <- if (ncol(wMat) == 1L) drop(wMat) else wMat - entries[[length(entries) + 1L]] <- TwasWeightsEntry( - variantIds = variantIds, - weights = wPayload, - fits = if (retainFits) fitVal else NULL, - cvResult = NULL, - standardized = isTRUE(standardized), - dataType = dataType) - } - } - list(study = studies, context = contexts, trait = traits, - method = methodsV, entry = entries) - } - rows <- buildEntries() + rows <- .buildTwasWeightEntries(weightsList, variantIds, ctx) TwasWeights( study = rows$study, context = rows$context, diff --git a/R/twasWeightsPipeline.R b/R/twasWeightsPipeline.R index d338cf27..a3e04780 100644 --- a/R/twasWeightsPipeline.R +++ b/R/twasWeightsPipeline.R @@ -829,6 +829,33 @@ combineTwasWeights <- function(..., ldSketch = NULL) { setGeneric("twasWeightsPipeline", function(data, ...) standardGeneric("twasWeightsPipeline")) +# Run the multivariate joint TWAS-weight fit over the (context, trait) grid for +# `traits`: dispatch each cis-region through the joint engine and merge per-region +# results. `marker` = the TwasJointPipeline config; `ctx` bundles the shared state +# (xRegions, data, norm, useCtx, fineMappingResult, dataDrivenPriorMatricesCv, +# cisWindow, verbose). +# @noRd +.twasRunMultivariateGrid <- function(traits, marker, ctx) { + synthSpec <- list(list(axes = c("context", "trait"), scope = NULL)) + labs <- vapply(ctx$xRegions, .twasRegionLabel, character(1)) + perRegion <- lapply(seq_along(ctx$xRegions), function(bi) { + .runJointSpecs(synthSpec, ctx$data, dataForm = "individual", pipeline = marker, + jointMethods = ctx$norm$tokens, contexts = ctx$useCtx, + traitIds = traits, + args = list(methodList = ctx$norm$methodList, + fineMappingResult = ctx$fineMappingResult, + dataDrivenPriorMatricesCv = ctx$dataDrivenPriorMatricesCv, + cisWindow = ctx$cisWindow, region = ctx$xRegions[[bi]], + regionIndex = bi, nRegions = length(ctx$xRegions), + verbose = ctx$verbose)) + }) + keep <- !vapply(perRegion, is.null, logical(1)) + perRegion <- perRegion[keep]; labs <- labs[keep] + if (length(perRegion) == 0L) return(NULL) + if (length(perRegion) == 1L) return(perRegion[[1L]]) + .twasMergeResultsByKey(perRegion, labs) +} + #' @rdname twasWeightsPipeline #' @export setMethod("twasWeightsPipeline", "QtlDataset", @@ -984,53 +1011,32 @@ setMethod("twasWeightsPipeline", "QtlDataset", # composed group per region -> per-method fit (the engine twas fitter) + # SR-TWAS ensemble layer, merged across regions. The SAME engine + fitter + # ensemble as every other multivariate path -- no separate fitting code. - runMultivariate <- function(traits) { - marker <- new("TwasJointPipeline", config = list( - cvFolds = cvFolds, samplePartition = samplePartition, - fitFullData = fitFullData, dataType = dataType, - retainFitDetail = retainFitDetail, standardized = FALSE, - ensemble = ensemble, ensembleR2Threshold = ensembleR2Threshold, - ensembleSolver = ensembleSolver, ensembleAlpha = ensembleAlpha, - maxCvVariants = maxCvVariants, cvThreads = cvThreads, - estimatePi = estimatePi, verbose = verbose, ldSketch = NULL)) - synthSpec <- list(list(axes = c("context", "trait"), scope = NULL)) - labs <- vapply(xRegions, .twasRegionLabel, character(1)) - perRegion <- lapply(seq_along(xRegions), function(bi) { - .runJointSpecs(synthSpec, data, dataForm = "individual", pipeline = marker, - jointMethods = norm$tokens, contexts = useCtx, - traitIds = traits, - args = list(methodList = norm$methodList, - fineMappingResult = fineMappingResult, - dataDrivenPriorMatricesCv = dataDrivenPriorMatricesCv, - cisWindow = cisWindow, region = xRegions[[bi]], - regionIndex = bi, nRegions = length(xRegions), - verbose = verbose)) - }) - keep <- !vapply(perRegion, is.null, logical(1)) - perRegion <- perRegion[keep]; labs <- labs[keep] - if (length(perRegion) == 0L) return(NULL) - if (length(perRegion) == 1L) return(perRegion[[1L]]) - .twasMergeResultsByKey(perRegion, labs) - } + # The joint-pipeline marker (config) + shared grid context are built once and + # used by both the multivariate grid (.twasRunMultivariateGrid) and the + # univariate engine path below. + marker <- new("TwasJointPipeline", config = list( + cvFolds = cvFolds, samplePartition = samplePartition, + fitFullData = fitFullData, dataType = dataType, + retainFitDetail = retainFitDetail, standardized = FALSE, + ensemble = ensemble, ensembleR2Threshold = ensembleR2Threshold, + ensembleSolver = ensembleSolver, ensembleAlpha = ensembleAlpha, + maxCvVariants = maxCvVariants, cvThreads = cvThreads, + estimatePi = estimatePi, verbose = verbose, ldSketch = NULL)) + twasGridCtx <- list(xRegions = xRegions, data = data, norm = norm, + useCtx = useCtx, fineMappingResult = fineMappingResult, + dataDrivenPriorMatricesCv = dataDrivenPriorMatricesCv, + cisWindow = cisWindow, verbose = verbose) # Top-level dispatch within the QtlDataset method body. if (multivariate) { # mvsusie / mr.mash: joint fit. If both nCtx == 1 and nTraits == 1 # we already rejected above via .twasCheckMultivariateY. - tw <- runMultivariate(allTraits) + tw <- .twasRunMultivariateGrid(allTraits, marker, twasGridCtx) } else { # Univariate methods ROUTED THROUGH THE ENGINE: one 1-condition group per # (context, trait), per region -> the SAME per-method fitter (+ ensemble # layer for >= 2 methods + resume cache) as the joint paths, merged across # regions. No separate per-(context, trait) fitting loop. - marker <- new("TwasJointPipeline", config = list( - cvFolds = cvFolds, samplePartition = samplePartition, - fitFullData = fitFullData, dataType = dataType, - retainFitDetail = retainFitDetail, standardized = FALSE, - ensemble = ensemble, ensembleR2Threshold = ensembleR2Threshold, - ensembleSolver = ensembleSolver, ensembleAlpha = ensembleAlpha, - maxCvVariants = maxCvVariants, cvThreads = cvThreads, - estimatePi = estimatePi, verbose = verbose, ldSketch = NULL)) univCell <- .lookupJointCell("univariate", "individual") scope <- list(studies = study, contexts = setNames(list(useCtx), study), @@ -1368,6 +1374,23 @@ setMethod("twasWeightsPipeline", "QtlSumStats", # rbind'd; the joint columns (when populated by either phase) are carried # through .rbindTwasWeights. +# Per-embedded-study TWAS-weights worker for .multiStudyPipelineDriver: recurse +# twasWeightsPipeline on one QtlDataset. `cfg` bundles the parent's forwarded args. +# @noRd +.twasPerStudy <- function(qd, cfg) do.call(twasWeightsPipeline, c(list( + data = qd, methods = cfg$methods, contexts = cfg$contexts, traitId = cfg$traitId, + region = cfg$region, cisWindow = cfg$cisWindow, jointRegions = cfg$jointRegions, + jointSpecification = NULL, fineMappingResult = cfg$fineMappingResult, + twasWeights = cfg$twasWeights, naAction = cfg$naAction, verbose = cfg$verbose), + cfg$dotArgs)) + +# Embedded-sumstats TWAS-weights worker for .multiStudyPipelineDriver. +# @noRd +.twasSumStats <- function(ss, cfg) do.call(twasWeightsPipeline, c(list( + data = ss, methods = cfg$methods, contexts = cfg$contexts, traitId = cfg$traitId, + jointSpecification = NULL, fineMappingResult = cfg$fineMappingResult, + twasWeights = cfg$twasWeights, verbose = cfg$verbose), cfg$dotArgs)) + #' @rdname twasWeightsPipeline #' @export setMethod("twasWeightsPipeline", "MultiStudyQtlDataset", @@ -1436,18 +1459,12 @@ setMethod("twasWeightsPipeline", "MultiStudyQtlDataset", } dotArgs <- list(...) - perStudyFn <- function(qd) do.call(twasWeightsPipeline, c(list( - data = qd, methods = methods, contexts = contexts, traitId = traitId, - region = region, cisWindow = cisWindow, jointRegions = jointRegions, - jointSpecification = NULL, fineMappingResult = fineMappingResult, - twasWeights = twasWeights, naAction = naAction, verbose = verbose), - dotArgs)) - sumStatsFn <- function(ss) do.call(twasWeightsPipeline, c(list( - data = ss, methods = methods, contexts = contexts, traitId = traitId, - jointSpecification = NULL, fineMappingResult = fineMappingResult, - twasWeights = twasWeights, verbose = verbose), dotArgs)) + cfg <- list(methods = methods, contexts = contexts, traitId = traitId, + region = region, cisWindow = cisWindow, jointRegions = jointRegions, + fineMappingResult = fineMappingResult, twasWeights = twasWeights, + naAction = naAction, verbose = verbose, dotArgs = dotArgs) .multiStudyPipelineDriver( - data, jointResult, perStudyFn, sumStatsFn, + data, jointResult, .twasPerStudy, .twasSumStats, cfg, .rbindTwasWeights, TwasWeights, "twasWeightsPipeline", noun = "weights") }) @@ -1552,6 +1569,16 @@ setMethod("twasWeightsPipeline", "ANY", zetaValid / zetaSum } +# Ensemble stacking objective (sum of squared residuals). `...` absorbs the +# gradient's extra optim args (PtP, Pty). +# @noRd +.ensembleObj <- function(z, Pvalid, yObs, ...) sum((yObs - Pvalid %*% z)^2) + +# Gradient of the ensemble stacking objective. `...` absorbs the objective's +# extra optim args (Pvalid, yObs). +# @noRd +.ensembleGrad <- function(z, PtP, Pty, ...) as.vector(2 * (PtP %*% z - Pty)) + # Solve ensemble stacking via L-BFGS-B (box-constrained optimization, then normalize). # Uses base R optim() with analytical gradient. No extra dependencies. # @param Pvalid Matrix of CV predictions for valid methods (n x Kvalid). @@ -1563,13 +1590,11 @@ setMethod("twasWeightsPipeline", "ANY", PtP <- crossprod(Pvalid) Pty <- as.vector(crossprod(Pvalid, yObs)) - fn <- function(z) sum((yObs - Pvalid %*% z)^2) - gr <- function(z) as.vector(2 * (PtP %*% z - Pty)) - fit <- tryCatch( optim( par = rep(1 / Kvalid, Kvalid), - fn = fn, gr = gr, + fn = .ensembleObj, gr = .ensembleGrad, + Pvalid = Pvalid, yObs = yObs, PtP = PtP, Pty = Pty, method = "L-BFGS-B", lower = rep(0, Kvalid) ), diff --git a/R/variantId.R b/R/variantId.R index 74f3c4a5..2ed007c6 100644 --- a/R/variantId.R +++ b/R/variantId.R @@ -288,6 +288,23 @@ variantIdToDf <- function(variantId) { parseVariantId(variantId) } +# Complement a DNA allele string (A<->T, C<->G) for strand flipping. +# @noRd +.strandFlip <- function(ref) chartr("ATCG", "TAGC", ref) + +# Ensure a data.frame has unique, non-empty column names (blanks become +# `unnamed_`, duplicates de-duplicated with make.unique). +# @noRd +.sanitizeNames <- function(df) { + nm <- colnames(df) + if (is.null(nm)) nm <- rep("unnamed", ncol(df)) + emptyIdx <- is.na(nm) | nm == "" + if (any(emptyIdx)) + nm[emptyIdx] <- paste0("unnamed_", seq_len(sum(emptyIdx))) + colnames(df) <- make.unique(nm, sep = "_") + df +} + #' Harmonize variant alleles against a reference #' #' The allele-harmonization engine for the package (used by summary-statistics @@ -343,18 +360,6 @@ harmonizeAlleles <- function(targetData, refVariants, colToFlip = NULL, removeStrandAmbiguous = TRUE, removeDups = FALSE, colToComplement = character(), ...) { - strandFlip <- function(ref) chartr("ATCG", "TAGC", ref) - - sanitizeNames <- function(df) { - nm <- colnames(df) - if (is.null(nm)) nm <- rep("unnamed", ncol(df)) - emptyIdx <- is.na(nm) | nm == "" - if (any(emptyIdx)) - nm[emptyIdx] <- paste0("unnamed_", seq_len(sum(emptyIdx))) - colnames(df) <- make.unique(nm, sep = "_") - df - } - if (is.data.frame(targetData)) { if (ncol(targetData) > 4 && all(c("chrom", "pos", "A2", "A1") %in% names(targetData))) { @@ -385,7 +390,7 @@ harmonizeAlleles <- function(targetData, refVariants, colToFlip = NULL, by = c("chrom", "pos"), suffix = c(".target", ".ref")) %>% as.data.frame() %>% - sanitizeNames() + .sanitizeNames() if (nrow(matchResult) == 0) { warning("No matching variants found between target data and reference variants.") @@ -401,8 +406,8 @@ harmonizeAlleles <- function(targetData, refVariants, colToFlip = NULL, mutate(variants_id_original = formatVariantId(chrom, pos, A2.target, A1.target), variants_id_qced = formatVariantId(chrom, pos, A2.ref, A1.ref)) %>% mutate(across(c(A1.target, A2.target, A1.ref, A2.ref), toupper)) %>% - mutate(flip1.ref = strandFlip(A1.ref), - flip2.ref = strandFlip(A2.ref)) %>% + mutate(flip1.ref = .strandFlip(A1.ref), + flip2.ref = .strandFlip(A2.ref)) %>% # AT / CG pairs cannot be distinguished from strand-flip without external # context; the keep rule below relies on this flag as a safety guard for # callers that may not have removed strand-ambiguous variants upstream. @@ -467,8 +472,8 @@ harmonizeAlleles <- function(targetData, refVariants, colToFlip = NULL, } if (flipStrand) { sIdx <- which(matchResult$strand_flip) - matchResult[sIdx, "A1.target"] <- strandFlip(matchResult[sIdx, "A1.target"]) - matchResult[sIdx, "A2.target"] <- strandFlip(matchResult[sIdx, "A2.target"]) + matchResult[sIdx, "A1.target"] <- .strandFlip(matchResult[sIdx, "A1.target"]) + matchResult[sIdx, "A2.target"] <- .strandFlip(matchResult[sIdx, "A2.target"]) } # Per-step QC counts (used by .runEntrySummaryStatsQc for "kept N of M @@ -544,6 +549,15 @@ harmonizeAlleles <- function(targetData, refVariants, colToFlip = NULL, out } +# Canonical per-id key for allele-aware matching: data.frame -> chrom/pos/A2/A1 +# formatted id (keeps allele order, so a ref/alt swap does NOT match); character +# vector -> normalized id. +# @noRd +.matchVariantKeyOf <- function(x) if (is.data.frame(x)) { + p <- parseVariantId(x) + formatVariantId(p$chrom, p$pos, p$A2, p$A1) +} else normalizeVariantId(x) + #' Allele-aware variant matcher (match by chrom/pos/ref/alt, not id string) #' #' The single matching primitive for the package: given two sets of variant @@ -584,11 +598,7 @@ matchVariants <- function(idsA, idsB, allowFlip = TRUE, # Exact-allele matching: canonicalize the id FORMAT (chr-prefix + separator, # rsID-safe) but keep allele order, then match by exact string identity, so # a ref/alt swap does NOT match (sign is always +1). - keyOf <- function(x) if (is.data.frame(x)) { - p <- parseVariantId(x) - formatVariantId(p$chrom, p$pos, p$A2, p$A1) - } else normalizeVariantId(x) - idxB <- match(keyOf(idsA), keyOf(idsB)) + idxB <- match(.matchVariantKeyOf(idsA), .matchVariantKeyOf(idsB)) matched <- which(!is.na(idxB)) return(list(idxA = matched, idxB = idxB[matched], sign = rep(1, length(matched)))) diff --git a/R/vcfWriter.R b/R/vcfWriter.R index 8941b4a1..aeaf5a2e 100644 --- a/R/vcfWriter.R +++ b/R/vcfWriter.R @@ -129,6 +129,21 @@ setMethod("writeSumstatsVcf", signature("FineMappingResultBase"), method = as.character(x$method)[r]) } +# Extract column `nm` from `df`, or an all-NA vector of length nSnps when absent. +# @noRd +.vcfCol <- function(df, nm, nSnps) { + if (!is.null(df) && nm %in% names(df)) df[[nm]] else rep(NA, nSnps) +} + +# Append one FORMAT field to the VCF geno accumulator `acc` (an environment +# holding geno + the header columns, plus nSnps for the matrix reshape). +# @noRd +.vcfAddGeno <- function(acc, name, vec, type, desc) { + acc$geno[[name]] <- matrix(vec, acc$nSnps) + acc$hdrRows <- c(acc$hdrRows, name); acc$hdrNum <- c(acc$hdrNum, "A") + acc$hdrType <- c(acc$hdrType, type); acc$hdrDesc <- c(acc$hdrDesc, desc) +} + # Internal worker: write one (study, context, trait, method) tuple to a # single VCF. When `splitByContext` / `splitByTrait` is in play the # output path is decorated with the corresponding tag(s) so multiple @@ -167,48 +182,43 @@ setMethod("writeSumstatsVcf", signature("FineMappingResultBase"), stop("writeSumstatsVcf: entry [", sn, "] has no variants to write") } nSnps <- nrow(base) - col <- function(df, nm) if (!is.null(df) && nm %in% names(df)) df[[nm]] - else rep(NA, nSnps) - geno <- list() - hdrRows <- character(0); hdrNum <- character(0) - hdrType <- character(0); hdrDesc <- character(0) - addGeno <- function(name, vec, type, desc) { - geno[[name]] <<- matrix(vec, nSnps) - hdrRows <<- c(hdrRows, name); hdrNum <<- c(hdrNum, "A") - hdrType <<- c(hdrType, type); hdrDesc <<- c(hdrDesc, desc) - } + acc <- new.env(parent = emptyenv()) + acc$nSnps <- nSnps + acc$geno <- list() + acc$hdrRows <- character(0); acc$hdrNum <- character(0) + acc$hdrType <- character(0); acc$hdrDesc <- character(0) # ES: posterior conditional effect (mvSuSiE/fSuSiE) when present, else the # marginal univariate beta (univariate susie). - es <- col(base, "conditional_effect") - if (all(is.na(es))) es <- col(m, "beta") + es <- .vcfCol(base, "conditional_effect", nSnps) + if (all(is.na(es))) es <- .vcfCol(m, "beta", nSnps) if (any(!is.na(es))) - addGeno("ES", es, "Float", + .vcfAddGeno(acc, "ES", es, "Float", "Effect size (posterior conditional effect, else marginal beta), effect allele") - se <- col(m, "se") + se <- .vcfCol(m, "se", nSnps) if (any(!is.na(se))) - addGeno("SE", se, "Float", "Standard error of the marginal effect-size estimate") - p <- col(m, "p") + .vcfAddGeno(acc, "SE", se, "Float", "Standard error of the marginal effect-size estimate") + p <- .vcfCol(m, "p", nSnps) if (any(!is.na(p))) { lp <- ifelse(is.na(p) | p <= 0, NA_real_, -log10(p)) - addGeno("LP", lp, "Float", "-log10 p-value of the marginal univariate effect") + .vcfAddGeno(acc, "LP", lp, "Float", "-log10 p-value of the marginal univariate effect") } - if (any(!is.na(col(m, "N")))) - addGeno("SS", as.integer(col(m, "N")), "Integer", "Sample size") - af <- col(base, "af"); if (all(is.na(af))) af <- col(m, "af") + if (any(!is.na(.vcfCol(m, "N", nSnps)))) + .vcfAddGeno(acc, "SS", as.integer(.vcfCol(m, "N", nSnps)), "Integer", "Sample size") + af <- .vcfCol(base, "af", nSnps); if (all(is.na(af))) af <- .vcfCol(m, "af", nSnps) if (any(!is.na(af))) - addGeno("AF", af, "Float", "Allele frequency (effect allele)") + .vcfAddGeno(acc, "AF", af, "Float", "Allele frequency (effect allele)") # Posterior fields, only when a posterior table is available. if (hasPost) { - pip <- col(base, "pip") + pip <- .vcfCol(base, "pip", nSnps) if (any(!is.na(pip))) - addGeno("PIP", pip, "Float", "Posterior inclusion probability") - lbf <- col(base, "logBF") + .vcfAddGeno(acc, "PIP", pip, "Float", "Posterior inclusion probability") + lbf <- .vcfCol(base, "logBF", nSnps) if (any(!is.na(lbf))) - addGeno("LBF", lbf, "Float", "Per-variant log Bayes factor (max single effect)") - lfsr <- col(base, "lfsr") + .vcfAddGeno(acc, "LBF", lbf, "Float", "Per-variant log Bayes factor (max single effect)") + lfsr <- .vcfCol(base, "lfsr", nSnps) if (any(!is.na(lfsr))) - addGeno("LFSR", lfsr, "Float", "Local false sign rate (per-condition posterior)") + .vcfAddGeno(acc, "LFSR", lfsr, "Float", "Local false sign rate (per-condition posterior)") # Credible sets are DYNAMIC: pecotmr does not assume any fixed coverage, so # we emit a CS (+ PUR) field for every cs_ # column the pipeline actually produced (e.g. cs_95 -> CS95 / PUR95). The @@ -218,13 +228,13 @@ setMethod("writeSumstatsVcf", signature("FineMappingResultBase"), cov <- sub("^cs_", "", cc) idx <- suppressWarnings(as.integer(sub(".*_", "", as.character(base[[cc]])))) idx[is.na(idx)] <- 0L - addGeno(paste0("CS", cov), idx, "Integer", + .vcfAddGeno(acc, paste0("CS", cov), idx, "Integer", sprintf("Credible-set index at %s%% coverage (0 = not captured)", cov)) pc <- paste0(cc, "_purity") if (pc %in% names(base)) { pur <- suppressWarnings(as.numeric(base[[pc]])) if (any(!is.na(pur))) - addGeno(paste0("PUR", cov), pur, "Float", + .vcfAddGeno(acc, paste0("PUR", cov), pur, "Float", sprintf("Purity (min abs corr) of the %s%% credible set", cov)) } } @@ -235,22 +245,22 @@ setMethod("writeSumstatsVcf", signature("FineMappingResultBase"), for (cc in grep("^(within_cs_pip|cs_logbf_|cs_effect_)", names(base), value = TRUE)) { v <- suppressWarnings(as.numeric(base[[cc]])) if (any(!is.na(v))) - addGeno(toupper(cc), v, "Float", + .vcfAddGeno(acc, toupper(cc), v, "Float", sprintf("Per-credible-set variant statistic (%s)", cc)) } } genoHeader <- DataFrame( - Number = hdrNum, Type = hdrType, Description = hdrDesc, - row.names = hdrRows) + Number = acc$hdrNum, Type = acc$hdrType, Description = acc$hdrDesc, + row.names = acc$hdrRows) .writeVcfImpl( - chrom = col(base, "chrom"), - pos = col(base, "pos"), - ref = col(base, "A2"), - alt = col(base, "A1"), - snpIds = col(base, "variant_id"), - geno = geno, + chrom = .vcfCol(base, "chrom", nSnps), + pos = .vcfCol(base, "pos", nSnps), + ref = .vcfCol(base, "A2", nSnps), + alt = .vcfCol(base, "A1", nSnps), + snpIds = .vcfCol(base, "variant_id", nSnps), + geno = acc$geno, genoHeader = genoHeader, sampleName = sn, outputPath = finalPath) diff --git a/man/overlapTopLoci.Rd b/man/dot-overlapPrefixNonKey.Rd similarity index 83% rename from man/overlapTopLoci.Rd rename to man/dot-overlapPrefixNonKey.Rd index 7e336750..2ee01590 100644 --- a/man/overlapTopLoci.Rd +++ b/man/dot-overlapPrefixNonKey.Rd @@ -1,31 +1,22 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/overlapTopLoci.R -\name{overlapTopLoci} -\alias{overlapTopLoci} -\alias{overlapTopLoci,QtlFineMappingResult,GwasFineMappingResult-method} +\name{.overlapPrefixNonKey} +\alias{.overlapPrefixNonKey} \title{Overlap QTL and GWAS top loci by allele-aware variant matching} \usage{ -overlapTopLoci(qtl, gwas, ...) - -\S4method{overlapTopLoci}{QtlFineMappingResult,GwasFineMappingResult}( - qtl, - gwas, - signalCutoff = 0.025, - type = c("data.frame", "GRanges"), - ... -) +.overlapPrefixNonKey(df, pfx, keyCols) } \arguments{ \item{qtl}{A \code{QtlFineMappingResult}.} \item{gwas}{A \code{GwasFineMappingResult}.} -\item{...}{Ignored.} - \item{signalCutoff}{PIP cutoff forwarded to \code{\link{getTopLoci}} for both inputs. Default 0.025.} \item{type}{\code{"data.frame"} (default) or \code{"GRanges"}.} + +\item{...}{Ignored.} } \value{ A \code{data.frame} (or \code{GRanges}) keyed on the QTL variant diff --git a/man/fineMappingPipeline.Rd b/man/fineMappingPipeline.Rd index 362c0d0b..a65185f7 100644 --- a/man/fineMappingPipeline.Rd +++ b/man/fineMappingPipeline.Rd @@ -499,6 +499,11 @@ deliberately not ported here: carries the bare token (\code{"susie"}, \code{"susieInf"}, \code{"mvsusie"}, ...) only. QC provenance is recorded on the sumstats' \code{qcInfo}. + \item An entry that \code{summaryStatsQc(pipCutoffToSkip = ...)} screened + out (recorded as \code{qcInfo$entryAudit[[i]]$pipScreenSkipped}, and + emptied to 0 variants) is \strong{skipped}, not fit: it produces no + row and a message with the screen reason. An all-screened collection + yields a valid empty result rather than an error. } } diff --git a/man/getP.Rd b/man/getP.Rd index c2b94cf5..561623d1 100644 --- a/man/getP.Rd +++ b/man/getP.Rd @@ -21,5 +21,5 @@ Numeric vector of p-values, or \code{NULL} if not available. Extract the association p-value vector from a \code{GwasSumStats} or \code{QtlSumStats} entry, selected by its identity tuple. Part of the first-class summary-statistic column set alongside - \code{\link{getZ}} / \code{\link{getBeta}} / \code{\link{getSE}}. + \code{\link{getZ}} / \code{\link{getBeta}} / \code{\link{getSe}}. } diff --git a/man/getSE.Rd b/man/getSe.Rd similarity index 81% rename from man/getSE.Rd rename to man/getSe.Rd index 0b374581..9b2c9e30 100644 --- a/man/getSE.Rd +++ b/man/getSe.Rd @@ -1,13 +1,13 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/AllGenerics.R, R/AllClasses.R -\name{getSE} -\alias{getSE} -\alias{getSE,SumStatsBase-method} +\name{getSe} +\alias{getSe} +\alias{getSe,SumStatsBase-method} \title{Get Effect-Size Standard Errors} \usage{ -getSE(x, ...) +getSe(x, ...) -\S4method{getSE}{SumStatsBase}(x, ...) +\S4method{getSe}{SumStatsBase}(x, ...) } \arguments{ \item{x}{A \code{GwasSumStats} or \code{QtlSumStats} object.} diff --git a/man/ldLoader.Rd b/man/ldLoader.Rd index 2b7146d4..003bbc71 100644 --- a/man/ldLoader.Rd +++ b/man/ldLoader.Rd @@ -41,8 +41,8 @@ default).} blocks larger than this to control memory usage.} } \value{ -A function \code{loader(g)} that, given a block index \code{g}, - returns the corresponding LD matrix or genotype matrix. +An \code{ldLoaderSpec} object (an opaque list describing the source). + Pass it with a block index to \code{\link{loadLdBlock}} to load one block. } \description{ Constructs a loader function that retrieves per-block LD matrices on @@ -69,8 +69,11 @@ Four modes are supported: # List mode with pre-computed LD R1 <- diag(10) R2 <- diag(15) -loader <- ldLoader(rList = list(R1, R2)) -loader(1) # returns R1 -loader(2) # returns R2 +spec <- ldLoader(rList = list(R1, R2)) +loadLdBlock(spec, 1) # returns R1 +loadLdBlock(spec, 2) # returns R2 } +\seealso{ +\code{\link{loadLdBlock}} +} diff --git a/man/loadLdBlock.Rd b/man/loadLdBlock.Rd new file mode 100644 index 00000000..1dbd5130 --- /dev/null +++ b/man/loadLdBlock.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ld.R +\name{loadLdBlock} +\alias{loadLdBlock} +\title{Load one LD block from an ldLoader spec} +\usage{ +loadLdBlock(spec, g) +} +\arguments{ +\item{spec}{An \code{ldLoaderSpec} object returned by \code{\link{ldLoader}}.} + +\item{g}{Integer block index (1-based).} +} +\value{ +The LD correlation matrix or genotype matrix for block \code{g}. +} +\description{ +Given an \code{ldLoaderSpec} (from \code{\link{ldLoader}}) and a block index +\code{g}, load the corresponding LD correlation matrix (or the genotype +matrix, in region mode with \code{returnGenotype = TRUE}). +} +\examples{ +spec <- ldLoader(rList = list(diag(10), diag(15))) +loadLdBlock(spec, 1) +} +\seealso{ +\code{\link{ldLoader}} +} diff --git a/pixi.toml b/pixi.toml index 1739ffb2..7461a2ed 100644 --- a/pixi.toml +++ b/pixi.toml @@ -1,10 +1,11 @@ [workspace] name = "r-pecotmr" channels = ["dnachun", "conda-forge", "bioconda"] -platforms = ["linux-64", "osx-arm64"] - -[system-requirements] -libc = { family="glibc", version="2.17" } +platforms = [ + { platform = "linux-64", glibc = "2.17" }, + { platform = "linux-aarch64", glibc = "2.17" }, + "osx-arm64" +] [tasks] devtools_document = "cd $GITHUB_WORKSPACE; R -e 'devtools::document()'" diff --git a/tests/testthat/test_QtlDataset.R b/tests/testthat/test_QtlDataset.R index 48909ed3..0de5f387 100644 --- a/tests/testthat/test_QtlDataset.R +++ b/tests/testthat/test_QtlDataset.R @@ -61,18 +61,18 @@ context("QtlDataset internal helpers") } # =========================================================================== -# .qtlResidualizeQR — pure linear algebra +# .qtlResidualizeQr — pure linear algebra # =========================================================================== -test_that(".qtlResidualizeQR: intercept-only residualization centers Y", { +test_that(".qtlResidualizeQr: intercept-only residualization centers Y", { set.seed(0) Y <- matrix(rnorm(20) + 5, nrow = 10, ncol = 2) - res <- pecotmr:::.qtlResidualizeQR(Y, C = NULL, scaleResiduals = FALSE) + res <- pecotmr:::.qtlResidualizeQr(Y, C = NULL, scaleResiduals = FALSE) # After removing the intercept, columns should have zero mean. expect_equal(unname(colMeans(res)), c(0, 0), tolerance = 1e-10) }) -test_that(".qtlResidualizeQR: covariate residualization removes the covariate signal", { +test_that(".qtlResidualizeQr: covariate residualization removes the covariate signal", { set.seed(1) n <- 50 C <- matrix(rnorm(n * 2), nrow = n, ncol = 2, @@ -80,7 +80,7 @@ test_that(".qtlResidualizeQR: covariate residualization removes the covariate si # Y = 0.5 * c1 - 0.3 * c2 + noise Y <- matrix(0.5 * C[, 1] - 0.3 * C[, 2] + rnorm(n, sd = 0.1), nrow = n, ncol = 1) - res <- pecotmr:::.qtlResidualizeQR(Y, C = C, scaleResiduals = FALSE) + res <- pecotmr:::.qtlResidualizeQr(Y, C = C, scaleResiduals = FALSE) # Residuals should be near-zero (only contain the noise). expect_lt(max(abs(res)), 0.5) # And uncorrelated with the covariates. @@ -88,22 +88,22 @@ test_that(".qtlResidualizeQR: covariate residualization removes the covariate si expect_lt(abs(cor(res[, 1], C[, 2])), 1e-8) }) -test_that(".qtlResidualizeQR: scaleResiduals = TRUE gives unit variance per column", { +test_that(".qtlResidualizeQr: scaleResiduals = TRUE gives unit variance per column", { set.seed(2) Y <- matrix(rnorm(30), nrow = 10, ncol = 3) - res <- pecotmr:::.qtlResidualizeQR(Y, C = NULL, scaleResiduals = TRUE) + res <- pecotmr:::.qtlResidualizeQr(Y, C = NULL, scaleResiduals = TRUE) sds <- apply(res, 2, sd) expect_equal(sds, c(1, 1, 1), tolerance = 1e-10) }) -test_that(".qtlResidualizeQR: constant residual columns survive the rescale step", { +test_that(".qtlResidualizeQr: constant residual columns survive the rescale step", { # Y is exactly its own mean -> residuals are 0, sd is 0 (and clamped to 1). Y <- matrix(5, nrow = 5, ncol = 1) - res <- pecotmr:::.qtlResidualizeQR(Y, C = NULL, scaleResiduals = TRUE) + res <- pecotmr:::.qtlResidualizeQr(Y, C = NULL, scaleResiduals = TRUE) expect_true(all(abs(res) < 1e-10)) }) -test_that(".qtlResidualizeQR: rank-deficient covariates are dropped by pivoted QR", { +test_that(".qtlResidualizeQr: rank-deficient covariates are dropped by pivoted QR", { set.seed(3) n <- 30 c1 <- rnorm(n) @@ -111,7 +111,7 @@ test_that(".qtlResidualizeQR: rank-deficient covariates are dropped by pivoted Q dimnames = list(NULL, c("a", "b", "c"))) Y <- matrix(rnorm(n), nrow = n, ncol = 1) # Even though `a` and `b` are collinear, the QR should not error. - expect_no_error(pecotmr:::.qtlResidualizeQR(Y, C = C, scaleResiduals = FALSE)) + expect_no_error(pecotmr:::.qtlResidualizeQr(Y, C = C, scaleResiduals = FALSE)) }) # =========================================================================== @@ -931,7 +931,7 @@ test_that("getResidualizedGenotypes: produces residualized matrix shape", { expect_equal(nrow(G), 12L) expect_equal(ncol(G), 6L) # When scaleResiduals = TRUE (the default), kept columns should have unit sd - # (constant columns are clamped to zero in .qtlResidualizeQR). + # (constant columns are clamped to zero in .qtlResidualizeQr). sds <- apply(G, 2L, sd) nonZero <- sds > 1e-6 expect_true(all(abs(sds[nonZero] - 1) < 1e-6)) diff --git a/tests/testthat/test_crossValidation.R b/tests/testthat/test_crossValidation.R index dcf37d11..474aef25 100644 --- a/tests/testthat/test_crossValidation.R +++ b/tests/testthat/test_crossValidation.R @@ -10,7 +10,7 @@ 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) { +mock_fit_fold <- function(Xtr, Ytr, j, ...) { list(weights = list(mock = matrix(1, ncol(Xtr), ncol(Ytr), dimnames = list(colnames(Xtr), NULL))), fits = list()) @@ -111,7 +111,7 @@ test_that("a degenerate fold (empty train/test) is skipped, not errored", { test_that("zero-variance predictions yield NA metrics with a message", { d <- mk_xy() - zero_fit <- function(Xtr, Ytr, j) { + zero_fit <- function(Xtr, Ytr, j, ...) { list(weights = list(mock = matrix(0, ncol(Xtr), ncol(Ytr), dimnames = list(colnames(Xtr), NULL))), fits = list()) @@ -134,7 +134,7 @@ test_that("the parallel fold path matches the serial one", { test_that("retainFits collects per-fold fits only when requested", { d <- mk_xy() - fit_with_model <- function(Xtr, Ytr, j) { + 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))) @@ -164,7 +164,7 @@ test_that("maxNumVariants subsamples from variantsToKeep when it already exceeds 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) + fit_with_null <- function(Xtr, Ytr, j, ...) list(weights = list( mock = matrix(1, ncol(Xtr), ncol(Ytr), dimnames = list(colnames(Xtr), NULL)), diff --git a/tests/testthat/test_fineMappingPipeline.R b/tests/testthat/test_fineMappingPipeline.R index 8f84346e..2504c735 100644 --- a/tests/testthat/test_fineMappingPipeline.R +++ b/tests/testthat/test_fineMappingPipeline.R @@ -477,16 +477,16 @@ test_that("combineFineMappingResults: row-binds same-class collections; rejects }) # =========================================================================== -# .fmExtractZN +# .fmExtractZn # =========================================================================== -test_that(".fmExtractZN: errors on missing SNP / Z / N columns", { +test_that(".fmExtractZn: errors on missing SNP / Z / N columns", { gr <- GenomicRanges::GRanges("chr1", IRanges::IRanges(100, 100)) - expect_error(pecotmr:::.fmExtractZN(gr, "x"), "no SNP mcol") + expect_error(pecotmr:::.fmExtractZn(gr, "x"), "no SNP mcol") S4Vectors::mcols(gr)$SNP <- "v1" - expect_error(pecotmr:::.fmExtractZN(gr, "x"), "no Z mcol") + expect_error(pecotmr:::.fmExtractZn(gr, "x"), "no Z mcol") S4Vectors::mcols(gr)$Z <- 1.0 - expect_error(pecotmr:::.fmExtractZN(gr, "x"), "no N mcol") + expect_error(pecotmr:::.fmExtractZn(gr, "x"), "no N mcol") }) # =========================================================================== diff --git a/tests/testthat/test_genotypeIo.R b/tests/testthat/test_genotypeIo.R index 97d36294..2a787c5c 100644 --- a/tests/testthat/test_genotypeIo.R +++ b/tests/testthat/test_genotypeIo.R @@ -392,12 +392,12 @@ test_that("readBim returns correct columns and types", { # =========================================================================== -# NoSNPsError / NoPhenotypeError custom conditions +# NoSnpsError / NoPhenotypeError custom conditions # =========================================================================== -test_that("NoSNPsError creates proper error condition", { - err <- NoSNPsError("test message") - expect_true(inherits(err, "NoSNPsError")) +test_that("NoSnpsError creates proper error condition", { + err <- NoSnpsError("test message") + expect_true(inherits(err, "NoSnpsError")) expect_true(inherits(err, "error")) expect_true(inherits(err, "condition")) expect_equal(err$message, "test message") @@ -1750,4 +1750,36 @@ test_that("loadGenotypeRegion warns on non-integer dosages without a sidecar", { "Non-integer genotype values detected") }) +# ---- fileIdx + snpInfo subsetting (memory-safe LD-sketch trimming) ----------- + +test_that(".withFileIdx + readers attach a sequential fileIdx to snpInfo", { + h <- readGenotypes(test_path("test_data/test_variants"), format = "plink2") + si <- getSnpInfo(h) + expect_true("fileIdx" %in% names(si)) + expect_identical(si$fileIdx, seq_len(nrow(si))) +}) + +test_that(".subsetGenotypeHandle keeps PLINK2 reads correct for kept variants", { + h <- readGenotypes(test_path("test_data/test_variants"), format = "plink2") + si <- getSnpInfo(h) + want <- c(2L, 4L, 5L, 9L) + dosFull <- pecotmr:::.dosageMatrix(h, want, meanImpute = TRUE) + # Keep a reordered superset so the wanted variants move to NEW row positions; + # fileIdx must rescue the (now-shifted) positional PLINK2 read. + keep <- c(9L, 5L, 4L, 2L, 7L, 1L) + hSub <- pecotmr:::.subsetGenotypeHandle(h, keep) + expect_equal(nrow(getSnpInfo(hSub)), length(keep)) + expect_identical(getSnpInfo(hSub)$fileIdx, keep) # original file order retained + wantSub <- match(want, keep) + dosSub <- pecotmr:::.dosageMatrix(hSub, wantSub, meanImpute = TRUE) + expect_equal(unname(dosFull), unname(dosSub)) +}) + +test_that(".subsetGenotypeHandle is NULL-safe and a no-op when nothing dropped", { + expect_null(pecotmr:::.subsetGenotypeHandle(NULL, TRUE)) + h <- readGenotypes(test_path("test_data/test_variants"), format = "plink2") + expect_identical( + pecotmr:::.subsetGenotypeHandle(h, seq_len(nrow(getSnpInfo(h)))), h) +}) + diff --git a/tests/testthat/test_jointEngine.R b/tests/testthat/test_jointEngine.R index c0cd5419..284e9974 100644 --- a/tests/testthat/test_jointEngine.R +++ b/tests/testthat/test_jointEngine.R @@ -644,7 +644,7 @@ test_that(".enumCrossContextIndividual: one group per trait in >= 2 contexts", { traits = list(S = c("G1", "G2"))) local_mocked_bindings( getStudy = function(data) "S", - .buildIndividualCrossContextXY = function(data, tid, scopedContexts, + .buildIndividualCrossContextXy = function(data, tid, scopedContexts, cisWindow, verbose, label, region = NULL) { if (tid == "G2") return(NULL) # skip branch (461) @@ -704,7 +704,7 @@ test_that(".enumCrossTraitIndividual: one group per context with >= 2 traits + p traits = list(S = c("G1", "G2"))) local_mocked_bindings( getStudy = function(data) "S", - .buildIndividualCrossTraitXY = function(data, cx, scopedTraits, cisWindow, + .buildIndividualCrossTraitXy = function(data, cx, scopedTraits, cisWindow, verbose, label, study, region = NULL) { if (cx == "c2") return(NULL) # skip branch (511) @@ -766,7 +766,7 @@ test_that(".enumComposedIndividual: one group joining every (context, trait) tup traits = list(S = c("gA", "gB"))) local_mocked_bindings( getStudy = function(data) "S", - .buildComposedIndividualXY = function(data, scope, study, cisWindow, + .buildComposedIndividualXy = function(data, scope, study, cisWindow, verbose, label, region = NULL) { Y <- matrix(0, 4, 3, dimnames = list(paste0("s", 1:4), c("c1:gA", "c1:gB", "c2:gA"))) @@ -782,7 +782,7 @@ test_that(".enumComposedIndividual: one group joining every (context, trait) tup test_that(".enumComposedIndividual: study not in scope / NULL xy -> empty", { local_mocked_bindings(getStudy = function(data) "S", - .buildComposedIndividualXY = function(...) NULL, .package = "pecotmr") + .buildComposedIndividualXy = function(...) NULL, .package = "pecotmr") expect_length(pecotmr:::.enumComposedIndividual( NULL, list(studies = "X")), 0L) # study not in scope expect_length(pecotmr:::.enumComposedIndividual( @@ -1109,8 +1109,8 @@ test_that(".twasFmHandoffCv: a token absent from the FM CV predictions -> NULL", expect_null(pecotmr:::.twasFmHandoffCv(NULL, "mvsusie")) }) -test_that("construct: empty rows -> NULL for both pipelines", { - empty <- pecotmr:::.jointRows() +test_that("construct: empty records -> NULL for both pipelines", { + empty <- list() expect_null(pecotmr:::construct(new("FmJointPipeline", config = list()), empty)) expect_null(pecotmr:::construct(new("TwasJointPipeline", config = list()), empty)) }) diff --git a/tests/testthat/test_jointSpecification.R b/tests/testthat/test_jointSpecification.R index 007b1437..ab00f9c3 100644 --- a/tests/testthat/test_jointSpecification.R +++ b/tests/testthat/test_jointSpecification.R @@ -1239,19 +1239,19 @@ test_that(".buildJointSumstatZMatrix: a mismatched SNP order across entries erro # Individual X/Y builders: the skip paths (return NULL / message branches) # ----------------------------------------------------------------------------- -test_that(".buildIndividualCrossContextXY: skips when a trait spans < 2 contexts", { +test_that(".buildIndividualCrossContextXy: skips when a trait spans < 2 contexts", { se1 <- .js_makeSe(traits = "g1") # g1 present se0 <- .js_makeSe(traits = "other") # g1 absent local_mocked_bindings( getPhenotypes = function(data, contexts) if (identical(contexts, "c1")) se1 else se0, .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildIndividualCrossContextXY( + expect_null(suppressMessages(pecotmr:::.buildIndividualCrossContextXy( NULL, "g1", c("c1", "c2"), cisWindow = 1000L, verbose = 1, label = "X"))) }) -test_that(".buildIndividualCrossContextXY: region path + complete-case skip", { +test_that(".buildIndividualCrossContextXy: region path + complete-case skip", { se <- .js_makeSe(traits = "g1", samples = paste0("s", 1:6)) samp <- paste0("s", 1:6) region <- GenomicRanges::GRanges("chr1", IRanges::IRanges(1, 10000)) @@ -1266,12 +1266,12 @@ test_that(".buildIndividualCrossContextXY: region path + complete-case skip", { c2 = ym(rnorm(6))) }, .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildIndividualCrossContextXY( + expect_null(suppressMessages(pecotmr:::.buildIndividualCrossContextXy( NULL, "g1", c("c1", "c2"), cisWindow = NULL, verbose = 1, label = "X", region = region))) }) -test_that(".buildIndividualCrossContextXY: too few shared samples skips", { +test_that(".buildIndividualCrossContextXy: too few shared samples skips", { se <- .js_makeSe(traits = "g1", samples = paste0("s", 1:6)) local_mocked_bindings( getPhenotypes = function(data, contexts) se, @@ -1282,7 +1282,7 @@ test_that(".buildIndividualCrossContextXY: too few shared samples skips", { list(c1 = matrix(0, 6, 1, dimnames = list(paste0("s", 1:6), "g1")), c2 = matrix(0, 6, 1, dimnames = list(paste0("s", 1:6), "g1"))), .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildIndividualCrossContextXY( + expect_null(suppressMessages(pecotmr:::.buildIndividualCrossContextXy( NULL, "g1", c("c1", "c2"), cisWindow = 1000L, verbose = 1, label = "X"))) }) @@ -1294,13 +1294,13 @@ test_that(".fmTraitsInRegion: filters traits by phenotype overlap with the regio c("g1", "g2")) # NULL region -> unchanged }) -test_that(".buildIndividualCrossTraitXY: skip branches (< 2 traits, region, complete)", { +test_that(".buildIndividualCrossTraitXy: skip branches (< 2 traits, region, complete)", { se2 <- .js_makeSe(traits = c("g1", "g2"), samples = paste0("s", 1:6)) samp <- paste0("s", 1:6) # < 2 scoped traits in the context -> NULL. local_mocked_bindings(getPhenotypes = function(data, contexts) se2, .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildIndividualCrossTraitXY( + expect_null(suppressMessages(pecotmr:::.buildIndividualCrossTraitXy( NULL, "cx", "g1", cisWindow = 1000L, verbose = 1, label = "X", study = "S"))) # region path + < 2 complete cases. @@ -1313,13 +1313,13 @@ test_that(".buildIndividualCrossTraitXY: skip branches (< 2 traits, region, comp cbind(g1 = c(NA, NA, NA, NA, NA, 1), g2 = rnorm(6)) |> `rownames<-`(samp), .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildIndividualCrossTraitXY( + expect_null(suppressMessages(pecotmr:::.buildIndividualCrossTraitXy( NULL, "cx", c("g1", "g2"), cisWindow = NULL, verbose = 1, label = "X", study = "S", region = GenomicRanges::GRanges("chr1", IRanges::IRanges(1, 9999))))) }) -test_that(".buildComposedIndividualXY: skip branches and single-context wrap", { +test_that(".buildComposedIndividualXy: skip branches and single-context wrap", { se <- .js_makeSe(traits = c("g1", "g2"), samples = paste0("s", 1:6)) samp <- paste0("s", 1:6) scope <- list(contexts = list(S = "c1"), traits = list(S = c("g1", "g2"))) @@ -1332,7 +1332,7 @@ test_that(".buildComposedIndividualXY: skip branches and single-context wrap", { .fmResidPheno = function(x, contexts, traitId = NULL, ...) matrix(rnorm(12), 6, 2, dimnames = list(samp, c("g1", "g2"))), .package = "pecotmr") - out <- suppressMessages(pecotmr:::.buildComposedIndividualXY( + out <- suppressMessages(pecotmr:::.buildComposedIndividualXy( NULL, scope, "S", cisWindow = 1000L, verbose = 1, label = "X")) expect_equal(ncol(out$Y), 2L) expect_setequal(colnames(out$Y), c("c1:g1", "c1:g2")) @@ -1340,7 +1340,7 @@ test_that(".buildComposedIndividualXY: skip branches and single-context wrap", { scope1 <- list(contexts = list(S = "c1"), traits = list(S = "g1")) local_mocked_bindings(getPhenotypes = function(data, contexts) .js_makeSe(traits = "g1"), .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXY( + expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXy( NULL, scope1, "S", cisWindow = 1000L, verbose = 1, label = "X"))) }) @@ -1538,7 +1538,7 @@ test_that("validateMethodsVsJointSpec: per-study methods with a context joint pa mp, list(list(axes = "context")))) # 542 }) -test_that(".buildIndividualCrossTraitXY: disjoint X/Y samples skip the context", { +test_that(".buildIndividualCrossTraitXy: disjoint X/Y samples skip the context", { se <- .js_makeSe(traits = c("g1", "g2"), samples = paste0("s", 1:6)) local_mocked_bindings( getPhenotypes = function(data, contexts) se, @@ -1548,12 +1548,12 @@ test_that(".buildIndividualCrossTraitXY: disjoint X/Y samples skip the context", .fmResidPheno = function(x, contexts, traitId = NULL, ...) matrix(0, 6, 2, dimnames = list(paste0("s", 1:6), c("g1", "g2"))), .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildIndividualCrossTraitXY( + expect_null(suppressMessages(pecotmr:::.buildIndividualCrossTraitXy( NULL, "cx", c("g1", "g2"), cisWindow = 1000L, verbose = 1, label = "X", study = "S"))) # 730 }) -test_that(".buildComposedIndividualXY: disjoint samples / missing trait col / NA rows skip", { +test_that(".buildComposedIndividualXy: disjoint samples / missing trait col / NA rows skip", { samp <- paste0("s", 1:6) se <- .js_makeSe(traits = c("g1", "g2"), samples = samp) scope <- list(contexts = list(S = c("c1", "c2")), @@ -1569,7 +1569,7 @@ test_that(".buildComposedIndividualXY: disjoint samples / missing trait col / NA matrix(rnorm(12), 6, 2, dimnames = list(samp, c("g1", "g2")))), c("c1", "c2")), .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXY( + expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXy( NULL, scope, "S", cisWindow = 1000L, verbose = 1, label = "X"))) # 774 # (b) one context's Y lacks the trait column -> tuple skipped -> < 2 yCols (779/785). local_mocked_bindings( @@ -1581,7 +1581,7 @@ test_that(".buildComposedIndividualXY: disjoint samples / missing trait col / NA list(c1 = matrix(0, 6, 1, dimnames = list(samp, "g1")), c2 = matrix(0, 6, 1, dimnames = list(samp, "zzz"))), # no g1/g2 .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXY( + expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXy( NULL, list(contexts = list(S = c("c1", "c2")), traits = list(S = "g1")), "S", cisWindow = 1000L, verbose = 1, label = "X"))) # 779/785 @@ -1597,7 +1597,7 @@ test_that(".buildComposedIndividualXY: disjoint samples / missing trait col / NA c2 = matrix(c(NA, NA, NA, NA, NA, 1), 6, 1, dimnames = list(samp, "g1"))), .package = "pecotmr") - expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXY( + expect_null(suppressMessages(pecotmr:::.buildComposedIndividualXy( NULL, list(contexts = list(S = c("c1", "c2")), traits = list(S = "g1")), "S", cisWindow = 1000L, verbose = 1, label = "X"))) # 788 diff --git a/tests/testthat/test_ld.R b/tests/testthat/test_ld.R index f4883698..5f6aaa38 100644 --- a/tests/testthat/test_ld.R +++ b/tests/testthat/test_ld.R @@ -1644,18 +1644,18 @@ test_that("ldLoader errors when multiple sources are provided", { # ldLoader: R_list branch # =========================================================================== -test_that("ldLoader with R_list returns a function", { +test_that("ldLoader with rList returns an ldLoaderSpec", { R <- list(matrix(c(1, 0.5, 0.5, 1), 2, 2)) - loader <- ldLoader(rList = R) - expect_type(loader, "closure") + spec <- ldLoader(rList = R) + expect_s3_class(spec, "ldLoaderSpec") }) test_that("ldLoader R_list returns correct matrix", { R1 <- matrix(c(1, 0.3, 0.3, 1), 2, 2) R2 <- matrix(c(1, 0.8, 0.8, 1), 2, 2) loader <- ldLoader(rList = list(R1, R2)) - expect_equal(loader(1), R1) - expect_equal(loader(2), R2) + expect_equal(loadLdBlock(loader, 1), R1) + expect_equal(loadLdBlock(loader, 2), R2) }) test_that("ldLoader R_list with max_variants downsamples", { @@ -1663,7 +1663,7 @@ test_that("ldLoader R_list with max_variants downsamples", { R <- matrix(0.1, 10, 10) diag(R) <- 1 loader <- ldLoader(rList = list(R), maxVariants = 5) - result <- loader(1) + result <- loadLdBlock(loader, 1) expect_equal(nrow(result), 5) expect_equal(ncol(result), 5) }) @@ -1672,7 +1672,7 @@ test_that("ldLoader R_list without max_variants returns full matrix", { R <- matrix(0.1, 10, 10) diag(R) <- 1 loader <- ldLoader(rList = list(R)) - result <- loader(1) + result <- loadLdBlock(loader, 1) expect_equal(nrow(result), 10) }) @@ -1680,7 +1680,7 @@ test_that("ldLoader R_list max_variants larger than matrix returns full matrix", R <- matrix(0.1, 3, 3) diag(R) <- 1 loader <- ldLoader(rList = list(R), maxVariants = 100) - result <- loader(1) + result <- loadLdBlock(loader, 1) expect_equal(nrow(result), 3) }) @@ -1688,25 +1688,25 @@ test_that("ldLoader R_list max_variants larger than matrix returns full matrix", # ldLoader: X_list branch # =========================================================================== -test_that("ldLoader with X_list returns a function", { +test_that("ldLoader with xList returns an ldLoaderSpec", { X <- list(matrix(rnorm(30), 10, 3)) - loader <- ldLoader(xList = X) - expect_type(loader, "closure") + spec <- ldLoader(xList = X) + expect_s3_class(spec, "ldLoaderSpec") }) test_that("ldLoader X_list returns correct matrix", { X1 <- matrix(1:12, 4, 3) X2 <- matrix(1:8, 4, 2) loader <- ldLoader(xList = list(X1, X2)) - expect_equal(loader(1), X1) - expect_equal(loader(2), X2) + expect_equal(loadLdBlock(loader, 1), X1) + expect_equal(loadLdBlock(loader, 2), X2) }) test_that("ldLoader X_list with max_variants downsamples columns", { set.seed(42) X <- matrix(rnorm(50), 10, 5) loader <- ldLoader(xList = list(X), maxVariants = 3) - result <- loader(1) + result <- loadLdBlock(loader, 1) expect_equal(nrow(result), 10) expect_equal(ncol(result), 3) }) @@ -1714,7 +1714,7 @@ test_that("ldLoader X_list with max_variants downsamples columns", { test_that("ldLoader X_list max_variants larger than ncol returns full matrix", { X <- matrix(rnorm(12), 4, 3) loader <- ldLoader(xList = list(X), maxVariants = 100) - result <- loader(1) + result <- loadLdBlock(loader, 1) expect_equal(ncol(result), 3) }) @@ -1757,7 +1757,7 @@ test_that("ldLoader ldInfo loads LD from PLINK2 files", { skip_if_not_installed("pgenlibr") plink_prefix <- file.path(test_data_dir, "test_variants") loader <- ldLoader(ldInfo = data.frame(LD_file = plink_prefix)) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_true(is.matrix(mat)) expect_equal(nrow(mat), 349L) expect_equal(ncol(mat), 349L) @@ -1769,7 +1769,7 @@ test_that("ldLoader ldInfo loads LD from VCF file", { skip_if_not_installed("VariantAnnotation") vcf_path <- file.path(test_data_dir, "test_variants.vcf.gz") loader <- ldLoader(ldInfo = data.frame(LD_file = vcf_path)) - mat <- suppressWarnings(loader(1)) + mat <- suppressWarnings(loadLdBlock(loader, 1)) expect_true(is.matrix(mat)) expect_equal(nrow(mat), 349L) expect_true(isSymmetric(mat)) @@ -1780,7 +1780,7 @@ test_that("ldLoader ldInfo loads LD from GDS file", { skip_if_not_installed("gdsfmt") gds_path <- file.path(test_data_dir, "test_variants.gds") loader <- ldLoader(ldInfo = data.frame(LD_file = gds_path)) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_true(is.matrix(mat)) expect_equal(nrow(mat), 349L) expect_true(isSymmetric(mat)) @@ -1790,7 +1790,7 @@ test_that("ldLoader ldInfo loads LD from PLINK1 files", { skip_if_not_installed("snpStats") plink1_prefix <- file.path(test_data_dir, "protocol_example.genotype") loader <- ldLoader(ldInfo = data.frame(LD_file = plink1_prefix)) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_true(is.matrix(mat)) expect_true(isSymmetric(mat)) }) @@ -1822,7 +1822,7 @@ test_that("ldLoader ldInfo loads pre-computed .cor.xz blocks", { ) loader <- ldLoader(ldInfo = data.frame(LD_file = ld_file, SNP_file = bim_file)) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_true(is.matrix(mat)) expect_true(isSymmetric(mat)) expect_true(nrow(mat) > 0) @@ -1833,7 +1833,7 @@ test_that("ldLoader ldInfo with max_variants subsamples", { plink_prefix <- file.path(test_data_dir, "test_variants") set.seed(42) loader <- ldLoader(ldInfo = data.frame(LD_file = plink_prefix), maxVariants = 20) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_equal(nrow(mat), 20L) expect_equal(ncol(mat), 20L) }) @@ -1846,8 +1846,8 @@ test_that("ldLoader ldInfo returns consistent LD across formats", { gds_path <- file.path(test_data_dir, "test_variants.gds") loader_plink <- ldLoader(ldInfo = data.frame(LD_file = plink_prefix)) loader_gds <- ldLoader(ldInfo = data.frame(LD_file = gds_path)) - mat_plink <- loader_plink(1) - mat_gds <- loader_gds(1) + mat_plink <- loadLdBlock(loader_plink, 1) + mat_gds <- loadLdBlock(loader_gds, 1) expect_equal(dim(mat_plink), dim(mat_gds)) }) @@ -1867,7 +1867,7 @@ test_that("ldLoader ld_meta_path loads LD from PLINK2 metadata", { file = meta_file, append = TRUE) region <- "chr21:17513228-17592874" loader <- ldLoader(ldMetaPath = meta_file, regions = region) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_true(is.matrix(mat)) expect_equal(nrow(mat), 349L) expect_equal(ncol(mat), 349L) @@ -1885,7 +1885,7 @@ test_that("ldLoader ld_meta_path loads LD from VCF metadata", { file = meta_file, append = TRUE) region <- "chr21:17513228-17592874" loader <- ldLoader(ldMetaPath = meta_file, regions = region) - mat <- suppressWarnings(loader(1)) + mat <- suppressWarnings(loadLdBlock(loader, 1)) expect_true(is.matrix(mat)) expect_equal(nrow(mat), 349L) }) @@ -3318,7 +3318,7 @@ test_that("ldLoader region mode subsamples and scales genotype matrices", { set.seed(1) loader <- ldLoader(ldMetaPath = meta_file, regions = geno_region_all, returnGenotype = TRUE, maxVariants = 20) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_equal(ncol(mat), 20L) expect_equal(nrow(mat), 100L) # samples expect_true(all(is.finite(mat))) # scaled, NAs replaced with 0 @@ -3334,7 +3334,7 @@ test_that("ldLoader region mode subsamples a correlation matrix", { set.seed(2) loader <- ldLoader(ldMetaPath = meta_file, regions = geno_region_all, returnGenotype = FALSE, maxVariants = 15) - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_equal(dim(mat), c(15L, 15L)) }) @@ -3361,7 +3361,7 @@ test_that("ldLoader ldInfo auto-detects companion file when SNP_file is absent", .package = "pecotmr" ) loader <- ldLoader(ldInfo = data.frame(LD_file = ld_file)) # no SNP_file column - mat <- loader(1) + mat <- loadLdBlock(loader, 1) expect_true(is.matrix(mat)) expect_true(isSymmetric(mat)) expect_true(nrow(mat) > 0) diff --git a/tests/testthat/test_qtlSumStats.R b/tests/testthat/test_qtlSumStats.R index 7de3aeff..918f0f6e 100644 --- a/tests/testthat/test_qtlSumStats.R +++ b/tests/testthat/test_qtlSumStats.R @@ -403,10 +403,10 @@ test_that("QtlSumStats: a non-GenotypeHandle ldSketch is rejected", { }) # =========================================================================== -# First-class summary-statistic accessors: getP / getBeta / getSE +# First-class summary-statistic accessors: getP / getBeta / getSe # =========================================================================== -test_that("getP / getBeta / getSE read the optional P/BETA/SE mcols", { +test_that("getP / getBeta / getSe read the optional P/BETA/SE mcols", { gr <- .qtlMakeEntryGr(4) S4Vectors::mcols(gr)$P <- c(0.1, 0.01, 1e-4, 0.5) S4Vectors::mcols(gr)$BETA <- c(0.2, -0.3, 0.4, 0.05) @@ -415,14 +415,14 @@ test_that("getP / getBeta / getSE read the optional P/BETA/SE mcols", { entry = list(gr), genome = "hg19") expect_equal(getP(x), c(0.1, 0.01, 1e-4, 0.5)) expect_equal(getBeta(x), c(0.2, -0.3, 0.4, 0.05)) - expect_equal(getSE(x), rep(0.1, 4)) + expect_equal(getSe(x), rep(0.1, 4)) }) -test_that("getP / getBeta / getSE return NULL when the entry omits them", { +test_that("getP / getBeta / getSe return NULL when the entry omits them", { y <- .qtlMakeOne(n = 4) # Z/N-only entry expect_null(getP(y)) expect_null(getBeta(y)) - expect_null(getSE(y)) + expect_null(getSe(y)) expect_equal(getZ(y), seq(1.0, by = 0.5, length.out = 4)) # Z still works }) diff --git a/tests/testthat/test_sumstatsQc.R b/tests/testthat/test_sumstatsQc.R index f31746ff..93dc2c4d 100644 --- a/tests/testthat/test_sumstatsQc.R +++ b/tests/testthat/test_sumstatsQc.R @@ -5502,3 +5502,38 @@ test_that("summaryStatsQc kriging QC sign-flips an LD-inconsistent variant and r # the audit's flip count equals the diagnostics rows marked flipped. expect_identical(sum(ea$krigingDiagnostics$flipped), ea$krigingFlipped) }) + +# ---- LD-sketch trimming to the summary-stats range / final variants ---------- + +test_that(".subsetSketchToRange keeps only panel variants in the entries' per-chrom span", { + h <- readGenotypes(test_path("test_data/test_variants"), format = "plink2") + si <- getSnpInfo(h) + ch <- names(sort(table(si$CHR), decreasing = TRUE))[1] # busiest chromosome + chIdx <- which(si$CHR == ch) + slice <- si[chIdx[seq_len(min(30L, length(chIdx)))], , drop = FALSE] + gr <- GenomicRanges::GRanges(paste0("chr", slice$CHR), + IRanges::IRanges(as.integer(slice$BP), width = 1L)) + S4Vectors::mcols(gr)$SNP <- slice$SNP + sub <- getSnpInfo(pecotmr:::.subsetSketchToRange(h, list(gr))) + lo <- min(slice$BP); hi <- max(slice$BP) + expected <- si$SNP[pecotmr:::canonChrom(si$CHR) == pecotmr:::canonChrom(ch) & + as.integer(si$BP) >= lo & as.integer(si$BP) <= hi] + expect_setequal(sub$SNP, expected) + expect_lt(nrow(sub), nrow(si)) +}) + +test_that(".subsetSketchToIds keeps exactly the entries' variants", { + h <- readGenotypes(test_path("test_data/test_variants"), format = "plink2") + si <- getSnpInfo(h) + sel <- c(3L, 10L, 40L, 200L) + gr <- GenomicRanges::GRanges(paste0("chr", si$CHR[sel]), + IRanges::IRanges(as.integer(si$BP[sel]), width = 1L)) + S4Vectors::mcols(gr)$SNP <- si$SNP[sel] + sub <- getSnpInfo(pecotmr:::.subsetSketchToIds(h, list(gr))) + expect_setequal(normalizeVariantId(sub$SNP), normalizeVariantId(si$SNP[sel])) +}) + +test_that(".subsetSketchToRange / .subsetSketchToIds are NULL-safe", { + expect_null(pecotmr:::.subsetSketchToRange(NULL, list())) + expect_null(pecotmr:::.subsetSketchToIds(NULL, list())) +})