add fine-mapping CV support and other features

NAMESPACE

@@ -83,7 +83,9 @@ export(fitSusieInfThenSusieRss)
 export(fixupExampleGenotypePaths)
 export(formatFinemappingOutput)
 export(fsusieGetCs)
+export(fsusieWeights)
 export(fsusieWrapper)
+export(getAf)
 export(getAnnotationMeta)
 export(getAnnotations)
 export(getBaseline)
@@ -95,6 +97,7 @@ export(getCorrelation)
 export(getCs)
 export(getCtwasMetaData)
 export(getCvPerformance)
+export(getCvResult)
 export(getDataType)
 export(getEigenList)
 export(getEnrichment)
@@ -273,6 +276,7 @@ exportMethods(colocboostPipeline)
 exportMethods(computeLdScores)
 exportMethods(estimateH2)
 exportMethods(fineMappingPipeline)
+exportMethods(getAf)
 exportMethods(getAnnotationMeta)
 exportMethods(getAnnotations)
 exportMethods(getBlockMetadata)
@@ -281,6 +285,7 @@ exportMethods(getContexts)
 exportMethods(getCorrelation)
 exportMethods(getCs)
 exportMethods(getCvPerformance)
+exportMethods(getCvResult)
 exportMethods(getDataType)
 exportMethods(getEigenList)
 exportMethods(getEnrichment)

R/AllGenerics.R

@@ -142,6 +142,20 @@ setGeneric("getN", function(x, ...) standardGeneric("getN"))
 #' @export
 setGeneric("getMaf", function(x, ...) standardGeneric("getMaf"))
 
+#' @title Get Effect-Allele Frequencies
+#' @description Extract the directional effect-allele (A1) frequency vector
+#'   for a \code{QtlDataset}. Unlike \code{\link{getMaf}}, the value is
+#'   \emph{not} folded to the minor allele: it is the frequency of the
+#'   dosage-counted allele (A1, the effect allele), matching the allele the
+#'   marginal effect sizes and the fine-mapping \code{af} column report.
+#' @param x A \code{QtlDataset} object.
+#' @param ... Class-specific selection arguments (e.g., \code{traitId},
+#'   \code{region}, \code{cisWindow}, \code{samples} for \code{QtlDataset}).
+#' @return Named numeric vector of effect-allele frequencies (names are
+#'   variant IDs), or an empty vector when no variants are selected.
+#' @export
+setGeneric("getAf", function(x, ...) standardGeneric("getAf"))
+
 #' @title Get Number of SNPs
 #' @description Number of SNPs in a \code{GwasSumStats} or
 #'   \code{QtlSumStats} entry, selected by its identity tuple.
@@ -397,6 +411,21 @@ setGeneric("getPip", function(x, ...) standardGeneric("getPip"))
 #' @export
 setGeneric("getSusieFit", function(x, ...) standardGeneric("getSusieFit"))
 
+#' @title Get Cross-Validation Result
+#' @description Extract the cross-validation payload stored on a
+#'   \code{FineMappingEntry} (or the matching entry of a
+#'   \code{FineMappingResult}). The payload is a list with components
+#'   \code{samplePartition} (a \code{data.frame} of \code{Sample}/\code{Fold}
+#'   assignments), \code{predictions} (a named list of per-method out-of-fold
+#'   prediction matrices), and \code{performance} (a named list of per-method
+#'   metric matrices). \code{NULL} when fine-mapping was run without
+#'   cross-validation (\code{cvFolds <= 1}).
+#' @param x A \code{FineMappingEntry} or \code{FineMappingResult}.
+#' @param ... Class-specific selection arguments.
+#' @return A list (the CV payload) or \code{NULL}.
+#' @export
+setGeneric("getCvResult", function(x, ...) standardGeneric("getCvResult"))
+
 #' @title Get Marginal Effects
 #' @description Extract per-variant marginal univariate effects from a
 #'   fine-mapping entry or result. Returns a \code{data.frame} with

R/FineMappingEntry.R

