fix snpInfo bug and unnest

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)

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

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

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") {

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 = <named list: methodKey -> (variants x outcomes) weight matrix,
 #                   rownames indexing colnames(Xtrain)>,
 #        fits    = <named list: methodKey -> 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")

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
 }