@@ -29,11 +29,15 @@ setClass("FineMappingEntry",
   representation(
     variantIds = "character",
     susieFit   = "ANY",
-    topLoci    = "data.frame"
+    topLoci    = "data.frame",
+    cvResult   = "ANY"
   ),
   validity = function(object) {
     errors <- character()
     n <- length(object@variantIds)
+    if (!is.null(object@cvResult) && !is.list(object@cvResult))
+      errors <- c(errors,
+        "cvResult must be NULL or a list (samplePartition/predictions/performance)")
     if (nrow(object@topLoci) > 0L) {
       # Minimal contract: variant_id + pip. Canonical projector columns
       # (marginal_*, posterior_*, etc.) are pipeline-populated; tests
@@ -85,13 +89,17 @@ setClass("FineMappingEntry",
 #'   (\code{pip, posterior_mean, posterior_sd, cs_*, cs_*_purity}),
 #'   pipeline stamps (\code{method, gene, event, grange_start,
 #'   grange_end}). Unfiltered: one row per variant in the fit.
+#' @param cvResult Optional cross-validation payload (list with
+#'   \code{samplePartition}, \code{predictions}, \code{performance}) recorded
+#'   when fine-mapping is run with \code{cvFolds > 1}. \code{NULL} otherwise.
 #' @return A \code{FineMappingEntry} object.
 #' @export
-FineMappingEntry <- function(variantIds, susieFit, topLoci) {
+FineMappingEntry <- function(variantIds, susieFit, topLoci, cvResult = NULL) {
   obj <- new("FineMappingEntry",
              variantIds = as.character(variantIds),
              susieFit   = susieFit,
-             topLoci    = as.data.frame(topLoci))
+             topLoci    = as.data.frame(topLoci),
+             cvResult   = cvResult)
   validObject(obj)
   obj
 }
@@ -110,6 +118,11 @@ setMethod("getVariantIds", "FineMappingEntry",
 setMethod("getSusieFit", "FineMappingEntry",
           function(x, ...) x@susieFit)
 
+#' @rdname getCvResult
+#' @export
+setMethod("getCvResult", "FineMappingEntry",
+          function(x, ...) x@cvResult)
+
 #' @rdname getTopLoci
 #' @export
 setMethod("getTopLoci", "FineMappingEntry",
@@ -223,10 +236,13 @@ setMethod("adjustPips", "FineMappingEntry",
         }
       }
     }
+    # cvResult is sample-indexed (partition + held-out predictions), not
+    # variant-indexed, so it carries through a variant-subsetting unchanged.
     new("FineMappingEntry",
         variantIds = common,
         susieFit   = fit,
-        topLoci    = newTopLoci)
+        topLoci    = newTopLoci,
+        cvResult   = x@cvResult)
   })
 
 #' @export

R/QtlDataset.R

@@ -26,14 +26,18 @@ setClass("QtlDataset",
     xvarCutoff         = "numeric",
     imissCutoff        = "numeric",
     keepSamples        = "character",
-    keepVariants       = "character"
+    keepVariants       = "character",
+    keepIndel          = "logical"
   ),
+  prototype = prototype(keepIndel = TRUE),
   validity = function(object) {
     errors <- character()
     if (length(object@study) != 1L || !nzchar(object@study))
       errors <- c(errors, "'study' must be a single non-empty character string")
     if (length(object@scaleResiduals) != 1L)
       errors <- c(errors, "'scaleResiduals' must be a single logical value")
+    if (length(object@keepIndel) != 1L || is.na(object@keepIndel))
+      errors <- c(errors, "'keepIndel' must be a single logical value")
     for (nm in c("mafCutoff", "macCutoff", "xvarCutoff", "imissCutoff")) {
       v <- methods::slot(object, nm)
       if (length(v) != 1L || is.na(v) || !is.finite(v) || v < 0)
@@ -132,6 +136,9 @@ setClass("QtlDataset",
 #' @param genotypeCovariates Numeric matrix of genotype-derived covariates
 #'   (e.g., ancestry PCs); rows are samples.
 #' @param scaleResiduals Logical (length 1). Default \code{TRUE}.
+#' @param keepIndel Logical (length 1). When \code{FALSE}, variants whose
+#'   alleles are not single nucleotides (indels) are dropped at extraction.
+#'   Default \code{TRUE} (keep all variants).
 #' @return A \code{QtlDataset} object.
 #' @export
 QtlDataset <- function(study, genotypes, phenotypes,
@@ -142,7 +149,8 @@ QtlDataset <- function(study, genotypes, phenotypes,
                        xvarCutoff = 0,
                        imissCutoff = 0,
                        keepSamples = character(0),
-                       keepVariants = character(0)) {
+                       keepVariants = character(0),
+                       keepIndel = TRUE) {
   obj <- new("QtlDataset",
              study              = as.character(study),
              genotypes          = genotypes,
@@ -154,7 +162,8 @@ QtlDataset <- function(study, genotypes, phenotypes,
              xvarCutoff         = as.numeric(xvarCutoff),
              imissCutoff        = as.numeric(imissCutoff),
              keepSamples        = as.character(keepSamples),
-             keepVariants       = as.character(keepVariants))
+             keepVariants       = as.character(keepVariants),
+             keepIndel          = isTRUE(keepIndel))
   validObject(obj)
   obj
 }
@@ -271,6 +280,12 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
   unique(idx)
 }
 
+# Internal: keepIndel slot read, tolerant of QtlDataset objects serialized
+# before the slot existed (treat a missing slot as TRUE = keep indels).
+.qtlKeepIndel <- function(x) {
+  isTRUE(tryCatch(x@keepIndel, error = function(e) TRUE))
+}
+
 # Internal: extract the panel dosage block (samples x variants) for the
 # requested region, narrow to the requested sample set, and apply lazy QC
 # (per-sample imiss filter, then per-variant max(mafCutoff,
@@ -283,6 +298,11 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
 #   variantIds : character vector of kept variant IDs (= colnames(geno))
 #   sampleIds  : character vector of kept sample IDs (= rownames(geno))
 #   maf        : numeric vector of per-variant MAF for kept variants
+#   af         : numeric vector of per-variant effect-allele (A1) frequency
+#                for kept variants. Directional (NOT folded to the minor
+#                allele): the frequency of the dosage-counted allele, which
+#                is A1 by the same convention the marginal betas use. `maf`
+#                is `pmin(af, 1 - af)`.
 .qtlExtractBlock <- function(x, traitId = NULL, region = NULL,
                              cisWindow = NULL, samples = NULL) {
   gr <- .qtlResolveVariantRegion(x, traitId = traitId, region = region,
@@ -294,7 +314,8 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
                           dimnames = list(x@genotypes@sampleIds, character(0))),
       variantIds = character(0),
       sampleIds  = x@genotypes@sampleIds,
-      maf        = numeric(0)
+      maf        = numeric(0),
+      af         = numeric(0)
     ))
   }
 
@@ -310,7 +331,27 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
                             dimnames = list(character(0), character(0))),
         variantIds = character(0),
         sampleIds  = character(0),
-        maf        = numeric(0)
+        maf        = numeric(0),
+        af         = numeric(0)
+      ))
+    }
+  }
+
+  # Drop indels (variants whose A1/A2 are not single nucleotides) unless kept,
+  # before materialization so we do not extract dosage we will immediately drop.
+  if (!.qtlKeepIndel(x)) {
+    si <- x@genotypes@snpInfo
+    snpMask <- nchar(as.character(si$A1[snpIdx])) == 1L &
+               nchar(as.character(si$A2[snpIdx])) == 1L
+    snpIdx <- snpIdx[snpMask]
+    if (length(snpIdx) == 0L) {
+      return(list(
+        geno       = matrix(numeric(0), nrow = 0L, ncol = 0L,
+                            dimnames = list(character(0), character(0))),
+        variantIds = character(0),
+        sampleIds  = character(0),
+        maf        = numeric(0),
+        af         = numeric(0)
       ))
     }
   }
@@ -335,7 +376,8 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
       geno       = dosage[integer(0), , drop = FALSE],
       variantIds = colnames(dosage),
       sampleIds  = character(0),
-      maf        = rep(NA_real_, ncol(dosage))
+      maf        = rep(NA_real_, ncol(dosage)),
+      af         = rep(NA_real_, ncol(dosage))
     ))
   }
   dosage <- dosage[keep, , drop = FALSE]
@@ -356,6 +398,10 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
     nObs <- colSums(!is.na(dosage))
     sumD <- colSums(dosage, na.rm = TRUE)
     p <- ifelse(nObs > 0L, sumD / (2 * nObs), NA_real_)
+    # `p` is the frequency of the dosage-counted allele (A1, the effect
+    # allele) -> exported directly as the directional `af`. `mafVec` folds
+    # it to the minor allele for the QC threshold and the `maf` accessor.
+    afVec <- p
     mafVec <- pmin(p, 1 - p)
     effectiveMaf <- max(x@mafCutoff, if (nSamp > 0L)
       x@macCutoff / (2 * nSamp) else 0)
@@ -371,8 +417,10 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
     }
     dosage <- dosage[, keepVarMask, drop = FALSE]
     mafVec <- mafVec[keepVarMask]
+    afVec  <- afVec[keepVarMask]
   } else {
     mafVec <- numeric(0)
+    afVec  <- numeric(0)
   }
 
   # Mean-impute remaining missing dosage cells so downstream linear
@@ -392,7 +440,8 @@ setMethod("getScaleResiduals", "QtlDataset", function(x) x@scaleResiduals)
     geno       = dosage,
     variantIds = colnames(dosage),
     sampleIds  = rownames(dosage),
-    maf        = mafVec
+    maf        = mafVec,
+    af         = afVec
   )
 }
 
@@ -416,6 +465,18 @@ setMethod("getMaf", "QtlDataset",
     out
   })
 
+#' @rdname getAf
+#' @export
+setMethod("getAf", "QtlDataset",
+  function(x, traitId = NULL, region = NULL, cisWindow = NULL,
+           samples = NULL, ...) {
+    block <- .qtlExtractBlock(x, traitId = traitId, region = region,
+                              cisWindow = cisWindow, samples = samples)
+    out <- block$af
+    names(out) <- block$variantIds
+    out
+  })
+
 #' @rdname getPhenotypes
 #' @export
 setMethod("getPhenotypes", "QtlDataset",

R/QtlFineMappingResult.R

@@ -165,6 +165,14 @@ setMethod("getCs", "QtlFineMappingResult",
     getCs(entry, coverage = coverage)
   })
 
+#' @rdname getCvResult
+#' @export
+setMethod("getCvResult", "QtlFineMappingResult",
+  function(x, study = NULL, context = NULL, trait = NULL, method = NULL, ...) {
+    entry <- getFineMappingResult(x, study, context, trait, method)
+    getCvResult(entry)
+  })
+
 #' @rdname getTopLoci
 #' @export
 setMethod("getTopLoci", "QtlFineMappingResult",